|
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: 20736 (0x5100)
Types: TextFile
Names: »gonimfreq22«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦this⟧ »gonimfreq22«
o gopr
nimfreq22=algol xref.no index.no list.yes
begin
integer i,j,k,l,m,iev,nev,it,pm,N,J,J1,J2,JJ,T,nrc,m0,k0,k01,k02,
Kmin,Kmax,dM,M1,M2,N1,Mmin1,Mmin,Mmax,gamma,lev,nfit,spinf;
boolean Jeven, meven, plus, closeres, disc, spind;
array khi,hf(1:3);
zone res(128,1,stderror);
nrc:= 28; read(in,dM,T); disc:= dM<0; spind:= false;
if disc then begin
readhead(in,hf,1); i:= 1;
open(res,4,string hf(increase(i)),0); inrec(res,128); dM:= -dM;
nfit:= nrc-(res(2) extract 12);
end else nfit:= 1;
begin
real X,Y,Z,A,B,C,E,F,freq,obs,w,LS,b1,b2,b3,b4,b5,b6,b7,b9,
kb,kc,kab,r1,r2,r3,r4,r5,r6,r7,r8,r9;
integer array no(1:nrc), tail(1:10);
array RC,dy(1:nrc), tau,D(1:14), error(1:nfit,1:nfit);
zone ev(128,1,stderror);
open(ev,4,<:nimev:>,0);
if monitor(42,ev,0,tail)>0 then create(ev,1,tail);
nev:= tail(10);
if disc then begin
i:= res(2) shift (-12) extract 12;
if i<>nrc then begin write(out,<:
***This version requires nrc = :>,<<dd>,nrc,<: not: :>,i);
goto stop end;
k:= (nrc+3)//4; j:= k+16; m:= k*4-nrc;
if m<>0 then m:= -m*12;
for i:=nrc step -1 until 1 do begin
no(i):= res(j) shift m extract 12; m:= m-12;
if m<-36 then begin j:= j-1; m:= 0 end;
end;
for i:=1 step 1 until 12 do D(i):= res(2+i); i:= 1;
write(out,nl,3,string D(increase(i)),nl,2);
spind:= D(1) shift (-32) = real<:D3:> shift (-32);
l:= res(1) extract 24;
k:= (nrc+nrc+l+5)//4+l+18;
for i:=1 step 1 until nrc do RC(i):= res(k+i);
for i:=1 step 1 until nfit do write(out,<<ddd>,no(i));
if spind then write(out,<: D3-spin:>);
for i:=1 step 1 until nrc do
write(out,nl,1,<<-dddddd.ddddd000000>,RC(i));
k:= k+nrc+nrc-nfit; w:= res(15);
for i:=1 step 1 until nfit do
for j:=1 step 1 until nfit do begin
if k=128 then begin k:= 0; inrec(res,128) end;
k:= k+1; error(i,j):= res(k)*w;
if i>j and error(i,j)-error(j,i)<>0 then begin
write(out,<:
***error matrix: i,j,k,l =:>,<<-dddd>,i,j,k,l);
goto stop end;
end;
write(out,nl,2); forceout(out);
close(res,true)
end else read(in,RC);
read(in,khi);
closeres:= outmedium(res);
pm:= dM//3; dM:= pm*3;
if disc then begin
i:= 1; write(res,<:<12>:>,nl,3,string D(increase(i)))
end;
write(res,nl,2,<:A :>,<<-ddddddd.ddd>,RC(1),
nl,1,<:B :>,RC(2),nl,1,<:C :>,RC(3),
nl,1,<:F :>,RC(4),nl,1,<:Af :>,RC(5),
nl,1,<:V6 :>,RC(6),nl,1);
for i:=1 step 1 until 9 do
write(res,<:T:>,<<d>,i,<< -dddddd.ddd>,RC(i+6),nl,1);
for i:=10 step 1 until 13 do
write(res,<:T:>,<<dd>,i,<< -dddddd.ddd>,RC(i+6),nl,1);
write(res,<< -ddddddd.ddd>,<:kb :>,RC(20),nl,1,<:kc :>,RC(21),nl,1,
<:kab:>,RC(22),nl,1);
for i:=4 step 1 until 7 do
write(res,<:b:>,<<d>,i,<< -ddddddd.ddd>,RC(19+i),nl,1);
write(res,<< -ddddddd.ddd>,<:b9 :>,RC(27),nl,1,<:hmmJ:>,RC(28),nl,1,
<:khi-aa:>,<<-ddddddd.ddd>,khi(1),nl,1,<:khi-bb:>,khi(2),
nl,1,<:khi-cc:>,khi(3),nl,1,<:deltaM:>,<<ddddd>,dM,nl,1);
rep:
read(in,m0,J1,k01,J2,k02,obs,w);
if m0<0 then goto stop;
J:= if J1>J2 then J1 else J2; k:= dM//3+2;
begin array u1(1:(J+J+1)*k+J+1);
A:= RC(1); B:= RC(2); C:= RC(3);
F:= RC(4); Z:=RC(5); X:= Z*Z; E:= (A-B)*(C-A)+9*X;
kb:= RC(20); kc:= RC(21); kab:= RC(22);
comment The standard constants are calculated in tau
using the Dowling relations and assuming tau-caca=tau-bcbc=0;
r4:= A/B; r8:= B/C; r3:= C/A;
r4:= -r4*r4; r8:= r8*r8; r3:= r3*r3;
r1:= 1/r4; r5:= 1/r8; r9:= 1/r3;
r2:= r1*r3; r6:= r4*r5; r7:= r8*r9;
for i:=1 step 1 until 3 do tau(i):= RC(i+6);
tau(4):= (RC(7)*r1 + RC(8)*r4 + RC(9)*r7)*0.5 + RC(10)*2;
tau(5):= (RC(7)*r2 + RC(8)*r5 + RC(9)*r8)*0.5 + RC(11)*2;
tau(6):= (RC(7)*r3 + RC(8)*r6 + RC(9)*r9)*0.5;
tau(8):= tau(4)+(RC(13)-RC(10))*2;
tau(9):= tau(6)+RC(14)*2;
tau(10):= tau(1)*3+RC(15); tau(11):= RC(16);
tau(12):= (RC(10)*2+RC(16)-4*RC(13))*r5+tau(1)*r3+RC(17);
tau(13):= RC(7)+RC(15)+RC(18);
comment Calculation of D-constants from standard constants,
and correction of rotational constants;
freq:= -1/16000;
D(1):= ((tau(2)+tau(3))*1.5+tau(5))*freq;
D(2):= -D(1)*2 + (tau(6)+tau(4))*4*freq;
D(3):= -D(1)-D(2)+tau(1)*4*freq;
D(4):= (tau(2)-tau(3))*freq;
D(5):= -D(4)-(tau(6)-tau(4))*2*freq;
D(6):= (tau(2)+tau(3)-tau(5)*2)*freq;
X:= (B+C)/2+D(6); Y:= (B-C)/4; Z:= A-X-D(6)*1.5; D(6):= D(6)/4;
freq:= 1/4000;
D(7):= RC(19)/4;
D(8):= (tau(11)+tau(12))*freq;
D(9):= (tau(10)*2-tau(11)-tau(12))*freq;
D(10):= (tau(11)-tau(12))*freq; freq:= freq*2;
D(11):= (tau(8)+tau(9))*freq;
D(12):= ((tau(1)+RC(12))*2-tau(8)-tau(9))*freq;
D(13):= (tau(8)-tau(9))*freq; D(14):= tau(13)/1000;
b1:= (kb+kc)/4; b2:= (kb-kc)/4; b3:= kab/4;
b4:= RC(23); b6:= RC(25); b7:= (r9+r4)*b4-(r9-r4)*b6; b5:= b7-b4;
b5:= b5+RC(24); b7:= b7+RC(26); b9:= RC(27);
freq:= 0; J:= J1; k0:= k01;
for i:=1 step 1 until nrc do dy(i):= 0;
for i:=1 step 1 until 3 do hf(i):= 0;
if closeres then write(out,<<dd>,m0,<:,:>);
for lev:=-1,1 do begin
Jeven:= J//2*2=J; JJ:=(J+1)*J;
Mmin:= m0-dM; Mmax:= m0+dM;
meven:= Mmin//2*2=Mmin;
if m0//3*3=m0 then begin
gamma:= abs k0 mod 10; k0:= k0//10;
plus:= gamma=0 == m0//2*2=m0;
if closeres then write(out,<<-ddd>,J,k0*10+sgn(k0)*gamma);
write(res,nl,2,<:m , J , K , gamma =:>,<<-dd>,m0,J,k0,gamma);
gamma:= if plus then 0 else 1; spinf:= if m0>6 then 22 else 11;
if Mmin<=0 then begin
Mmin:= 0; meven:= true; M1:= J; M2:= 1-gamma;
end else if meven then begin
M1:= J; M2:= (J+2)//2;
end else begin
M1:= J+1; M2:= M1//2;
end end else begin
if closeres then write(out,<<-ddd>,J,k0);
write(res,nl,2,<:m , J , K =:>,<<-dd>,m0,J,k0); gamma:= 0;
spinf:= if spind then 16 else 11;
if meven then begin
M1:= J; M2:= (J+2)//2;
end else begin
M1:= J+1; M2:= M1//2
end end;
w:= w shift 8 add m0 shift 8 add J shift 8 add
(k0+(if k0<0 then 256 else 0)) shift 8 add gamma shift 16;
iev:= -1; setposition(ev,0,0); inrec(ev,128);
for iev:=iev+2 while iev<=nev and ev(iev mod 128) shift 16<>w do
if iev mod 128 = 127 then inrec(ev,128);
if iev<=nev and m0//3*3<>m0 then begin
Mmin:= Mmin - ev(iev mod 128) shift (-32) extract 8;
if -,meven==Mmin//2*2=Mmin then begin
meven:= -,meven;
if meven then begin M1:= J; M2:= (J+2)//2
end else begin M1:= J+1; M2:= M1//2
end end end;
N:= ((J*2+1)*(Mmax-Mmin)//3+M1)//2+M2;
l:= if N<10 then 10 else N;
if system(2,i,tau) < (N+3)*N*4+160 then begin
write(res,<: stack:>,nl,1);
if closeres then write(out,nl,1); goto rep end;
begin
integer k2,m2;
real go,gu,p,q,y,lambda,XJ,YJ,ZJ,FJ,AfJ,V,Vk,psum;
array h(1:l,1:N), a,b,u(1:l), c(0:(Mmax-Mmin)//3);
integer procedure index(m,k);
value m,k; integer m,k;
begin integer kmax;
kmax:= if Jeven == m//2*2=m then J else J-1;
index:= if abs k>kmax or k<gamma and m=0 then 0 else
((J*2+1)*(m-Mmin)//3+M1)//2+(k-kmax)//2+M2;
end;
for i:=1 step 1 until N do
for j:=1 step 1 until i do h(i,j):= 0;
XJ:= (X-D(1)*JJ)*JJ; YJ:= Y-D(4)*JJ; ZJ:= Z-D(2)*JJ;
FJ:= F+D(8)*JJ; AfJ:= RC(5)+D(11)*JJ/2;
V:= RC(6)/4-b4*JJ/2; E:= -(D(7)*m0*m0+F)*m0*m0;
for m:=Mmin step 3 until Mmax do begin
Kmax:= if Jeven==meven then J else J-1;
Kmin:= if m<>0 then -Kmax else gamma*2;
q:= if m=0 and Kmin=2 then sqrt(J*(J-1)*(J+1)*(J+2)*2) else 0;
p:= if m=0 and Kmin=0 then sqrt(2) else 1;
for k:=Kmin step 2 until Kmax do begin
i:= index(m,k); k2:=k*k; m2:= m*m; Vk:= V-b5*k2/2-b7*(m+3)*k;
if i<1 or i>N then begin
write(res,<<ddd>,m,k); goto stop end;
h(i,i):= E+k2*ZJ+m2*FJ+k*m*2*AfJ-D(3)*k2*k2
+((RC(28)*'-6*JJ+D(7))*m2+D(9)*k2)*m2+(D(12)*k2+D(14)*m2)*k*m;
if m+6>Mmax then
h(i,i):= h(i,i)-V*V/((m*12+36)*FJ);
if m-6<Mmin and Mmin<0 then
h(i,i):= h(i,i)+V*V/((m*12-36)*FJ);
if k<Kmax then begin
y:= sqrt((J-k)*(J-k-1)*(J+k+1)*(J+k+2))*p;
h(i+1,i):= (YJ-(k2+k+k+2)*D(5)
+D(10)*m2/2+D(13)*(k+1)*m/2)*y;
if k>Kmin then h(i+1,i-1):= -q*y*D(6);
end;
if m=3 and k>0 then begin
l:= index(m,-k); go:= Vk+b7*k*6;
h(i,l):= h(i,l)+(if plus then go else -go);
l:= index(m,-k+2);
if l>0 then h(i,l):= h(i,l)+(if plus then -q else q)*(b6+b9)/4;
if k=1 then begin
l:= i-1; h(l,l):= h(l,l)+(if plus then -q else q)*(b6-b9)/4 end;
l:= index(m,-k-2);
if l>0 then h(i,l):= h(i,l)+(if plus then -y else y)*(b6-b9)/4;
end;
m2:= m+3;
if m2<=Mmax then begin
j:= index(m2,k-1);
if j>0 then begin
go:= sqrt((J+k)*(J-k+1));
h(j,i):= ((m+m2)*b2+(k+k-1)*b3)*go*p;
if m=0 and k>0 then
h(index(3,-k+1),i):= (if plus then go else -go)
*(b1*3-(k+k-1)*b3);
end;
j:= index(m2,k+1);
if j>0 then begin
go:= sqrt((J-k)*(J+k+1));
h(j,i):= ((m+m2)*b1+(k+k+1)*b3)*go*p;
if m=0 and k>0 then
h(index(3,-k-1),i):= (if plus then go else -go)
*(b2*3-(k+k+1)*b3);
end;
m2:= m+6;
if m2<=Mmax then begin
h(index(m2,k),i):= -Vk*p;
if m=0 and k>0 then begin
Vk:= Vk+b7*k*6;
h(index(m2,-k),i):= if plus then -Vk else Vk;
end;
j:= index(m2,k-2);
if j>0 then begin
if k=0 and m=0 then q:= y;
h(j,i):= q*(b6-b9)/4; j:= index(m2,-k+2);
if m=0 and k>0 then
h(j,i):= h(j,i)+(if plus then q else -q)*(b6+b9)/4;
if m=0 and k=2 then h(j,i):= h(j,i)/sqrt(2);
end;
j:= index(m2,k+2);
if j>0 then begin
h(j,i):= y*(b6+b9)/4; if m=0 and k>0 and k<Kmax then
h(index(m2,-k-2),i):= (if plus then y else -y)*(b6-b9)/4;
end end end;
p:= 1; q:= y;
end k;
if m=0 then begin
if -,plus then h(1,1):= h(1,1)+D(6)*(JJ-2)*JJ else
if Kmax>0 then h(2,2):= h(2,2)-D(6)*(JJ-2)*JJ
end;
meven:= -,meven;
end m;
if m0=0 then begin
m:= J//2+1-gamma; k:= 0; pm:= (J-k0)//2+1; l:= m+1-pm;
end else begin
i:= m0+(if k0<-3*J-1 then -6 else if k0<-J then -3 else 0);
j:= if i=m0 then k0 else if i=m0-3 then k0+2*J+1 else k0+4*J+2;
m:= if Jeven==i//2*2=i then J+1 else J;
k:= index(i,-m+1)-1; pm:= (m+1-j)//2; l:= index(i,j);
end;
if iev<=nev then begin
lambda:= ev(iev mod 128 + 1);
i:= ev(iev mod 128) shift (-40) extract 8;
l:= l+(if i<128 then i else i-256);
end else begin
for i:=k+2 step 1 until k+m do begin
a(i-k):= h(i,i); b(i-k):= h(i,i-1)**2 end;
a(1):= h(k+1,k+1); b(1):= 0;
go:= a(m)*2; gu:= -go;
for lambda:=(go+gu)/2 while go-lambda>10 do begin
p:= 0; q:= 1; it:= 0;
for i:=1 step 1 until m do begin
y:= (a(i)-lambda)*q-b(i)*p;
p:= q; q:= y;
if p>=0 == q>=0 then it:= it+1;
if it=pm then i:= m;
end i;
if it=pm then gu:= lambda else go:= lambda;
end bisection;
write(res,nl,1,<<ddd>,iev,<< -dddddddd.ddd>,lambda);
lambda:= lambda+(if m0=3 then -V/(72*FJ)+V/(12*AfJ*k0)
else V/(FJ*2*(m0*m0-9)))*V;
end;
Vk:= h(l,l);
for i:=l-1 step -1 until 1 do b(i):= h(l,i);
for i:=l+1 step 1 until N do begin
b(i-1):= h(i,l);
for j:=l-1 step -1 until 1 do h(i-1,j):= h(i,j);
for j:=l+1 step 1 until i do h(i-1,j-1):= h(i,j)
end;
it:= 0; N:=N-1; y:= lambda;
repeat1:
for i:=N step -1 until 1 do begin
for j:=i step -1 until 1 do begin
p:= h(i,j);
for k:=i+1 step 1 until N do
p:= p-h(i,k)*h(j,k)*a(k);
if i=j then a(i):= p-lambda
else h(j,i):= p/a(i);
end j;
p:= -b(i);
for k:=i+1 step 1 until N do
p:= p-h(i,k)*u(k);
u(i):= p;
end i;
q:= 1;
for i:=1 step 1 until N do begin
p:= if u(i)=0 then 0 else u(i)/a(i);
for k:=i-1 step -1 until 1 do
p:= p-h(k,i)*u(k);
u(i):= p; q:= q+p*p;
end i;
go:= Vk-lambda;
for i:=1 step 1 until N do
go:= go + b(i)*u(i);
go:= go/q;
if abs go>'-3 and lambda+go<>lambda
and (it<2 or abs go < abs y) then begin
write(res,nl,1,<<ddd>,it,<< -dddddddd.ddd>,lambda);
lambda:= lambda+go; it:= it+1; y:= go; goto repeat1;
end;
q:= 1/sqrt(q);
for i:=N step -1 until l do u(i+1):= u(i)*q;
for i:=l-1 step -1 until 1 do u(i):= u(i)*q;
u(l):= q; N:= N+1;
q:= 0; j:= 0;
for i:=1 step 1 until N do
if abs u(i)>q then begin q:= abs u(i); j:= i end;
write(res,<< ddd>,l,j); forceout(res);
getposition(ev,0,i); setposition(ev,0,i); outrec(ev,128);
ev(iev mod 128 + 1):= lambda; j:= l:= (j-l) shift 8;
if j<>0 or iev>nev then begin
if iev>nev then nev:= iev;
if j<>0 and iev<nev then
l:= l+ev(iev mod 128) shift (-32) extract 16;
m:= w shift (-16) extract 24;
i:= w shift(-40) extract 8;
ev(iev mod 128):= 0.0 shift 6 add l shift 8 add i shift 24 add m;
end;
lambda:= lambda+XJ;
if lev=1 then begin
LS:= 0; j:= if J1>J2 then J2 else J1; i:= l:= 1;
k:= Mmin1-Mmin;
if k<0 then begin k:= -k;
i:= i+k//6*(2*J1+1); Mmin1:= Mmin;
if k mod 6=3 then i:= i+N1
end else if k>0 then begin
l:= l+k//6*(2*J+1);
if k mod 6=3 then l:= l+M1//2+M2
end;
meven:= Mmin1//2*2=Mmin1; Jeven:= j//2*2=j;
for m:= Mmin1 step 3 until Mmax do begin
Kmax:= if Jeven==meven then j else j-1;
Kmin:= if m<>0 then -Kmax else if J1=J2 then 2 else
if gamma=0 then 0 else 2;
if m<>0 then begin
if J1-Kmax=2 then i:= i+1;
if J2-Kmax=2 then l:= l+1
end else if J1=J2 then begin
if gamma=0 then l:=2 else i:= 2 end;
for k:=Kmin step 2 until Kmax do begin
LS:= LS+u1(i)*u(l)*(if J1=J2 then k else sqrt((j+1+k)*(j+1-k)));
i:=i+1; l:= l+1 end;
if J1-Kmax=2 then i:= i+1;
if J2-Kmax=2 then l:= l+1;
meven:= -,meven;
end;
LS:= (if J1=J2 then LS*(J+J+1)/JJ else LS/(j+1))*LS
*exp((E+freq)*4.7993'-5/T);
Jeven:= J//2*2=J
end else for i:=1 step 1 until N do u1(i):= u(i);
freq:= freq+lambda*lev;
for i:=1 step 1 until N do
for j:=4 step 1 until 10 do h(j,i):= 0;
for j:=1 step 1 until 10 do a(j):= 0;
for j:=1 step 1 until 14 do tau(j):= 0;
it:= (Mmax-Mmin)//3; for i:=0 step 1 until it do c(i):= 0; it:= -1;
meven:= Mmin//2*2=Mmin; psum:= 0;
for m:=Mmin step 3 until Mmax do begin
Kmax:= if Jeven==meven then J else J-1;
k:= if m<>0 then -Kmax else gamma*2;
gu:= if m=0 and k=2 then sqrt(J*(J-1)*(J+1)*(J+2))/4 else 0;
p:= if m=0 and k=0 then sqrt(2) else 1; it:= it+1;
for k:=k step 2 until Kmax do begin
i:= index(m,k); y:= u(i); k2:=k*k; m2:= m*m; c(it):= c(it)+y*y;
h(1,i):= y*m2; h(2,i):= y*m*k; h(3,i):= y*k2; h(8,i):= y*m2*m2;
if k<Kmax then begin
go:= sqrt((J-k)*(J-k-1)*(J+k+1)*(J+k+2))*p/2;
h(4,i+1):= go*y; h(4,i):= h(4,i)+go*u(i+1); go:= go/2
end;
if m+6>Mmax then psum:= psum-y*y*V/((12*m+36)*FJ);
if m-6<Mmin and Mmin<0 then psum:= psum+y*y*V/((12*m-36)*FJ);
if m=3 and k>0 then begin
l:= index(m,-k); q:= if plus then -0.5 else 0.5;
h(5,l):= h(5,l)+q*y;
h(5,i):= h(5,i)+q*u(l);
l:= index(m,-k+2);
if l>0 then begin q:= if plus then -gu else gu;
h(9,l):= h(9,l)+q*y; if k>1 then
h(9,i):= h(9,i)+q*u(l);
if k=1 then begin l:= i-1; h(9,l):= h(9,l)-q*u(l) end;
end;
l:= index(m,-k-2);
if l>0 then begin q:= if plus then go else -go;
h(9,l):= h(9,l)+q*y;
h(9,i):= h(9,i)+q*u(l);
end end;
m2:= m+3;
if m2<=Mmax then begin
j:= index(m2,k-1);
if j>0 then begin
w:= sqrt((J+k)*(J-k+1))*p; q:= (m+m2)*w; w:= (k+k-1)*w;
h(7,j):= h(7,j)+q*y; h(7,i):= h(7,i)+q*u(j);
h(10,j):= h(10,j)+w*y; h(10,i):= h(10,i)+w*u(j);
if m=0 and k>0 then begin
j:= index(m2,-k+1);
if -,plus then begin q:= -q; w:= -w end;
h(6,j):= h(6,j)+q*y; h(6,i):= h(6,i)+q*u(j);
h(10,j):= h(10,j)-w*y; h(10,i):= h(10,i)-w*u(j);
end end;
j:= index(m2,k+1);
if j>0 then begin
w:= sqrt((J-k)*(J+k+1))*p; q:= (m+m2)*w; w:= (k+k+1)*w;
h(6,j):= h(6,j)+q*y; h(6,i):= h(6,i)+q*u(j);
h(10,j):= h(10,j)+w*y; h(10,i):= h(10,i)+w*u(j);
if m=0 and k>0 then begin
j:= index(m2,-k-1);
if -,plus then begin q:= -q; w:= -w end;
h(7,j):= h(7,j)+q*y; h(7,i):= h(7,i)+q*u(j);
h(10,j):= h(10,j)-w*y; h(10,i):= h(10,i)-w*u(j);
end end;
m2:= m+6;
if m2<=Mmax then begin
j:= index(m2,k); q:= p/2;
h(5,j):= h(5,j)+q*y;
h(5,i):= h(5,i)+q*u(j);
if m=0 and k>0 then begin
j:= index(m2,-k); q:= if plus then 0.5 else -0.5;
h(5,j):= h(5,j)+q*y;
h(5,i):= h(5,i)+q*u(j);
end;
j:= index(m2,k-2);
if j>0 then begin if k=0 and m=0 then gu:= go;
h(9,j):= h(9,j)-gu*y;
h(9,i):= h(9,i)-gu*u(j);
if m=0 and k>0 then begin
j:= index(m2,-k+2); q:= if plus then gu else -gu;
h(9,j):= h(9,j)+q*y;
h(9,i):= h(9,i)+q*u(j);
end end;
j:= index(m2,k+2);
if j>0 then begin
h(9,j):= h(9,j)+go*y;
h(9,i):= h(9,i)+go*u(j);
if m=0 and k>0 then begin
j:= index(m2,-k-2); q:= if plus then -go else go;
h(9,j):= h(9,j)+q*y;
h(9,i):= h(9,i)+q*u(j);
end end end end;
p:= 1; gu:= go;
end k; meven:= -,meven
end m;
for i:=1 step 1 until N do begin
for j:=1 step 1 until 10 do a(j):= a(j)+h(j,i)*u(i);
tau(1) := tau(1)+h(3,i)**2; tau(2) := tau(2)+h(3,i)*h(4,i);
tau(3) := tau(3)+h(4,i)**2; tau(4) := tau(4)+h(8,i)*h(4,i);
tau(5) := tau(5)+h(1,i)*h(3,i); tau(6) := tau(6)+h(1,i)*h(4,i);
tau(7) := tau(7)+h(2,i)*h(3,i); tau(8) := tau(8)+h(2,i)*h(4,i);
tau(9) := tau(9)+h(1,i)*h(2,i); tau(10):= tau(10)+h(3,i)*h(5,i);
tau(11):= tau(11)+h(4,i)*h(5,i);tau(12):= tau(12)+h(2,i)*h(5,i);
end;
p:= (JJ-a(3))/2; q:= a(4)/2;
ZJ:= a(3); YJ:= p+q; XJ:= p-q;
go:= (khi(1)*ZJ+khi(2)*YJ+khi(3)*XJ)/2;
p:= if J=0 then 0 else go/(J*(2*J-1));
q:= if J=0 then 0 else -go/JJ;
y:= go/((J+1)*(2*J+3));
hf(1):= hf(1)+p*lev; hf(2):= hf(2)+q*lev; hf(3):= hf(3)+y*lev;
write(res,sp,3,<<-dddddddd.ddd>,lambda,nl,2);
comment eigen vector:
for i:=1 step 1 until N do
write(res,<< -d.dddddd>,u(i),nl,if i mod 5 = 0 then 1 else 0);
for i:=0 step 1 until it do
write(res,nl,if i mod 5=0 then 1 else 0,<< d.dddddd>,c(i));
b(1):= (JJ*a(1)-tau(5)+tau(6))/4;
b(2):= (JJ*a(1)-tau(5)-tau(6))/4;
b(3):= (JJ*a(2)-tau(7)+tau(8))/2;
b(4):= (JJ*a(2)-tau(7)-tau(8))/2;
p:= (JJ-a(3)*2)*JJ+tau(1); q:= (JJ*a(4)-tau(2))*2;
u(1):= tau(1)/4+b(2)*r3+tau(5)*1.5+tau(7)+tau(9);
u(2):= (p+q+tau(3))/16;
u(3):= (p-q+tau(3))/16;
u(5):= (p-tau(3))/8; p:= JJ*a(3)-tau(1);
u(6):= (p-tau(2))/4+b(4);
u(4):= (p+tau(2))/4+b(3);
go:= (JJ-a(3))/2; q:= a(4)/2; p:= go+q; q:= go-q;
dy(1):= dy(1)+a(3)*lev; dy(2):= dy(2)+p*lev;
dy(3):= dy(3)+q*lev;
gu:= tau(10)+2*tau(12);
dy(4):= dy(4)+(a(1)-psum*V/RC(4))*lev;
dy(5):= dy(5)+a(2)*2*lev; dy(6):= dy(6)+(psum/2-a(5)/2)*lev;
dy(7):= dy(7)+(u(1)+(u(4)*r1+u(5)*r2+u(6)*r3)*0.5)*lev;
dy(8):= dy(8)+(u(2)+(u(4)*r4+u(5)*r5+u(6)*r6)*0.5)*lev;
dy(9):= dy(9)+(u(3)+(u(4)*r7+u(5)*r8+u(6)*r9)*0.5)*lev;
dy(10):= dy(10)+(u(4)+b(2)*r5-b(3))*2*lev;
dy(11):= dy(11)+u(5)*2*lev; dy(12):= dy(12)+tau(7)*lev;
dy(13):= dy(13)+(b(3)-b(2)*r5*2)*2*lev;
dy(14):= dy(14)+b(4)*2*lev;
dy(15):= dy(15)+(tau(5)/2+tau(9))*lev;
dy(16):= dy(16)+(b(1)+b(2)*r5)*lev;
dy(17):= dy(17)+b(2)*lev; dy(18):= dy(18)+tau(9)*lev;
dy(19):= dy(19)+a(8)*lev/4;
dy(20):= dy(20)+(a(6)+a(7))/4*lev;
dy(21):= dy(21)+(a(6)-a(7))/4*lev;
dy(22):= dy(22)+a(10)/4*lev;
dy(23):= dy(23)+(a(5)*JJ+(r9+r4)*gu-tau(10))*lev;
dy(24):= dy(24)+tau(10)*lev;
dy(25):= dy(25)+(tau(11)-(r9-r4)*gu)*lev;
dy(26):= dy(26)+tau(12)*2*lev;
dy(27):= dy(27)+a(9)*lev;
dy(28):= dy(28)+JJ*a(8)*'-6*lev;
J:= J2; k0:= k02; Mmin1:= Mmin; N1:= M1//2+M2;
end end;
for i:=7 step 1 until 18 do dy(i):= dy(i)/1000;
A:= sum(dy(j)*RC(j),j,1,nrc);
if disc then begin array x(1:nfit);
for i:=1 step 1 until nfit do
x(i):= sum(error(i,j)*dy(no(j)),j,1,nfit);
B:= sum(x(i)*dy(no(i)),i,1,nfit);
B:= if B>=0 then sqrt(B) else -sqrt(-B);
end else B:= 0;
write(res,nl,2);
for j:=1 step 1 until nrc do begin
write(res,<< -dddddddd.ddd>,dy(j)*RC(j));
if j mod 5=0 then write(res,nl,1)
end;
LS:= LS*freq*freq*spinf*'-9;
write(res,nl,2,<< -dddddd.ddd>,<:Transition-freq::>,freq,
hf(1),hf(2),hf(3),nl,1,sp,16,
B,if obs=0 then 0 else obs-freq,A-freq,LS,nl,3);
forceout(res);
if closeres then begin write(out,<< -ddddddd.dd>,freq,
<<-dddd.dd>,B,if obs=0 then 0 else obs-freq,LS,
dy(4)*RC(4)+dy(6)*RC(6)+sum(dy(i)*RC(i),i,20,22),
dy(28)*RC(28),nl,1);
forceout(out);
end; goto rep end;
stop:
if tail(10)<nev then begin
tail(10):= nev; monitor(44,ev,0,tail) end;
close(ev,true); write(res,<:<25>:>); close(res,closeres) end end;
▶EOF◀