|
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: »algvibrot«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦7e928b248⟧ »algbib«
└─⟦this⟧
;gosav time.180
lookup vibrot
if ok.no
(lookup xbmatrix
if ok.no
(r=algol extxbmatrix index.no
rename r.xbmatrix
permanent xbmatrix.17)
vibrot=set 148
permanent vibrot.17
vibrot=algol index.no list.no xref.no)
\f
VIBROT 8-4-1975. GOS
begin
boolean ha,hb,nl,sp,centrifugal,zeta,displ,closeres,fclear,msa;
real T;
integer linecount,page,space,Na,N3,BL,G,ibl,ig,nt,i,j,k,l,m,n,bpl;
integer array tail(1:10);
array head(1:12), species(1:5);
zone res(128,1,stderror);
procedure cr;
begin linecount:= linecount+1;
write(res,nl,1); if linecount>60 then writepage
end cr;
procedure writepage;
begin integer i;
page:= page+1; linecount:= 3; i:= 1;
write(res,<:<12>:>,nl,3,string head(increase(i)),
sp,space,<<dd>,page, nl, 2);
writespecies
end writepage;
procedure writespecies;
begin integer i;
i:= 1; write(res,string species(increase(i)));
cr; cr end;
closeres:=outmedium(res); page:= -1;
nl:= false add 10; sp:= false add 32;
space:= 67 - readhead(in,head,1); species(1):= 0;
writepage; read(in,m);
hb := m >= 1000;
displ:= m mod 1000 >= 100;
zeta := m mod 100 >= 10; m:= m mod 10;
msa := m mod 4 >= 2;
centrifugal:= m mod 2 = 1;
if msa then read(in,T);
read(in,G,Na,BL); N3:= abs Na*3;
begin
integer array S(1:BL,1:2);
for ibl:=1 step 1 until BL do read(in,S(ibl,1));
if xbmatrix(Na,BL,S,hb,ha,res) <> 0 then goto stop;
nt:= if -,centrifugal then 1 else if ha then 4 else 6;
if G<>0 then begin
zone f(128*2,2,stderror);
l:= m:= 0; for ibl:=1 step 1 until BL do begin
l:= l+(S(ibl,1)*N3+127)//128;
m:= m+(S(ibl,1)*(S(ibl,1)+1)+254)//256 end;
ha:= zeta or centrifugal;
if ha then begin
open(f,4,<:kbtrack:>,0); tail(1):= l;
monitor(48,f,0,tail); monitor(40,f,0,tail);
close(f,true) end;
open(f,4,<:ftrack:>,0); i:= monitor(42,f,0,tail); fclear:= i<>0;
write(res,false add 10,3,<:
SYMMETRY FORCE CONSTANTS:
:>);
if i=0 then begin
for ibl:=1 step 1 until BL do begin
array fn(1:S(ibl,1));
write(res,<:
Species:>,<<ddd>,ibl,<:<10>:>); n:= S(ibl,1); k:= 129;
for i:= 1 step 1 until n do begin
m:= n+1-i; write(res,<:<10>:>);
for j:=1 step 1 until m do begin
if k=129 then begin k:= 1; inrec(f,128) end;
fn(j):= f(k); k:= k+1 end j;
m:= m+1; for m:= m-1 while fn(m)=0 do;
for j:=1 step 1 until m do begin
write(res,<< -dd.dddddd>,fn(j));
if j<>m and j mod 6=0 then write(res,<:,<10>:>)
end j;
write(res,<:<10>:>)
end end ibl end else begin
tail(1):= m; monitor(40,f,0,tail);
for ibl:=1 step 1 until BL do begin
write(res,<:
Species:>,<<ddd>,ibl,<:<10>:>); n:= S(ibl,1); k:= 129;
for i:= 1 step 1 until n do begin
l:= 0; m:= n+1-i;
for j:=1 step 1 until m do begin
if k=129 then begin k:= 1; outrec(f,128) end;
if l<>10 and l<>115 then begin read(in,f(k));
repeatchar(in); readchar(in,l);
if (j-1) mod 6 = 0 then write(res,<:<10>:>);
write(res,<< -dd.dddddd>,f(k)) end else f(k):= 0; k:= k+1
end j;
if l<>10 and l<>115 then begin write(res,<:
***delimiter in F-input.:>); goto stop end; write(res,<:<10>:>)
end end ibl end;
close(f,true);
if centrifugal then begin
m:= 0;
for ibl:=1 step 1 until BL do
m:= m + (S(ibl,1)*nt+127)//128;
tail(1):= m; open(f,4,<:jtrack:>,0);
monitor(48,f,0,tail); monitor(40,f,0,tail);
close(f,true);
end;
if msa then begin
open(f,4,<:msatrack:>,0); monitor(48,f,0,tail);
tail(1):= ((N3+1)*N3+254)//256;
monitor(40,f,0,tail); monitor(50,f,17,tail);
close(f,true)
end end G<>0;
if displ then begin integer r;
writepage;
write(res,<:B-matrix (S = B * X):>); bpl:= 0;
for ibl:=1 step 1 until BL do begin
zone B((S(ibl,1)*N3+127)//128*128,1,stderror);
open(B,4,<:btrack:>,0); setposition(B,0,bpl);
n:= S(ibl,1); inrec(B,n*N3); bpl:= bpl+(n*N3+127)//128;
if linecount+n>54 then begin
writepage; write(res,<:B-matrix (S = B * X):>) end;
cr; cr; cr; write(res,<:Species:>,<<ddd>,ibl); k:= -1;
rep:
k:= k+2; l:= if Na>k then k+1 else Na; cr; cr;
for j:=k step 1 until l do begin
for i:=1 step 1 until 3 do
write(res,<: :>,false add (87+i),1,<<dd>,j,<: :>);
write(res,<: :>) end;
cr; m:= (n-1)*N3;
for r:=0 step N3 until m do begin
cr;
for j:=k step 1 until l do begin
for i:=(j-1)*3+1 step 1 until 3*j do if abs B(r+i)<5'-7 then
write(res,<: 0 :>) else
write(res,<<-ddd.dddddd>,B(r+i));
write(res,<: :>) end
end r;
if l<Na then begin
if linecount+n>57 then begin
writepage; write(res,<:B-matrix (S = B * X):>);
cr; cr; cr; write(res,<:Species:>); write(res,<<ddd>,ibl) end;
goto rep
end;
close(B,true) end ibl
end displ;
for ig:=1 step 1 until G do begin
real Ia, Ib, Ic; boolean planar;
array f2(1:if centrifugal then N3-6 else 3), Ma(1:Na),
U(1:3,1:3), t(1:if centrifugal then 9 else 1);
readhead(in,species,1);
if atomic(Ma,Na) then goto stop;
if linecount+S(1,1)>52 or ig=1 then writepage else
begin cr; cr; cr; cr; writespecies end;
write(res,<:Atomic masses::>); cr;
for i:=1 step 1 until Na do begin
write(res,<<dddd.dddddd>,Ma(i)); if i mod 6 = 0 then cr end;
comment If wanted the tensor of inertia is calculated, and a
transformation to the principal axes system is performed;
if ha then begin
array IP(1:3);
zone X, UX((N3+127)//128*128,1,stderror);
open(X,4,<:xtrack:>,0); inrec(X,N3); Ib:= 1/sum(Ma(i),i,1,Na);
for i:=-2 step 1 until 0 do begin
Ia:= sum(X(i+j*3)*Ma(j),j,1,Na)*Ib;
for k:=i+3 step 3 until N3 do X(k):= X(k)-Ia
end i, center of mass corr.;
Ia:= sum(Ma(i)*X(3*i-2)**2,i,1,Na);
Ib:= sum(Ma(i)*X(3*i-1)**2,i,1,Na);
Ic:= sum(Ma(i)*X(3*i)**2,i,1,Na);
U(1,1):= Ib+Ic; U(2,2):= Ia+Ic; U(3,3):= Ia+Ib;
U(2,1):= -sum(Ma(i)*X(3*i-2)*X(3*i-1),i,1,Na);
U(3,2):= -sum(Ma(i)*X(3*i-1)*X(3*i),i,1,Na);
U(3,1):= -sum(Ma(i)*X(3*i)*X(3*i-2),i,1,Na);
tridql(3,IP,U);
Ia:= IP(1); Ib:= IP(2); Ic:= IP(3);
open(UX,4,<:xtrack:>,0); setposition(UX,0,(N3+127)//128);
outrec(UX,N3);
for k:=N3-3 step -3 until 0 do
for i:=1 step 1 until 3 do UX(k+i):= sum(U(i,j)*X(k+j),j,1,3);
close(X,true); close(UX,true);
planar:= abs (Ic-Ia-Ib) < (Ia+Ib+Ic)*'-8;
if centrifugal then begin for i:=1 step 1 until 9 do t(i):= 0 end
end ha;
for i:=1 step 1 until Na do Ma(i):= 1/Ma(i); bpl:= 0;
for ibl:=1 step 1 until BL do begin
n:= S(ibl,1);
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 L(1:n,1:n), f(1:n), V,K(1:n*(n+1)//2);
k:= N3*n;
begin
zone B((k+127)//128*128,1,stderror);
open(B,4,<:btrack:>,0); setposition(B,0,bpl);
inrec(B,k); ri:= rl:= 0;
for i:=1 step 1 until n do begin rj:= 0;
for j:=1 step 1 until i do begin
V(rl+j):= sum(B(ri+k)*B(rj+k)*Ma((k+2)//3),k,1,N3);
rj:= rj+N3
end j;
ri:= ri+N3; rl:= rl+i
end i; close(B,true) end B;
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*i 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 ha or displ or msa 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. Now the matrix product
Vtr * F * V is formed in L;
end ha or displ;
S(ibl,2):= red:= n-red; k:= n*(n+1)//2;
begin array FS(1:k); zone f(128,1,stderror);
m:= 0; for i:=ibl-1 step -1 until 1 do
m:= m+ (S(i,1)*(S(i,1)+1)+254)//256;
open(f,4,<:ftrack:>,0); setposition(f,0,m);
j:= 129; ri:= 0;
for i:=1 step 1 until k do begin
if j=129 then begin j:= 1; inrec(f,128) end;
FS(i):= f(j); j:= j+1
end;
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+ FS(rl)*V(rk); rk:= rk+k; rl:= rl+n-k end;
for k:=k step 1 until n do begin
W:= W+ FS(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 end i;
close(f,true) end FS, f;
begin array U(1:red,1:red);
comment L is diagonalized by the transformation U * L * Utr using
Householder's method. The eigenvalues are produced in in-
creasing 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, and the frequencies
and L are printed;
ri:= 0;
for i:=1 step 1 until n do begin
for j:=1 step 1 until red do
L(i,j):= sum(V(ri+k)*U(j,k),k,1,if i<red then i else red);
ri:= ri+i
end i;
if linecount+n>56 then writepage else cr; cr;
write(res,<:Species:>,<<ddd>,ibl);
k:= red+5; l:= red+1; m:= 0;
for i:=ibl-1 step -1 until 1 do m:= m+S(i,1);
B: cr; write(res,<: :>); k:= k-5;
l:= if l>5 then l-5 else 1;
for j:=k step -1 until l do
write(res,<<___-ddddd.dd>,
sign(f(j))*sqrt(abs f(j)*1.697372'6)); cr;
for i:=1 step 1 until n do begin
hb:= false;
for j:=k step -1 until l do hb:= hb or abs(L(i,j))>'-6;
if hb then begin
write(res,<:<10> S:>,<<dd>,m+i,<: :>);
for j:=k step -1 until l do if abs(L(i,j))<5'-7 then
write(res,<: 0 :>) else
write(res,<<-dddd.dddddd>,L(i,j))
end else linecount:= linecount-1
end i;
linecount:= linecount+n;
if l>1 then begin
if linecount+n>59 then begin writepage; cr;
write(res,<:Species:>,<<ddd>,ibl) end; cr;
goto B end;
if ha or displ or msa then begin
comment The left inverse of L ( = U * K ) is formed in L;
m:= n-red; rj:= m*(m+1)//2; l:= m+1;
for j:=1 step 1 until n do begin
if j>l then begin l:= l+1; rj:= rj+j end else rj:= rj+1;
for i:=1 step 1 until red do begin
W:= 0; rk:= rj;
for k:=l step 1 until n do begin
W:= W+ U(i,k-m)*K(rk); rk:= rk+k end;
L(i,j):= W
end i end j end ha
end U;
k:= N3*n;
if ha or displ or msa then begin
zone B((k+127)//128*128,1,stderror),
X((N3+127)//128*128,1,stderror);
open(B,4,<:btrack:>,0); setposition(B,0,bpl); inrec(B,k);
if ha then begin close(B,true);
open(B,4,<:kbtrack:>,0); setposition(B,0,bpl); outrec(B,k) end;
comment The matrix product U * K * B is calculated in B. The
transformation X = A * Q is Atr = U * K * B * my;
for j:=1 step 1 until N3 do begin
for i:=1 step 1 until red do begin
W:= 0; rj:= j;
for k:=1 step 1 until n do begin
W:= W+ L(i,k)*B(rj); rj:= rj+ N3 end;
X(i):= W
end; rj:= j;
for i:=1 step 1 until red do begin
B(rj):= X(i); rj:= rj+ N3 end
end j;
if msa then begin
comment Mean square amplitudes are calculated;
zone sigma(128,1,stderror);
array delta(1:red);
open(sigma,4,<:msatrack:>,0);
for k:=1 step 1 until red do begin
W:= sqrt(f(k)); VV:= exp(-1874.500134*W/T);
delta(k):= 0.012939112*(1+VV)/(1-VV)/W
end;
m:= 128;
for i:=1 step 1 until N3 do
for j:=1 step 1 until i do begin
m:= m+1; W:= 0; ri:= i; rj:= j;
for k:=1 step 1 until red do begin
W:= W + delta(k)*B(ri)*B(rj);
ri:= ri+N3; rj:= rj+N3
end;
if m=129 then begin
m:=1; swoprec(sigma,128);
if ibl=1 then begin
for k:=1 step 1 until 128 do sigma(k):= 0
end end;
sigma(m):= sigma(m) + W*Ma((i-1)//3+1)*Ma((j-1)//3+1)
end i,j;
close(sigma,true)
end msa;
if displ then begin
if linecount+red>54 then begin
writepage; cr; write(res,<:Species:>,<<ddd>,ibl) end;
cr; cr; cr; write(res,<:A-matrix ( X = A * Q ) ::>); k:= -1;
rep:
k:= k+2; l:= if Na>k then k+1 else Na; cr; cr;
for j:=k step 1 until l do begin
for i:=1 step 1 until 3 do
write(res,<: :>,false add (87+i),1,<<dd>,j,<: :>);
write(res,<: :>) end;
cr;
for ri:=(red-1)*N3 step -N3 until 0 do begin
cr;
for j:=k step 1 until l do begin
for i:=(j-1)*3+1 step 1 until 3*j do
if abs B(ri+i)*Ma(j)<5'-7 then
write(res,<: 0 :>) else
write(res,<<-ddd.dddddd>,B(ri+i)*Ma(j));
write(res,<: :>) end
end ri;
if l<Na then begin
if linecount+red>57 then begin
writepage; cr; write(res,<:Species:>,<<ddd>,ibl) end;
goto rep
end end displ;
if ha then begin
comment The KB-matrix is transformed to the principal axes system;
ri:= 0;
for i:=1 step 1 until red do begin
for j:=N3-3 step -3 until 0 do
for k:=1 step 1 until 3 do
X(j+k):= sum(B(ri+j+l)*U(k,l),l,1,3);
for j:=1 step 1 until N3 do B(ri+j):= X(j);
ri:= ri+N3
end i;
m:= 0;
for i:=ibl-1 step -1 until 1 do m:= m+S(i,2);
for i:=1 step 1 until red do f2(m+i):= f(i);
if centrifugal then begin
array J(1:red,1:nt);
open(X,4,<:xtrack:>,0); setposition(X,0,(N3+127)//128);
inrec(X,N3);
for j:=1 step 1 until 3 do begin
k:= j mod 3 + 1; l:= k mod 3 + 1; rk:= k; rl:= l;
for i:=1 step 1 until red do begin W:= 0;
for m:=N3-3 step -3 until 0 do
W:= W+ B(m+rk)*X(m+k) + B(m+rl)*X(m+l);
J(i,j):= W; rk:= rk+N3; rl:= rl+N3 end i
end j;
for j:=4 step 1 until nt do begin
k:= (j-1) mod 3 + 1; l:= k mod 3 + 1; rk:= k;
for i:=1 step 1 until red do begin W:= 0;
for m:=N3-3 step -3 until 0 do W:= W- B(m+rk)*X(m+l);
J(i,j):= W; rk:= rk+N3 end i
end j;
close(X,true);
if zeta then begin
if linecount+red>56 then begin
writepage; cr; write(res,<:Species:>,<<ddd>,ibl) end;
cr; cr; cr;
m:= 0; for i:=ibl-1 step -1 until 1 do m:= m+S(i,2);
write(res,<: Jaa Jbb Jcc Jab:>);
if nt=6 then write(res,<: Jbc Jca:>); cr;
for i:=red step -1 until 1 do begin
cr; write(res,<:Q:>,<<dd>,m+red+1-i);
for j:=1 step 1 until nt do if abs J(i,j)<2.5'-7 then
write(res,<: 0 :>) else
write(res,<<-ddd.dddddd>,J(i,j)*2) end;
end zeta;
comment Now the t-sums are formed;
for i:=1 step 1 until red do f(i):= 1/f(i);
for k:=1 step 1 until nt do
t(k):= t(k) + sum(f(i)*J(i,k)**2,i,1,red);
t(7):= t(7) + sum(f(i)*J(i,1)*J(i,2),i,1,red);
t(8):= t(8) + sum(f(i)*J(i,2)*J(i,3),i,1,red);
t(9):= t(9) + sum(f(i)*J(i,3)*J(i,1),i,1,red);
comment J contains 1/2*Jgg',Q (sqrt(u)*Å). my' elements are
calculated in MHz corresponding to dimensionless
normal coordinates. 1/lambda (uÅ/mdyn) was stored
in f(i);
k:= 0;
for i:=ibl-1 step -1 until 1 do
k:= k+1+(S(i,2)-1)//(128//nt);
open(X,4,<:jtrack:>,0); setposition(X,0,k);
W:= -3.251942'5;
for i:=red step -1 until 1 do begin
VV:= W*f(i)**0.25; outrec(X,nt);
for j:=1 step 1 until nt do X(j):= J(i,j)*VV/
(case j of (Ia*Ia,Ib*Ib,Ic*Ic,Ia*Ib,Ib*Ic,Ic*Ia))
end;
close(X,true)
end centrifugal;
end ha;
close(B,true) end ha or displ
end L,f,V,K;
bpl:= bpl + (N3*n+127)//128 end ibl;
comment Now the moments of inertia, rotational constants and
centrifugal distortion constants are printed;
if centrifugal then begin
k:= 0; l:= 3; m:= 2;
for ibl:=1 step 1 until BL do k:= k+ S(ibl,2);
end else k:= l:= m:= 0;
if ha then begin
integer jbl,ni,nj,noi,noj,pi,pj,ri1,ri2,rj1,rj2;
real h,V,W,Z2,Ia2,Ib2,Ic2;
array my(1:k,1:nt), alfa(1:k,1:l,1:m), f(1:k);
if centrifugal then begin
for i:=1 step 1 until k do begin
f(i):= sqrt(f2(i));
for j:=1,2,3 do alfa(i,j,2):= 0
end;
if linecount>54 then writepage else begin cr; cr end; cr;
write(res,<:Moments of Inertia ::>,<<___ddd.ddddd>,Ia,Ib,Ic);
cr; h:= 505376.42; if Ia>0 then Ia:= h/Ia; Ib:= h/Ib; Ic:= h/Ic;
write(res,<:Rotational constants ::>,<<___dddddd.dd>,Ia,Ib,Ic);
hb:= false add 10; linecount:= linecount+4;
for i:=1 step 1 until 9 do t(i):= t(i)*2.0753268'-20;
Ia2:= Ia*Ia; Ib2:= Ib*Ib; Ic2:= Ic*Ic;
write(res,hb,2,<:tau-aaaa = :>,<<-dd.ddddd00>,-Ia2*Ia2*t(1),
<: tau-aabb = :>,-Ia2*Ib2*t(7),
<: tau-abab = :>,-Ia2*Ib2*t(4));
write(res,hb,1,<:tau-bbbb = :>,<<-dd.ddddd00>,-Ib2*Ib2*t(2),
<: tau-bbcc = :>,-Ib2*Ic2*t(8),
<: tau-bcbc = :>,-Ib2*Ic2*t(5));
write(res,hb,1,<:tau-cccc = :>,<<-dd.ddddd00>,-Ic2*Ic2*t(3),
<: tau-ccaa = :>,-Ic2*Ia2*t(9),
<: tau-caca = :>,-Ic2*Ia2*t(6));
begin zone z(128,1,stderror);
open(z,4,<:jtrack:>,0); l:= n:= 0;
for ibl:=1 step 1 until BL do begin
setposition(z,0,l); m:= S(ibl,2);
for i:=1 step 1 until m do begin inrec(z,nt);
for j:=1 step 1 until nt do my(n+i,j):= z(j) end;
n:= n+m; l:= l+1+(m-1)//(128//nt)
end end;
if zeta then begin
if linecount+k>50 then writepage else begin cr; cr end;
cr; write(res,<:First Derivatives of my-Tensor (MHz)::>);
cr; cr;
write(res,<: my-aa my-bb my-cc my-ab:>);
if nt=6 then write(res,<: my-bc my-ca:>); cr; cr;
for i:=1 step 1 until k do begin
write(res,<:q:>,<<dd>,i,<: :>);
for j:=1 step 1 until nt do begin
h:= my(i,j);
if abs h<'-4 then
write(res,<: 0 :>) else
write(res,if abs h<'5 then <<-ddddd.dddd>
else <<-d.dddddd'd>,h);
end; cr;
end;
for m:=1 step 1 until nt do begin noi:= 1;
for ibl:=1 step 1 until BL do begin
nj:= n:= S(ibl,2)+noi-1; noj:= noi;
for jbl:=ibl step 1 until BL do begin
array my2(noi:n,noj:nj);
h:= 0;
for i:=noi step 1 until n do
for j:=if ibl=jbl then i else noj step 1 until nj do
begin my2(i,j):= (case m of
((my(i,1)*my(j,1)/Ia+my(i,4)*my(j,4)/Ib
+(if nt=6 then my(i,6)*my(j,6)/Ic else 0))*2,
(my(i,4)*my(j,4)/Ia+my(i,2)*my(j,2)/Ib
+(if nt=6 then my(i,5)*my(j,5)/Ic else 0))*2,
((if nt=6 then my(i,6)*my(j,6)/Ia+my(i,5)*my(j,5)/Ib
else 0)+my(i,3)*my(j,3)/Ic)*2,
(my(i,1)*my(j,4)+my(i,4)*my(j,1))/Ia
+(my(i,2)*my(j,4)+my(i,4)*my(j,2))/Ib+(if nt=6
then (my(i,5)*my(j,6)+my(i,6)*my(j,5))/Ic else 0),
(my(i,4)*my(j,6)+my(i,6)*my(j,4))/Ia
+(my(i,2)*my(j,5)+my(i,5)*my(j,2))/Ib
+(my(i,3)*my(j,5)+my(i,5)*my(j,3))/Ic,
(my(i,1)*my(j,6)+my(i,6)*my(j,1))/Ia
+(my(i,4)*my(j,5)+my(i,5)*my(j,4))/Ib
+(my(i,3)*my(j,6)+my(i,6)*my(j,3))/Ic))*0.375;
if abs my2(i,j)>h then h:= abs my2(i,j)
end;
if m<4 and ibl=jbl then begin
for i:= noi step 1 until n do alfa(i,m,1):= -my2(i,i)/4
end;
if h>'-4 then begin
k:= noj-4;
rep:
k:= k+4; l:= k+4; if l>nj then l:= nj;
ni:= if ibl=jbl then l else n;
if linecount+ni-noi>54 then writepage else begin cr; cr end;
cr; if k=noj or linecount=6 then
write(res,<:Second Derivatives of my-:>,
case m of (<:aa:>,<:bb:>,<:cc:>,<:ab:>,<:bc:>,<:ca:>),
<: (MHz), Species:>,<<dd>,ibl,<: * Species:>,jbl);
cr; cr; write(res,<: :>);
for j:=k step 1 until l do
write(res,<: q:>,<<dd>,j);
cr;
for i:=noi step 1 until ni do begin
cr; write(res,<: q:>,<<dd>,i,false add 32,
if i>k then (i-k)*12+2 else 2);
for j:=if i>k then i else k step 1 until l do
if abs my2(i,j)<'-4
then write(res,<: 0 :>)
else write(res,<<-dddddd.dddd>,my2(i,j))
end;
if l<nj then goto rep end;
noj:= nj+1; if jbl<BL then nj:= S(jbl+1,2)+noj-1
end jbl;
noi:= n+1
end ibl end m end zeta else
for m:=1 step 1 until 3 do
for i:=1 step 1 until n do alfa(i,m,1):= -(case m of
(my(i,1)**2/Ia+my(i,4)**2/Ib+
(if nt=6 then my(i,6)**2/Ic else 0),
my(i,4)**2/Ia+my(i,2)**2/Ib+
(if nt=6 then my(i,5)**2/Ic else 0),
(if nt=6 then my(i,6)**2/Ia+my(i,5)**2/Ib else 0)+
my(i,3)**2/Ic))*0.1875;
end centrifugal;
if msa then begin
comment Mean square amplitudes are printed;
zone sigma(128,1,stderror);
open(sigma,4,<:msatrack:>,0);
if linecount+N3>50 then writepage else begin cr; cr end; cr;
write(res,<:Mean Binary Products of Cartesian Displacements::>);
cr; cr; cr; k:= -5;
rep:
if linecount>53 then writepage else begin cr; cr end;
k:= k+6; l:= k+5; if l>N3 then l:= N3; write(res,sp,1);
for j:=k step 1 until l do begin
write(res,sp,if j mod 3=1 then 9 else 7);
write(res,case j mod 3 + 1 of (<:Z:>,<:X:>,<:Y:>));
noj:= (j-1)//3+1;
if noj<10 then write(res,<<d>,noj,sp,1)
else write(res,<<dd>,noj)
end;
m:= (k+1)*k//2;
setposition(sigma,0,(m-1)//128); m:= m mod 128 + 128;
for i:=k step 1 until N3 do begin
if i mod 3 = 1 then cr; cr;
write(res,case i mod 3 + 1 of (<:Z:>,<:X:>,<:Y:>));
noi:= (i-1)//3 + 1;
if noi<10 then write(res,<<d>,noi,sp,1)
else write(res,<<dd>,noi);
for j:=k step 1 until l do begin
if m>128 then begin inrec(sigma,128); m:= m-128 end;
write(res,sp,if j mod 3 = 1 then 3 else 1);
if abs sigma(m)<5'-6 then write(res,<: 0 :>)
else write(res,<<-d.dddddd>,sigma(m));
m:= m+1; if j=i then j:= l
end;
m:= m+k-1; if i>l then m:= m+i-l
end i;
if l<N3 then goto rep;
close(sigma,true)
end msa;
noi:= 0; n:= S(1,2); k:= S(1,1)*N3; l:= (k+127)//128; pi:= 0;
for ibl:=1 step 1 until BL do begin
zone BI(l*128,1,stderror); open(BI,4,<:kbtrack:>,0);
setposition(BI,0,pi); inrec(BI,k);
noj:= noi; nj:= n; pj:= pi; pi:= pi+l;
for jbl:=ibl step 1 until BL do begin
array Z(1:n,1:nj);
zone BJ(l*128,1,stderror); open(BJ,4,<:kbtrack:>,0);
setposition(BJ,0,pj); inrec(BJ,k); pj:= pj+l;
for m:=1 step 1 until 3 do begin
ri1:= m mod 3+ 1; l:= 0;
ri2:= ri1 mod 3+ 1; W:= 0;
for i:=1 step 1 until n do begin rj1:= ri1-l; rj2:= ri2-l;
for j:=1 step 1 until nj do begin
V:= 0; for k:=N3-3 step -3 until 0 do
V:= V+ (BI(ri1+k)*BJ(rj2+k)-BI(ri2+k)*BJ(rj1+k))*Ma(k//3+1);
Z(i,j):= V; W:= W+ abs V;
rj1:= rj1+N3; rj2:= rj2+N3 end j;
ri1:= ri1+N3; ri2:= ri2+N3; l:= l+N3
end i;
if centrifugal then begin
k:= noi+n; l:= noj+nj;
for i:=noi+1 step 1 until k do begin
V:= f2(i);
for j:=(if ibl=jbl then i else noj)+1
step 1 until l do begin
h:= V-f2(j); if abs h >'-7 then begin
Z2:= Z(i-noi,j-noj)**2
*(case m of (Ia2,Ib2,Ic2))*5.1205987'-8;
h:= 4*V/h;
alfa(i,m,2):= alfa(i,m,2)-Z2*(h-1)/f(i);
alfa(j,m,2):= alfa(j,m,2)-Z2*(3-h)/f(j);
end end end i end;
if W>'-3 and zeta then begin
if linecount+n>54 then writepage else begin cr; cr end; cr;
write(res,<:Zeta-:>,false add (64+m),1,<: matrix, Species:>,
<<dd>,ibl,<: * Species:>,jbl);
k:= nj+5; l:= nj+1;
rep: cr; cr; write(res,<: :>);
k:= k-5; l:= if l>5 then l-5 else 1;
for j:=k step -1 until l do begin
write(res,<: Q:>,<<dd>,noj+nj+1-j) end; cr;
for i:=n step -1 until 1 do begin
cr; write(res,<: Q:>,<<dd>,noi+n+1-i,<: :>);
for j:=k step -1 until l do if abs Z(i,j)<5'-7 then
write(res,<: 0 :>) else
write(res,<<-dddd.dddddd>,Z(i,j))
end i;
if l>1 then begin
if linecount+n>57 then begin writepage; cr;
write(res,<:Zeta-:>,false add (64+m),1,<: matrix, Species:>,
<<dd>,ibl,<: * Species:>,jbl) end;
goto rep end end zeta-output;
end m;
noj:= noj+nj; if jbl<BL then begin
nj:= S(jbl+1,2); k:= S(jbl+1,1)*N3; l:= (k+127)//128 end;
close(BJ,true) end jbl;
noi:= noi+n; if ibl<BL then begin
n:= S(ibl+1,2); k:= S(ibl+1,1)*N3; l:= (k+127)//128 end;
if ibl=BL and ig=G then monitor(48,BI,0,tail);
close(BI,true) end ibl;
begin array GS(1:3,1:2);
for j:=1,2,3 do for k:=1,2 do GS(j,k):= 0;
m:= 0; W:= real << -dddd.dddddd>; h:= real << -dddddd.dddd>;
if linecount+S(1,2)*6>50 then writepage; cr; cr; write(res,<:
Corrections to Rotational Constants and Moments of Inertia.
gamma(anh.), (Mhz*cm), alfa(harm.), alfa(Cor.),
alfa(harm.+Cor.), (MHz), eps(harm.+Cor.), (u*Å**2).
:>); linecount:= linecount+4; cr;
write(res,<: Vib.:>);
for j:=1,2,3 do write(res,sp,13,false add (64+j),1);
write(res,sp,13,<:ID:>); cr; cr;
for ibl:=1 step 1 until BL do begin
n:= S(ibl,2); noi:= m+1; m:= m+n;
if linecount+n*6>59 and ibl>1 then writepage;
cr; write(res,<:Species:>,<<ddd>,ibl); cr;
for i:=noi step 1 until m do begin
l:= noi+m-i; linecount:= linecount+4; cr;
write(res,<:q:>,<<dd>,i,sp,7);
for j:=1,2,3 do
write(res,string h,my(i,j)*5.7566890'-4/f(l));
write(res,nl,1,sp,10);
for j:=1,2,3 do write(res,string h,alfa(i,j,1));
write(res,nl,1,<< dddd.dd>,1302.832*f(l));
for j:=1,2,3 do write(res,string h,alfa(l,j,2));
write(res,nl,1,sp,10);
for j:=1,2,3 do begin
V:= alfa(i,j,1):= alfa(i,j,1)+alfa(l,j,2);
GS(j,1):= GS(j,1)+V; write(res,string h,V)
end;
write(res,nl,1,sp,10); Z2:= 0;
for j:=1,2,3 do begin
V:= alfa(i,j,1)/(case j of(Ia2,Ib2,Ic2))*505376.42;
Z2:= Z2+(if j=3 then V else -V);
GS(j,2):= GS(j,2)+V; write(res,string W,V)
end;
if planar then write(res,string W,Z2); cr;
end i end ibl;
if linecount>56 then writepage; cr; cr;
write(res,<:Ground State Corrections (alfa/2,eps/2)::>);
cr; cr; write(res,sp,10);
for j:=1,2,3 do write(res,string h,GS(j,1)/2);
cr; write(res,sp,10);
for j:=1,2,3 do write(res,string W,GS(j,2)/2);
if planar then
write(res,string W,(GS(3,2)-GS(1,2)-GS(2,2))/2);
cr; cr;
end GS end ha end ig;
stop:
write(res,<:<10>:>,<:<25>:>); close(res,closeres);
if fclear then begin
open(res,4,<:ftrack:>,0); monitor(48,res,0,tail); close(res,true)
end;
open(res,4,<:xtrack:>,0); monitor(48,res,0,tail); close(res,true);
open(res,4,<:btrack:>,0); monitor(48,res,0,tail); close(res,true);
end end
▶EOF◀