|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 23040 (0x5a00)
Types: TextFile
Names: »ch2txt«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦d3633872f⟧ »cnch«
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
└─⟦4334b4c0b⟧
└─⟦d3633872f⟧ »cnch«
└─⟦this⟧
(mode list.yes
lookup ch2txt
lookup rydterm
if ok.yes
mode 11.yes
clear user ch2
lookup ch2det
if ok.no
contract from.cnch ch2det
lookup ryproc alutproc rydseg rydstruct ryglobal statevar stateloop
if ok.no
contract from.crypr ryproc alutproc rydseg rydstruct ryglobal statevar stateloop
ch2=set 290 disc3
scope user ch2
lookup rydrun
if ok.yes
mode 10.yes
lookup rydlist
if ok.yes
mode 15.yes
if 11.yes 15.no
o rydterm
clear rydlist rydgtest
head cpu
ch2=algol xref.no connect.no fp.yes
head cpu
;ch2 mg12 l1.0 l2.1 x1.0.52 x2.0.45 print.no
;ch2 mg12 mix.yes l1.1
;ch2 mg12 print.no l1.1 x1.0.0.5 x2.0.67
;ch2 mg12 mix.yes l1.2 l2.1
;ch2 mg12 print.no l1.2 l2.1 x1.0.62 x2.0.42
;ch2 mg12 min.yes iterations.25 l1.0 l2.1 x1.0.52 x2.0.45
;ch2 mg12 min.yes l1.1 l2.1 x1.0.0.5 x2.0.67
;ch2 mg12 min.yes l1.2 l2.1 x1.0.62 x2.0.42
mode 10.no 15.no list.no)
1980-07-30
2-channelprogram
begin
boolean xu,print,endstate,mix,min,plot,weight,lufanoplot,
type1;
real Ip1,Ip2,I1au,I2au,ind,rel,Ecm,nstar1,nstar2,sum1,sum2,
wsum,digam1,digam2,
maxns,my1,my2,nsmin,eps,eps1,x,y,
ny1,ny2,pymax;
integer l1,l2,L,n,l,J,states,i,j,maxiter,maxiter1,c,char,par,Jc,si;
array PSI(1:6),pname(1:3);
real procedure ny1ny2(PSI,ny2);
value ny2; real ny2;
array PSI;
begin
real x,my1,my2,my10,my20,my11,my21,dny,ny,ny1,c2,c1,
s22,s21,s1,s2,st2,ct2;
my10:=PSI(1); my11:=PSI(4);
my20:=PSI(2); my21:=PSI(5);
my2:=my20-my21/ny2/ny2/2;
ct2:=cos(PSI(3))**2;
st2:=sin(PSI(3))**2;
s22:=sin(pi*(ny2+my2));
s2:=sin(pi*my2);
c2:=cos(pi*my2);
ny1:=ny:=0;
for my1:=my10,my10-my11/2/ny1/ny1 while abs dny>'-7 do
begin
ny:=ny1;
s21:=sin(pi*(ny2+my1));
c1:=cos(pi*my1);
s1:=sin(pi*my1);
x:=(st2*c2*s21+ct2*c1*s22)/
(st2*s2*s21+ct2*s1*s22);
ny1:=arctan(-x/pi);
dny:=ny1-ny;
end loop;
ny1ny2:=ny1;
end ny1ny2;
<*algol copy.platsymtxt*>
algol list.off copy.ryproc;
algol list.off copy.rydstruct;
algol list.off copy.alutproc;
algol list.off copy.rydseg;
real procedure cot(x);
value x; real x;
cot:=if x=0.0 then maxreal else cos(x)/sin(x);
real procedure B(ny,l);
value ny,l;
real ny; integer l;
begin
real t,ny2;
ny2:=ny*ny;
t:=1/ny2/ny2**(l+1);
for l:=l step -1 until 1 do t:=t*(ny2-l*l);
B:=t;
end B;
real procedure g(x,l);
value x,l;
real x; integer l;
begin
real y,z;
z:=y:=digamma(x-l);
for i:=2*l+1 step -1 until 1 do z:=z+1/(x+i-l-1);
g:=B(x,l)*(y+z-2*ln(x))/2/pi;
end g;
real procedure compnstar(par,ny,nyc,states);
value states; integer states;
array par,ny,nyc;
begin
integer i,j,k;
real theta,ns,nsmin,nsmax,y0,y1,nstar1,
my10,my11,my20,my21,ns1,dns1;
theta:=par(3);
my10:=par(1); my11:=par(4);
my20:=par(2); my21:=par(5);
j:=1;
nsmin:=ny(1,1)-.2;
nsmax:=ny(1,states)+.3;
for i:=1 step 1 until states do nyc(1,i):=0;
y0:=det(nsmin,my10,my11,my20,my21,theta);
dns1:=.1;
for ns:=nsmin+dns1 step dns1 until nsmax do
begin
y1:=det(ns,my10,my11,my20,my21,theta);
if sign(y1)<>sign(y0) then
begin
ns1:=nstar1:=ns-.5*dns1;
if -,zero1(nstar1,det(nstar1,my10,my11,my20,my21,theta),
nstar1-dns1*.5,dns1*.5+nstar1,'-6) then
write(out,"nl",1,"*",2,<:ch2 bisec :>,
det(ns1-dns1*.5,my10,my11,my20,my21,theta),
det(ns1+dns1*.5,my10,my11,my20,my21,theta),
det(ns1,my10,my11,my20,my21,theta),ns1,nstar2,"nl",1) else
nyc(1,j):=nstar1;
j:=j+1;
end y0<>y1;
y0:=y1;
end ns loop;
end compnstar;
real procedure det(nstar1,my10,my11,my20,my21,theta);
value nstar1,my10,my11,my20,my21,theta;
real nstar1,my10,my11,my20,my21,theta;
begin
real mu1,mu2,nyp,nym,nup,num,cnyp,snyp,cnym,snym,
cnup,snup,cnum,snum,e1,e2,ct2;
ct2:=cos(2*theta);
algol list.on copy.ch2det;
end det;
real procedure ddet(d,nstar1,my10,my11,my20,my21,theta);
value nstar1,my10,my11,my20,my21,theta;
real nstar1,my10,my11,my20,my21,theta;
array d;
begin
real mu1,mu2,e1,e2,t2,st2,ct2,nyp,nym,nup,num,
cnyp,snyp,cnym,snym,cnup,snup,cnum,snum,d1,d2;
e1:=.5/nstar1/nstar1;
mu1:=my10-my11*e1;
e2:=I2au-I1au+.5/nstar1/nstar1;
nstar2:=sqrt(.5/e2);
mu2:=my20-my21*e1;
t2:=2*theta;
st2:=sin(t2); ct2:=cos(t2);
nyp:=pi*(nstar1+nstar2);
nym:=pi*(nstar1-nstar2);
cnyp:=cos(nyp); snyp:=sin(nyp);
cnym:=cos(nym); snym:=sin(nym);
nup:=pi*(mu1+mu2);
num:=pi*(mu1-mu2);
cnup:=cos(nup); snup:=sin(nup);
cnum:=cos(num); snum:=sin(num);
d1:=.5*(cnyp*snup+cnup*snyp);
d2:=.5*(cnym*snum+ct2*cnum*snym);
d(1):=d1-d2;
d(2):=d1+d2;
if my11<>0 then d(4):=-e1*d(1);
if my21<>0 then d(5):=-e1*d(2);
d(3):=st2*snum*snym;
end ddet;
real procedure d2det(d2,nstar1,my10,my11,my20,my21,theta);
value nstar1,my10,my11,my20,my21,theta;
real nstar1,my10,my11,my20,my21,theta;
array d2;
begin
real mu1,mu2,e1,e2,t2,st2,ct2,nyp,nym,nup,num,
cnyp,snyp,cnym,snym,cnup,snup,cnum,snum,da,db;
e1:=.5/nstar1/nstar1;
mu1:=my10-my11*e1;
e2:=I2au-I1au+.5/nstar1/nstar1;
nstar2:=sqrt(.5/e2);
mu2:=my20-my21*e1;
t2:=2*theta;
st2:=sin(t2); ct2:=cos(t2);
nyp:=pi*(nstar1+nstar2);
nym:=pi*(nstar1-nstar2);
cnyp:=cos(nyp); snyp:=sin(nyp);
cnym:=cos(nym); snym:=sin(nym);
nup:=pi*(mu1+mu2);
num:=pi*(mu1-mu2);
cnup:=cos(nup); snup:=sin(nup);
cnum:=cos(num); snum:=sin(num);
da:=.5*(cnyp*cnup-snup*snyp);
db:=.5*(cnym*cnum-ct2*snum*snym);
d2(1):=da-db;
d2(3):=da+db;
d2(6):=2*ct2*snym*snum;
d2(4):=st2*cnum*snym;
d2(5):=-d2(4);
d2(2):=.5*(cnym*cnum+cnyp*cnup
-snyp*snup-ct2*snym*snum);
if my11<>0 then
begin
d2(10):=e1*e1*d2(1);
d2(7):=-e1*d2(1);
d2(8):=-e1*d2(2);
d2(9):=-e1*d2(4);
end;
if my21<>0 then
begin
d2(15):=e1*e1*d2(3);
d2(11):=-e1*d2(2);
d2(12):=-e1*d2(3);
d2(13):=-e1*d2(5);
d2(14):=e1*e1*d2(2);
end;
end d2det;
procedure printsym(M,par,fac);
value par,fac; integer par;
real fac;
array M;
begin
integer i,j;
for i:=1 step 1 until par do
for j:=1 step 1 until i do
write(out,"nl",if j=1 then 1 else 0,
<< -d.dddddddd>,fac*M(i*(i-1)//2+j));
end printsym;
procedure moment(M,M2,ny,nyc,PSI,par,states);
value par,states;
integer par,states;
array M,M2,ny,nyc,PSI;
begin
integer i,j,k,idx;
real my10,my11,my20,my21,theta,detv;
array d(1:6),corr,d2,dQ2(1:15),ds(1:6);
my10:=PSI(1); my11:=PSI(4);
my20:=PSI(2); my21:=PSI(5);
theta:=PSI(3);
for i:=1 step 1 until par do
for j:=1 step 1 until i do
begin
idx:=i*(i-1)//2+j;
dQ2(idx):=M2(idx):=M(idx):=0;
end;
d(4):=d(5):=d(6):=0;
for i:=1 step 1 until 6 do ds(i):=0;
for k:=1 step 1 until states do
begin
detv:=det(ny(1,k),my10,my11,my20,my21,theta);
ddet(d,ny(1,k),my10,my11,my20,my21,theta);
d2det(d2,ny(1,k),my10,my11,my20,my21,theta);
for i:=1 step 1 until 6 do ds(i):=ds(i)+2*detv*d(i);
for i:=1 step 1 until par do
for j:=1 step 1 until i do
begin
idx:=i*(i-1)//2+j;
M(idx):=M(idx)+d(i)*d(j);
M2(idx):=M2(idx)+d2(idx);
dQ2(idx):=dQ2(idx)+d(i)*d(j)-detv*d2(idx);
end M loop;
end states;
if true then
begin
write(out,"nl",2,<:Moment matrix :>,my10,my11,my20,my21,theta);
write(out,"nl",1,<: :>,ds(1),ds(4),ds(2),ds(5),ds(3));
printsym(M,par,1);
end M write;
for i:=1 step 1 until par do
for j:=1 step 1 until i do
begin
idx:=i*(i-1)//2+j;
corr(idx):=if i=j then 1 else M(idx)/sqrt(M(i*(i+1)//2)**2+M(j*(j+1)//2)**2);
end corr;
if false then
begin
write(out,"nl",2,<:correlation matrix:>);
printsym(corr,par,1);
end;
if true then
begin
write(out,"nl",1,<:2. derivative det F :>);
printsym(M2,par,1);
end M2;
if true then
begin
write(out,"nl",2,<:2. derivative Q:>);
printsym(dQ2,par,1);
end;
overflows:=0;
if -,symin(par,M) then alarm("nl",1,"*",1,<:ch2 symin :>) else
if false then
begin
write(out,"nl",2,<:invert M matrix:>);
if overflows>0 then write(out,"sp",2,"*",3);
printsym(M,par,1.0);
write(out,"nl",2);
end inv;
overflows:=-1;
end moment;
<*procedure plotlufano(states,ny,param);
value states; integer states;
array ny,param;
begin
real x,xmax;
integer i;
xmax:=ny(1,states);
plotform(0,(if xmax<12 then 2 else 1)*xmax+2,13);
setmargin(3,plotyform-1);
writeplot(ff,1,<:det<22>F(n*)<22> :>,false add ryalf(l1+256),1,<:-series:>,"sp",2);
plotatsym(S,atno,Z);
plotadmini(0,xmax+1,-pymax,pymax,0);
plotgraph(x,det(x,PSI(1),PSI(4),PSI(2),PSI(5),PSI(3)),ny(1,1)-.1,xmax,.05);
for i:=states step -1 until 1 do
plotpoint(ny(1,i),det(ny(1,i),PSI(1),PSI(4),PSI(2),PSI(5),PSI(3)),1);
plotclose;
end plotlufano;*>
real procedure detsum(states,par,ny,my,Bv,G,W);
value states; integer states;
array par,ny,my,Bv,G,W;
if states<=0 then detsum:=0 else
if type1 then
begin
integer i;
real d,det,my10,my11,my20,my21,theta,ct2,e1,e2,
mu1,mu2,nyp,nym,nup,num,cnyp,snyp,cnym,snym,
cnup,snup,cnum,snum;
d:=0;
ct2:=cos(2*par(3));
my10:=par(1); my11:=par(4);
my20:=par(2); my21:=par(5);
for i:=1 step 1 until states do
begin
nstar1:=ny(1,i);
algol copy.ch2det;
d:=d+det*det*W(i);
end states;
detsum:=d;
end else
begin
own integer eval;
integer i;
real D,p,q,r,s,mu1,mu2,ny1,ny2,theta,ctheta,stheta,
my10,my11,my20,my21,cot1,cot2,s11,s12,s21,s22;
D:=0;
theta:=par(3);
my10:=par(1); my11:=par(4);
my20:=par(2); my21:=par(5);
ctheta:=cos(theta); stheta:=sin(theta);
for i:=1 step 1 until states do
begin
ny1:=ny(1,i); ny2:=ny(2,i);
mu1:=my10-my11/ny(1,i)**2/2;
mu2:=my20-my21/ny(1,i)**2/2;
if xu then
begin
cot1:=cot(pi*mu1); cot2:=cot(pi*mu2);
p:=(ctheta*ctheta*cot1+
stheta*stheta*cot2+
G(1,i))/Bv(1,i);
q:=stheta*stheta*(cot2-cot1)/
sqrt(Bv(1,i)*Bv(2,i));
r:=(stheta*stheta*cot1+
ctheta*ctheta*cot2+
G(2,i))/Bv(2,i);
s:=sqrt((p-r)**2+4*q*q);
mu1:=(p+r+s)/2; mu2:=(p+r-s)/2;
theta:=if abs theta<smallreal then theta else (p-mu1)/q;
mu1:=arctan(1/mu1/pi);
mu2:=arctan(1/mu2/pi);
theta:=arctan(theta);
stheta:=sin(theta); ctheta:=cos(theta);
end xu;
s11:=sin(pi*(ny1+mu1));
s12:=sin(pi*(ny1+mu2));
s21:=sin(pi*(ny2+mu1));
s22:=sin(pi*(ny2+mu2));
D:=D+(stheta**2*s21*
s12+
ctheta**2*s22*
s11)**2*W(i);
end summing over states;
eval:=eval+1;
if print and min and eval<=maxiter1 then
begin
write(out,"nl",2,<:evaluation det:>,eval,<: minimum :>,D,
"nl",1,<:my10 :>,my10,<: my11 :>,my11,
"nl",1,<:my20 :>,my20,<: my21 :>,my21,
"nl",1,<:theta :>,theta,"nl",1);
end;
detsum:=D;
end detsum;
real procedure ddetsum(states,no,par,ny);
value states,no; integer states,no;
array par,ny;
if states<=0 then ddetsum:=0 else
begin
integer i;
real my10,my11,my20,my21,theta,dd;
array d(1:5);
d(4):=d(5):=0;
theta:=par(3);
my10:=par(1); my11:=par(4);
my20:=par(2); my21:=par(5);
dd:=0;
for i:=1 step 1 until states do
begin
ddet(d,ny(1,i),my10,my11,my20,my21,theta);
dd:=dd+d(no)*2*det(ny(1,i),my10,my11,my20,my21,theta);
end summing over states;
ddetsum:=dd;
end ddetsum;
procedure mixcalc(PSI,ny,nyc,nn,ll,JJ,states);
value states; integer states;
array PSI,ny,nyc;
integer array nn,ll,JJ;
begin
real theta,stheta,ctheta,ny1,ny2,my1,my2,Na;
integer i,j,char;
array cm(1:2,1:2),A1,A2,N1,N2,Nn,D1,D2(1:2),z1,z2(1:2,1:states);
theta:=PSI(3);
stheta:=sin(theta); ctheta:=cos(theta);
if mix then write(out,"nl",2,<:calculation of mixing coefficients only:>);
write(out,"nl",2);
writeatsym(out,S,atno,Z);
write(out,"nl",1,<:state:>,"sp",if finestruct then 6 else 2,
<:z1:>,"sp",5,<:z2:>,"sp",5,<:norm:>,
"sp",4,<:Nn:>,"sp",4,<:A1:>,"sp",5,<:A2:>);
if false then write(out,"sp",5,<:D1:>,"sp",5,<:D2:>);
write(out,"sp",5,<:n*1:>,"sp",4,<:n*1(c):>);
for j:=1 step 1 until states do
begin
ny1:=ny(1,j); ny2:=ny(2,j);
my1:=PSI(1)-PSI(4)/(ny1**2)/2;
my2:=PSI(2)-PSI(5)/(ny1**2)/2;
cm(1,1):=cos(theta)*sin(pi*(ny2+my2));
cm(1,2):=-si*sin(theta)*sin(pi*(ny2+my1));
cm(2,1):=si*sin(theta)*sin(pi*(ny1+my2));
cm(2,2):=cos(theta)*sin(pi*(ny1+my1));
for i:=1,2 do
begin
Na:=cm(i,1)*cm(i,1)+cm(i,2)*cm(i,2);
A1(i):=cm(i,1)/sqrt(Na);
A2(i):=cm(i,2)/sqrt(Na);
N1(i):=(ctheta*cos(pi*(ny1+my1))*A1(i)-
stheta*cos(pi*(ny1+my2))*A2(i))**2;
N2(i):=(sin(theta)*cos(pi*(ny2+my1))*A1(i)+
ctheta*cos(pi*(ny2+my2))*A2(i))**2;
Nn(i):=sqrt(ny1**3*N1(i)+ny2**3*N2(i)
+PSI(4)*A1(i)*A1(i)+PSI(5)*A2(i)*A2(i));
z1(i,j):=(-1)**(l1//2+1)*ny1**(3/2)*(
ctheta*cos(pi*(ny1+my1))*A1(i)/Nn(i)-
stheta*cos(pi*(ny1+my2))*A2(i)/Nn(i));
z2(i,j):=(-1)**(l2//2+1)*ny2**(3/2)*(
stheta*cos(pi*(ny2+my1))*A1(i)/Nn(i)+
ctheta*cos(pi*(ny2+my2))*A2(i)/Nn(i));
D1(i):=A1(i)/Nn(i)*(ny1**(3/2));
D2(i):=A2(i)/Nn(i)*(ny2**(3/2));
<*D1 and D1 are the coefficients in the calculation of f-values
see Lu and Lee 2.13*>
end i;
write(out,"nl",1);
char:=writestate(out,l1,nn(j),ll(j),JJ(j));
for i:=1,2 do
begin
write(out,"nl",if i=1 then 0 else 1,"sp",if i=1 then 1 else char+1,
<<-dd.ddd>,
z1(i,j),z2(i,j),z1(i,j)*z1(i,j)+z2(i,j)*z2(i,j),
Nn(i),A1(i),A2(i));
if false then write(out,<<-dd.ddd>,D1(i),D2(i));
if i=1 then write(out,<< -dd.ddd>,ny1,nyc(1,j));
end;
end for states;
end mixcalc;
real procedure hkena(n,PSI,PHI,S,rho,sig,fn,maxeval);
value n,rho,sig,maxeval;
integer n,maxeval;
real rho,sig;
array PSI,PHI,S;
real fn;
comment Hooke and Jeeves direct search
n = number of variables
PSI(1:n) = starting points and subsequent base points
PHI(1:n) = points for call of fn
S(1:n) = step lengths
spsi = function value at base point
rho = reduction factor for step lengths
sig = significant figures required in PSI
fn = the function
maxeval = maximum no of function values permitted
;
begin
boolean cont;
integer l,k,eval;
array theta(1:n);
real sphi,ss,stheta,spsi;
real procedure E;
begin
for k:=1 step 1 until n do
begin
PHI(k):=PHI(k)+S(k);
sphi:=fn;
eval:=eval+1;
if ind=0 or sphi>=ss then
begin
S(k):=-S(k);
PHI(k):=PHI(k)+2*S(k);
sphi:=fn;
eval:=eval+1;
if ind=0 or sphi>=ss then PHI(k):=PHI(k)-S(k) else
ss:=sphi;
end else ss:=sphi;
end k;
printpar;
end E;
procedure printpar;
if print then
begin
write(out,"nl",3);
for l:=1 step 1 until n do write(out,PSI(l));
write(out,"nl",1,<:spsi =:>,spsi);
write(out,"nl",1,<:iteration =:>,eval);
end print;
for k:=1 step 1 until n do PHI(k):=PSI(k);
spsi:=fn;
eval:=1;
printpar;
label1:
cont:=true;
for k:=k while eval<=maxeval and cont do
begin
ss:=spsi;
for k:=1 step 1 until n do PHI(k):=PSI(k);
E;
if ss<spsi then
label2:
begin
for k:=1 step 1 until n do
begin
if PHI(k)>PSI(k) and S(k)<0 then S(k):=-S(k);
if PHI(k)<PSI(k) and S(k)>0 then S(k):=-S(k);
theta(k):=PHI(k);
PHI(k):=2*PHI(k)-PSI(k);
end for k;
stheta:=ss;
ss:=sphi:=fn;
eval:=eval+1;
printpar;
if ind=0 or ss>=stheta then
begin
for k:=1 step 1 until n do
PSI(k):=theta(k);
spsi:=stheta;
goto label1;
end else
begin
for k:=1 step 1 until n do if abs(PHI(k)-PSI(k))>.5*abs S(k) then
goto label2;
end;
end ss<sphi;
cont:=false;
for k:=1 step 1 until n do cont:=cont or abs S(k)>=sig*abs PSI(k);
for k:=1 step 1 until n do S(k):=rho*S(k);
if print then write(out,"nl",2,<:step lengths reduced:>);
end main loop;
maxiter:=eval;
hkena:=fn;
end hkena;
algol list.off copy.ryglobal;
comment START OF PROGRAM;
readbfp(<:type1:>,type1,true);
readbfp(<:plot:>,plot,false);
readbfp(<:plf:>,lufanoplot,false);
readbfp(<:min:>,min,false);
readbfp(<:mix:>,mix,false);
readbfp(<:list:>,list,true);
readbfp(<:xu:>,xu,false);
readbfp(<:print:>,print,false);
readbfp(<:weight:>,weight,true);
readifp(<:si:>,si,-1);
readifp(<:par:>,par,5);
readifp(<:iterations:>,maxiter,1000);
maxiter1:=maxiter;
readrfp(<:ymax:>,pymax,.5);
readrfp(<:nsmin:>,nsmin,1.5);
readrfp(<:nsmax:>,nsmax,35);
<*if plot or lufanoplot then
begin
if readsfp(<:plotter:>,pname) then setplotname(string inc(pname),0) else
setplotname(<:houstona:>,0);
end plot*>;
if mix then
begin
read(in,PSI(1),PSI(4),
PSI(2),PSI(5),
PSI(3));
write(out,"nl",1,<:parameters read::>,
"nl",1,PSI(1),PSI(4),"nl",1,
PSI(2),PSI(5),"nl",1,PSI(3));
end;
connectinp(1);
connectlso;
readchar(in,i); <*dummy for a repeatchar in readatsym*>
readatsym(in,S,atno,Z);
write(out,"ff",0,"nl",4);
writeatsym(out,S,atno,Z);
read(in,i,i);
if readifp(<:l1:>,l1,0) then l1:=2*l1; Ip1:=readr(<:Ip1:>);
write(out,"sp",4,false add ryalf(l1+256),1,<:-series:>);
I1au:=Ip1*100/ryinf/Z/Z/2;
write(out,"sp",4,I1au);
if readifp(<:l2:>,l2,2) then l2:=2*l2; Ip2:=readr(<:Ip2:>);
write(out,"sp",4,false add ryalf(l2+256),1,<:-limit:>);
I2au:=Ip2*100/ryinf/Z/Z/2;
write(out,"sp",5,I2au);
readifp(<:j:>,Jc,-1);
if Jc>=0 then
begin
finestruct:=true;
write(out,<: J = :>,Jc,if Jc mod 2=1 then <:/2:> else <::>);
end;
readrfp(<:eps:>,eps,8);
eps1:=eps:=10**(-eps);
readrfp(<:steplengths:>,rel,6);
rel:=10**(-rel);
states:=nsmax-nsmin+10;
begin
array PHI,ACCP,ACCP1(1:6),Bv,ny,nyc,my,G(1:2,1:states),W(1:states);
integer array nn,ll,JJ(1:states);
procedure calcparam(par,states,PSI,ny,my,W);
value par,states; integer par,states;
array PSI,ny,my,W;
begin
array PHI(1:6),M,M2(1:par*(par+1)//2),dPSI(1:6);
integer i;
real minval,s2val;
ind:=1;
for i:=1 step 1 until 6 do ACCP(i):=ACCP1(i);
for i:=1 step 1 until par do PHI(i):=PSI(i);
for i:=par+1 step 1 until 6 do PSI(i):=PHI(i):=0;
eps:=eps1;
minval:=if mix then detsum(states,PSI,ny,my,Bv,G,W) else
if min then minimum(par,i,PHI,FN(PHI),delta(i,PHI),eps,PSI) else
hkena(par,PSI,PHI,ACCP,.5,rel,detsum(states,PHI,ny,my,Bv,G,W),maxiter1);
compnstar(PSI,ny,nyc,states);
s2val:=sqrt(detsum(states,PSI,ny,my,Bv,G,W)/(states-par));
if weight then write(out,"nl",1,<:weight function 1/n*/n*:>) else
write(out,"nl",1,<:weight function 1:>);
write(out,"nl",2,<<d>,par,<:-parameters:>,
"nl",1,<:deviation s :>,<< d.dddddddd>,s2val);
moment(M,M2,ny,nyc,PSI,par,states);
write(out,"nl",2,<:covariance matrix:>);
printsym(M,par,s2val**2);
for i:=1 step 1 until par do
dPSI(i):=s2val*sqrt(M(i*(i+1)//2));
write(out,"nl",1,<< d.dddddd>,<:my1(0):>,PSI(1),
<: +- :>,<< d.d'-d>,dPSI(1));
if par>3 then write(out,<: my1(1):>,<< d.ddddd>,PSI(4),
<: +- :>,<< d.d'-d>,dPSI(4));
write(out,"nl",1,<:my2(0):>,<< d.dddddd>,PSI(2),
<: +- :>,<< d.d'd>,dPSI(2));
if par>4 then write(out,<: my2(1):>,<< d.ddddd>,PSI(5),
<: +- :>,<< d.d'd>,dPSI(5));
write(out,"nl",1,<:theta :>,<< d.dddddd>,PSI(3),
<: +- :>,<< d.d'd>,dPSI(3),"nl",1,
<:cos :>,<< d.dddddd>,cos(PSI(3)),<: sin :>,sin(PSI(3)));
if -,mix and -,min then write(out,"nl",1,<:iterations:>,<< ddddd>,maxiter,"sp",2);
write(out,"nl",if mix or min then 1 else 0,<:minimum :>,
<< d.dddd'-dd>,minval);
if min then write(out,"sp",4,<:gradient norm:>,<<d.d'-dd>,eps);
if -,mix and -,min then write(out,<: steplength :>,<< -d.ddd'-dd>,ACCP(1));
end calcparam;
real procedure FN(x);
array x;
FN:=detsum(states,x,ny,my,Bv,G,W);
real procedure delta(i,x);
value i; integer i;
array x;
delta:=ddetsum(states,i,x,ny);
for i:=1 step 1 until 6 do
begin
if -,mix then PSI(i):=0;
ACCP(i):=ACCP1(i):=if par<=3 then .05 else .01;
end;
if -,mix then
begin
if readrfp(<:x1:>,PSI(1),.5) then
write(out,"nl",1,<:x1 = :>,PSI(1));
if readrfp(<:x2:>,PSI(2),.5) then
write(out,"nl",1,<:x2 = :>,PSI(2));
if readrfp(<:theta:>,PSI(3),.6) then
write(out,"nl",1,<:theta = :>,PSI(3),cos(PSI(3)),sin(PSI(3)));
end;
states:=0;
endstate:=false;
wsum:=0;
if list then write(out,"nl",2,<:state:>,"sp",if finestruct then 10 else 6,<:Ecm:>,
"sp",6,<:n*1:>,"sp",8,<:n*2:>,"sp",7,<:my1:>,"sp",8,<:my2:>);
maxns:=0;
for states:=states while -,endstate do
begin
readstate(in,L,n,l,J);
if L<0 then L:=l;
read(in,Ecm);
repeatchar(in);
for c:=readchar(in,char) while c>=7 and char<>25 do;
if char=101 or char=105 or char=60 then endstate:=true;
repeatchar(in);
if L=l1 and (Jc>=0 =>Jc=J) then
begin
nstar1:=if Ecm<Ip1 then Z*sqrt(ryinf/100/(Ip1-Ecm)) else 0;
nstar2:=if Ecm<Ip2 then Z*sqrt(ryinf/100/(Ip2-Ecm)) else 0;
if nstar1>=nsmin and nstar1<=nsmax then
begin
j:=states:=states+1;
nn(j):=n;
ll(j):=l;
JJ(j):=J;
ny(1,j):=nstar1;
my(1,j):=entier nstar1 - nstar1+1;
ny(2,j):=nstar2;
my(2,j):=nstar2 - entier nstar2;
W(j):=if weight then 1/nstar1/nstar1 else 1;
wsum:=wsum+W(j);
if nstar1>maxns then maxns:=nstar1;
if print or list then
begin
write(out,"nl",1);
writestate(out,L,n,l,J);
write(out,"sp",2,<<ddddddddd.d>,Ecm,<< dd.ddddd>,nstar1,nstar2,
my(1,j),my(2,j));
end;
if xu then
begin
digam1:=sum1:=digamma(nstar1-.5*l1);
for i:=1 step 1 until l1+1 do sum1:=sum1+1/ny(1,j)+i-l1/2;
digam2:=sum2:=digamma(nstar2-.5*l2);
for i:=1 step 1 until l2+1 do sum2:=sum2+1/ny(2,j)+i-l2/2;
Bv(1,j):=B(ny(1,j),l1/2);
G(1,j):=(digam1+sum1-2*ln(ny(1,j)))*Bv(1,j)/2/pi;
Bv(2,j):=B(ny(2,j),l2/2);
G(2,j):=(digam2+sum2-2*ln(ny(2,j)))*Bv(2,j)/2/pi;
end xu;
end nsmin and nsmax;
end l=l1;
end states;
wsum:=states/wsum;
for j:=1 step 1 until states do W(j):=W(j)*wsum;
for j:=if mix then par else 3 step 1 until par do
begin
write(out,"ff",1,"nl",1);
write(out,"nl",1,"sp",30);
writedate(out,0.0);
write(out,"nl",1,"sp",20);
writeatsym(out,S,atno,Z);
write(out,"nl",2,<:Ip1 = :>,<<dd ddd ddd.dd>,Ip1,"sp",4,
false add ryalf(l1+256),1,<:-series:>,
<< d.ddddd>,I1au);
write(out,"nl",1,<:Ip2 = :>,<<dd ddd ddd.dd>,Ip2,"sp",4,
false add ryalf(l2+256),1,<:-limit:>,
<< d.ddddd>,I2au);
write(out,"nl",2);
if abs PSI(3)<'-2 and -,mix then PSI(3):=.2;
calcparam(j,states,PSI,ny,my,W);
<*if lufanoplot then plotlufano(states,ny,PSI);*>
mixcalc(PSI,ny,nyc,nn,ll,JJ,states);
write(out,"nl",4,"ff",1);
end parameters;
<*if plot then
begin
plotform(0,12,14);
setmargin(1,13);
writeplot(ff,1,<:z2**2 :>);
plotatsym(S,atno,Z);
writeplot(sp,3,false add ryalf(l1+256),1,<:-series:>);
if Jc>=0 then writeplot(sp,2,<:J =:>,Jc,if Jc mod 2=1 then <:/2:> else <::>);
plotsubform(0,12,0,12,false);
plotscale(0,maxns,0,1.1);
plotframe(0.0,0.0);
for j:=1 step 1 until states do plotpoint(ny(1,j),z2(1,j)**2,2);
plotform(0,12,14);
setmargin(1,13);
writeplot(ff,1,<:Lu-Fano plot :>);
plotatsym(S,atno,Z);
writeplot(sp,3,false add ryalf(l1+256),1,<:-series:>);
plotsubform(0,12,0,12,false);
plotscale(-.1,1.1,-.1,1.1);
plotframe(0.0,0.0);
for j:=1 step 1 until states do
begin
x:=-my(1,j); y:=ny(2,j);
plotpoint(x,y,1);
if list then write(out,j,x,y);
end;
plotgraph(x,ny1ny2(PSI,x),'-2,1,1/100);
end plot*>;
end block for array
if fpout then closeout;
end;
;0.049633,0.252385,.300763,-1.73794,.008038
;.6153,.27465,2.616,-1.846,-0.0482
▶EOF◀