DataMuseum.dk

top - metrics - download

    Length: 4608 (0x1200)
    Types: TextFile
    Names: »algsymosc«

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

TextFile

;kemlab5 1
lookup banddiag
if ok.no
(banddiag= set 8
banddiag=algol extbanddiag index.no)
symosc=set 50
symosc=algol list.no index.no
\f





SYMMETRIC OSCILLATOR af 4-6-73. GOS

begin
comment
     Program til beregning af egenværdier og egenvektorer
svarende til Hamiltonfunktionen
   1/2 (G0 + G2 X**2 + G4 X**4) P**2
           + V2 X**2 + V4 X**4 + V6 X**6.
Som input kræves:
    1) Matrixordenen for harmonisk basis, N.
    2) Scaling Constant in 1/cm.
    3) G0, G2 og G4 i enhedssystemet u,Å,rad.
    4) V2, V4 og V6 i 1/cm.
    5) Numrene på laveste og højeste egenværdi, der
       ønskes beregnet, m1 og m2.
    6) Et tal = 0, hvis forventninsværdier ønskes beregnet,
       og = 1, hvis egenvektorer ønskes udskrevet.
Egenværdierne beregnes i 1/cm, forventninsværdier i Å;

zone res(128,1,stderror);
real G,G2,G4,A,B,C,D,L,Q2,Q4,Q6,p,q,r,s,sq1,sq2,sq3;
integer N,m1,m2,outvector,i,j,k,k1,k2,m,n;
array head(1:12);
boolean even,sp,nl,closeres,pqq;
even:= true; sp:= false add 32; nl:= false add 10;
readhead(in,head,1);
read(in,N,D,G,G2,G4,A,B,C,m1,m2,outvector); i:= 1;
pqq:= N<0; N:= abs N;
closeres:= outmedium(res); p:= 1493.236445;
write(res,<:<12>:>,nl,3,string head(increase(i)),
   nl,1,<:Matrix dimension:>,<<  ddd>,N,
   <:    Scaling: g =:>,<<dddd>,D,nl,1,
   <:G0,G2,G4: :>,<<  -d.ddddddd000>,G,G2,G4,nl,1,
   <:V2,V4,V6: :>,<<  -dddd.dddd000>,A,B,C);
L:= G*33.71498443/D; G2:= G2/G*L; G4:= G4/G*L*L; G:= D;
A:= A/D*L; B:= B/D*L*L; C:= C/D*L**3;
if G4<>0 then A:= A+G2*G2*0.125;
m:= if G4<>0 or C<>0 then 3 else 2;
begin
comment oprettelse af H-matrix.;
array H(1:N,0:m),ev(m1:m2),x(m1:m2,1:N),
      xe(m1:m2,1:if pqq then N else 1);
matrix:
sq1:= sq2:= sq3:= 0;  p:= A+0.5; q:= (A-0.5)/2; 
r:= (G2+B*6)/4; s:= (B+B-G2)/8;
D:= (G4+C*10)/4; Q2:= (G4+C*30)/16;
Q4:= (-G4+C*6)/8; Q6:= (-G4+C*2)/16;
n:= if even then -2 else -1;
for i:=1 step 1 until N do begin 
   n:=n+2; k:= n*(n+1);
   H(i,0):= (n+0.5)*(p+(k+1.5)*D) + (k+0.5)*r;
   H(i,1):= (q+(n-0.5)*B + ((n-1)*n+1)*Q2)*sq1;
   H(i,2):= (s+(n-1.5)*Q4)*sq1*sq2;
   if Q6<>0 then
   H(i,3):= Q6*sq1*sq2*sq3;
   sq3:= sq2;  sq2:= sq1;  sq1:= sqrt((n+1)*(n+2))
end i;
i:=1;  if -,even and (outvector>0 or (m2-m1)>6) then
write(res,<:<12>:>,nl,2,string head(increase(i)),nl,1)
else write(res,nl,3);
if even then write(res,<:Even eigenfunctions and -values:>) else
write(res,<:Odd eigenfunctions and -values:>);
write(res,nl,1,<:g,g2,g4: :>,<< -dddd.ddddd>,G,
         <<  -d.dddddddd>,G2,G4,nl,1,<:A, B, C: :>,A,B,C);
if m1=0 then goto ud;
banddiag(N,if Q6<>0 then 3 else 2,m1,m2,ev,H,x);
if outvector>0 then begin
k1:= m1-4;  n:= if even then -2 else -1;
rep:  write(res,nl,2);
k1:= k1+4; k2:= k1+3; if k2>m2 then k2:= m2;
for i:=k1 step 1 until k2 do
   write(res,sp,6,<:n = :>,<<dd>,n+i+i,sp,6);
write(res,nl,2);
for i:=k1 step 1 until k2 do
   write(res,<<  -d.dddddddddd'-d>,ev(i));
for j:=1 step 1 until N do  begin
   write(res,nl,1);
   for i:=k1 step 1 until k2 do
      write(res,<<    -d.dddddddd>,x(i,j),sp,3)
end;
if k2<m2 then goto rep
end else begin
for i:= m1 step 1 until m2 do begin
write(res,nl,2,<<-dddd.ddddddd000>,G*ev(i));
Q2:= 0.5; Q4:= 0.75; Q6:= 1.875;
p:= q:= r:= sq1:= sq2:= sq3:= 0;
n:= if even then -2 else -1;
for j:=1 step 1 until N do begin
   s:= r; r:= q; q:= p; p:= x(i,j); n:= n+2;
   Q2:= Q2+n*p*p+sq1*p*q;
   Q4:= Q4+1.5*(n+1)*n*p*p+(n+n-1)*sq1*p*q+0.5*sq1*sq2*p*r;
   Q6:= Q6+((n+1.5)*n*2.5+5)*n*p*p+((n-1)*n+1)*3.75*sq1*p*q
          +(n-1.5)*1.5*sq1*sq2*p*r+0.25*sq1*sq2*sq3*p*s;
   sq3:= sq2; sq2:= sq1; sq1:= sqrt((n+1)*(n+2))
end;
write(res,sp,14,<<      d.dddddddd000>,Q2*L,nl,1,sp,30,
          Q4*L**2,nl,1,sp,30,Q6*L**3);
end end;
ud:  if even then begin 
   if pqq then begin
      for i:=m1 step 1 until m2 do
      for j:=1  step 1 until N  do xe(i,j):= x(i,j)
   end;
   even:= false; goto matrix end;
if pqq then begin
   A:= sqrt(L); write(res,nl,2,<:
Matricen <v(even)! pqq+qqp !v(odd)>*(2pi/ih) , Å::>,nl,1);
   for m:=m1 step 1 until m2 do begin
      q:= 0; r:= xe(m,1); s:= xe(m,2);
      for n:=1 step 1 until N do begin
         p:= q; q:= r; r:= s; s:= if n+2<=N then xe(m,n+2) else 0;
         xe(m,n):= -p*sqrt((n+n-1)*(n-1)*(n+n-3))
                   -q*(n+n-1)*sqrt(n-0.5) + r*(n+n)*sqrt(n)
                   +s*sqrt((n+1)*(n+n+1)*(n+n));
      end n;
      write(res,nl,2);
      for n:=m1 step 1 until m2 do
         write(res,<<  -d.ddddddddd>,sum(xe(m,i)*x(n,i),i,1,N)*A);
   end m
end pqq;
write(res,<:<25>:>); close(res,closeres);
end end;
permanent symosc.15
▶EOF◀