|
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: 26880 (0x6900)
Types: TextFile
Names: »algrotfit0«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦7e928b248⟧ »algbib«
└─⟦this⟧
;gosav lines.10000
lookup rotfit
if ok.no
(rotfit=set 95
permanent rotfit.17
rotfit=algol list.no index.no xref.no)
\f
FITTING OF ROTATIONAL CONSTANTS. 3-3-1975. GOS.
begin
integer i,j,k,l,m,m1,N,N1,J1,J2,Q,Q1,linecount,page,space,version;
real s, so, starttime, stoptime, cpu, time;
boolean slut,aires,delres,iterer,oblat,
sp,closeres,quartic,sextic,hb,ID;
array Rot, dRot, Int, dInt(1:4), D, head(1:12);
zone res(128,1,stderror);
procedure writecr;
begin linecount:= linecount+1;
if linecount>60 then writepage;
write(res,<:<10>:>)
end writecr;
procedure writepage;
begin integer i;
page:= page+1; linecount:= 3; i:= 1;
write(res,<:<10><12>:>,false add 10,3,string head(increase(i)),
false add 32, space,<<dd>, page, false add 10, 2)
end writepage;
closeres:=outmedium(res); sp:= false add 32;
dataind:
cpu:= systime(1,0,starttime);
space:= 68 - readhead(in,head,1); read(in,m1,stoptime,Rot);
stoptime:= stoptime*60; N:= Rot(4);
slut:= m1>0; m1:= abs m1; m:= m1//1000;
delres:= m1 mod 1000 >= 100; version:= m1 mod 10;
aires := m1 mod 100 >= 10; iterer := N<>0;
ID:= version > 4;
page:= 0; writepage; write(res,<:
Calculation by ROTFIT of 3-3-1975
:>); if ID then begin
Rot(1):= Rot(1)/505376; version:= version-4;
write(res,<:with inertial defect as first parameter.:>);
writecr end;
case version of begin
begin Q1:= 8; Q:= 15;
linecount:= linecount+6; write(res,<:
D1 = Dj, D2 = Djk, D3 = Dk, D4 = dj, D5 = djk, (kHz)
D6 = Hj, D7 = Hjk, D8 = Hkj, D9 = Hk,
D10= hj, D11= hjk, D12= hk, (Hz).
:>) end;
begin Q1:= 8; Q:= 15;
linecount:= linecount+8; write(res,<:
D1 = tau-aaaa, D2 = tau-bbbb, D3 = tau-cccc,
D4 = tau-abab + tau-caca*(A-B)/(A-C),
D5 = tau-bcbc + tau-caca*(B-C)/(A-C),
D6 = Hj, D7 = Hjk, D8 = Hkj, D9 = Hk
D10= hj, D11= hjk, D12= hk, (Hz).
:>) end;
begin Q1:= 7; Q:= 8;
linecount:= linecount+7; write(res,<:
D1 = tau-abab + 2*tau-cccc,
D2 = tau-aaaa + tau-bbbb - 2*tau-abab - 8*tau-cccc,
D3 = tau-aaaa - tau-bbbb, D4 = tau-cccc (kHz).
D5 = Hjk (Hz).
:>) end;
begin Q1:= Q:= 8;
linecount:= linecount+5; write(res,<:
Version for analysis of transitions between K-doublet levels
(see J.Mol.Structure 8 (1971) 217). Constants in kHz.
:>) end;
end;
s:= if ID then 1/(1/Rot(3)-1/Rot(2)-Rot(1)) else Rot(1);
oblat:= 2*Rot(2)>(s+Rot(3));
N1:= if iterer then N else 1;
l:= Q-m;
begin
integer k11,k12,k21,k22,type1,type2,min,xx,yy,zz,xy,yz,zx;
array DR(3:Q), tr(1:N1), an(1:l,1:l+1),
col(1:Q,1:Q+1), AO,lo(1:Q);
integer array nf(0:m),qn(1:N1,1:3),MATRIX(1:8),lw,cif,exp(1:Q);
min:= 1;
for i:=4,i+1 while k<>10 do begin
if i>Q then begin write(res,<:
***delimiter error in C.D.constant input.:>); goto stop end;
read(in,DR(i)); repeatchar(in); readchar(in,k) end;
for j:=i step 1 until Q do DR(j):= 0;
for j:=1 step 1 until 12 do D(j) := 0;
for j:=1 step 1 until m do read(in,nf(j));
for i:=1 step 1 until Q do cif(i):= 0;
writecr;
if m=1 then write(res,<:Fixed Constant::>) else
if m>1 then write(res,<:Fixed Constants::>);
for i:=1 step 1 until m do begin
if i-1 = i//2*2 and i<>1 then begin writecr;
write(res,sp, 16) end;
l:= nf(i); AO(l):= 0; write(res,<: :>);
if l<4 then begin
so:= Rot(l);
if ID and l=1 then begin
so:= so*505376; write(res,<:ID = :>) end else
write(res, case l of (<:A = :>,<:B = :>,<:C = :>));
write(res,<<dd ddd.dddd>, so)
end else begin
so:= DR(l);
write(res,<:D:>,<<d>,l-3,<: = :>, if abs so=0 then <<d>
else <<-ddd.ddd0000>,so)
end
end;
if m<>0 then writecr;
if iterer then begin
integer q,max,my;
boolean b11,b12,b21,b22;
max:=1; writecr;writecr;write(res,<:Input transitions::>);writecr;
for k:=1 step 1 until N do begin
read(in, J1, k11, k12, J2, k21, k22, my, tr(k));
writecr; write(res,<<ddd>,J1,k11,k12,J2,k21,k22,my,
<<___ddd ddd.ddd>,tr(k));
repeatchar(in); for i:=readchar(in,l) while l<>10 and l<>9 do;
b11:=k11=k11//2*2; b12:=k12=k12//2*2;
b21:=k21=k21//2*2; b22:=k22=k22//2*2;
type1:=if b11 and b12 then 1 else
if b11 and-,b12 then 2 else
if -,b11 and b12 then 3 else 4;
type2:=if b21 and b22 then 1 else
if b21 and-,b22 then 2 else
if -,b21 and b22 then 3 else 4;
b11:= b11==b21; b22:= b12==b22;
if my<>(if b11 and-,b22 then 1 else
if b22 and-,b11 then 3 else 2)or abs(J1-J2)>1 then begin
N1:=N1-1; write(res,<: ***deleted:>);
if k=min then max:=min:=min+1;
goto linieslut end;
comment The reading of a transition have been performed, and
the quantum numbers have been controlled for correct tran-
sition type, my, and will be packed in the array qn. To
obtain print out arranged in order of increasing J, a number
giving the index of the line to follow is stored in qn(k,3);
qn(k,1):= J1 shift 7 add J2 shift 7 add k11 shift 3 add type1;
qn(k,2):= k12 shift 7 add k21 shift 7 add k22 shift 3 add type2;
qn(k,3):= min; q:=max;
for i:=min,qn(i,3) while i<>min do
if qn(i,1)<=qn(k,1) then q:=i else i:= max;
if q=max then begin
if qn(k,1)>=qn(max,1) then max:=k else
if qn(k,1)<qn(min,1) then min:=k end;
qn(k,3):= qn(q,3); qn(q,3):= k;
linieslut: end k;
writecr; writecr; write(res,<:Total number: :>,<<ddd>,N1);
writecr
end reading and control of transitions;
if oblat then begin
for i:=1 step 1 until 8 do
MATRIX(i):= case i of (1,4,2,3,2,3,1,4) end
else begin
for i:=1 step 1 until 8 do
MATRIX(i):= case i of (1,2,3,4,2,1,4,3) end;
systime(1,starttime,time);
l:= Q-m; nf(0):= 0; s:= 0; so:= 2**64;
if oblat then begin
xx:= 1; yy:= 2; zz:= 3; xy:= 4; yz:= 5; zx:= 6
end else begin
zz:= 1; xx:= 2; yy:= 3; zx:= 4; xy:= 5; yz:= 6 end;
igen:
m1:= if -,iterer then 0 else if N1=l then 1 else N1-l;
quartic:= sextic:= false;
for i:= 4 step 1 until Q1 do quartic:= quartic or DR(i)<>0;
for i:= Q1+1 step 1 until Q do sextic:= sextic or DR(i)<>0;
begin
real p,qd,y,a,b,c,a2,b2,c2,X,Y,Z,r1,r2,r3,r4,r5,r6,r7,r8,r9,DW;
integer j0,k0,K2;
boolean ha;
array ai(1:Q+1);
b:= Rot(2); c:= Rot(3); Int(2):= 1/b; Int(3):= 1/c;
if ID then begin
Int(4):= Rot(1); Int(1):= Int(3)-Int(4)-Int(2);
a:= Rot(1):= 1/Int(1);
a2:= a*a; b2:= a2/(b*b); c2:= a2/(c*c)
end else begin
a:= Rot(1); Int(1):= 1/a;
Int(4):= Int(3)-Int(1)-Int(2)
end;
p:= 1/(a-c); Rot(4):= (b+b-a-c)*p;
ha:= s=so; hb:= false;
X:= Rot(xx); Y:= Rot(yy); Z:= Rot(zz);
if ha then begin
ai(1):= -(1+Rot(4))*p; ai(2):= p+p;
ai(3):= -(1-Rot(4))*p; p:= 0;
if ID then begin
ai(3):= ai(3)+ai(1)*c2;
ai(2):= ai(2)-ai(1)*b2; ai(1):= ai(1)*a2 end;
for i:=1 step 1 until 3 do begin
dRot(i):= y:= sqrt(col(i,i)*s); dInt(i):= y/(Rot(i)**2);
for j:=i step 1 until 3 do p:= p + col(i,j)*ai(i)*ai(j)
end;
dRot(4):= sqrt(p*s); p:= 0;
if ID then begin
ai(1):= a2; ai(2):= -b2; ai(3):= c2
end else begin
ai(1):= 1/a**2; ai(2):= 1/b**2; ai(3):= -1/c**2 end;
for i:=1 step 1 until 3 do
for j:=i step 1 until 3 do p:= p + col(i,j)*ai(i)*ai(j);
p:= sqrt(p*s);
if ID then begin
dInt(4):= dRot(1); dRot(1):= p; dInt(1):= p/a2
end else dInt(4):= p
end ha else if iterer then writepage;
case version of begin
begin
for i:=1 step 1 until 5 do D(i):= DR(i+3)*'-3;
for i:=6 step 1 until 12 do D(i):= DR(i+3)*'-6;
end vers1;
begin
comment The standard constants are calculated in ai (in kHz)
using the Dowling relations and assuming tau-caca=0;
r4:= a*Int(2); r8:= b*Int(3); r3:= c*Int(1);
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 j:=1 step 1 until 3 do ai(j):= DR(j+3)*'3;
ai(4):= (DR(4)*r1 + DR(5)*r4 + DR(6)*r7)*0.5 + DR(7)*2;
ai(5):= (DR(4)*r2 + DR(5)*r5 + DR(6)*r8)*0.5 + DR(8)*2;
ai(6):= (DR(4)*r3 + DR(5)*r6 + DR(6)*r9)*0.5;
for j:=4 step 1 until 6 do ai(j):= ai(j)*'3;
for j:=6 step 1 until 12 do D(j):= DR(j+3)*'-6;
hb:= true end vers2;
begin
comment As in version 2 the standard constants are calcu-
lated, but in this case it is assumed that both constants
tau-bcbc = 0 and tau-caca = 0;
r4:= a*Int(2); r8:= b*Int(3); r3:= c*Int(1);
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;
p:= r1; r1:= (p+r4)*0.5; r4:= (p-r4)*0.25; r7:= 2*r1+r7*0.5-4;
p:= r2; r2:= (p+r5)*0.5; r5:= (p-r5)*0.25; r8:= 2*r2+r8*0.5;
p:= r3; r3:= (p+r6)*0.5; r6:= (p-r6)*0.25; r9:= 2*r3+r9*0.5;
ai(1):= DR(4) + (DR(5)+DR(6))*0.5 + DR(7)*2;
ai(2):= DR(4) + (DR(5)-DR(6))*0.5 + DR(7)*2;
ai(3):= DR(7);
ai(4):= (DR(4)+DR(5)*0.5)*r1 + DR(6)*r4 + DR(7)*r7 + DR(4)*2;
ai(5):= (DR(4)+DR(5)*0.5)*r2 + DR(6)*r5 + DR(7)*r8;
ai(6):= (DR(4)+DR(5)*0.5)*r3 + DR(6)*r6 + DR(7)*r9;
D(7):= DR(8)*'-6;
hb:= true end vers3;
begin
c:= 1/(X-Y); a:= (Y-Z)*c; b:= (Z-X)*c;
r1:= (X-Y)/(Z+Z-X-Y);
a:= DR(4)+DR(6)*r1*r1*2;
b:= DR(5)-DR(7)*r1;
ai(xx):= (a+b+DR(8))*4; ai(yz):= (-b-DR(7))*2+DR(8)*4;
ai(yy):= (a-b+DR(8))*4; ai(zx):= ( b-DR(7))*2+DR(8)*4;
ai(zz):= (DR(6)-DR(7)+DR(8))*4; ai(xy):= (-a+DR(8))*4;
hb:= true
end vers4;
end;
if hb then begin
comment Calculation of D-constants from standard constants,
and correction of rotational constants;
c:= -1/16000;
D(1):= ((ai(xx)+ai(yy))*1.5+ai(xy))*c;
D(2):= -D(1)*2 + (ai(yz)+ai(zx))*4*c;
D(3):= -D(1)-D(2)+ai(zz)*4*c;
D(4):= (ai(xx)-ai(yy))*c;
D(5):= -D(4)-(ai(yz)-ai(zx))*2*c;
DW := (ai(xx)+ai(yy)-ai(xy)*2)*c;
X:= X+DW; Y:= Y+DW; Z:= Z-DW*1.5; DW:= DW/4
end else begin
comment Calculation of standard constants from D-constants
assuming tau-xxxx + tau-yyyy - 2*tau'-xxyy = 0;
p:= -4000; DW:= 0;
ai(xx):= (D(1) + D(4)*2)*p;
ai(yy):= (D(1) - D(4)*2)*p;
ai(zz):= (D(1) + D(2) + D(3))*p;
ai(xy):= D(1)*p;
ai(yz):= (D(1) + D(2)*0.5 - D(4) - D(5))*p;
ai(zx):= (D(1) + D(2)*0.5 + D(4) + D(5))*p
end;
j:= 6; p:= 505376;
for i:=1 step 1 until Q do
lw(i):= if i>3 then 14 else 16;
exp(1):= k0:= if ID and ha then
entier(ln(abs Rot(1))/ln10)-entier(ln(dRot(1)/2)/ln10)+1
else cif(1);
for i:=2,3 do exp(i):= cif(i);
exp(4):= if ha then
entier(ln(abs Rot(4))/ln10)-entier(ln(dRot(4)/2)/ln10)+1
else 7;
write(res,<:
Rotational constants (MHz) and kappa:
:>); for i:=1 step 1 until 4 do begin y:= Rot(i);
write(res,string format(lw(i),exp(i),y,lo(i)),y) end;
if ha then begin write(res,<:
+-:>); for i:=1 step 1 until 4 do begin y:= dRot(i);
if y=0 then write(res,sp,lw(i)-5,<:fixed:>)
else write(res,string lo(i),y) end end;
exp(1):= k0;
exp(4):= if ID then cif(1) else if ha then
entier(ln(abs Int(4))/ln10)-entier(ln(dInt(4)/2)/ln10)+1
else 4;
for i:=2,3 do exp(i):= cif(i);
write(res,<:
Inertial constants and defect (uÅ**2):
:>); for i:=1 step 1 until 4 do begin y:= Int(i)*p;
write(res,string format(lw(i),exp(i),y,lo(i)),y) end;
if ha then begin write(res,<:
+-:>); for i:=1 step 1 until 4 do begin y:= dInt(i)*p;
if y=0 then write(res,sp,lw(i)-5,<:fixed:>)
else write(res,string lo(i),y) end end;
for i:=4 step 1 until Q do exp(i):= cif(i);
if quartic then begin write(res,<:
Quartic distortion constants:
:>); for i:=4 step 1 until Q1 do begin y:= DR(i);
write(res,string format(lw(i),exp(i),y,lo(i)),y) end;
j:= j+5;
if ha then begin write(res,<:
+-:>); for i:=4 step 1 until Q1 do begin y:= sqrt(col(i,i)*s);
if y=0 then write(res,sp,lw(i)-5,<:fixed:>)
else write(res,string lo(i),y) end end;
write(res,<:
Standard distortion constants (kHz):
:>); for i:=1 step 1 until 6 do
write(res,<<-dddd.ddd000>,ai(i)) end;
if sextic then begin write(res,<:
Sextic distortion constants:
:>); k0:= if Q>Q1+4 then Q1+4 else Q;
for i:=Q1+1 step 1 until k0 do begin y:= DR(i);
write(res,string format(lw(i),exp(i),y,lo(i)),y) end; j:= j+3;
if ha then begin write(res,<:
+-:>); for i:=Q1+1 step 1 until k0 do begin y:= sqrt(col(i,i)*s);
if y=0 then write(res,sp,lw(i)-5,<:fixed:>)
else write(res,string lo(i),y) end end;
ha:= false; for i:=k0+1 step 1 until Q do ha:= ha or DR(i)<>0;
if ha then begin writecr; writecr;
write(res,<: :>); ha:= s=so;
for i:=k0+1 step 1 until Q do begin y:= DR(i);
write(res,string format(lw(i),exp(i),y,lo(i)),y) end;
if ha then begin writecr; write(res,<:+-:>);
for i:=k0+1 step 1 until Q do begin y:= sqrt(col(i,i)*s);
if y=0 then write(res,sp,lw(i)-5,<:fixed:>)
else write(res,string lo(i),y) end
end end else ha:= s=so end;
if ha then j:= j//3 + j;
linecount:= linecount+j+4;
if ha and -,iterer then begin
writecr; writecr; writecr;
write(res,<:Significant Digits and Correlation Matrix::>);
writecr; linecount:=linecount - 4;
for i:=1 step 1 until l do begin p:= an(i,i);
AO(i):= if p>0 then 1/sqrt(p) else 0 end;
N:= (l-1) mod 6 + 1; j0:= 0; k:= 1;
Cor: for i:=1 step 1 until N do begin writecr;
for j:=k step 1 until i-1 do write(res,sp,8);
if j=i then begin j:= j+1;
for j0:=j0+1 while cif(j0)= 0 do;
write(res,<<__dd>,cif(j0),<: :>) end;
for j:=j step 1 until N do
write(res,<<__-d.ddd>,an(i,j)*AO(i)*AO(j))
end i;
writecr; k:= N+1; N:= N+6; if k<l then goto Cor;
linecount:= linecount+4 end;
write(res,false add 10,3);
if iterer then begin
write(res,<: Transition obs. frequency obs.- t:>);
if quartic then write(res,<: distortion corrections:>);
write(res,<:<10>:>,sp,33,<:calc. :>)
end else begin
read(in,N1); if N1<0 then goto stop;
if linecount > 50 then writepage;
write(res,<: Transition calc. frequency sigma inten-:>);
if quartic then write(res,<: distortion corrections:>);
write(res,<:<10>:>,sp,38,<: sitet:>)
end;
if quartic then begin
write(res,<: total:>);
if sextic then write(res,<: sextic:>);
if m>0 then write(res,<: fixed:>)
end; writecr;
p:= Y; Y:= (X-p)*0.25; X:= (X+p)*0.5; Z:= Z-X;
if iterer then s:= 0; Q:= Q+1;
begin
integer line, lc, matr, nr, type, JJ, it;
real norm,lambda,Li,dE,F,G,G1,G2,H,H1,H2,
Pxy,Pz,Pzxy,Pxyxy,Pzzd,Pz4,Pz6;
array atot(1:if iterer then N1 else 1,1:l+1), imp(1:2,1:5);
lc:= 0;
for line:= min, qn(line,3) while line<>min do begin
if iterer then begin
N:= qn(line,1); lc:= lc+1;
J1 := N shift (-17) extract 7;
J2 := N shift (-10) extract 7;
k11:= N shift (-3) extract 7;
type1:= N extract 3;
N:= qn(line,2);
k12:= N shift (-17) extract 7;
k21:= N shift (-10) extract 7;
k22:= N shift (-3) extract 7;
type2:= N extract 3;
Li:= tr(line)
end else begin
boolean b11,b12,b21,b22;
nylinie: J1:= N1;
if J1<0 then goto stop;
read(in, k11, k12, J2, k21, k22, j, Li);
repeatchar(in);
for N:= readchar(in,k) while k<>10 and k<>9 do;
read(in,N1); qn(line,3):= 1; min:= 0;
b11:=k11=k11//2*2; b12:=k12=k12//2*2;
b21:=k21=k21//2*2; b22:=k22=k22//2*2;
type1:=if b11 and b12 then 1 else
if b11 and-,b12 then 2 else
if -,b11 and b12 then 3 else 4;
type2:=if b21 and b22 then 1 else
if b21 and-,b22 then 2 else
if -,b21 and b22 then 3 else 4;
b11:= b11==b21; b22:= b12==b22;
if j<>(if b11 and-,b22 then 1 else
if b22 and-,b11 then 3 else 2)
or abs(J1-J2)>1 then goto nylinie
end;
for j:=1,2 do imp(j,4):= 0;
for j:=1 step 1 until Q do ai(j):= 0;
for j:= 1,2,3 do imp(1,j):= imp(2,j):= 0;
dE:= 0;
for k:=-1,1 do begin
if k=-1 then begin
j0:= J1; type:= type1; nr:= (J1-k11)//2+1
end else begin
j0:= J2; type:= type2; nr:= (J2-k21)//2+1
end;
JJ:= j0*(j0+1);
matr:= MATRIX(if j0=j0//2*2 then type else (type+4));
k0:= if matr=1 then -2 else if matr=2 then 0 else -1;
N:= (j0-k0)//2; p:= if matr=1 then 2 else 1;
F:= (X - (D(1) - D(6)*JJ)*JJ)*JJ;
G:= Z - (D(2) - D(7)*JJ)*JJ;
G1:= -D(3) + D(8)*JJ; G2:= D(9);
H:= Y - (D(4) - D(10)*JJ)*JJ;
H1:= (-D(5) + D(11)*JJ)*0.5; H2:= D(12)*0.5;
if j0<>0 then begin
array alfa, beta, gamm, S, Pdif(1:N), PG(1:15), PS(1:6);
comment The matrix elements are calculated. Diagonal elements
in alfa and off diagonal elements in beta and gamm(a).
Pdif contains
Pdif(j) = <j-1 , Px**2 - Py**2 , j> * 2 ,
however Pdif(1) = < 1 , Px**2 - Py**2 , 1>. ;
beta(1):= Pdif(1):= gamm(1):= 0;
k0:= k0+2; K2:= k0*k0; Pz4:= K2*K2;
for j:=1 step 1 until N do begin
alfa(j):= F + (G + (G1 + G2*K2)*K2)*K2;
if j<>N then begin
p:= p*(j0-k0)*(j0-k0-1)*(j0+k0+1)*(j0+k0+2);
y:= H2*Pz4+H1*K2+H;
k0:= k0+2; K2:= k0*k0; Pz4:= K2; Pz4:= Pz4*K2;
beta(j+1):= (H2*Pz4+H1*K2+y)**2*p;
Pdif(j+1):= -sqrt(p);
gamm(j+1):= -DW*Pdif(j)*Pdif(j+1); p:= 1
end end j;
if matr>=3 then begin
j:= if matr=3 then -1 else 1;
alfa(1):= alfa(1) + (H+(H1+H2)*2)*j*JJ;
if N>1 then begin
p:= -DW*Pdif(2)*JJ*j; qd:= beta(2);
beta(2):= qd + p*p - sqrt(qd)*p*2
end;
Pdif(1):= j*JJ//2
end else if matr=2 then alfa(1):= alfa(1)+DW*(JJ-2)*JJ
else if N>1 then alfa(2):= alfa(2)-DW*(JJ-2)*JJ;
comment Eigenvalue and eigenvector are evaluated;
j:= if oblat then nr else N+1-nr;
H:= (if j>1 then sqrt(beta(j)) else 0) +
(if j<N then sqrt(beta(j+1)) else 0);
if j=1 or j=N then H:= H*3;
G:= alfa(j) + H; H:= alfa(j) - H;
for lambda:= (G+H)*0.5 while G<>lambda do begin
p:= 0; qd:= 1; it:= 0;
for i:=1 step 1 until N do begin
y:= (alfa(i)-lambda)*qd - beta(i)*p;
p:= qd; qd:= y;
if p>=0 == qd>=0 then it:=it+1;
if it=nr then i:= N
end sturm;
if it=nr then H:= lambda else G:= lambda
end bisection;
begin
real a1,a2,b1,c1,c2,s1,s2,ab1,ac2;
array a, b(1:N), d(-2:1);
for i:= 2 step 1 until N do beta(i):= - sqrt(beta(i));
d(-2):= gamm(j); d(-1):= beta(j); G:= alfa(j);
for i:= j+1 step 1 until N do begin
alfa(i-1):= alfa(i); beta(i-1):= beta(i); gamm(i-1):= gamm(i)
end;
alfa(N):= beta(N):= gamm(N):= 0;
d(0):= beta(j); beta(j):= gamm(j); gamm(j):= 0;
if j<N then begin d(1):= gamm(j+1); gamm(j+1):= 0 end
else d(1):= 0;
it:= 0; N:= N-1; F:= lambda;
repeat:
for i:=1 step 1 until N do S(i):= 0;
for i:=-2 step 1 until 1 do if d(i)<>0 then S(i+j):= -d(i);
a1:= a2:= b1:= c1:= c2:= s1:= s2:= ab1:= ac2:= 0;
for i:=N step -1 until 1 do begin
a(i):= 1/(alfa(i)-lambda-b1*ab1-c2*ac2);
b(i):= beta(i)-c1*ab1;
if i<=j+1 then S(i):= S(i)-s1*ab1-s2*ac2;
s2:= s1; s1:= S(i); a2:= a1; a1:= a(i); b1:= b(i);
c2:= c1; c1:= gamm(i); ab1:= a1*b1; ac2:= a2*c2
end;
s1:= s2:= 0; norm:= 1;
for i:=1 step 1 until N do begin
S(i):= (S(i)-s1*b(i)-s2*gamm(i))*a(i);
s2:= s1; s1:= S(i); norm:= norm + s1*s1
end;
norm:= 1/norm;
if N>1 and DW<>0 then begin
H:= G-lambda;
for i:=-2 step 1 until 1 do
if d(i)<>0 then H:= H + d(i)*S(j+i);
H:= H*norm;
if it<2 or (abs H<abs F and lambda+H<>lambda) then begin
lambda:= lambda+H; it:= it+1; F:= H; goto repeat
end end end;
for i:=N step -1 until j do S(i+1):= S(i);
S(j):= 1; N:= N+1;
Pz:= Pz4:= Pz6:= Pxy:= Pzxy:= Pxyxy:= Pzzd:=
G:= F:= H:= H1:= 0; j:= 0;
for i:=N step -1 until 1 do begin
G:= S(i); H:= H*G+H1; H1:= F; F:= F*G; G2:= G*G*K2;
Pz := Pz + G2; G2:= G2*K2;
Pz4 := Pz4 + G2; G2:= G2*K2;
Pz6 := Pz6 + G2;
Pxy := Pxy + F;
Pzxy:= Pzxy + (j+K2)*F;
Pzzd:= Pzzd + F*j*j + F*K2*K2;
Pxyxy:= Pxyxy + H*H;
H:= Pdif(i); F:= H*G; k0:= k0-2; j:= K2; K2:= k0*k0
end;
H:= F+H1/2; F:= F*G;
Pxy := (Pxy + F)*norm;
Pzxy:= (Pzxy + F + F)*norm;
Pzzd:= (Pzzd + F + F)*norm;
Pxyxy:= (Pxyxy/4+H*H)*norm;
PG(zz):= Pz:= Pz*norm; F:= JJ-Pz;
PG(xx):= (F+Pxy)*0.5;
PG(yy):= (F-Pxy)*0.5;
Pz4:= Pz4*norm;
for i:=1,2,3 do imp((3+k)//2,i):= PG(i);
imp(1,4):= imp(1,4) + Pz4;
imp(2,4):= imp(2,4) + Pz4*k;
imp((3+k)//2,5):= lambda;
if hb then begin
G:= (JJ-Pz-Pz)*JJ+Pz4; F:= (JJ+JJ)*Pxy-Pzxy;
PS(xx):= (G+F+Pxyxy)/16;
PS(yy):= (G-F+Pxyxy)/16;
PS(zz):= Pz4/4;
PS(xy):= (G-Pxyxy)/8; G:= JJ*Pz-Pz4;
PS(yz):= (G+G-Pzxy)/8;
PS(zx):= (G+G+Pzxy)/8;
end calculation of standard angular momenta;
case version of begin
begin p:= 0.001;
PG(4):= -p*JJ*JJ; PG(5):= -JJ*Pz*p;
PG(6):= -Pz4*p; PG(7):= -JJ*Pxy*(p+p);
PG(8):= -Pzxy*p; PG(12):= Pz6*norm*'-6;
PG(15):= Pzzd*'-6;
for i:=7,8 do PG(i+6):= -PG(i)*JJ*p;
for i:=4,5,6 do PG(i+5):= -PG(i)*JJ*p
end vers1;
begin p:= 0.000001;
PG(4):= PS(1)+(PS(4)*r1+PS(5)*r2+PS(6)*r3)*0.5;
PG(5):= PS(2)+(PS(4)*r4+PS(5)*r5+PS(6)*r6)*0.5;
PG(6):= PS(3)+(PS(4)*r7+PS(5)*r8+PS(6)*r9)*0.5;
PG(7):= PS(4)*2; PG(8):= PS(5)*2;
PG(9):= p*JJ*JJ*JJ; PG(10):= p*JJ*JJ*Pz;
PG(11):= JJ*Pz4*p; PG(12):= Pz6*norm*p;
PG(13):= Pxy*JJ*JJ*(p+p); PG(14):= Pzxy*JJ*p;
PG(15):= Pzzd*p
end vers2;
begin p:= 0.001; F:= PS(1)+PS(2);
G:= F+PS(4)*r1+PS(5)*r2+PS(6)*r3;
PG(4):= (G + PS(4)*2)*p; PG(5):= G*0.5*p;
PG(6):= ((PS(1)-PS(2))*0.5+PS(4)*r4+PS(5)*r5+PS(6)*r6)*p;
PG(7):= (F+F+PS(3)+PS(4)*r7+PS(5)*r8+PS(6)*r9)*p;
PG(8):= p*p*JJ*JJ*Pz
end vers3;
begin
PG(4):= Pxyxy;
PG(5):= Pxy*JJ;
PG(6):= Pz4 + Pxyxy*r1*r1*2;
PG(7):= -(Pz + Pxy*r1)*JJ;
PG(8):= JJ*JJ;
for i:=4 step 1 until 8 do PG(i):= PG(i)*0.001
end vers4;
end;
dE:= dE + k*lambda;
if ID then begin
PG(3):= PG(3)+PG(1)*c2;
PG(2):= PG(2)-PG(1)*b2;
PG(1):= PG(1)*a2
end;
for j:=Q-1 step -1 until 1 do ai(j):= ai(j) + PG(j)*k
end matrix element block end k;
if iterer then begin
dE:= ai(Q):= Li - dE; it:= 0; s:= s + dE*dE
end else begin
it:= Li; Li:= dE; dE:= 0 end;
if ha then begin H:= 0;
for j:=Q-1 step -1 until 1 do
for k:=Q-1 step -1 until j do H:= H + ai(j)*ai(k)*col(j,k);
if iterer then H:= if 1-H<=0 or m1=1 then 9.9
else dE/sqrt(so*(1-H)) else dE:= sqrt(H*s)
end;
F:= G:= 0;
if quartic then
for j:=4 step 1 until Q1 do F:= F + ai(j)*DR(j);
if sextic then
for j:=Q1+1 step 1 until Q-1 do G:= G + ai(j)*DR(j);
writecr; write(res,<<ddd>,J1,k11,k12,<: :>,J2,k21,k22,<: :>,
<<dddddd.ddd>,Li,<<-ddd.ddd>,dE);
if -,ha then write(res,<: :>) else
if iterer and m1>1 then begin if abs H<10 or dE*dE<so*100
then write(res,<< -d.d>,H) else begin
for j:=l+1 step -1 until 1 do atot(lc,j):= 0;
write(res,<< -ddd>,H,<: ***excluded:>);
m1:= m1-1; s:= s-dE*dE; goto endline
end end else
if it<'5 and it>0 then write(res,<<_ddddd>,it)
else write(res,<: :>);
if F+G<>0 then begin
write(res,<<-ddddd.ddd>,F+G);
if delres then begin
J1:= -1; k:= if sextic then Q-1 else Q1;
flere: J1:= J1+5; J2:= J1+4; if J2>k then J2:= k;
writecr; write(res,sp,25);
for j:=J1 step 1 until J2 do begin
p:= ai(j)*DR(j);
if p=0 then write(res,<: 0 :>)
else write(res,<<-dddd.ddd>,p) end;
if J2<>k then goto flere;
writecr
end else begin
if G<>0 then begin
if abs(G)<'-2 then write(res,<<__-d.d'-d>,G)
else write(res,<<-ddd.ddd>,G,<: :>) end;
p:= 0;
for j:=1 step 1 until m do begin
k:= nf(j); if k>3 then p:= p+ai(k)*DR(k) end;
if p<>0 then begin
if abs(p)<'-2 then write(res,<<__-d.d'-d>,p)
else write(res,<<-ddd.ddd>,p,<: :>) end end
end ;
if aires then begin writecr;
writecr; write(res,<:dE:dA l.::>);
for j:= 1,2,3 do
write(res,<< -d.ddd ddd'-d>,imp(1,j));
writecr; write(res,<:dE:dA u.::>);
for j:= 1,2,3 do
write(res,<< -d.ddd ddd'-d>,imp(2,j));
writecr; write(res,<: u. - l. ::>);
for j:= 1,2,3 do
write(res,<< -d.ddd ddd'-d>,imp(2,j)-imp(1,j));
writecr; write(res,<: u. + l. ::>);
for j:= 1,2,3 do
write(res,<< -d.ddd ddd'-d>,imp(2,j)+imp(1,j));
writecr; write(res,<:Pz4, u. + l. , u. - l. : :>,
<< -d.ddd ddd'-d>,imp(1,4),imp(2,4));
writecr; write(res,<:Energies: :>,
<< d.dddddddd'-dd>,imp(1,5),imp(2,5));
writecr end;
if iterer then begin
for j:=1 step 1 until m do
for k:=nf(j)+1-j step 1 until Q-j do ai(k):= ai(k+1);
for j:=1 step 1 until l+1 do atot(lc,j):= ai(j)
end;
endline: end step over observed transitions, line;
if ID then Rot(1):= Int(4);
if -,ha then begin
for i:=1 step 1 until N1 do begin
p:= atot(i,l+1)**2;
if m1>1 and (s-p)/m1 < p*'-4 then begin
for j:=l+1 step -1 until 1 do atot(i,j):= 0;
s:= s-p; m1:= m1-1
end end end;
s:= s/m1; hb:= so = 2**64; Q:= Q-1;
if so>=s then begin for j:=1 step 1 until l do begin
p:= 0;
for i:=1 step 1 until N1 do begin
qd:= abs atot(i,j); if qd>p then p:= qd
end;
lo(j):= p
end end;
Z:= ortho(atot,an,N1,l)
end atot;
comment The proceduren ortho solves the equation a x = b,
where the rectangular matrix a as well as b is stored in
atot. an is a l*(l+1) array containing by return the solu-
tion x (in column no. l+1) and the covariance matrix;
ha:= so<s; write(res,<:<10>:>);
if ha then begin
write(res,<:<10>Stop on divergency:>);
for i:=1 step 1 until 3 do Rot(i):= Rot(i)-AO(i);
for i:=4 step 1 until Q do DR(i):= DR(i)-AO(i)
end else begin
k:= Q; p:= l/(if N1=l then 0.1 else sqrt(s));
for i:=1 step 1 until l do lo(i):= lo(i)*p;
for j:=m step -1 until 0 do begin
for i:=nf(j)+1 step 1 until k do begin
p:= AO(i):= an(i-j,l+1);
hb:= hb or lo(i-j)*abs p > 1;
cif(i):= -entier(-ln(lo(i-j)*2)/ln10)+1;
if i<=3 then begin
Rot(i):= AO(i)+Rot(i);
cif(i):= entier(ln(abs Rot(i))/ln10)+cif(i)
end else begin
DR(i):= AO(i)+DR(i);
cif(i):= entier(ln(abs DR(i))/ln10)+cif(i)
end end;
k:= nf(j)-1
end;
j0:= Q;
for k0:=m step -1 until 0 do begin
for i:=nf(k0)+1 step 1 until j0 do begin
J2:= Q;
for K2:=m step -1 until k0 do begin
k:= if i>nf(K2) then i else nf(K2)+1;
for j:=k step 1 until J2 do
col(i,j):= an(i-k0,j-K2)*(if i=j then 1 else 2);
J2:= nf(K2)-1
end K2 end i;
i:= nf(k0); for j:=i step -1 until 1 do col(j,i):= 0;
if i<>0 then begin cif(i):= 0;
for j:=i+1 step 1 until Q do col(i,j):= 0
end; j0:= i-1
end k0 end;
y:= systime(1,starttime,p); qd:= p-time; time:= p;
iterer:= (hb or so>s+s) and p+qd<stoptime;
if -,iterer and N1=l then s:= 1; so:= s;
write(res,<:<10>rms deviation: :>,<<dd.dddd>,sqrt(s),
<: time, cpu: :>,(y-cpu)/60,<: , real: :>,p/60);
if -,ha then write(res,<:
Gramdet.: :>,<<d.d'-dd>,Z);
if closeres then begin
write(out,<:rms::>,<<__dd.dddd>,sqrt(s),<: tid::>,
time/60,<:<10>:>); forceout(out)
end;
writepage end;
goto igen;
stop:
if -,slut then goto dataind end;
write(res,<:<10><12><25>:>); close(res,closeres) end
▶EOF◀