|
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: 45312 (0xb100)
Types: TextFile
Names: »asm«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦7e928b248⟧ »algbib«
└─⟦this⟧
;algasmixx
;gosav
asmixx=set 91
permanent asmixx.17
asmixx=algol list.no index.no
\f
begin
integer i,j,k,l,m,m1,N,N1,J1,J2,Q,linecount,page,space,version;
real s, so, starttime, stoptime, cpu, time, Z, steplength;
boolean slut,deltarot,deltacd,iterer,oblat,atype,
closeres,quartic,couple,fixed,hb;
array Rot(1:9), D(1:10), head(1:12);
zone res(128,1,stderror);
procedure writecr;
begin linecount:= linecount+1;
write(res,<:<10>:>); if linecount>59 then writepage
end writecr;
procedure writepage;
begin integer i;
page:= page+1; linecount:= 3; i:= 1;
write(res,<:<12>:>,false add 10,3,string head(increase(i)),
false add 32, space,<<dd>, page, false add 10, 2)
end writepage;
closeres:=outmedium(res);
dataind:
cpu:= systime(1,0,starttime);
space:= 68 - readhead(in,head,1);
read(in,m1); repeatchar(in); readchar(in,j);
atype:= false add j;
read(in,stoptime,Rot,N,steplength);
stoptime:= stoptime*60;
slut:= m1>0; m1:= abs m1; m:= m1//1000; fixed:=m>0;
deltarot:= m1 mod 1000 >= 100; version:= m1 mod 10;
deltacd := m1 mod 100 >= 10; iterer := N<>0;
page:= 0; writepage; write(res,<:
Calculation by ASMIX-X,YZ of 23-3-1976
:>);
case version of begin
begin
linecount:= linecount+4; Q:=19; write(res,<:
D1 = Dj, D2 = Djk, D3 = Dk, D4 = dj, D5 = djk, (kHz)
:>) end vers1;
begin
linecount:= linecount+6; Q:=19; 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), (MHz).
:>) end vers2;
begin
Q:=9
end vers3;
begin
Q:=19; 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 vers4;
end;
oblat:= 2*Rot(2)>(Rot(1)+Rot(3));
j:= atype extract 12;
atype:= j<>98;
if atype and j<>97 then atype:= oblat;
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,v1,v2,
xxt,yyt,zzt,xyt,yzt,zxt;
array DR(1:10),tr(1:N1),an(1:l,1:l+1),AO(1:Q),cif(1:l),DW(0:1);
integer array qn(1:N1,1:3), MATRIX(1:16);
boolean array nf(1:Q);
min:=i:= 1;
if version<>3 then for i:=1,i+1 while k<>10 do begin
if i>10 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 10 do DR(j):= 0;
for j:=1 step 1 until 10 do D(j) := 0;
for j:=1 step 1 until Q do nf(j):=false;
for j:=1 step 1 until m do begin
read(in,k); nf(k):=true; AO(k):=0
end;
if iterer then begin
integer q,max,my;
boolean b11,b12,b21,b22;
max:=1; writecr;writecr;write(res,<:Input transitions:
_(v->v'),(J,k-,k+)->(J'k-'k+')_dip______freq.:>); writecr;
for k:=1 step 1 until N do begin
read(in, v1, v2, J1, k11, k12, J2, k21, k22, my, tr(k));
writecr; write(res,<<ddd>,v1,v2,<:__:>,J1,k11,k12,<:__:>,J2,k21,
k22,<:__:>,my,<<___ddd ddd.ddd>,tr(k));
repeatchar(in); for i:=readchar(in,j) while j<>10 and j<>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,<:***fault***:>);
if k=min then max:=min:=min+1;
goto linieslut end;
qn(k,1):= v1 shift 2 add v2 shift 6 add J1 shift 6
add k11 shift 6 add k12;
qn(k,2):= J2 shift 6 add k21 shift 6 add k22 shift 3
add type1 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; writecr
end reading and control of transitions;
write(res,<:Steplength::>,<< d.dd>,steplength);
writecr; writecr;
if m=1 then write(res,<:Fixed Constant::>) else
if m>1 then write(res,<:Fixed Constants::>);
if fixed then for i:=1 step 1 until Q do
if nf(i) then begin
writecr; AO(i):= 0;
if i<10 then begin
so:= Rot(i);
if deltarot or i>6 then write(res,<:_____:>,case i of
(<:_(A0+A1)/2:>,<:_(B0+B1)/2:>,<:_(C0+C1)/2:>,
<:___(A0-A1):>,<:___(B0-B1):>,<:___(C0-C1):>,
<:delta Evib:>,<:__Gvib-rot:>,<:__Fvib-rot:>),<: = :>)
else write(res,<:_____________:>,case i of
(<:A0 = :>,<:B0 = :>,<:C0 = :>,
<:A1 = :>,<:B1 = :>,<:C1 = :>));
write(res,if abs so=0 then <<d> else <<-ddd ddd.dddd>,so);
end else begin
k:=i-9; so:= DR(k);
if deltacd then write(res,case k of
(<:(D1(0)+D1(1))/2 = :>,<:(D2(0)+D2(1))/2 = :>,
<:(D3(0)+D3(1))/2 = :>,<:(D4(0)+D4(1))/2 = :>,
<:(D5(0)+D5(1))/2 = :>,<: D1(0)-D1(1) = :>,
<:____D2(0)-D2(1) = :>,<:____D3(0)-D3(1) = :>,
<:____D4(0)-D4(1) = :>,<:____D5(0)-D5(1) = :>))
else begin
j:=if k>5 then k-5 else k;
write(res,<:__________D:>,<<d>,j,if k>5 then <:(1) = :>
else <:(0) = :>)
end;
write(res,if abs so=0 then <<d> else <<-ddd.ddd0000>,so);
end;
end;
if fixed then writecr; systime(1,starttime,time);
s:= 0; so:= 2**64;
m1:= if -,iterer then 0 else if N1=l then 1 else N1-l;
if atype then begin
for i:=1 step 1 until 16 do
MATRIX(i):= case i of (1,4,2,3,2,3,1,4,4,1,3,2,3,2,4,1);
xx:= 1; yy:= 2; zz:= 3; xy:= 4; yz:= 5; zx:= 6
end else begin
for i:=1 step 1 until 16 do
MATRIX(i):= case i of (1,2,3,4,2,1,4,3,4,3,2,1,3,4,1,2);
zz:= 1; xx:= 2; yy:= 3; zx:= 4; xy:= 5; yz:= 6
end;
if oblat then begin
xxt:=1; yyt:=2; zzt:=3; xyt:=4; yzt:=5; zxt:=6;
end else begin
zzt:=1; xxt:=2; yyt:=3; zxt:=4; xyt:=5; yzt:=6;
end;
igen:
quartic:= false; couple:=Rot(7)<>0 or Rot(8)<>0 or Rot(9)<>0;
for i:= 1 step 1 until 10 do quartic:= quartic or DR(i)<>0;
begin
real p,pp,y,a0,b0,c0,a1,b1,c1,r1,r2,r3,r4,r5,r6,r7,r8,r9,
s1,s2,s3,s4,s5,s6,s7,s8,s9,X0,Y0,Z0,X1,Y1,Z1;
boolean ha,hm;
array ai(1:Q+1);
ha:= s=so; hb:=false;
if deltarot then begin
p:=Rot(1);pp:=Rot(4);a0:=ai(1):=0.5*(p+p+pp);a1:=ai(4):=a0-pp;
p:=Rot(2);pp:=Rot(5);b0:=ai(2):=0.5*(p+p+pp);b1:=ai(5):=b0-pp;
p:=Rot(3);pp:=Rot(6);c0:=ai(3):=0.5*(p+p+pp);c1:=ai(6):=c0-pp;
X0:=ai(xx); Y0:=ai(yy); Z0:=ai(zz);
X1:=ai(3+xx); Y1:=ai(3+yy); Z1:=ai(3+zz);
end else begin
X0:=Rot(xx); Y0:=Rot(yy); Z0:=Rot(zz);
X1:=Rot(3+xx); Y1:=Rot(3+yy); Z1:=Rot(3+zz);
a0:=Rot(1); b0:=Rot(2); c0:=Rot(3);
a1:=Rot(4); b1:=Rot(5); c1:=Rot(6);
end;
if -,ha and iterer then writepage;
if deltacd then begin
for i:=1 step 1 until 5 do begin
j:=i+5; p:=DR(i); pp:=DR(j); y:=D(i):=0.5*(p+p+pp); D(j):=y-pp
end i;
end else for i:=1 step 1 until 10 do D(i):=DR(i);
case version of begin
begin
if -,oblat==atype then begin
ai(xxt):= -(D(1) + D(4)*2)*4;
ai(yyt):= -(D(1) - D(4)*2)*4;
ai(zzt):= -(D(1) + D(2) + D(3))*4;
ai(xyt):= -D(1)*4;
ai(yzt):= -(D(1) + D(2)*0.5 - D(4) - D(5))*4;
ai(zxt):= -(D(1) + D(2)*0.5 + D(4) + D(5))*4;
ai(6+xxt):= -(D(6) + D(9)*2)*4;
ai(6+yyt):= -(D(6) - D(9)*2)*4;
ai(6+zzt):= -(D(6) + D(7) + D(8))*4;
ai(6+xyt):= -D(6)*4;
ai(6+yzt):= -(D(6) + D(7)*0.5 - D(9) - D(10))*4;
ai(6+zxt):= -(D(6) + D(7)*0.5 + D(9) + D(10))*4;
hb:= true
end else begin
for i:=1 step 1 until 10 do D(i):= D(i)*'-3;
DW(0):=DW(1):=0
end end vers1;
begin
comment The standard constants are calculated in ai (in kHz)
using the Dowling relations and assuming tau-caca=0;
r4:= a0/b0; r8:= b0/c0; r3:= c0/a0;
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;
s4:= a1/b1; s8:= b1/c1; s3:= c1/a1;
s4:= -s4*s4; s8:= s8*s8; s3:= s3*s3;
s1:= 1/s4; s5:= 1/s8; s9:= 1/s3;
s2:= s1*s3; s6:= s4*s5; s7:= s8*s9;
for j:=1,2,3 do ai(j):= D(j)*'3;
for j:=7,8,9 do ai(j):=D(j-1)*'3;
ai(4) :=(D(1)*r1 + D(2)*r4 + D(3)*r7)*0.5 + D(4)*2;
ai(5) :=(D(1)*r2 + D(2)*r5 + D(3)*r8)*0.5 + D(5)*2;
ai(6) :=(D(1)*r3 + D(2)*r6 + D(3)*r9)*0.5;
ai(10):=(D(6)*s1 + D(7)*s4 + D(8)*s7)*0.5 + D(9)*2;
ai(11):=(D(6)*s2 + D(7)*s5 + D(8)*s8)*0.5 + D(10)*2;
ai(12):=(D(6)*s3 + D(7)*s6 + D(8)*s9)*0.5;
for j:=4,5,6,10,11,12 do ai(j):= ai(j)*'3;
hb:=true
end vers2;
begin
for i:=1 step 1 until 10 do D(i):=0; DW(0):=DW(1):=0
end vers3;
begin
if -,oblat and atype then begin
r1:= (Y0-Z0)/(X0+X0-Y0-Z0);
s1:= (Y1-Z1)/(X1+X1-Y1-Z1)
end else begin
r1:= (X0-Y0)/(Z0+Z0-X0-Y0);
s1:= (X1-Y1)/(Z1+Z1-X1-Y1)
end;
a0:=D(1)+D(3)*r1*r1*2;
b0:=D(2)-D(4)*r1;
ai(xxt):=(a0+b0+D(5))*4;
ai(yyt):=(a0-b0+D(5))*4;
ai(zzt):=(D(3)-D(4)+D(5))*4;
ai(yzt):=(-b0-D(4))*2+D(5)*4;
ai(zxt):=( b0-D(4))*2+D(5)*4;
ai(xyt):=(-a0+D(5))*4;
a1:=D(6)+D(8)*s1*s1*2;
b1:=D(7)-D(9)*s1;
ai(6+xxt):=(a1+b1+D(10))*4;
ai(6+yyt):=(a1-b1+D(10))*4;
ai(6+zzt):=(D(8)-D(9)+D(10))*4;
ai(6+yzt):=(-b1-D(9))*2+D(10)*4;
ai(6+zxt):=( b1-D(9))*2+D(10)*4;
ai(6+xyt):=(-a1+D(10))*4;
hb:=true
end vers4;
end case version;
if hb then begin
comment Calculation of D-constants from standard constants,
and correction of rotational constants;
p:= -1/16000;
D(1) := ((ai(xx)+ai(yy))*1.5+ai(xy))*p;
D(2) := -D(1)*2 + (ai(yz)+ai(zx))*4*p;
D(3) := -D(1)-D(2)+ai(zz)*4*p;
D(4) := (ai(xx)-ai(yy))*p;
D(5) := -D(4)-(ai(yz)-ai(zx))*2*p;
DW(0):= (ai(xx)+ai(yy)-ai(xy)*2)*p;
D(6) := ((ai(6+xx)+ai(6+yy))*1.5+ai(6+xy))*p;
D(7) := -D(6)*2 + (ai(6+yz)+ai(6+zx))*4*p;
D(8) := -D(6)-D(7)+ai(6+zz)*4*p;
D(9) := (ai(6+xx)-ai(6+yy))*p;
D(10):= -D(9)-(ai(6+yz)-ai(6+zx))*2*p;
DW(1):= (ai(6+xx)+ai(6+yy)-ai(6+xy)*2)*p;
X0:=X0+DW(0); Y0:=Y0+DW(0); Z0:=Z0-DW(0)*1.5;
X1:=X1+DW(1); Y1:=Y1+DW(1); Z1:=Z1-DW(1)*1.5;
DW(0):=DW(0)/4; DW(1):=DW(1)/4;
end hb;
write(res,<:<10><10>:>,if deltarot then
<:Mean- and differences of :> else <::>,
<:Rotational constants (MHz):<10> :>); k:=0; m:=5;
for j:=1,4 do begin
for i:=j step 1 until j+2 do
write(res,<<__-ddd_ddd.dddddd>,Rot(i));
if ha then begin
write(res,<:<10>+-:>);
for i:=j step 1 until j+2 do
if nf(i) then write(res,<:_______Fixed_____:>)
else begin
k:=k+1; write(res,<<___dd0_000.000000>,sqrt(an(k,k)*s))
end nf;
m:=m+1
end ha;
if j=1 then write(res,<:<10><10> :>);
end j;
if couple then begin
m:=m+4;
write(res,<:<10>
Energy difference between coupling states
and :>,false add (if atype then 97 else 98),1,
<:-type coupling constant (MHz):<10> :>);
for i:=7,8,9 do write(res,<<_-d_ddd_ddd.dddddd>,Rot(i));
if ha then begin
write(res,<:<10>+-:>);
for i:=7,8,9 do
if nf(i) then write(res,<:________Fixed_____:>)
else begin
k:=k+1; write(res,<<__d_d00_000.000000>,sqrt(an(k,k)*s))
end nf;
m:=m+1
end ha;
end couple;
if quartic then begin
m:=m+5;
write(res,<:<10><10>:>, if deltacd then
<:Mean- and differences of :> else <::>,
<:Quartic distortion constants:<10> :>);
for j:=1,6 do begin
for i:=j step 1 until j+4 do
write(res,<<-dddd.dddd000>,DR(i));
if ha then begin
write(res,<:<10>+-:>);
for i:=j+9 step 1 until j+13 do
if nf(i) then write(res,<:____Fixed____:>)
else begin
k:=k+1; write(res,<<_dd00.0000000>,sqrt(an(k,k)*s))
end nf;
m:=m+1
end ha;
if j=1 then write(res,<:<10><10> :>);
end j;
end quartic;
linecount:=linecount+m;
if ha and -,iterer then begin
writecr; writecr; writecr;
write(res,<:Significant Digits and Correlation Matrix::>);
pp:= ln(10);
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; k:= 1;
Cor: for i:=1 step 1 until N do begin writecr;
for j:=k step 1 until i-1 do write(res,false add 32,8);
if j=i then begin j:= j+1;
write(res,<<__dd>,entier(ln(cif(i))/pp)+2,<:____:>) end;
for j:=j step 1 until N do
write(res,<<__-d.ddd>,an(i,j)*AO(i)*AO(j)*0.5)
end i;
writecr; k:= N+1; N:= N+6; if k<l then goto Cor;
linecount:= linecount+4 end ha;
write(res,false add 10,3);
if iterer then begin
write(res,false add 32,30,<:obs______obs.-_______:>);
if couple then write(res,<:coupling__:>);
if quartic then write(res,<:centrifugal:>); writecr;
write(res,<:_______Transition__________frequency___calc.___t:>,
<:___:>);
end else begin
read(in,N1); if N1<0 then goto stop;
if linecount > 50 then writepage;
write(res,false add 32,29,<:calc.____________:>);
if couple then write(res,<:coupling__:>);
if quartic then write(res,<:centrifugal:>);
write(res,<:
_______Transition_________frequency____sigma__:>);
end;
if couple then write(res,<:contrib.__:>);
if quartic then write(res,<:distortion:>);
writecr; if iterer then s:=0;
p:=Y0; Y0:=(X0-p)*0.25; X0:=(X0+p)*0.5; Z0:=Z0-X0;
p:=Y1; Y1:=(X1-p)*0.25; X1:=(X1+p)*0.5; Z1:=Z1-X1;
begin
integer line, matr, nr, type, JJ, J0, vk, i1, i2, k1, k2;
real lambda,Li,dE,Px,Pyz,S,S2;
array atot(1:if iterer then N1 else 1,1:l+1),PG(1:Q),
PS(1:if hb then 12 else 1),
Pxy,Pz,Pzxy,Pxyxy,Pz4,F,G,H,G1,H1,G2,norm(0:1);
integer array kk(0:1);
for line:= min, qn(line,3) while line<>min do begin
if iterer then begin
N:= qn(line,1);
v1:=N shift (-20) extract 2;
v2:=N shift (-18) extract 2;
J1:=N shift (-12) extract 6;
k11:=N shift (-6) extract 6;
k12:=N extract 6;
N:= qn(line,2);
J2:=N shift (-18) extract 6;
k21:=N shift (-12) extract 6;
k22:=N shift (-6) extract 6;
type1:=N shift (-3) extract 3;
type2:=N extract 3;
Li:= tr(line)
end else begin
boolean b11,b12,b21,b22;
nylinie: v1:= N1;
if v1<0 then goto stop;
read(in, v2, J1, k11, k12, J2, k21, k22, j);
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 step 1 until Q+1 do ai(j):= 0;
dE:= 0;
for k:=-1,1 do begin
if k=-1 then begin
J0:=J1; type:=type1; nr:=(J1-k11)//2; vk:=v1
end else begin
J0:=J2; type:=type2; nr:=(J2-k21)//2; vk:=v2
end;
if J0=0 then begin
if vk=1 then begin dE:=dE+Rot(7)*k; ai(7):=ai(7)+k end
end else begin
if J0//2*2<>J0 then type:=type+4; matr:=MATRIX(type);
j:=if matr=1 then 2 else if matr=2 then 0 else 1;
nr:=(J0+j)//2-nr; JJ:=J0*(J0+1);
if vk=1 then matr:=MATRIX(type+8);
hm:=matr=1 or matr=3;
if hm then begin
i:=1; j:=0; i1:=1
end else begin
i:=0; j:=1; i1:=0
end;
F(i) :=(X0-D(1)*JJ)*JJ; F(j) :=(X1-D(6)*JJ)*JJ+Rot(7);
G(i) :=Z0-D(2)*JJ; G(j) :=Z1-D(7)*JJ;
H(i) :=Y0-D(4)*JJ; H(j) :=Y1-D(9)*JJ;
G1(i):=-D(3); G1(j):=-D(8);
H1(i):=-0.5*D(5); H1(j):=-0.5*D(10);
G2(i):=DW(0); G2(j):=DW(1);
Px:=0.5*Rot(8); Pyz:=0.5*Rot(9);
if matr=1 or matr=4 then begin
p:=2; k1:=-1; N:=J0+1
end else begin
p:=1; k1:=0; N:=J0
end;
begin
array u(1:N,1:N),Pdif,Pdifx,ev(1:N);
integer array sym(1:N);
for i:=1 step 1 until N do
for j:=1 step 1 until i do u(i,j):=0;
Pdif(1):=Pdifx(1):=0; if N>1 then Pdif(2):=0;
for i:=1 step 1 until N do begin
i2:=i mod 2; k1:=k1+1; k2:=k1*k1;
u(i,i):=F(i2)+(G(i2)+G1(i2)*k2)*k2;
if i<N then begin
p:=p*(J0-k1)*(J0+k1+1); j:=i+1;
Pdifx(j):=sqrt(p);
u(j,i):=Pdifx(j)*(Pyz*(k1+k1+1)
+(if i1=i2 then Px else -Px));
if j<N then begin
p:=p*(J0-k1-1)*(J0+k1+2); j:=j+1;
Pdif(j):=-sqrt(p);
u(j,i):=(H(i2)+H1(i2)*(k2+(k1+2)**2))*Pdif(j);
end;
end;
if i>4 then u(i,i-4):=-G2(i2)*Pdif(i)*Pdif(i-2);
p:=1
end;
if matr=1 or matr=4 then begin
u(2,2):=u(2,2)+(H(0)+H1(0)*2)*JJ;
Pdif(2):=JJ//2;
if N>3 then begin
u(4,2):=u(4,2)-G2(0)*Pdif(4)*JJ;
u(3,3):=u(3,3)-G2(1)*(JJ-2)*JJ
end;
end else begin
u(1,1):=u(1,1)-(H(1)+H1(1)*2)*JJ;
Pdif(1):=-JJ//2;
if N>2 then u(3,1):=u(3,1)+G2(1)*Pdif(3)*JJ;
if N>1 then u(2,2):=u(2,2)+G2(0)*(JJ-2)*JJ
end;
tridql(N,ev,u);
k1:=(N+1)//2; k2:=N//2;
for i:=1 step 1 until N do begin
p:=pp:=0;
for j:=1 step 2 until N do
p:=p+u(i,j)**2;
for j:=2 step 2 until N do
pp:=pp+u(i,j)**2; if k2=0 then pp:=0;
sym(i):=if (p>pp)==hm then 0 else 1;
end;
k1:=0;
for i:=1,i+1 while k1<>nr do begin
if vk=sym(i) then k1:=k1+1;
end;
nr:=i-1; lambda:=ev(nr);
for i:=0,1 do norm(i):=F(i):=G(i):=H(i):=H1(i):=Pz(i):=
Pz4(i):=Pxy(i):=Pzxy(i):=Pxyxy(i):=0;
kk(0):=kk(1):=0; p:=pp:=Px:=Pyz:=0; k1:=J0; k2:=k1*k1;
for i:=N step -1 until 1 do begin
i2:=i mod 2; S:=u(nr,i); S2:=S*S;
norm(i2):=norm(i2)+S2;
H(i2):=H(i2)*S+H1(i2); H1(i2):=F(i2);
F(i2):=F(i2)*S; p:=p*S; pp:=pp*S;
Px:=Px+p; Pyz:=Pyz+pp; S2:=S2*k2;
p:=Pdifx(i)*S*(if i1=i2 then -1 else 1);
pp:=Pdifx(i)*S*(k1+k1-1);
Pz(i2):=Pz(i2)+S2;
Pz4(i2):=Pz4(i2)+S2*k2;
Pxy(i2):=Pxy(i2)+F(i2);
Pzxy(i2):=Pzxy(i2)+(kk(i2)+k2)*F(i2);
Pxyxy(i2):=Pxyxy(i2)+H(i2)*H(i2); H(i2):=Pdif(i);
F(i2):=H(i2)*S; k1:=k1-1; kk(i2):=k2; k2:=k1*k1;
end;
G(1):=u(nr,1); if N>1 then G(0):=u(nr,2);
for i:=0,1 do begin
H(i):=F(i)+H1(i)/2; F(i):=F(i)*G(i);
Pxy(i):=Pxy(i)+F(i); Pzxy(i):=Pzxy(i)+F(i)*2;
Pxyxy(i):=Pxyxy(i)/4+H(i)*H(i);
F(i):=norm(i)*JJ-Pz(i);
end i;
if hm then begin
i:=1; j:=0
end else begin
i:=0; j:=1
end;
PG(zz):=Pz(i); PG(3+zz):=Pz(j);
PG(xx):=0.5*(F(i)+Pxy(i)); PG(3+xx):=0.5*(F(j)+Pxy(j));
PG(yy):=0.5*(F(i)-Pxy(i)); PG(3+yy):=0.5*(F(j)-Pxy(j));
PG(7):=norm(j); PG(8):=Px; PG(9):=Pyz;
if hb then begin
y:=(norm(i)*JJ-Pz(i)*2)*JJ+Pz4(i);
pp:=2*JJ*Pxy(i)-Pzxy(i);
PS(xx):=(y+pp+Pxyxy(i))/16;
PS(yy):=(y-pp+Pxyxy(i))/16;
PS(zz):=Pz4(i)/4;
PS(xy):=(y-Pxyxy(i))/8;
y:=JJ*Pz(i)-Pz4(i);
PS(yz):=(y+y-Pzxy(i))/8;
PS(zx):=(y+y+Pzxy(i))/8;
y:=(norm(j)*JJ-Pz(j)*2)*JJ+Pz4(j);
pp:=2*JJ*Pxy(j)-Pzxy(j);
PS(6+xx):=(y+pp+Pxyxy(j))/16;
PS(6+yy):=(y-pp+Pxyxy(j))/16;
PS(6+zz):=Pz4(j)/4;
PS(6+xy):=(y-Pxyxy(j))/8;
y:=JJ*Pz(j)-Pz4(j);
PS(6+yz):=(y+y-Pzxy(j))/8;
PS(6+zx):=(y+y+Pzxy(j))/8;
if version=1 or version=4 then begin
for k1:=0,1 do begin
m:= if k1=i then 0 else 3;
Pz(k1):= PG(zzt+m);
Pxy(k1):= PG(xxt+m)-PG(yyt+m); m:= m+m;
Pz4(k1):= PS(zzt+m)*4;
Pzxy(k1):= (PS(zxt+m)-PS(yzt+m))*4;
Pxyxy(k1):= (PS(xxt+m)+PS(yyt+m)-PS(xyt+m))*4
end end end hb;
case version of begin
begin
p:=-0.001;
PG(10):=p*JJ*JJ*norm(i);
PG(11):=p*JJ*Pz(i);
PG(12):=p*Pz4(i);
PG(13):=(p+p)*JJ*Pxy(i);
PG(14):=p*Pzxy(i);
PG(15):=p*JJ*JJ*norm(j);
PG(16):=p*JJ*Pz(j);
PG(17):=p*Pz4(j);
PG(18):=(p+p)*JJ*Pxy(j);
PG(19):=p*Pzxy(j)
end vers1;
begin
p:=0.000001;
PG(10):= PS(1)+(PS(4)*r1+PS(5)*r2+PS(6)*r3)*0.5;
PG(11):= PS(2)+(PS(4)*r4+PS(5)*r5+PS(6)*r6)*0.5;
PG(12):= PS(3)+(PS(4)*r7+PS(5)*r8+PS(6)*r9)*0.5;
PG(13):= PS(4)*2; PG(14):= PS(5)*2;
PG(15):=PS(7)+(PS(10)*s1+PS(11)*s2+PS(12)*s3)*0.5;
PG(16):=PS(8)+(PS(10)*s4+PS(11)*s5+PS(12)*s6)*0.5;
PG(17):=PS(9)+(PS(10)*s7+PS(11)*s8+PS(12)*s9)*0.5;
PG(18):=PS(10)*2; PG(19):=PS(11)*2;
end vers2;
begin
end vers3;
begin
PG(10):=Pxyxy(i);
PG(11):=Pxy(i)*JJ;
PG(12):=Pz4(i)+Pxyxy(i)*r1*r1*2;
PG(13):=-(Pz(i)+Pxy(i)*r1)*JJ;
PG(14):=norm(i)*JJ*JJ;
PG(15):=Pxyxy(j);
PG(16):=Pxy(j)*JJ;
PG(17):=Pz4(j)+Pxyxy(j)*s1*s1*2;
PG(18):=-(Pz(j)+Pxy(j)*s1)*JJ;
PG(19):=norm(j)*JJ*JJ;
for i:=10 step 1 until 19 do PG(i):=PG(i)*0.001
end vers4;
end case version;
dE:= dE + k*lambda;
for j:=1 step 1 until Q do ai(j):= ai(j) + PG(j)*k
end u;
end J0<>0;
end k;
if iterer then begin
dE:= ai(Q+1):= Li - dE; s:= s + dE*dE
end else begin
Li:= dE; dE:= 0
end;
if deltarot then for i:=1,2,3 do begin
p:=ai(i); j:=i+3; ai(i):=p+ai(j); ai(j):=0.5*(p-ai(j))
end;
if quartic then begin
if deltacd then for i:=10 step 1 until 14 do begin
p:=ai(i); j:=i+5; ai(i):=p+ai(j); ai(j):=0.5*(p-ai(j))
end deltacd;
pp:=0;
for i:=10 step 1 until 19 do pp:=pp+ai(i)*DR(i-9)
end quartic;
if couple then p:=ai(8)*Rot(8)+ai(9)*Rot(9);
k:=Q+1; k1:=-1;
if fixed then
for i:=1 step 1 until Q do
if nf(i) then begin
k:=k-1; k1:=k1+1; for j:=i-k1 step 1 until k do ai(j):=ai(j+1);
end nf;
if ha then begin y:= 0;
for j:=l step -1 until 1 do
for k:=l step -1 until j do y:= y + ai(j)*ai(k)*an(j,k);
if iterer then y:= if 1-y<=0 then 9.9 else dE/sqrt(so*(1-y))
else dE:= sqrt(y*s)
end;
writecr; write(res,<<dd>,v1,v2,<: :>,<<ddd>,J1,k11,k12,<: :>,
J2,k21,k22,<: :>,<<ddd ddd.ddd>,Li,<<-ddd.ddd>,dE);
if iterer then write(res,<< -d.d>,if ha then y else 0);
if couple then write(res,<< -dddd.ddd>,p);
if quartic then write(res,<< -dddd.ddd>,pp);
if iterer then
for j:=1 step 1 until l+1 do atot(line,j):= ai(j)
end step over observed transitions, line;
s:= s/m1; hb:= so = 2**64;
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
pp:= abs atot(i,j); if pp>p then p:= pp
end;
cif(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;
write(res,<:<10>:>);
if so<s then begin
write(res,<:<10>Stop on divergency:>);
for i:=1 step 1 until 9 do Rot(i):= Rot(i)-AO(i);
if quartic then for i:=1 step 1 until 10 do DR(i):=DR(i)-AO(i+9)
end else begin
k:= 0; p:= l/sqrt(s);
for i:=1 step 1 until l do cif(i):= cif(i)*p;
for j:=1 step 1 until Q do begin
if -,nf(j) then begin
k:=k+1; p:=AO(j):=an(k,l+1)*steplength;
hb:=hb or cif(k)*p>steplength;
if j<10 then begin
Rot(j):=AO(j)+Rot(j);
cif(k):=cif(k)*abs Rot(j)
end else begin
DR(j-9):=AO(j)+DR(j-9);
cif(k):=cif(k)*abs DR(j-9)
end;
end nf;
end;
end so;
for i:=1 step 1 until l do
for j:=i+1 step 1 until l do an(i,j):=an(i,j)*2;
y:= systime(1,starttime,p); pp:= p-time; time:= p;
iterer:= (hb or so>s+s) and p+pp<stoptime; so:= s;
write(res,<:<10>rms deviation: :>,<<dd.dddd>,sqrt(s),
<: time, cpu: :>,(y-cpu)/60,<: , real: :>,p/60,<:
Gramdet.: :>,<<d.d'-dd>,Z);
if closeres then begin
write(out,<:rms::>,<<__dd.dddd>,sqrt(s),<:__tid::>,
p/60,<:<10>:>); setposition(out,0,0)
end;
writepage end;
goto igen;
stop:
if -,slut then goto dataind end;
write(res,<:<10><12><25>:>); close(res,closeres) end
;goasm
;gosav time.300 lines.5000
gtxt=set 50
permanent gtxt.15
gtxt=edit algasmixx c.50
l1,r/91/110/,
l2,r/no/yes/,
l./MATRIX/,r/;/, AUXMAT(1:8);/
l./if oblat/,l3,i/
for i:=1 step 1 until 8 do
AUXMAT(i):= case i of (1,2,3,4,2,1,4,3);
/,p1,l./integer line/,r/J0/J0,N0/
l1,r/dE/dE,dE0/,l./goto nylinie/,
l3,r/d/dE0:=d/
l./nr:=/,g2?//2?//2+1?
l./nr:=(J0+j)/,r/nr/N0/,r/-nr//
l./array u/,r/u(1:N,1:N)/e(1:N,0:4)/,r/ev/b/,l1,d2
l./u(i,i)/,r/u(i,i)/e(i,0)/,r/;/; e(i,3):= 0;/
l4,r/u(j,i)/e(j,1)/
l5,r/u(j,i)/e(j,2)/
l3,r/u(i,i-4)/e(i,4)/,l4,g/u(2,2)/e(2,0)/
l3,g/u(4,2)/e(4,2)/,l1,g/u(3,3)/e(3,0)/
l3,g/u(1,1)/e(1,0)/,l2,g/u(3,1)/e(3,2)/,l1,g/u(2,2)/e(2,0)/
l2,d./lambda/,l-1,r/;/; k1:=k;/,l1,i?
begin
integer it,i0,k; real g,h,y,p,q;
array a(1:N),d(-4:3),u(1:N,1:4);
j:=i2:= if i1=vk then 0 else -1;
if N0>2 and DW(0)<>0 then begin
comment The pentadiagonal matrix of the uncoupled
state is reduced to tridiagonal form by succesive
Jacobi rotations, Swartz: Numer. Math. 12, 231 (1968),
modified by replacing a(i,j) and b by v(i+j,j) and p;
real c,s,c2,s2,cs,u,u1;
array v(1:N0,0:2);
for i:=1 step 1 until N0 do begin
j:=j+2;
for k:=0,1,2 do v(i,k):=e(j,k+k)
end;
comment v(1,1):=v(1,2):=v(2,2):=0;
for k:=1 step 1 until N0-2 do begin
for j:=k+2 step 2 until N0 do begin
if j=k+2 then begin
if v(k+2,2)=0 then goto endk;
p:=-v(k+1,1)/v(k+2,2);
end else begin
if g=0 then goto endk;
p:=-v(j-1,2)/g;
end;
s2:=1/(p*p+1); s:=sqrt(s2); c:=p*s;
c2:=c*c; cs:=c*s;
u:=c2*v(j-1,0)-2*cs*v(j,1)+s2*v(j,0);
u1:=s2*v(j-1,0)+2*cs*v(j,1)+c2*v(j,0);
v(j,1):=cs*(v(j-1,0)-v(j,0))+(c2-s2)*v(j,1);
v(j-1,0):=u; v(j,0):=u1;
u:=c*v(j-1,1)-s*v(j,2);
v(j,2):=s*v(j-1,1)+c*v(j,2);
v(j-1,1):=u;
if j<>k+2 then
v(j-1,2):=c*v(j-1,2)-s*g;
if j<N0 then begin
u:=c*v(j+1,2)-s*v(j+1,1);
v(j+1,1):=s*v(j+1,2)+c*v(j+1,1);
v(j+1,2):=u
end;
if j+2<=N0 then begin
g:=-s*v(j+2,2);
v(j+2,2):=c*v(j+2,2)
end end j;
endk: end k;
for i:=1 step 1 until N0 do begin
a(i):=v(i,0); b(i):=v(i,1) end
end else begin
for i:=1 step 1 until N0 do begin
j:=j+2; a(i):=e(j,0); b(i):=e(j,2) end
end;
b(1):= 0;
i0:=j:= if atype then nr else N0+1-nr;
h:=(if j>1 then abs b(j) else 0) +
(if j<N0 then abs b(j+1) else 0);
if j=1 or j=N0 then h:=h+h;
g:=a(j)+h; h:=a(j)-h;
for i:=2 step 1 until N0 do b(i):=b(i)*b(i);
for lambda:=(g+h)*0.5 while g<>lambda do begin
p:=0; q:=1; it:=0;
for i:=1 step 1 until N0 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=nr then i:=N0
end;
if it=nr then h:=lambda else g:=lambda
end bisection;
dE0:=dE0+k1*lambda;
if Rot(7)=0 and couple and atype and -,oblat then begin
real procedure fac(m,n);
value m,n; integer m,n;
begin real f; integer i;
f:=1;
for i:=n+1 step 1 until m do f:=f*i;
fac:=f;
end fac;
k:= AUXMAT(type);
k2:= if k1=1 then k21 else k11;
if k2>0 then begin
g:=H(0)*3-G(0)/2; y:=(-H(0)-G(0)/2)/g;
h:=fac(J0+k2,J0-k2)/((fac(k2-1,1)**2)*8**(k2-1))*g*y**k2;
p:= k2*Rot(8); p:= p*p;
p:= p/(sqrt(p+h*h/4)+h/2);
lambda:= lambda+(if k=1 or k=4 then p else -p);
end end;
i0:=i2+i0+i0; g:=e(i0,0);
for i:=1,2,3,4 do begin
d(-i):= if i0-i>0 then e(i0,i) else 0;
d(i-1):= if i0+i<=N then e(i0+i,i) else 0
end i;
for i:=i0+1 step 1 until N do
for j:=0 step 1 until 4 do
e(i-1,j):= if j<i-i0 then e(i,j) else
if j<4 then e(i,j+1) else 0;
it:=0; N:=N-1; y:=lambda;
repeat:
for i:=N step -1 until 1 do b(i):=0;
for i:=-4 step 1 until 3 do
if d(i)<>0 then b(i0+i):=-d(i);
for i:=N step -1 until 1 do begin
for j:=0,j+1 while i-j>0 and j<=4 do begin
p:=e(i,j); k:=i;
for k:=k+1 while k<=N and k-i+j<=4 do
p:=p-a(k)*u(k,k-i)*u(k,k-i+j);
if j=0 then a(i):=p-lambda
else u(i,j):=p/a(i)
end j;
p:=b(i); k:=i;
for k:=k+1 while k<=N and k-i<=4 do
p:=p-u(k,k-i)*b(k);
b(i):=p
end i;
q:=1;
for i:=1 step 1 until N do begin
p:=if b(i)=0 then 0 else b(i)/a(i); k:=i;
for k:=k-1 while k>0 and i-k<=4 do
p:=p-u(i,i-k)*b(k);
b(i):=p; q:=q+p*p
end i;
if N>0 then begin
h:=g-lambda;
for i:=-4 step 1 until 3 do
if d(i)<>0 then h:=h+d(i)*b(i0+i);
h:=h/q;
if it<2 or (abs h<abs y and lambda+h<>lambda) then begin
lambda:=lambda+h; it:=it+1; y:=h; goto repeat
end
end;
q:=1/sqrt(q);
for i:=N step -1 until i0 do b(i+1):=b(i)*q;
for i:=i0-1 step -1 until 1 do b(i):=b(i)*q;
b(i0):=q; N:=N+1
end;
?
l./u(nr,i)/,r/u(nr,i)/b(i)/
l./u(nr,1)/,g/u(nr,/b(/
l./p:=ai(8/,r/p:=/p:=Li-dE-dE0;
/,d,f
i gtxt
;goasm2
;gosav time.300 lines.5000
gtxt=set 50
permanent gtxt.15
gtxt=edit algasmixx c.50
l1,r/91/110/,
l2,r/no/yes/,
l./Rot/,r/9/10/,l./MATRIX/,r/;/, AUXMAT(1:8);/
l./if oblat/,l3,i/
for i:=1 step 1 until 8 do
AUXMAT(i):= case i of (1,2,3,4,2,1,4,3);
/,l./couple/,r/; /;
/,r/;/ or Rot(10)<>0;/,l./type coup/,l1,r/9/9,10/,r/_ddd_d/dd/,
l./integer line/,r/J0/J0,N0/
l1,r/dE/dE,dE0,dP3/,l./goto nylinie/,
l3,r/d/dE0:=dP3:=d/
l./nr:=/,g2?//2?//2+1?
l./nr:=(J0+j)/,r/nr/N0/,r/-nr//
l./array u/,r/u(1:N,1:N)/e(1:N,0:4)/,r/ev/b/,l1,d2
l./u(i,i)/,r/u(i,i)/e(i,0)/,r/;/; e(i,3):= 0;/
l4,r/u(j,i)/e(j,1)/
l5,r/u(j,i)/e(j,2)/
l3,r/u(i,i-4)/e(i,4)/,l4,g/u(2,2)/e(2,0)/
l3,g/u(4,2)/e(4,2)/,l1,g/u(3,3)/e(3,0)/
l3,g/u(1,1)/e(1,0)/,l2,g/u(3,1)/e(3,2)/,l1,g/u(2,2)/e(2,0)/
l2,d./lambda/,l-1,r/;/; k1:=k;/,l1,i?
begin
integer it,i0,k; real g,h,y,p,q;
array a(1:N),d(-4:3),u(1:N,1:4);
j:=i2:= if i1=vk then 0 else -1;
if N0>2 and DW(0)<>0 then begin
comment The pentadiagonal matrix of the uncoupled
state is reduced to tridiagonal form by succesive
Jacobi rotations, Swartz: Numer. Math. 12, 231 (1968),
modified by replacing a(i,j) and b by v(i+j,j) and p;
real c,s,c2,s2,cs,u,u1;
array v(1:N0,0:2);
for i:=1 step 1 until N0 do begin
j:=j+2;
for k:=0,1,2 do v(i,k):=e(j,k+k)
end;
comment v(1,1):=v(1,2):=v(2,2):=0;
for k:=1 step 1 until N0-2 do begin
for j:=k+2 step 2 until N0 do begin
if j=k+2 then begin
if v(k+2,2)=0 then goto endk;
p:=-v(k+1,1)/v(k+2,2);
end else begin
if g=0 then goto endk;
p:=-v(j-1,2)/g;
end;
s2:=1/(p*p+1); s:=sqrt(s2); c:=p*s;
c2:=c*c; cs:=c*s;
u:=c2*v(j-1,0)-2*cs*v(j,1)+s2*v(j,0);
u1:=s2*v(j-1,0)+2*cs*v(j,1)+c2*v(j,0);
v(j,1):=cs*(v(j-1,0)-v(j,0))+(c2-s2)*v(j,1);
v(j-1,0):=u; v(j,0):=u1;
u:=c*v(j-1,1)-s*v(j,2);
v(j,2):=s*v(j-1,1)+c*v(j,2);
v(j-1,1):=u;
if j<>k+2 then
v(j-1,2):=c*v(j-1,2)-s*g;
if j<N0 then begin
u:=c*v(j+1,2)-s*v(j+1,1);
v(j+1,1):=s*v(j+1,2)+c*v(j+1,1);
v(j+1,2):=u
end;
if j+2<=N0 then begin
g:=-s*v(j+2,2);
v(j+2,2):=c*v(j+2,2)
end end j;
endk: end k;
for i:=1 step 1 until N0 do begin
a(i):=v(i,0); b(i):=v(i,1) end
end else begin
for i:=1 step 1 until N0 do begin
j:=j+2; a(i):=e(j,0); b(i):=e(j,2) end
end;
b(1):= 0;
i0:=j:= if atype then nr else N0+1-nr;
h:=(if j>1 then abs b(j) else 0) +
(if j<N0 then abs b(j+1) else 0);
if j=1 or j=N0 then h:=h+h;
g:=a(j)+h; h:=a(j)-h;
for i:=2 step 1 until N0 do b(i):=b(i)*b(i);
for lambda:=(g+h)*0.5 while g<>lambda do begin
p:=0; q:=1; it:=0;
for i:=1 step 1 until N0 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=nr then i:=N0
end;
if it=nr then h:=lambda else g:=lambda
end bisection;
dE0:=dE0+k1*lambda;
if Rot(7)=0 and couple and atype and -,oblat then begin
real procedure fac(m,n);
value m,n; integer m,n;
begin real f; integer i;
f:=1;
for i:=n+1 step 1 until m do f:=f*i;
fac:=f;
end fac;
k:= AUXMAT(type);
k2:= if k1=1 then k21 else k11;
if k2>0 then begin
g:=H(0)*3-G(0)/2; y:=(-H(0)-G(0)/2)/g;
h:=fac(J0+k2,J0-k2)/((fac(k2-1,1)**2)*8**(k2-1))*g*y**k2;
h:= h/2; p:= k2*Rot(8); q:= k2**3*Rot(10);
k2:= if k=1 or k=4 then 1 else -1;
g:=p*p; y:= sqrt(g+h*h);
lambda:= lambda+g/(y+h)*k2; g:=p+q;
dP3:= dP3+(g+p)*q/(sqrt(g*g+h*h)+y)*k1*k2;
comment dP3 is the contribution from a perturbation G*Pa**3;
end end;
i0:=i2+i0+i0; g:=e(i0,0);
for i:=1,2,3,4 do begin
d(-i):= if i0-i>0 then e(i0,i) else 0;
d(i-1):= if i0+i<=N then e(i0+i,i) else 0
end i;
for i:=i0+1 step 1 until N do
for j:=0 step 1 until 4 do
e(i-1,j):= if j<i-i0 then e(i,j) else
if j<4 then e(i,j+1) else 0;
it:=0; N:=N-1; y:=lambda;
repeat:
for i:=N step -1 until 1 do b(i):=0;
for i:=-4 step 1 until 3 do
if d(i)<>0 then b(i0+i):=-d(i);
for i:=N step -1 until 1 do begin
for j:=0,j+1 while i-j>0 and j<=4 do begin
p:=e(i,j); k:=i;
for k:=k+1 while k<=N and k-i+j<=4 do
p:=p-a(k)*u(k,k-i)*u(k,k-i+j);
if j=0 then a(i):=p-lambda
else u(i,j):=p/a(i)
end j;
p:=b(i); k:=i;
for k:=k+1 while k<=N and k-i<=4 do
p:=p-u(k,k-i)*b(k);
b(i):=p
end i;
q:=1;
for i:=1 step 1 until N do begin
p:=if b(i)=0 then 0 else b(i)/a(i); k:=i;
for k:=k-1 while k>0 and i-k<=4 do
p:=p-u(i,i-k)*b(k);
b(i):=p; q:=q+p*p
end i;
if N>0 then begin
h:=g-lambda;
for i:=-4 step 1 until 3 do
if d(i)<>0 then h:=h+d(i)*b(i0+i);
h:=h/q;
if it<2 or (abs h<abs y and lambda+h<>lambda) then begin
lambda:=lambda+h; it:=it+1; y:=h; goto repeat
end
end;
q:=1/sqrt(q);
for i:=N step -1 until i0 do b(i+1):=b(i)*q;
for i:=i0-1 step -1 until 1 do b(i):=b(i)*q;
b(i0):=q; N:=N+1
end;
?
l./u(nr,i)/,r/u(nr,i)/b(i)/
l./u(nr,1)/,g/u(nr,/b(/
l./end k;/,l1,r/if/dE:= dE+dP3;
if/,l./p:=ai(8/,r/p:=/p:=Li-dE-dE0;
/,d,l./>,Li/,r/dE/dP3/,f
i gtxt
goasm3
;gosav time.900 lines.5000
gtxt=set 50
permanent gtxt.15
gtxt=edit algasmixx c.50
l1,r/91/110/,
l./MATRIX/,r/;/, AUXMAT(1:8);/
l./if oblat/,l3,i/
for i:=1 step 1 until 8 do
AUXMAT(i):= case i of (1,2,3,4,2,1,4,3);
/,l./integer line/,r/J0/J0,N0/
l1,r/dE/dE,dE0,dP3/,l./goto nylinie/,
l3,r/d/dE0:=dP3:=d/
l./nr:=/,g2?//2?//2+1?
l./nr:=(J0+j)/,r/nr/N0/,r/-nr//
l./array u/,r/u(1:N,1:N)/e(1:N,0:4)/,r/ev/b/,l1,d2
l./u(i,i)/,r/u(i,i)/e(i,0)/,r/;/; e(i,3):= 0;/
l4,r/u(j,i)/e(j,1)/
l5,r/u(j,i)/e(j,2)/
l3,r/u(i,i-4)/e(i,4)/,l4,g/u(2,2)/e(2,0)/
l3,g/u(4,2)/e(4,2)/,l1,g/u(3,3)/e(3,0)/
l3,g/u(1,1)/e(1,0)/,l2,g/u(3,1)/e(3,2)/,l1,g/u(2,2)/e(2,0)/
l2,d./lambda/,l-1,r/;/; k1:=k;/,l1,i?
begin
integer it,i0,k; real g,h,y,p,q;
array a(1:N),d(-4:3),u(1:N,1:4);
j:=i2:= if i1=vk then 0 else -1;
if N0>2 and DW(0)<>0 then begin
comment The pentadiagonal matrix of the uncoupled
state is reduced to tridiagonal form by succesive
Jacobi rotations, Swartz: Numer. Math. 12, 231 (1968),
modified by replacing a(i,j) and b by v(i+j,j) and p;
real c,s,c2,s2,cs,u,u1;
array v(1:N0,0:2);
for i:=1 step 1 until N0 do begin
j:=j+2;
for k:=0,1,2 do v(i,k):=e(j,k+k)
end;
comment v(1,1):=v(1,2):=v(2,2):=0;
for k:=1 step 1 until N0-2 do begin
for j:=k+2 step 2 until N0 do begin
if j=k+2 then begin
if v(k+2,2)=0 then goto endk;
p:=-v(k+1,1)/v(k+2,2);
end else begin
if g=0 then goto endk;
p:=-v(j-1,2)/g;
end;
s2:=1/(p*p+1); s:=sqrt(s2); c:=p*s;
c2:=c*c; cs:=c*s;
u:=c2*v(j-1,0)-2*cs*v(j,1)+s2*v(j,0);
u1:=s2*v(j-1,0)+2*cs*v(j,1)+c2*v(j,0);
v(j,1):=cs*(v(j-1,0)-v(j,0))+(c2-s2)*v(j,1);
v(j-1,0):=u; v(j,0):=u1;
u:=c*v(j-1,1)-s*v(j,2);
v(j,2):=s*v(j-1,1)+c*v(j,2);
v(j-1,1):=u;
if j<>k+2 then
v(j-1,2):=c*v(j-1,2)-s*g;
if j<N0 then begin
u:=c*v(j+1,2)-s*v(j+1,1);
v(j+1,1):=s*v(j+1,2)+c*v(j+1,1);
v(j+1,2):=u
end;
if j+2<=N0 then begin
g:=-s*v(j+2,2);
v(j+2,2):=c*v(j+2,2)
end end j;
endk: end k;
for i:=1 step 1 until N0 do begin
a(i):=v(i,0); b(i):=v(i,1) end
end else begin
for i:=1 step 1 until N0 do begin
j:=j+2; a(i):=e(j,0); b(i):=e(j,2) end
end;
b(1):= 0;
i0:=j:= if atype then nr else N0+1-nr;
h:=(if j>1 then abs b(j) else 0) +
(if j<N0 then abs b(j+1) else 0);
if j=1 or j=N0 then h:=h+h;
g:=a(j)+h; h:=a(j)-h;
for i:=2 step 1 until N0 do b(i):=b(i)*b(i);
for lambda:=(g+h)*0.5 while g<>lambda do begin
p:=0; q:=1; it:=0;
for i:=1 step 1 until N0 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=nr then i:=N0
end;
if it=nr then h:=lambda else g:=lambda
end bisection;
dE0:=dE0+k1*lambda;
if Rot(7)=0 and couple and atype and -,oblat then begin
real procedure fac(m,n);
value m,n; integer m,n;
begin real f; integer i;
f:=1;
for i:=n+1 step 1 until m do f:=f*i;
fac:=f;
end fac;
k:= AUXMAT(type);
k2:= if k1=1 then k21 else k11;
if k2>0 then begin
g:=H(0)*3-G(0)/2; y:=(-H(0)-G(0)/2)/g;
h:=fac(J0+k2,J0-k2)/((fac(k2-1,1)**2)*8**(k2-1))*g*y**k2;
h:= h/2; p:= k2*Rot(8);
q:= p*k2*k2*(-7.01667'-3+k2*k2*1.47696'-5);
k2:= if k=1 or k=4 then 1 else -1;
g:=p*p; y:= sqrt(g+h*h);
lambda:= lambda+g/(y+h)*k2; g:=p+q;
dP3:= dP3+(g+p)*q/(sqrt(g*g+h*h)+y)*k1*k2;
comment dP3 is the contribution from a perturbation G*Pa**3;
end end;
i0:=i2+i0+i0; g:=e(i0,0);
for i:=1,2,3,4 do begin
d(-i):= if i0-i>0 then e(i0,i) else 0;
d(i-1):= if i0+i<=N then e(i0+i,i) else 0
end i;
for i:=i0+1 step 1 until N do
for j:=0 step 1 until 4 do
e(i-1,j):= if j<i-i0 then e(i,j) else
if j<4 then e(i,j+1) else 0;
it:=0; N:=N-1; y:=lambda;
repeat:
for i:=N step -1 until 1 do b(i):=0;
for i:=-4 step 1 until 3 do
if d(i)<>0 then b(i0+i):=-d(i);
for i:=N step -1 until 1 do begin
for j:=0,j+1 while i-j>0 and j<=4 do begin
p:=e(i,j); k:=i;
for k:=k+1 while k<=N and k-i+j<=4 do
p:=p-a(k)*u(k,k-i)*u(k,k-i+j);
if j=0 then a(i):=p-lambda
else u(i,j):=p/a(i)
end j;
p:=b(i); k:=i;
for k:=k+1 while k<=N and k-i<=4 do
p:=p-u(k,k-i)*b(k);
b(i):=p
end i;
q:=1;
for i:=1 step 1 until N do begin
p:=if b(i)=0 then 0 else b(i)/a(i); k:=i;
for k:=k-1 while k>0 and i-k<=4 do
p:=p-u(i,i-k)*b(k);
b(i):=p; q:=q+p*p
end i;
if N>0 then begin
h:=g-lambda;
for i:=-4 step 1 until 3 do
if d(i)<>0 then h:=h+d(i)*b(i0+i);
h:=h/q;
if it<2 or (abs h<abs y and lambda+h<>lambda) then begin
lambda:=lambda+h; it:=it+1; y:=h; goto repeat
end
end;
q:=1/sqrt(q);
for i:=N step -1 until i0 do b(i+1):=b(i)*q;
for i:=i0-1 step -1 until 1 do b(i):=b(i)*q;
b(i0):=q; N:=N+1
end;
?
l./u(nr,i)/,r/u(nr,i)/b(i)/
l./u(nr,1)/,g/u(nr,/b(/
l./end k;/,l1,r/if/dE:= dE+dP3;
if/,l./p:=ai(8/,r/p:=/p:=Li-dE-dE0;
/,d,f
i gtxt
;gosav time.900 lines.10000
asmixx
<Ethyl vinyl ether, g.s. (E - comp., s = 81.75)>
18111a 10
16307.064653 2816.942907 2478.088434
0 0 0 0 -0.114246 0
41 1
0.477893 -2.21924 36.997862 0.0821998 0.488354
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
0 1a 34 7 27 33 8 26 2 22315.41
0 1a 34 7 27 33 8 25 3 22316.58
0 1a 34 7 28 33 8 26 3 22328.56
0 1a 34 7 28 33 8 25 2 22329.71
0 1a 32 7 25 31 8 24 2 33409.73
0 1a 32 7 25 31 8 23 3 33410.84
0 1a 32 7 26 31 8 24 3 33415.56
0 1a 32 7 26 31 8 23 2 33416.62
0 1a 31 7 24 30 8 23 2 38924.46
0 1a 31 7 24 30 8 22 3 38925.58
0 1a 31 7 25 30 8 23 3 38928.31
0 1a 31 7 25 30 8 22 2 38929.42
0 1a 36 8 28 35 9 27 2 39370.82
0 1a 36 8 28 35 9 26 3 39371.87
0 1a 36 8 29 35 9 27 3 39372.42
0 1a 36 8 29 35 9 26 2 39373.42
0 1a 37 8 29 36 9 28 2 33847.52
0 1a 37 8 29 36 9 27 3 33848.56
0 1a 37 8 30 36 9 28 3 33849.64
0 1a 37 8 30 36 9 27 2 33850.70
0 1a 39 8 31 38 9 30 2 22743.70
0 1a 39 8 31 38 9 29 3 22744.77
0 1a 39 8 32 38 9 30 3 22748.04
0 1a 39 8 32 38 9 29 2 22749.11
0 1a 44 9 35 43 10 34 2 23172.96
0 1a 44 9 36 43 10 34 3 23174.65
0 1a 44 9 36 43 10 33 2 23175.63
0 0 30 8 23 30 8 22 1 1.115
0 0 31 7 25 31 7 24 1 3.845
0 0 31 8 24 31 8 23 1 1.085
0 0 32 7 26 32 7 25 1 5.805
0 0 33 8 26 33 8 25 1 1.16
0 0 34 7 28 34 7 27 1 13.14
0 0 35 9 27 35 9 26 1 1.00
0 0 36 8 29 36 8 28 1 1.55
0 0 36 9 28 36 9 27 1 1.06
0 0 37 8 30 37 8 29 1 2.14
0 0 38 9 30 38 9 29 1 1.07
0 0 39 8 32 39 8 31 1 4.34
0 0 43 10 34 43 10 33 1 0.98
0 0 44 9 36 44 9 35 1 1.69
0 1a 44 9 35 43 10 33 3 0
-1
;
goasmtest
;kemlab5 1 time.300
gtxt2=set 50
permanent gtxt2.15
gtxt2=edit gtxt
l./goto nylinie/,r/goto nylinie/j:=j; s:=0/
l./d(-4:3)/,l1,i/
procedure writemat(res,e,n0,n,k);
value n0,n,k; integer n0,n,k; array e; zone res;
begin
real a; integer i,j;
for i:=1 step 1 until n do begin
write(res,<:<10>:>);
for j:=0,j+1 while j<i and j<=k do begin a:=e(n0+i,j);
write(res,if abs(a)<'6 then << -dddddd.ddd>
else << -d.dddddd'd>,a);
end end;
write(res,<:<10>:>);
end;
writemat(res,e,0,N,4);
/,l./repeat1/,i/
writemat(res,u,N0,N0,2);
/,l./q:=1;/,i/
write(out,<:<10>:>,<<ddddddd.ddd>,lambda,<< -d.ddddd'-d>,y,<:<10>:>);
for i:=1 step 1 until N0 do begin
for j:=0 step 1 until 3 do begin
q:=if j=0 then a(i) else if j=3 then b(i) else u(i,j);
if j>=i and j<3 then write(out,<: :>) else
write(out,if abs q<'6 then << -dddddd.dd>
else << -d.ddddd'd>,q);
end;
write(out,<:<10>:>);
end;
setposition(out,0,0);
/,l./i2+i0+i0/,i/
s:=s+k1*lambda;
/,l./repeat2/,i/
writemat(res,e,0,N,4);
/,l./q:=1;/,i/
write(out,<:<10>:>,<<ddddddd.ddd>,lambda,<< -d.ddddd'-d>,y,<:<10>:>);
for i:=1 step 1 until N do begin
for j:=0 step 1 until 5 do begin
q:=if j=0 then a(i) else if j=5 then b(i) else u(i,j);
if j>=i and j<5 then write(out,<: :>) else
write(out,if abs q<'6 then << -dddddd.dd>
else << -d.dddd'dd>,q);
end;
write(out,<:<10>:>);
end;
setposition(out,0,0);
/,l./>,Li/,r/d d/ddd/,r/Li/Li,<: :>,s);
/,d,f
i gtxt2
▶EOF◀