|
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: 9984 (0x2700)
Types: TextFile
Names: »algqrotfit«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦7e928b248⟧ »algbib«
└─⟦this⟧
;kemlab5 1
t=algol index.no
begin
comment INPUTREGLER:
1) En overskrift i < >.
2) Maximal regnetid i min.
3) ID, B og C (ID<>0).
4) Antal overgange.
5) Antal frekvenser.
7) Startnummer for kombinering af frekvenserne.
Første gang 1, ellers skriver programmet
startnummer for fortsat beregning.
8) Maximal rms for udskrift af rotationskonst.
9) Frekvenserne.
10) Overgangene (ingen dipol-kontrol).
For hver overgang angives Kvantetal, samt
antallet af frekvenser, der kan tænkes
tilordnet overgangen.
De mulige frekvenser specificeres efter
inputrækkefølgen;
integer i,m,N1,ncr,xx,yy,zz,k,prod1;
real smax,starttime,stoptime;
boolean oblat;
array ROT(1:3),head(1:12);
systime(1,0,starttime);
readhead(in,head,1); read(in,stoptime,ROT,N1,m,ncr,smax);
stoptime:=stoptime*60;
oblat:=2*ROT(2)>(1/(1/ROT(3)-1/ROT(2)-ROT(1)/505376)+ROT(3));
i:=1; write(out,<:<12><10><10><10>:>,string head(increase(i)),
<:<10><10>Calculation by Q-ROTFIT of 10-12-1975<10>:>);
begin
integer array qn1,qn2,ntr,qntr(1:N1),itr(1:40);
array tr(1:m);
begin
integer J1,k11,k12,J2,k21,k22,m1,type1,type2;
integer array MATRIX(1:8);
boolean b11,b12,b21,b22;
if oblat then begin
for i:=1 step 1 until 8 do
MATRIX(i):= case i of (1,4,2,3,2,3,1,4);
xx:=1; yy:=2; zz:=3
end else begin
for i:=1 step 1 until 8 do
MATRIX(i):= case i of (1,2,3,4,2,1,4,3);
zz:=1; xx:=2; yy:=3
end;
write(out,<:<10><10>Input frequencies:<10>:>);
for i:=1 step 1 until m do begin
read(in,tr(i));
write(out,<:<10>:>,<< dd>,i,<< dd ddd.ddd>,tr(i));
end;
m1:=0; write(out,<:<10><10>Input transitions:<10>:>);
for k:=1 step 1 until N1 do begin
read(in,J1,k11,k12,J2,k21,k22,m);
write(out,<:<10>:>,<<ddd>,J1,k11,k12,J2,k21,k22);
for i:=1 step 1 until m do begin
read(in,itr(m1+i));
write(out,<:<10>:>,false add 32,18,<<ddd>,itr(m1+i));
end i;
m1:=m1+m;
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;
qn1(k):=J1 shift 7 add ((J1-k11)//2+1) shift 7 add
MATRIX(if J1=J1//2*2 then type1 else (type1+4));
qn2(k):=J2 shift 7 add ((J2-k21)//2+1) shift 7 add
MATRIX(if J2=J2//2*2 then type2 else (type2+4));
qntr(k):=m;
end k;
end reading of transitions;
begin
integer line,matr,nr,JJ,it,j0,k0,K2,N,j,prod,
M,M0,cn,pr;
real norm,lambda,Li,dE,F,G,H,Pxy,Pz,y,s,so,p,qd,
X,Y,Z,time,a,b,c,an12;
array atot(1:N1,1:3),ai(1:2),rot(1:3),an(1:2),PG(1:3);
prod:=1; M:=if N1=2 then 1 else N1-2;
for i:=1 step 1 until N1-1 do prod:=prod*qntr(i);
systime(1,starttime,time); prod1:=prod*qntr(N1)+1;
write(out,<:<10><10>Start values:<10><10>ID=:>,<<-d.dddd>,
ROT(1),<:<10> B=:>,<< dd ddd.dddd>,ROT(2),<:<10> C=:>,ROT(3),
<:<10><10>Rotational Constants are printed when rms < :>,
<<dd.dd>,smax,<:<10>First combination::>,<< d>,ncr);
ROT(1):=ROT(1)/505376;
nycomb:
cn:=ncr-1; pr:=prod;
for i:=N1 step -1 until 2 do begin
k:=cn//pr; ntr(i):=k+1; cn:=cn-k*pr; pr:=pr/qntr(i-1)
end;
ntr(1):=cn+1; j:=0;
for i:=1 step 1 until N1 do begin
ntr(i):=itr(ntr(i)+j); j:=j+qntr(i)
end;
for i:=1 step 1 until N1 do begin
cn:=ntr(i);
for pr:=i+1 step 1 until N1 do
if cn=ntr(pr) and tr(cn)<>0 then begin
ncr:=ncr+1; goto nycomb
end;
end;
for i:=2,3 do rot(i):=ROT(i); ncr:=ncr+1;
if ncr>prod1 then goto stop; so:=1'612;
igen:
a:=rot(1):=1/(1/rot(3)-1/rot(2)-ROT(1));
b:=rot(2); c:=rot(3); b:=a*a/(b*b); c:=a*a/(c*c);
s:=0; X:=rot(xx); Y:=rot(yy); Z:=rot(zz);
y:=Y; Y:=(X-y)*0.25; X:=(X+y)*0.5; Z:=Z-X; M0:=M;
for line:=1 step 1 until N1 do begin
dE:= 0; Li:=tr(ntr(line));
if Li>0 then begin
for i:=1,2 do ai(i):=0;
for k:=-1,1 do begin
if k=-1 then begin
N:=qn1(line); j0:=N shift (-14) extract 7;
nr:=N shift (-7) extract 7; matr:=N extract 7
end else begin
N:=qn2(line); j0:=N shift (-14) extract 7;
nr:=N shift (-7) extract 7; matr:=N extract 7;
end;
JJ:= j0*(j0+1); G:=X*JJ;
k0:= if matr=1 then -2 else if matr=2 then 0 else -1;
N:= (j0-k0)//2; p:= if matr=1 then 2 else 1;
if j0<>0 then begin
array alfa,beta,S,Pdif(1:N);
beta(1):= Pdif(1):= 0;
k0:= k0+2; K2:= k0*k0;
for j:=1 step 1 until N do begin
alfa(j):= G + Z*K2;
if j<>N then begin
p:= p*(j0-k0)*(j0-k0-1)*(j0+k0+1)*(j0+k0+2);
k0:= k0+2; K2:= k0*k0;
beta(j+1):=Y**2*p;
Pdif(j+1):= -sqrt(p);
p:=1;
end j<>N;
end j;
if matr>=3 then begin
j:= if matr=3 then -1 else 1;
alfa(1):= alfa(1) + Y*j*JJ;
Pdif(1):= j*JJ//2
end;
j:=if oblat then nr else N+1-nr;
H:= (if j>1 then sqrt(beta(j)) else 0) +
(if j<N then sqrt(beta(j+1)) else 0);
if j=1 or j=N then H:= H+H;
G:= alfa(j) + H; H:= alfa(j) - H;
for lambda:= (G+H)*0.5 while G<>lambda do begin
p:= 0; qd:= 1; it:= 0;
for i:=1 step 1 until N do begin
y:= (alfa(i)-lambda)*qd - beta(i)*p;
p:= qd; qd:= y;
if p>=0 == qd>=0 then it:=it+1;
if it=nr then i:= N
end sturm;
if it=nr then H:= lambda else G:= lambda
end bisection;
for i:=1 step 1 until N do
alfa(i):=alfa(i)-lambda;
for i:=2 step 1 until N do
beta(i-1):=-sqrt(beta(i));
S(1):=S(N):=1;
if j>1 then S(2):=-alfa(1)/beta(1);
i:=1; for i:=i+1 while i<j do S(i+1):=
-(S(i)*alfa(i)+S(i-1)*beta(i-1))/beta(i);
y:=1/S(j);
for i:=j-1 step -1 until 1 do S(i):=S(i)*y;
if j<N then S(N-1):=-alfa(N)/beta(N-1);
i:=N; for i:=i-1 while i>j do S(i-1):=
-(S(i)*alfa(i)+S(i+1)*beta(i))/beta(i-1);
y:=1/S(j); S(j):=1;
for i:=j+1 step 1 until N do S(i):=S(i)*y;
norm:=Pz:=G:=F:=Pxy:=0;
for i:=N step -1 until 1 do begin
G:=S(i); H:=G*G; F:=F*G;
norm:=norm+H; H:=H*K2;
Pz:=Pz+H; Pxy:=Pxy+F; F:=Pdif(i)*G;
k0:=k0-2; K2:=k0*k0;
end i;
F:=F*G; norm:=1/norm; Pxy:=(Pxy+F)*norm;
PG(zz):=Pz:=Pz*norm; F:=JJ-Pz;
PG(xx):=(F+Pxy)*0.5; PG(yy):=(F-Pxy)*0.5;
dE:= dE + k*lambda;
PG(3):=PG(3)+PG(1)*c; PG(2):=PG(2)-PG(1)*b;
for i:=1,2 do ai(i):=ai(i)+PG(i+1)*k;
end j0<>0;
end k;
dE:=Li-dE; s:=s+dE*dE;
for i:=1,2 do atot(line,i):=ai(i); atot(line,3):=dE;
end else begin
M0:=M0-1;
for i:=1,2,3 do atot(line,i):=0;
end Li=0;
end line;
if M0<1 then goto nycomb;
for j:=1,2 do begin
if j=2 then begin
qd:= an(1);
p:= 0;
for i:=1 step 1 until N1 do
p:= p + atot(i,2)*atot(i,1);
an12:= p; p:= qd*p;
for i:=1 step 1 until N1 do
atot(i,2):= atot(i,2)-atot(i,1)*p
end j=2;
p:= qd:= 0;
for i:=1 step 1 until N1 do begin
p:= p + atot(i,j)**2;
qd:= qd + atot(i,j)*atot(i,3)
end i;
an(j):= 1/p; ai(j):= qd;
end j;
ai(2):=an(2)*ai(2);
ai(1):=an(1)*(ai(1)-ai(2)*an12);
y:=systime(1,starttime,p); qd:=p-time; time:=p;
if p+qd>stoptime then goto stop;
s:=sqrt(s/M0);
if so-s<1'-2 then begin
if so<smax then begin
write(out,<:<10>:>);
for i:=1 step 1 until N1 do write(out,<< dd>,ntr(i));
write(out,<:<10><10>:>,<<___dd.dd>,so);
for i:=1,2,3 do write(out,<< ddd ddd.ddd>,rot(i));
so:=(rot(2)*2-rot(1)-rot(3))/(rot(1)-rot(3));
write(out,<: kappa =:>,<< -d.ddd>,so,<:<10>:>);
end;
goto nycomb;
end else begin
for i:=2,3 do rot(i):=rot(i)+ai(i-1);
end;
so:=s;
goto igen;
end;
stop:
if ncr<prod1 then write(out,<:<10><10>Start new calculation:>,
<: with combination no.:>,<<ddd>,ncr) else
write(out,<:<10><10>Calculation complete.:>);
end;
end;
rename t.qrotfit
permanent qrotfit.17
▶EOF◀