DataMuseum.dk

top - metrics - download

    Length: 6912 (0x1b00)
    Types: TextFile
    Names: »algtdosc«

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦7e928b248⟧ »algbib« 
            └─⟦this⟧ 

TextFile

;gosav
lookup banddiag
if ok.no
(banddiag= set 8
permanent banddiag.15
banddiag=algol extbanddiag index.no)
tdosc=set 60
tdosc=algol list.no xref.no index.no
\f





TWODIMENSIONAL OSCILLATOR , 24-8-77. GOS

begin
comment
     Program for calculating eigenvalues, eigenvectors and
expectation values for the Hamiltonian
   1/2 P (G0 + G2 X**2 + G4 X**4 + G6 X**6) P
   + 1/2 G0 (Pz**2 - 1/4) / X**2
           + V2 X**2 + V4 X**4 + V6 X**6.
Input specifications:
    1) Text string in < >.
    2) Order of matrix in harmonic basis, N.
    3) Energy spacing of basis, 1/cm.
    4) G0, G2, G4 and G6.
    5) V2, V4 and V6.
    6) Option for units: 0 , if G- and V-constants in 1/cm,
       1 , if G- and V-constants in terms of u, Å, rad and mdyn.
    7) Two vibrational quantum numbers indicating range
       of levels to be calculated.
    8) The maximal change in vibrational quantum number, dv,
       of calculated transitions (delta(l) = 0).
Energies are calculated in 1/cm. Expectation values in Å or rad,
momenta divided by h-bar, however;

zone res(128,1,stderror);
real lo,G,G2,G4,G6,A,B,C,E0,LS,h00,h01,h02,h03,h0l,h2l,h20,h21,h22,
     h40,h41,h60,Q2,Q4,Q6,P2,P2Q2,P2Q4,P4,P4Q2,p,q,r,s,t,
     sq1,sq2,sq3,sq4,cpu,time;
integer N,m1,m2,dv,nw,i,j,k,l,k1,k2,m,u,v,v1,v2,lines;
array head(1:12);
boolean even,sp,nl,closeres;
even:= true; sp:= false add 32; nl:= false add 10;
readhead(in,head,1);  lo:= real <<-dddd.dddd000000>;
read(in,N,LS,G,G2,G4,G6,A,B,C,nw,v1,v2,dv); i:= 1;
closeres:= outmedium(res);
write(res,<:<12>:>,nl,2,string head(increase(i)),
   nl,1,<:Matrix dimension:>,<<  ddd>,N,
   <:   Energy Spacing : :>,<<ddddd>,LS,nl,1,
   <:G0,G2,G4,G6: :>,string lo,G,G2,G4,G6,
   nl,1,<:V2,V4,V6   : :>,A,B,C);
p:= G; G:= LS; LS:= p/LS*(if nw=0 then 1 else 33.71498443);
G2:= G2*LS/p; G4:= G4*LS*LS/p; G6:= G6*LS**3/p;
p:= (if nw=0 then 1 else 50344.44349)/G;
A:= A*LS*p; B:= B*LS**2*p; C:= C*LS**3*p;
h00:= (G2*3/4+B)/2;      h01:= (1+G4*15/4+C*7)/2+A;
h02:= (G2+B*6)/4;        h03:= (G4+C*10)/4;
h0l:= -(G2+B+B)/4;       h2l:= -(G4+C*6)/8;
h20:= -1/2+G4*13/8+A+C*3;h21:= B+B; 
h22:= (G4+C*30)/8;       h40:= -G2/2+B;
h41:= -G4/2+C*3;         h60:= G4/2-C;
m:= if G6<>0 then 4 else if h60<>0 then 3 else 2;
if dv>-1 then begin
  nw:= 0; v:= v2-v1+2; v:= v//2*(v+1)//2;
end else v:= 1;
E0:= 0; lines:= 0;
cpu:= systime(1,0,time);

begin array w(1:v,1:10);

matrix:
v:= if even == v1 mod 2 = 0 then v1 else v1+1;
for l:=v step 2 until v2 do begin
m1:= if v1<l then 1 else (v1-l+2)//2;
m2:= (v2-l+2)//2;
begin
comment oprettelse af H-matrix.;
array H(1:N,0:m),ev(m1:m2),x(m1:m2,1:N);
sq1:= sq2:= sq3:= sq4:= 0;  u:= l+1; p:= l*l;
for i:=1 step 1 until N do begin 
   H(i,0):= ((h03*u+h02)*u+h01)*u+h00+((u+u)*h2l+h0l)*p
            +G6*((u*u*5-p*6+109)*u*u+(p-1)*(p-34))/16;
   H(i,1):= -((h22*(u-1)+h21)*(u-1)+h20+h2l*p
                +G6*((u-1)**2-p+33)*(u-1)/4)*sq1;
   H(i,2):= (h41*(u-2)+h40-G6*((u-2)**2-25/4)/2)*sq2;
   if m>2 then
   H(i,3):= (h60+G6*(u-3))*sq3;
   if m>3 then
   H(i,4):= -G6/2*sq4;
   u:= u+1; q:= sqrt(u*u-p)/2; u:= u+1;
   sq4:= sq3*q; sq3:= sq2*q; sq2:= sq1*q; sq1:= q;
end i;
banddiag(N,m,m1,m2,ev,H,x);
if l=0 then E0:= if m1=1 and dv=0 then ev(1)*G else 0;
k1:= m1-4;  v:= l+m1*2-4;
rep: lines:= lines+1; i:= 1;
if (lines-1) mod 5=0 and lines>1 then
write(res,<:<12>:>,nl,2,string head(increase(i)));
write(res,nl,2,sp,5);
k1:= k1+4; k2:= k1+3; if k2>m2 then k2:= m2;
for i:=k1 step 1 until k2 do
   write(res,sp,4,<:v,l =:>,<<ddd>,v+i+i,l,sp,1);
write(res,nl,2,<:E      :>);
for i:=k1 step 1 until k2 do
   write(res,string lo,G*ev(i)-E0);
if dv=-1 then begin
for j:=1 step 1 until N do  begin
   write(res,nl,1,sp,5);
   for i:=k1 step 1 until k2 do
      write(res,string lo,x(i,j))
end;
end else begin
for i:= k1 step 1 until k2 do begin
u:= l+1; Q4:= P4:= 0.5; P2Q2:= 0.75; t:=l*l;
Q2:= Q6:= P2:= P2Q4:= P4Q2:= p:= q:= r:= sq1:= sq2:= sq3:= 0;
for j:=1 step 1 until N do begin
   s:= r; r:= q; q:= p+p; p:= x(i,j);
   Q2:= Q2+u*p*p-sq1*p*q;
   P2:= P2+u*p*p+sq1*p*q;
   Q4:= Q4+(u*u*3-t)/2*p*p-(u-1)*2*sq1*p*q+sq1*sq2*p*r;
   P2Q2:= P2Q2+(u*u-t)/2*p*p-sq1*sq2*p*r;
   P4:= P4+(u*u*3-t)/2*p*p+(u-1)*2*sq1*p*q+sq1*sq2*p*r;
   Q6:= Q6+(u*u*5-3*t+7)*u/2*p*p-((u-2)*u*5-t+9)*0.75*sq1*p*q
          +(u-2)*3*sq1*sq2*p*r-sq1*sq2*sq3*p*s;
   P2Q4:= P2Q4+(u*u-t+7.5)*u/2*p*p-((u-2)*u-t+14)/4*sq1*p*q
          -(u-2)*sq1*sq2*p*r+sq1*sq2*sq3*p*s;
   P4Q2:= P4Q2+(u*u-t+3.5)*u/2*p*p+((u-2)*u-t+6)/4*sq1*p*q
          -(u-2)*sq1*sq2*p*r-sq1*sq2*sq3*p*s;
   u:= u+1; s:= (u*u-t)/4; u:= u+1;
   sq3:= sq2; sq2:= sq1; sq1:= sqrt(s)
end;
for j:=1 step 1 until 10 do w(nw+i,j):= case j of
   (G*ev(i),Q2*LS,Q4*LS*LS,Q6*LS**3,P2/LS,P2Q2,P2Q4*LS,
   P4/(LS*LS),P4Q2/LS,0.0 shift 24 add (v+i+i) shift 24 add l);
end i;
for j:=2 step 1 until 4 do begin
   write(res,nl,1,<:X:>,<<d>,j+j-2,sp,5);
   for i:=k1 step 1 until k2 do
      write(res,string lo,w(nw+i,j));
end;
for j:=5 step 1 until 7 do begin
   write(res,nl,1,case j-4 of (<:P2  :>,<:P2X2:>,<:P2X4:>),sp,3);
   for i:=k1 step 1 until k2 do
      write(res,string lo,w(nw+i,j));
end;
for j:=8,9 do begin
   write(res,nl,1,if j=8 then <:P4     :> else <:P4X2   :>);
   for i:=k1 step 1 until k2 do
      write(res,string lo,w(nw+i,j));
end end;
if k2<m2 then goto rep
end H;  nw:= nw+m2-m1+1 end l;
if even then begin even:= false; goto matrix end;
cpu:= systime(1,time,time)-cpu;

if dv>0 then begin
array tr(1:4);    i:= 1;
write(res,<:<12>:>,nl,2,string head(increase(i)),
      nl,2,<:Transitions:>,nl,2);
for k:= 1 step 1 until dv do
for l:=v1 step 1 until v2-k do
for v:= l step 1 until v2-k do begin
  m:= k mod 2;
  r:= 0.0 shift 24 add    v   shift 24 add l;
  rep:  even:= true;
  s:= 0.0 shift 24 add (v+k) shift 24 add (l+m);
  for i:=1 step 1 until nw do begin
    q:= w(i,10);
    if q=r or q=s then begin
    if even then begin
      for j:=1 step 1 until 4 do tr(j):= w(i,j);
      even:= false
    end else begin
      p:= if q=r then -1 else 1;
      write(res,nl,1,<<dd>,r shift (-24) extract 24,r extract 24,
            sp,2,s shift (-24) extract 24,s extract 24,
            <<   dddd.dddd>,(w(i,1)-tr(1))*p);
      for j:=2 step 1 until 4 do
      write(res,<<   -d.ddddddd>, (w(i,j)-tr(j))*p);
      i:= nw
  end end end i;
  if m=1 and l>v1 then begin
    m:= -1; goto rep
end end v end tr end w;
write(res,nl,2,<<ddd.d>,<:cpu-time: :>,cpu,
      <:    real-time: :>,time,nl,2,<:<25>:>);
close(res,closeres) end
permanent tdosc.15
▶EOF◀