|
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: 36096 (0x8d00)
Types: TextFile
Names: »algflinda«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦7e928b248⟧ »algbib«
└─⟦this⟧
;gosav time.180
lookup flinda
if ok.no
(lookup xbmatrix
if ok.no
(r=algol extxbmatrix index.no
rename r.xbmatrix
permanent xbmatrix.17)
flinda=set 180
permanent flinda.17
flinda=algol index.no list.no xref.no)
\f
F-LINDA. 2-1-1975. GOS.
begin
comment The input rules are:
1) A text string in <> (max 59 characters) to appear
in top of each output page.
2) The allowed computing time (min.), and a number m:
If m=1 all inputs are required,
if m=0 only input 17 are performed, and
if m=-1 valence force data (15-17) are performed.
3) The number of molecules, NM, and for each of these
the following inputs 4-12. The individual molecules
are numbered by letters: a, b, c . . .
4) The name of the molecule in <>.
5) A number m : If m=1 centrifugal distortion constants
can be included in the calculations.
If m = 10: Angular internal coordinates are defined in
Angstrom units by multiplication with a function depen-
ding on the equilibrium bondlengths (see VIBROT spec.).
6) The number of isotopic species, G.
7) The number of atoms in the molecule, N.
8) The number of symmetry species, BL.
9) For each of these the number of symmetry coordinates
including redundants.
10) The coordinates of the N atoms (see VIBROT spec.).
11) Specifications for the internal coordinates and the sym-
metry coordinates (see VIBROT spec.).
12) For each symmetry species:
a) The substitution in <> (max 17 characters),
b) N atomic masses ordered in accordance with the
B-matrix,
c) The total number of normal frequencies, NV. If only
centrifugal distortion data are available put NV=0.
d) NV observed frequencies ordered in symmetry species
and decreasing. Following every frequency is given
an uncertainty, dnu . If dnu = 0 the appropiate fre-
quency is neglected in the calculations.
if m=1 :
e) The number of observed centrifugal distortion const.
f) For each of these the nine coefficients defining
the observed constant in terms of the basic constants
tau-vxyz, ordered as the index sequence:
aaaa, bbbb, cccc, abab, bcbc, caca, aabb, bbcc, ccaa.
g) The observed values, each followed by an uncertainty.
13) The number of shifts (i.e. isotopic differences), and
for each of these:
Two integers indicating the two observetions to be
subtracted, followed by the uncertainty of the shift.
14) The number of the last performed iteration, IT.
15) A maximal change in F, Xmax.
16) The valence force data:
a) The number of force constants (NKT) including fixed.
b) For each of these:
The number of depending symmetry force constants
(Fij). If the force constant should be fixed
the negative number is given. Next the start
value is given and for each symmetry force
constant is specified:
c) symmetry species no. (delimiter mol no., see 3),
i , j , factor. ( i=<j ).
17) A text string in <> (may be empty), giving the name of
a jobfile to be submitted to TABASCO;
integer NM,IT,NO,NF,NK,NKF,NS,On,Ol,i,j,k,l,m,n,space,page;
real Xmax, s, so, h, starttime, stoptime, time;
boolean ha, hb, newinp, closeres, centrifugal;
array head(1:10);
zone res(128,1,stderror);
integer array tail(1:10);
integer procedure split;
begin split:= h shift (-36) extract 12; h:= h shift 12 end;
space:= 30-readhead(in,head,1); read(in,stoptime,m);
newinp:= m>0;
for i:=2 step 1 until 10 do tail(i):= 0;
s:= 0; i:= k:= 0; j:= 48;
for j:=j+8 while i<4 do begin
if j>48 then begin
k:= k+1; j:= 8; h:= head(k)
end;
n:= h shift (-48+j) extract 8;
if n<91 then n:= n+32;
if n<=122 and n>=97 then begin
i:= i+1; s:= s add n shift 8
end end;
if newinp then begin read(in,NM); page:= 0 end else begin
s:= s add 49 shift 8;
open(res,4,string s,0); monitor(42,res,0,tail);
setposition(res,0,tail(1)-1); inrec(res,128);
h:= res(128);
NM:= split; IT:= split; NO:= split; On:= split;
page:= NM extract 9; NM:= NM shift (-9);
h:= res(127); Xmax:= res(126); so:= res(125);
NF:= split; NK:= split; NKF:= split; NS:= split;
s:= s shift (-16) shift 8; close(res,true);
end;
begin
integer linecount,date,G,Na,N3,BL,NV,NC,NT,
NKT,ig,ibl,red,inm;
array mol(1:5), track(1:(NM+1)*4);
procedure cr;
begin linecount:= linecount+1; write(res,<:<10>:>) end;
procedure writepage;
begin integer i,j; real lo; boolean nl,sp;
nl:= false add 10; sp:= false add 32;
page:= page+1; linecount:= 3; i:= (date//100) mod 100; j:= 1;
lo:= if i<=9 then real(<<-d>) else real(<<-dd>);
write(res,nl,1,<:<12>:>,nl,3,string head(increase(j)),sp,space
+(if i<=9 then 1 else 0),<<dd>,date//10000,string lo,
-i,<<-dd>,-(date mod 100),sp,19,<<ddd>,page,nl,2)
end writepage;
procedure set(z,a,tail);
value a; real a; zone z; integer array tail;
begin
monitor(48,z,0,tail);
if monitor(40,z,0,tail)>0 then begin
write(out,<:
***:>,string a,<: not created.:>); goto stop
end;
monitor(50,z,17,tail);
end set;
systime(1,0,starttime);
date:= round systime(2,starttime,time); stoptime:= stoptime*60;
closeres:= outmedium(res);
for i:=1 step 1 until 4 do
track(i):= s add (48+i) shift 8;
h:= s shift (-16) shift 8;
for i:=1 step 1 until 4 do begin
n:= case i of(109,102,118,100);
for j:=1 step 1 until NM do
track(4*j+i):= h add (48+j) shift 8 add n shift 8
end;
if newinp then begin
zone z, Obs(128,1,stderror);
integer procedure Ostep;
begin On:= On+1;
if On=129 then begin
On:= 1; Ol:= Ol+1; outrec(Obs,128) end;
Ostep:= On
end;
procedure pack(n);
value n; integer n;
begin z(m):= z(m) shift 12 add n end;
open(Obs,4,string track(3),0);
tail(1):= 5; set(Obs,track(3),tail);
On:= 128; Ol:= NO:= 0;
for inm:=1 step 1 until NM do begin
comment Input 4-12 for each molecule inm;
readhead(in,mol,1);
read(in,m,G,Na,BL);
centrifugal:= m mod 10= 1; hb:= m//10= 1; N3:= abs Na*3;
if centrifugal then begin
ig:= G; ibl:= BL end else ig:= ibl:= 1;
begin
integer array S(1:BL,1:2), Vpl(1:BL,1:G), Dpl(1:ibl,1:ig);
array species(1:G*3);
procedure writespecies;
begin integer i,j;
i:= (ig-1)*3 + 1; j:= 1;
write(res,<: of :>,string species(increase(i)),
string mol(increase(j)))
end writespecies;
for ibl:=1 step 1 until BL do read(in,S(ibl,1)); writepage;
i:= 1; write(res,string mol(increase(i))); cr; cr;
if xbmatrix(Na,BL,S,hb,ha,res) <> 0 then goto stop;
NT:= if -,centrifugal then 1 else if ha then 4 else 6;
comment Disctracks are dimensioned for the
matrices F, V and D;
Vpl(1,1):= 0;
for ig:=1 step 1 until G do
for ibl:=1 step 1 until BL do begin
n:= S(ibl,1); m:= ((n+1)*n+254)//256;
if ibl<BL then Vpl(ibl+1,ig):= Vpl(ibl,ig) + m else
if ig<G then Vpl( 1,ig+1):= Vpl(BL,ig) + m
end;
tail(1):= Vpl(BL,G)+m;
for m:=2,3 do begin
h:= track(4*inm+m); open(z,4,string h,0);
set(z,h,tail); close(z,true);
end;
if centrifugal then begin
Dpl(1,1):= 0;
for ig:=1 step 1 until G do
for ibl:=1 step 1 until BL do begin
n:= S(ibl,1); m:= (n*NT+127)//128;
if ibl<BL then Dpl(ibl+1,ig):= Dpl(ibl,ig) + m else
if ig<G then Dpl( 1,ig+1):= Dpl(BL,ig) + m
end;
open(z,4,string track(4*inm+4),0); tail(1):= Dpl(BL,G)+m;
set(z,track(4*inm+4),tail); close(z,true);
NC:= N3
end else NC:= 1;
begin
array X(1:Na,1:3);
open(z,4,<:xtrack:>,0); inrec(z,128); k:= 0;
for i:=1 step 1 until Na do
for j:=1 step 1 until 3 do begin
k:= k+1; X(i,j):= z(k) end;
monitor(48,z,0,tail); close(z,true);
hb:= false; NV:= 0;
for ig:=1 step 1 until G do begin
real Ia, Ib, Ic;
array M(1:Na), UX(1:NC,1:3), U(1:3,1:3);
readhead(in,species,(ig-1)*3+1);
if atomic(M,Na) then goto stop;
write(res,false add 10,3,<:Atomic masses:>);
writespecies;
for i:=1 step 1 until Na do begin
if (i-1) mod 6 = 0 then cr;
write(res,<<dddd.dddddd>,M(i))
end;
read(in,n);
NV:= NV+n; Obs(Ostep):= n; n:= n+n;
for i:=1 step 1 until n do read(in,Obs(Ostep));
if centrifugal then begin
array IP(1:3);
read(in,n); Obs(Ostep):= n; hb:= n<>0;
if hb then begin
for i:=1 step 1 until n do
for j:=1 step 1 until 9 do read(in,Obs(Ostep));
NO:= NO+n; n:= n+n;
for i:=1 step 1 until n do read(in,Obs(Ostep));
comment This completes the input. The tensor of inertia is cal-
culated and a transformation to the principal axes system
is performed;
Ib:= 1/sum(M(i),i,1,Na);
for j:=1 step 1 until 3 do begin
Ia:= sum(X(i,j)*M(i),i,1,Na)*Ib;
for i:=1 step 1 until Na do X(i,j):= X(i,j)-Ia
end j, center of mass corr.;
Ia:= sum(M(i)*X(i,1)**2,i,1,Na);
Ib:= sum(M(i)*X(i,2)**2,i,1,Na);
Ic:= sum(M(i)*X(i,3)**2,i,1,Na);
U(1,1):= Ib+Ic; U(2,2):= Ic+Ia; U(3,3):= Ia+Ib;
U(2,1):= -sum(M(i)*X(i,1)*X(i,2),i,1,Na);
U(3,2):= -sum(M(i)*X(i,2)*X(i,3),i,1,Na);
U(3,1):= -sum(M(i)*X(i,3)*X(i,1),i,1,Na);
tridql(3,IP,U); s:= 6.0657791;
Ia:= s/IP(1); Ib:= s/IP(2); Ic:= s/IP(3);
for i:=1 step 1 until Na do
for j:=1 step 1 until 3 do UX(i,j):= sum(U(j,k)*X(i,k),k,1,3)
end hb end centrifugal;
for i:=1 step 1 until Na do M(i):= 1/M(i); NF:= 0;
for ibl:=1 step 1 until BL do begin
n:= S(ibl,1); l:= n*(n+1)//2; k:= (l+127)//128*128;
comment The matrix product G = B * M * Btr is calculated and
stored in the linear array V, representing a lower tri-
angular matrix;
begin
integer red,ri,rj,rk,rl; real W,VV;
array K(1:l), B(1:n,1:N3); zone V(k,1,stderror);
open(z,4,<:btrack:>,0); setposition(z,0,NF); k:= 128;
for i:=1 step 1 until n do
for j:=1 step 1 until N3 do begin
k:= k+1; if k=129 then begin
k:= 1; NF:= NF+1; inrec(z,128) end;
B(i,j):= z(k)
end;
if ig=G and ibl=BL then monitor(48,z,0,tail); close(z,true);
open(V,4,string track(4*inm+3),0); setposition(V,0,Vpl(ibl,ig));
outrec(V,l); ri:= 0;
for i:=1 step 1 until n do
for j:=1 step 1 until i do begin
ri:= ri+1;
V(ri):= sum(B(i,k)*B(j,k)*M((k+2)//3),k,1,N3) end;
comment Now the triangular matrix V ( V * Vtr = G ) is
calculated in V;
red:= rl:= 0; rk:= n*(n-1)//2;
for j:=1 step 1 until n do K(rk+j):= 1;
for j:=1 step 1 until n-red do begin l:= j+red; ri:= rl;
for i:=l step 1 until n do begin
W:= V(ri+l);
for k:=j-1 step -1 until 1 do
W:= W- V(rl+k)*V(ri+k);
if i=l then begin
if W<'-6 then begin red:= red+1; rl:= rl+l; l:= j+red; rk:= ri;
for k:=i step 1 until n do
begin K(rk+i):= 0; rk:= rk+k end end
else begin W:= V(ri+j):= sqrt(W);
VV:= K(ri+i):= 1/W end
end else V(ri+j):= W*VV;
ri:= ri+i
end i;
ri:= rl;
for i:=l-1 step -1 until j do begin ri:= ri-i; V(ri+j):= 0 end;
rl:= rl+l
end j;
comment red is the number of redundancies. The matrix V occupies
the first n-red columns of the array V;
if hb then begin
l:= red; rl:= n*(n-1)//2; rj:= rl+n;
for j:=n step -1 until 1 do if K(rl+j)<>0 then begin
VV:= K(rj); ri:= rj:= rj-j; k:= j-1+red-l;
for i:=j step 1 until k do begin K(ri+j):= 0; ri:= ri+i end;
if i<>j then K(ri+j):= VV; ri:= ri+i;
for i:=k+2 step 1 until n do begin
W:= 0; rk:= rj+j-l;
for k:=j+1 step 1 until i do begin
W:= W- K(ri+k)*V(rk+j); rk:= rk+k end;
K(ri+j):= W*VV; ri:= ri+i
end i
end else begin rj:= rj-j; l:=l-1 end;
comment The (left) inverse of V has been calculated in the last
n-red rows of the array K, and V is transferred to the
drum;
end hb;
close(V,true); red:= n-red;
if ig=1 then S(ibl,2):= red else if S(ibl,2)<>red then begin
write(res,<:
Error in determination of redundancies.
:>); goto stop
end;
comment If necessary the Jab,s matrices are calculated in D
and stored on the drum;
if hb then begin
array D(1:n,1:NT), BU(1:N3), I(1:if red<NT then NT else red);
open(z,4,string track(4*inm+4),0); setposition(z,0,Dpl(ibl,ig));
for j:=1 step 1 until NT do
I(j):= case j of (Ia*Ia, Ib*Ib, Ic*Ic, Ia*Ib, Ib*Ic, Ic*Ia);
for i:=1 step 1 until n do begin
for k:=N3-3 step -3 until 0 do
for j:=1 step 1 until 3 do
BU(k+j):= sum(B(i,k+l)*U(j,l),l,1,3);
for j:=1 step 1 until 3 do begin
rk:= j mod 3+ 1; rl:= rk mod 3+ 1; W:= 0; ri:= 0;
for k:=1 step 1 until Na do begin
W:= W + BU(ri+rk)*UX(k,rk) + BU(ri+rl)*UX(k,rl);
ri:= ri+3 end;
D(i,j):= W*I(j)
end j;
for j:=4 step 1 until NT do begin
rk:= (j-1) mod 3+ 1; rl:= rk mod 3+ 1; W:= 0; ri:= 0;
for k:=1 step 1 until Na do begin
W:= W- BU(ri+rk)*UX(k,rl); ri:= ri+3 end;
D(i,j):= W*I(j)
end j
end i;
rj:= n-red; rl:= rj*(rj+1)//2;
for j:=1 step 1 until NT do begin
l:= rj; ri:= rl;
for i:=1 step 1 until red do begin
l:= l+1; I(i):= sum(K(ri+k)*D(k,j),k,1,l); ri:= ri+l end;
for i:=1 step 1 until red do D(i,j):= I(i)
end j;
k:= 128;
for i:=1 step 1 until red do
for j:=1 step 1 until NT do begin
k:= k+1; if k=129 then begin
k:= (red-i)*NT+NT-j+1;
outrec(z,if k>128 then 128 else k); k:= 1 end;
z(k):= D(i,j)
end;
close(z,true)
end hb end K, B, V;
end ibl end ig end X;
n:= 0; for ibl:=1 step 1 until BL do n:= n+S(ibl,2);
if NV//n*n<>NV then begin
write(res,<:
Error in number of observed frequencies. inm, n, NV = :>,
<<-ddd>,inm, n,NV); goto stop
end;
NO:= NO+NV;
comment moldata are stored on disc;
open(z,4,string track(4*inm+1),0);
tail(1):= 1; set(z,track(4*inm+1),tail); outrec(z,128);
m:= 1; z(m):= 0; pack( G); pack(BL); pack(NT);
z(m):= z(m) shift 12;
for i:=1 step 1 until 5 do z(1+i):= mol(i);
for i:=G*3 step -1 until 1 do z(6+i):= species(i);
m:= 6+G*3; n:= 2;
for ibl:=1 step 1 until BL do begin
if n=2 then begin m:= m+1; n:= 0; z(m):= 0 end;
n:= n+1; pack(S(ibl,1)); pack(S(ibl,2))
end;
z(m):= z(m) shift ((2-n)*24); n:= 4;
for ibl:=1 step 1 until BL do
for ig:=1 step 1 until G do begin
if n=4 then begin m:= m+1; n:= 0; z(m):= 0 end;
n:= n+1; pack(Vpl(ibl,ig))
end;
z(m):= z(m) shift ((4-n)*12);
if centrifugal then begin n:= 4;
for ibl:=1 step 1 until BL do
for ig:=1 step 1 until G do begin
if n=4 then begin m:= m+1; n:= 0; z(m):= 0 end;
n:= n+1; pack(Dpl(ibl,ig))
end;
z(m):= z(m) shift ((4-n)*12) end;
close(z,true)
end S, Vpl, Dpl
end inm;
read(in,NS); Obs(Ostep):= NS;
for i:=1 step 1 until NS do begin
read(in,k,l,s);
Obs(Ostep):= 0.0 shift 14 add k shift 24 add l;
Obs(Ostep):= s
end;
NO:= NO+NS;
tail(1):= Ol; monitor(44,Obs,0,tail); close(Obs,true);
On:= (Ol-1)*128+On; so:= 0;
read(in,IT)
end else begin newinp:= m=-1; NKT:= NK+NKF end;
comment The precalculations have been completed by the transfer
of the V- and D-matrices to the drum, and the final
inputs (15-16) are performed;
if newinp then begin
read(in,Xmax,NKT);
begin
array fv(1:NKT*3);
zone f(128,1,stderror);
open(f,4,string track(1),0); tail(1):= 10; set(f,track(1),tail);
if m=-1 then writepage;
write(res,false add 10,3,<:
VALENCE FORCE DATA:
:>); k:= 127; NF:= NKF:= 0;
for n:=1 step 1 until NKT do begin
read(in,m); repeatchar(in); i:= 32;
for i:=i while i=32 do readchar(in,i);
repeatchar(in); j:= NKT+n*2-1;
if i=60 then begin
i:= readhead(in,fv,j);
s:= fv(j); so:= fv(j+1);
if i>11 then i:= 11
end else i:= 0;
if i<6 then s := s shift (8*i-48)
else so:= so shift (8*i-96);
for i:=i+1 step 1 until 6 do s := s shift 8 add 32;
for i:=i step 1 until 11 do so:= so shift 8 add 32;
fv(j):= s; fv(j+1):= so shift 8;
read(in,fv(n)); write(res,<:
Param.:>,<<dd>,n,<:: :>,string fv(increase(j)),
<<__-dd.dddddd>,fv(n));
hb:= m<0; m:= abs m; NF:= NF+m+m;
s:= 0.0 shift 12 add (if hb then 2048+m else m);
if hb then begin
write(res,<: (fixed):>); NKF:= NKF+1 end;
for l:=1 step 1 until m do begin
if l mod 3 = 1 then write(res,<:<10> :>);
k:= k+2; if k=129 then begin
k:= 1; outrec(f,128) end;
read(in,ibl);
if NM>1 then begin
repeatchar(in); readchar(in,inm)
end else inm:= 97;
read(in,i,j,f(k+1));
repeatchar(in); readchar(in,ig);
if ig=47 then begin
read(in,so); f(k+1):= f(k+1)/so;
repeatchar(in); readchar(in,ig)
end;
if ig=115 then f(k+1):= sign(f(k+1))*sqrt(abs f(k+1));
write(res,<<ddd>,ibl,false add inm,1,i,j,<< -d.dddd>,f(k+1));
inm:= inm-96; if inm<1 or inm>NM then goto stop;
f(k):= s shift 12 add j shift 12 add i shift 6
add inm shift 6 add ibl;
s:= 0
end end n;
NK:= NKT-NKF;
tail(1):= (NF+131)//128;
monitor(44,f,0,tail);
k:= 129; l:= NKT*3;
close(f,true); open(f,4,string track(2),0);
tail(1):= (l+NK+127)//128; set(f,track(2),tail);
for i:=1 step 1 until l do begin
if k=129 then begin k:= 1; outrec(f,128) end;
f(k):= fv(i); k:= k+1
end;
close(f,true) end end kdata;
begin
integer Oc, p, q;
real nu, dnu, H, DF;
array Obs(1:On), fv(1:NKT), B(1:NK), K(1:NF);
zone z(128,1,stderror);
integer procedure kstep;
begin k:= k+1; if k>128 then begin
k:= k mod 128; inrec(z,128) end;
kstep:= k
end kstep;
open(z,4,string track(1),0); k:= 128;
for i:=1 step 1 until NF do K(i):= z(kstep);
close(z,true);
open(z,4,string track(2),0); k:= 128;
for i:=1 step 1 until NKT do fv(i):= z(kstep);
if -, newinp then begin
k:= k+NKT+NKT;
for i:= 1 step 1 until NK do B(i):= z(kstep)
end else so:= 0;
close(z,true);
open(z,4,string track(3),0); k:= 128;
for i:=1 step 1 until On do Obs(i):= z(kstep);
close(z,true);
open(z,4,string track(4),0);
if newinp then begin
tail(1):= (NO-1)//(128//(NK+1))+1; set(z,track(4),tail)
end else monitor(42,z,0,tail);
close(z,true); open(z,4,<:atrack:>,0);
set(z,real<:atrack:>,tail); close(z,true);
systime(1,starttime,time); ha:= false;
Ol:= On-NS-NS;
igen:
if NS>0 then begin p:= NS; q:= NK+1 end else p:= q:= 1;
begin
array ashifts(1:p,1:q);
for i:=1 step 1 until p do
for j:=1 step 1 until q do ashifts(i,j):= 0;
writepage; s:= 0; NO:= Na:= On:= 0; IT:= IT+1;
for inm:=1 step 1 until NM do begin
open(z,4,string track(4*inm+1),0); inrec(z,128);
h:= z(1); G:= split; BL:= split; NT:= split;
for i:=1 step 1 until 5 do mol(i):= z(1+i);
centrifugal:= NT>1;
if centrifugal then begin
ig:= G; ibl:= BL end else ig:= ibl:= 1;
begin
integer array S(1:BL,1:2), Vpl(1:BL,1:G), Dpl(1:ibl,1:ig);
array species(1:G*3);
procedure writespecies;
begin integer i,j;
i:= (ig-1)*3 + 1; j:= 1;
write(res,<: of :>,string species(increase(i)),
string mol(increase(j)))
end writespecies;
for i:=G*3 step -1 until 1 do species(i):= z(6+i);
m:= 6+G*3; n:= 2;
for ibl:=1 step 1 until BL do begin
if n=2 then begin m:= m+1; n:= 0; h:= z(m) end;
n:= n+1; S(ibl,1):= split; S(ibl,2):= split
end;
n:= 4;
for ibl:=1 step 1 until BL do
for ig:=1 step 1 until G do begin
if n=4 then begin m:= m+1; n:= 0; h:= z(m) end;
n:= n+1; Vpl(ibl,ig):= split
end;
if centrifugal then begin n:= 4;
for ibl:=1 step 1 until BL do
for ig:=1 step 1 until G do begin
if n=4 then begin m:= m+1; n:= 0; h:= z(m) end;
n:= n+1; Dpl(ibl,ig):= split
end end;
close(z,true);
for ig:=1 step 1 until G do begin
On:= On+1; NV:= Obs(On);
if centrifugal then begin
Oc:= NV*2+On+1; NC:= Obs(Oc) end else NC:= 0;
hb:= NC<>0;
if linecount+NV+NC>55 then writepage else begin cr;cr;cr end;
if NV<>0 then begin
write(res,<: -1:>); cr;
write(res,<:Frequencies in cm :>);
if hb then write(res,<: and :>)
end;
if hb then begin
write(res,<:CD-constants in MHz:>);
l:= 9; m:= NK; n:= NC
end else l:= m:= n:= 1;
writespecies; cr; cr;
write(res,
<: Obs. Calc. Obs.-calc. delta:>);
begin
array C(1:n,1:l), t(1:l), acd(1:n,1:m);
if hb then begin
for i:=1 step 1 until NC do begin
for j:=1 step 1 until 9 do begin
Oc:= Oc+1; C(i,j):= Obs(Oc) end;
for j:=1 step 1 until NK do acd(i,j):= 0 end;
for j:=1 step 1 until 9 do t(j):= 0
end;
for ibl:=1 step 1 until BL do begin
n:= S(ibl,1); red:= S(ibl,2); l:= (n+1)*n//2;
k:= if -,hb then 1 else n;
begin
integer ri, rj, rk, rl, inmbl;
array V(1:l), L(1:n,1:red), f(1:red), D(1:k,1:NT);
open(z,4,string track(4*inm+3),0);
setposition(z,0,Vpl(ibl,ig)); k:= 128;
for i:=1 step 1 until l do begin
k:= k+1; if k=129 then begin
k:= l-i+1; inrec(z,if k>128 then 128 else k); k:= 1 end;
V(i):= z(k)
end; close(z,true);
comment The matrix V has been transferred from the drum, and
the matrix product Vtr * F * V is formed in L;
begin
real W; array F(1:l);
for i:=1 step 1 until l do F(i):= 0;
k:= 1; inmbl:= inm shift 6 add ibl;
for rl:=1 step 1 until NKT do begin
m:= ((K(k) shift (-36) extract 12) mod 2048)*2+k-2;
for k:= k step 2 until m do
if K(k) extract 12 = inmbl then begin
i:= K(k) shift (-12) extract 12;
j:= K(k) shift (-24) extract 12;
i:= ((i-1)*(n+n-i))//2+j;
F(i):= F(i)+fv(rl)*K(k+1)
end end l;
open(z,4,string track(4*inm+2),0);
setposition(z,0,Vpl(ibl,1)); k:= 128;
for i:=1 step 1 until l do begin
k:= k+1; if k=129 then begin
k:= l-i+1; outrec(z,if k>128 then 128 else k); k:= 1 end;
z(k):= F(i)
end; close(z,true);
ri:= 0;
for i:=1 step 1 until n do begin rj:= 0; l:= i-1;
for j:=1 step 1 until red do begin
W:= 0; rl:= n*(j-1)-rj+i; rj:= rk:= rj+j;
for k:=j step 1 until l do begin
W:= W+ F(rl)*V(rk); rk:= rk+k; rl:= rl+n-k end;
for k:=k step 1 until n do begin
W:= W+ F(ri+k)*V(rk); rk:= rk+k end;
L(i,j):= W
end j; ri:= ri+n-i
end i; ri:= 0;
for i:=1 step 1 until red do begin ri:= ri+i;
for j:=i step 1 until red do begin
W:= 0; rk:= ri;
for k:=i step 1 until n do begin
W:= W+ V(rk)*L(k,j); rk:= rk+k end;
L(i,j):= W
end j;
for j:=i-1 step -1 until 1 do L(i,j):= L(j,i)
end i end F;
if NV<>0 then begin
begin array U(1:red,1:red);
comment L is diagonalized by the transformation U * L * Utr
using Householder's method. The eigenvalues are pro-
duced in ascending order;
if red>1 then begin
for j:=1 step 1 until red do
for i:=1 step 1 until j do U(j,i):= L(i,j);
tridql(red,f,U)
end else begin
f(1):= L(1,1); U(1,1):= 1 end;
comment The matrix L = V * Utr is calculated;
ri:= 0;
for i:=1 step 1 until n do begin
for j:=1 step 1 until red do L(i,red+1-j):=
sum(V(ri+k)*U(j,k),k,1,if i<red then i else red);
ri:= ri+i
end i;
comment The matrices Jab,Q = U * Jab,s are calculated in
J, and the t-sums are formed. Further the matrix pro-
ducts L * (LAMBDA)&-1 * Jab,Q are formed in D;
if hb then begin array J(1:red,1:NT);
open(z,4,string track(4*inm+4),0); setposition(z,0,Dpl(ibl,ig));
k:= 128;
for i:=1 step 1 until red do
for j:=1 step 1 until NT do D(i,j):= z(kstep);
close(z,true);
for j:=1 step 1 until NT do
for i:=1 step 1 until red do
J(red+1-i,j):= sum(U(i,k)*D(k,j),k,1,red);
for i:=1 step 1 until red do begin
h:= 1/f(red+1-i);
for j:=1 step 1 until NT do begin
H:= D(i,j):= J(i,j); J(i,j):= H*h end
end;
for k:=1 step 1 until NT do
t(k):= t(k)-sum(D(i,k)*J(i,k),i,1,red);
t(7):= t(7) - sum(D(i,1)*J(i,2),i,1,red);
t(8):= t(8) - sum(D(i,2)*J(i,3),i,1,red);
t(9):= t(9) - sum(D(i,3)*J(i,1),i,1,red);
for j:=1 step 1 until NT do
for i:=1 step 1 until n do
D(i,j):= sum(L(i,k)*J(k,j),k,1,red)
end end U;
comment The frequencies are printed and the corresponding
contributions to the normal equation are calculated;
begin array avib(1:red,1:NK);
k:= 1; q:= 0;
for l:=1 step 1 until NKT do begin
m:= K(k) shift (-36) extract 12;
if m>2048 then k:= (m mod 2048)*2+k else begin
q:= q+1; m:= m+m+k-2;
for p:=1 step 1 until red do avib(p,q):= 0;
for k:=k step 2 until m do
if K(k) extract 12 = inmbl then begin
i:= K(k) shift (-12) extract 12;
j:= K(k) shift (-24) extract 12;
h:= K(k+1)*(if i=j then 1 else 2);
for p:=1 step 1 until red do
avib(p,q):= avib(p,q) + L(i,p)*L(j,p)*h
end end end l;
comment The elements d(lambda)/dK have been calculated in avib;
k:= 128//(NK+1); open(z,4,<:atrack:>,0);
setposition(z,0,Na//k); swoprec(z,(Na mod k)*(NK+1));
cr;
for m:=1 step 1 until red do begin
NO:= NO+1; On:= On+1; nu:= Obs(On); On:= On+1; dnu:= Obs(On);
h:= f(red+1-m)*1.697372'6;
if h<0 then begin write(res,<:
***negative lambda, inm, ig, ibl, red, n = :>,
<<dddd>,inm,ig,ibl,red,n); goto stop end;
h:= sqrt(h); cr;
write(res,<<________-dddd.d>, nu, h, nu-h, dnu);
H:= 0.848686'6/h; h:= nu-h;
for i:=1 step 1 until NK do avib(m,i):= avib(m,i)*H;
for j:=1 step 1 until NS do begin
H:= Obs(Ol+j+j-1);
p:= H extract 24; q:= H shift (-24) extract 24;
l:= if NO=p then 1 else if NO=q then -1 else 0;
if l<>0 then begin
ashifts(j,NK+1):= ashifts(j,NK+1) + h*l;
for i:=1 step 1 until NK do
ashifts(j,i):= ashifts(j,i) + avib(m,i)*l
end end j;
if dnu<>0 then begin
swoprec(z,NK+1); z(NK+1):= h:= h/dnu;
for i:=1 step 1 until NK do z(i):= avib(m,i)/dnu;
Na:= Na+1; s:= s + h*h
end end k;
close(z,true)
end avib end NV<>0
else begin
array J(1:red,1:NT);
open(z,4,string track(4*inm+4),0);
setposition(z,0,Dpl(ibl,ig)); k:= 128;
for i:=1 step 1 until red do
for j:=1 step 1 until NT do D(i,j):= z(kstep);
close(z,true);
comment Fs in L is inverted (by the method INVRS from CACM
66-P1-0. The t-sums are formed, and the matrix products
D = V * Fs&-1 * Jab,s are calculated;
m:= red-1;
for k:=0 step 1 until m do begin
H:= 1/L(1,1);
for i:=2 step 1 until red do f(i-1):= L(1,i);
for i:=1 step 1 until m do begin
L(i,red):= h:= -f(i)*H;
for j:=i step 1 until m do L(i,j):= L(i+1,j+1) + f(j)*h
end i;
L(red,red):= -H
end k;
for i:=1 step 1 until red do
for j:=i step 1 until red do L(i,j):= L(j,i):= -L(i,j);
for j:=1 step 1 until NT do
for i:=1 step 1 until red do J(i,j):= sum(L(i,k)*D(k,j),k,1,red);
for k:=1 step 1 until NT do
t(k):= t(k) - sum(D(i,k)*J(i,k),i,1,red);
t(7):= t(7) - sum(D(i,1)*J(i,2),i,1,red);
t(8):= t(8) - sum(D(i,2)*J(i,3),i,1,red);
t(9):= t(9) - sum(D(i,3)*J(i,1),i,1,red);
ri:= 0;
for i:=1 step 1 until n do begin
for j:=1 step 1 until NT do
D(i,j):= sum(V(ri+k)*J(k,j),k,1,if i<red then i else red);
ri:= ri+i end i;
end NV=0;
comment Now the derivatives of the CD-constants are calculated;
if hb then begin array ac(1:9);
k:= 1; q:= 0;
for l:=1 step 1 until NKT do begin
m:= K(k) shift (-36) extract 12;
if m>2048 then k:= (m mod 2048)*2+k else begin
q:= q+1; m:= m+m+k-2;
for i:=1 step 1 until 9 do ac(i):= 0;
for k:=k step 2 until m do
if K(k) extract 12 = inmbl then begin
i:= K(k) shift (-12) extract 12;
j:= K(k) shift (-24) extract 12;
h:= K(k+1)*(if i=j then 1 else 2);
for p:=1 step 1 until NT do ac(p):= ac(p)+h*D(i,p)*D(j,p);
h:= h*0.5;
ac(7):= ac(7) + (D(i,1)*D(j,2)+D(i,2)*D(j,1))*h;
ac(8):= ac(8) + (D(i,2)*D(j,3)+D(i,3)*D(j,2))*h;
ac(9):= ac(9) + (D(i,3)*D(j,1)+D(i,1)*D(j,3))*h
end k;
for p:=1 step 1 until NC do
acd(p,q):= acd(p,q) + sum(C(p,j)*ac(j),j,1,9)
end end l end NC<>0
end V, L, D, f
end bl step;
comment The centrifugal distortion constants are printed and
their contributions to the A-matrix are calculated;
if hb then begin cr;
k:= 128//(NK+1); open(z,4,<:atrack:>,0);
setposition(z,0,Na//k); swoprec(z,(Na mod k)*(NK+1));
for k:=1 step 1 until NC do begin
NO:= NO+1; Oc:= Oc+1; nu:= Obs(Oc); Oc:= Oc+1; dnu:= Obs(Oc);
h:= sum(C(k,i)*t(i),i,1,9); cr;
write(res,<<____-dd.dddd000>, nu, h, nu-h, dnu);
h:= nu-h;
for j:=1 step 1 until NS do begin
H:= Obs(Ol+j+j-1);
p:= H extract 24; q:= H shift (-24) extract 24;
l:= if NO=p then 1 else if NO=q then -1 else 0;
if l<>0 then begin
ashifts(j,NK+1):= ashifts(j,NK+1) + h*l;
for i:=1 step 1 until NK do
ashifts(j,i):= ashifts(j,i) + acd(k,i)*l
end end j;
if dnu<>0 then begin
swoprec(z,NK+1); z(NK+1):= h:= h/dnu;
for i:=1 step 1 until NK do z(i):= acd(k,i)/dnu;
Na:= Na+1; s:= s + h*h
end end k;
close(z,true); On:= Oc;
end NC<>0 else if centrifugal then On:= On+1
end C, t, acd end g step end S, Vpl, Dpl end inm;
On:= On+1;
if NS>0 then begin
if linecount+NS>55 then writepage else
begin cr; cr; cr end; write(res,<:Shifts::>); cr;
k:= 128//(NK+1); open(z,4,<:atrack:>,0);
setposition(z,0,Na//k); swoprec(z,(Na mod k)*(NK+1));
for i:=1 step 1 until NS do begin
On:= On+1; h:= Obs(On); On:= On+1; dnu:= Obs(On);
k:= h shift (-24) extract 24;
l:= h extract 24;
cr; write(res,<: Obs. no.::>,<<ddd>,l,
<: - Obs. no.::>,k,<<___-ddd.dd00000>,ashifts(i,NK+1),dnu);
if dnu<>0 then begin
swoprec(z,NK+1);
for j:=NK+1 step -1 until 1 do z(j):= ashifts(i,j)/dnu;
Na:= Na+1; s:= s + z(NK+1)**2
end end;
close(z,true);
end end ashifts;
s:= s/(if Na=NK then 1 else Na-NK);
hb:= s>=so and so<>0; NC:= NK+1; writepage;
if -, hb then begin
zone z2(128*2,2,stderror);
open(z2,4,<:atrack:>,0); open(z,4,string track(4),0);
for k:=1 step 1 until Na do begin
outrec(z,NC); inrec(z2,NC);
for i:=1 step 1 until NC do z(i):= z2(i)
end k;
close(z,true); close(z2,true);
end z2;
begin array amax, X, TD(1:NK);
begin
real p,q,length,gramdet,amin,df,dfu,dfl,h0,hu,hl;
integer dfsteps;
array T(1:NK,1:NK), a(1:Na,1:NC);
if hb then begin
h:= 0;
for i:=1 step 1 until NK do
if h<abs B(i) then h:= abs B(i);
if h<'-5 then goto packstop;
Xmax:= h*0.5
end else begin
for i:=1 step 1 until NK do B(i):= 0; so:= s;
end;
open(z,4,string track(4),0);
DF:= dfu:= dfl:= 0; dfsteps:= -1;
newDF:
setposition(z,0,0);
for i:=1 step 1 until Na do begin inrec(z,NC);
for j:=1 step 1 until NC do a(i,j):= z(j) end;
gramdet:= 1; dfsteps:= dfsteps+1;
if dfsteps=0 then begin
comment For calculation of significant digits the largest
a(i,j) in each column are determined;
for j:=1 step 1 until NK do begin
p:= 0;
for i:=1 step 1 until Na do begin
q:= abs a(i,j); if q>p then p:= q
end;
amax(j):= p
end end;
for j:=1 step 1 until NK do begin
length:= DF; for i:=1 step 1 until Na do
length:= length + a(i,j)**2;
if j=1 or amin>length then amin:= length;
l:= j-1;
for k:=1 step 1 until l do begin h:= T(k,k);
p:= 0; for i:=1 step 1 until Na do
p:= p + a(i,j)*a(i,k);
T(k,j):= p; q:= p*h;
for i:=1 step 1 until Na do a(i,j):= a(i,j)-a(i,k)*q
end k, orthogonalisation;
p:= 1;
for k:=l step -1 until 1 do begin
h:= T(k,k); q:= T(k,j);
for i:=k+1 step 1 until l do
q:= q + T(j,i)*T(k,i);
q:= -q*h; T(j,k):= q; p:= p+q*q
end;
p:= p*DF; q:= 0;
for i:=1 step 1 until Na do begin
p:= p + a(i,j)**2;
q:= q + a(i,j)*a(i,NC)
end i;
T(j,j):= 1/p; X(j):= q; gramdet:= gramdet*p/length
end j;
h:= 0;
for j:=NK step -1 until 1 do begin q:= T(j,j);
p:= X(j); for k:=j+1 step 1 until NK do
p:= p - T(j,k)*X(k);
H:= X(j):= q*p; h:= h + H*H
end solution;
h:= 0;
for i:=1 step 1 until NK do
if h<abs X(i) then h:= abs X(i);
if closeres then begin
write(out,<< d.dddddd'-d>,DF,h,<:<10>:>);
setposition(out,0,0) end;
if h>Xmax or (h*1.25<Xmax and DF>0) then begin
if DF=0 then H:= amin else
if dfu=0 or dfl=0 then
H:= DF*ln(hu/(hu-h))/ln(hu/(hu-Xmax*0.9))
else begin
if h>Xmax then begin h0:= hl; df:= dfl end
else begin h0:= hu; df:= dfu end;
p:= ln(df*h0/(DF*h))/(h-h0);
H:= if p>0 then DF*h*exp(p*h)/(Xmax*0.9*exp(Xmax*0.9*p))
else DF + (DF-df)*(Xmax*0.9-h)/(h-h0)
end;
if h>Xmax then begin hu:= h; dfu:= DF end
else begin hl:= h; dfl:= DF end;
DF:= H; goto newDF
end;
for i:=1 step 1 until NK do B(i):= X(i)-B(i);
close(z,true);
for i:=1 step 1 until NK do
for j:=i step 1 until NK do begin
q:= if i=j then 1 else T(j,i); q:= q*T(j,j);
for k:=j+1 step 1 until NK do
q:= q + T(k,i)*T(k,k)*T(k,j);
T(i,j):= q
end;
comment The variances and correlation matrix elements are printed;
k:= 1; l:= NK mod 6; s:= sqrt(s);
if closeres then begin
write(out,<:<10>Iteration no.:>,<<ddd>,IT,
<: sigma = :>,<<d.dd'-d>,s,
<: DFsteps: :>,<<ddd>,dfsteps,
<:<10>Gramdet.= :>,<<d.d'-dd>,gramdet,
<: Xmax = :>,<<d.ddd>,Xmax,
<: DF = :>,<<d.d'-d>,DF,<:<10><10>:>);
setposition(out,0,0) end;
write(res,<:Iteration no.:>,<<ddd>,IT,
<: sigma = :>,<<d.dd'-d>,s,
<: DFsteps: :>,<<ddd>,dfsteps); cr;
write(res,<:Gramdet.= :>,<<d.d'-dd>,gramdet,
<: Xmax = :>,<<d.ddd>,Xmax,
<: DF = :>,<<d.d'-d>,DF);
cr; cr; write(res,<:Correlation and variance::>); cr;
for i:=1 step 1 until NK do begin
p:= sqrt(T(i,i)); T(i,i):= p*s; TD(i):= 1/p end;
for i:=NK-1 step -1 until 1 do
for j:= i+1 step 1 until NK do T(i,j):= T(i,j)*TD(i)*TD(j);
for i:=1 step 1 until NK do TD(i):= T(i,i);
Cor:
if linecount+l>62 then writepage;
for i:=1 step 1 until l do begin cr;
for j:=k step 1 until i-1 do write(res,<: :>);
for j:=j step 1 until l do write(res,<<__-d.dd000>,T(i,j))
end i;
cr; k:= l+1; l:= l+6; if k<NK then goto Cor
end T,a;
comment The new fv matrix is calculated and printed;
begin array fn(1:NKT*2);
open(z,4,string track(2),0);
k:= 128+NKT; l:= NKT*2; if k>255 then setposition(z,0,k//128-1);
for i:=1 step 1 until l do fn(i):= z(kstep);
if linecount+NKT>57 then writepage else begin cr; cr; cr end;
write(res,
<:Parameter. f new sigma delta digits:>);
cr; cr; j:= 1; k:=1; l:= 0;
hb:= DF<>0; s:= NK/s; H:= ln(10);
for i:=1 step 1 until NKT do begin
write(res,<<ddd>,i,<: :>,string fn(increase(j)));
m:= K(k) shift (-36) extract 12;
if m>2048 then begin
write(res,<<__-dd.dddddd>,fv(i),<: (fixed):>);
m:= m mod 2048
end else begin
l:= l+1; fv(i):= fv(i)+B(l);
n:= entier(ln(s*amax(l)*abs fv(i))/H)+2;
hb:= hb or abs B(l)*amax(l)*s > 1;
write(res,<<__-dd.dddddd>,fv(i),TD(l),B(l),<<_____dd>,n)
end;
cr; k:= m+m+k;
end i;
for i:=1 step 1 until NK do B(i):= X(i);
end fn end TD,X;
close(z,true);
H:= time; systime(1,starttime,time);
if hb and 2*time-H<stoptime then goto igen;
packstop:
begin
comment Packing of input for a continued calculation;
procedure pack(n);
value n; integer n;
begin z(m):= z(m) shift 12 add n end;
open(z,4,string track(1),0); monitor(42,z,0,tail);
setposition(z,0,tail(1)-1); swoprec(z,128);
NM:= NM shift 9 add page;
m:= 128; z(m):= 0; pack(NM); pack(IT); pack(NO); pack(On);
m:= 127; z(m):= 0; pack(NF); pack(NK); pack(NKF); pack(NS);
z(126):= Xmax; z(125):= so;
close(z,true); open(z,4,string track(2),0); k:= 128;
for i:=1 step 1 until NKT do begin
k:= k+1;
if k=129 then begin k:= NKT-i+1;
swoprec(z,if k>128 then 128 else k); k:= 1 end;
z(k):= fv(i);
end;
setposition(z,0,(3*NKT)//128); k:= m:= 3*NKT mod 128;
for i:=1 step 1 until NK do begin
k:= k+1;
if k=129 then begin k:= m+NK-i+1;
swoprec(z,if k>128 then 128 else k); k:= 1 end;
z(k):= B(i);
end;
close(z,true) end
end Obs, fv, B, K;
readhead(in,mol,1);
if hb then begin
i:= 1; l:= 0; k:=psubmit(string mol(increase(i)),l);
i:= 1; write(out,string mol(increase(i)));
if k=0 then write(out,<: job nr. :>,<<dddd>,l) else
if k>0 then begin
i:= (k shift (-12)) -2;
write(out,case i of (<: no jobfile:>,<: que full:>,
<: submitting stopped:>,<: unknown:>,
<: area process error:>))
end else write(out,<: wait answer error:>);
end else begin
cr; cr; cr; write(out,<:Iterations completed:>) end;
stop:
write(res,<:<10><12><25>:>); close(res,closeres);
open(res,4,<:atrack:>,0); monitor(48,res,0,tail); close(res,true);
end end
▶EOF◀