DataMuseum.dk

top - metrics - download

    Length: 5376 (0x1500)
    Types: TextFile
    Names: »algmathieu«

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

TextFile

;kemlab5 1
r=algol index.no
\f




MATHIEU-matrix diagonalisering.         8-3-1971.   GOS.

begin
comment input reglerne er:
   1) Barrierens multiplicitet,  N.
   2) En reduceret barrierehøjde, her defineret
         B = V*Ired/(2*hbar**2)= N*N*s/16.
   3) Symmetriparameteren sigma, 0<=sigma<=N//2.
   4) Mindste og største værdi for vibrationskvantetallene, vmin og vmax.
   5) Et heltal, t (normalt = 1), der bestemmer trunkeringen;
integer N, n, i, v, vmin, vmax, d, k, k1, l, m, s, t;
real A1, A2, B, b, y, q, q1, p, p1, r, a, a1, g, h, lambda, norm;
boolean even, wang, vlige, closeres;
zone res(128,1,stderror);
procedure writecr;
   begin write(res,<:<10>:>) end;
closeres:= outmedium(res);
newNB:
read(in,N); if N<1 then goto stop;
read(in,B); write(res,<:
N = :>,<<ddd>,N,<:    B = :>,<<ddd.ddddd00000>,B,
<:     s = :>, 16*B/(N*N));
newvd:
read(in,d); writecr; writecr;
if d<0 then begin writecr; writecr; writecr; goto newNB end;
read(in,vmin,vmax,t);
write(res,<:sigma, vmin, vmax og t = :>,<<dddd>,d,vmin,vmax,t);
even:= d=0;  wang:= even or (d+d)=N;
b:= 1/(if wang then 4*B else B);  lambda:= -2;
t:= t*(if wang then 2*N else N);

for v:=vmin step 1 until vmax do begin
h:= lambda;  g:= 4*v*v/B+1;
if wang then begin
   vlige:= v//2*2=v;
   k1:= if even and vlige then 0 else if -,even then N else 2*N;
   A1:= b*k1*k1+ (if even then 0 else if vlige then -1 else 1);
   A2:= if k1=0 then 2 else 1;  n:=  v//2;  s:= k1+t+N*4
end else begin
   s:= t+N+N;  n:= v end;

for lambda:= (g+h)*0.5 while g<>lambda do begin
if wang then begin
   q:= A1-lambda; p:= A2;  k:= k1;  m:= if q<0 then 1 else 0;
   for k:= k+N+N while m<=n and k<=s do begin
      y:= (b*k*k-lambda)*q- p;
      q1:= p; p:= q; q:= y;
      if q>=0 == p<0 then begin m:= m+1; if s<k+t then s:= k+t end;
      if q*q1< p*p and s<k+t then s:= k+t
end end else begin
   q:= d*d*b-lambda;  p:= p1:= 1;  r:= 0;
   m:= if q<0 then 1 else 0;     k:= 0;
   for k:= k+N while m<=n and k<=s do begin q1:= r;
   for l:= -k+d,k+d do begin
      a:= l*l*b-lambda;
      y:= q*a-p;  p:= p1;  p1:= q;  q:= y;
      y:= p*a-r;  r:= p;   p:= y;
      if q>=0 == p1<0 then begin m:= m+1; if s<k+t then s:= k+t end
   end;
   if (q*q1< r*r or (q>=0 == r<0)) and s<k+t then s:= k+t
end end;
if m>n then g:= lambda else h:= lambda
end bisektion;

y:= lambda*0.25+0.5;   writecr;
m:= if wang then ((s-k1)//(N+N)+1) else s//N;
write(res,<:v = :>,<<dd>,v,<:  orden: :>,<<dd>,m);
write(res,<:  lambda = :>,<<-dd.ddd ddd ddd 000>, y);
write(res,<:  b = :>,<<-ddd.ddd ddd dd0 00>, y*16*B/(N*N));
if wang then begin
array U(1:m);
U(m):= 1;  U(m-1):= b*s*s-lambda;  n:= 1;  norm:= 1;
k:= s-N-N;     a:= b*k*k-lambda;
k:= k-N-N;     a1:= b*k*k-lambda;
for i:=m-2 step -1 until n do begin
   U(i):= U(i+1)*a- U(i+2);
   a:= a1;  k:= k-N-N;  a1:= b*k*k-lambda;
   if abs a< abs a1 then n:= i
end;
if n=1 then begin
   U(1):= U(1)/sqrt(A2);
   for i:=m-1 step -1 until 1 do norm:= norm+U(i)**2
end else begin
   g:= 1/U(n);
   for i:=n+1 step 1 until m do begin
      y:= U(i):= U(i)*g;  norm:= norm+y*y end;
   U(1):= A2;  U(2):= A1-lambda;  k:= k1;
   for i:=3 step 1 until n do begin
      k:= k+N+N;  U(i):= (b*k*k-lambda)*U(i-1) - U(i-2) end;
   U(1):= sqrt(A2);  g:= 1/U(n);  U(n):= 1;
   for i:=n-1 step -1 until 1 do begin
      y:= U(i):= U(i)*g;  norm:= norm+y*y end
end;
A1:= U(1)*U(2); g:= 1/norm; 
if even then begin
   if vlige then A1:= A1*sqrt(2) end
else A1:= A1+(if vlige then 0.5 else -0.5)*U(1)**2; 
for i:= 3 step 1 until m do A1:= A1+U(i-1)*U(i); 
write(res,<:<10>         cosN  = :>,<<-dd.dddddddd000>,A1*g); 
A1:= U(1)*U(3); 
if even then begin
   if vlige then A1:= A1*sqrt(2) + U(2)**2*0.5
            else A1:= A1 - U(1)**2*0.5 end
else A1:= A1 + (if vlige then 1 else -1)*U(1)*U(2);
for i:=4 step 1 until m do A1:= A1 + U(i-2)*U(i);
write(res,<:<10>         cos2N = :>,<<-dd.dddddddd000>,A1*g);
g:= sqrt(g);  k:= -3;  l:= 0;  writecr;
rep:
writecr;  k:= k+4;  l:= l+4;  if l>m then l:= m;
for i:=k step 1 until l do write(res,<<___-d.ddd ddd'-dd>, U(i)*g);
if l<m then goto rep
end egenvektor_i_wang-basis else begin
array U(-m:m);
n:= -m;  k:= -s+d;  U(-m):= b*k*k-lambda;  y:= abs U(-m);
for i:=-m+1 step 1 until m do begin
   k:= k+N;  U(i):= b*k*k-lambda;  a:= abs U(i);
   if a<y then begin y:= a;  n:= i end
end;
a:= U(m-1);  U(m-1):= U(m);  U(m):= norm:= 1;
for i:= m-2 step -1 until n do begin
   a1:= U(i);  U(i):= U(i+1)*a- U(i+2);  a:= a1 end;
g:= 1/U(n);
for i:=n+1 step 1 until m do begin
   y:= U(i):= U(i)*g;  norm:= norm+y*y end;
a:= U(-m+1);  U(-m+1):= U(-m);  U(-m):= 1;
for i:= -m+2 step 1 until n do begin
   a1:= U(i);  U(i):= U(i-1)*a- U(i-2);  a:= a1 end;
g:= 1/U(n);  U(n):= 1;
for i:=n-1 step -1 until -m do begin
   y:= U(i):= U(i)*g;  norm:= norm+y*y end;
g:= 1/norm;  A1:= 0;
for i:= -m+1 step 1 until m do A1:= A1+U(i-1)*U(i); 
write(res,<:<10>         cosN  = :>,<<-dd.dddddddd000>,A1*g); 
A1:= 0; 
for i:=-m+2 step 1 until m do A1:= A1 + U(i-2)*U(i);
write(res,<:<10>         cos2N = :>,<<-dd.dddddddd000>,A1*g);
g:= sqrt(g);  k:= -m;  l:= k+ (m-1) mod 4;
rep:
writecr;  if l>m then l:= m;
for i:=k step 1 until l do write(res,<<___-d.ddd ddd'-dd>, U(i)*g);
k:= l+1;  l:= if k=0 then 0 else k+3;
if k<=m then goto rep end egenvektor_for_E-level;
writecr end v;
goto newvd;
stop: write(res,<:<25>:>); close(res,closeres)  end

rename r.mathieu
permanent mathieu.15
▶EOF◀