DataMuseum.dk

top - metrics - download

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

└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
    └─⟦b2ec5d50f⟧ 
        └─⟦7b6e66aaa⟧ »crypr« 
            └─⟦this⟧ 
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦84e44a383⟧ »crypr« 
            └─⟦this⟧ 

TextFile

<* Wigner 6J and 9J symbols and fractional parantage coefficients
1979-02-16
*>

procedure init_factorial(f);
array f;
begin
integer i;
f(0):=1;
for i:=1 step 1 until 25 do
f(i):=f(i-1)*i;
end init_factorial;


real procedure sjs(j1,j2,j3,l1,l2,l3);
value j1,j2,j3,l1,l2,l3; integer j1,j2,j3,l1,l2,l3;
begin
integer w,wmin,wmax;
real omega;

real procedure delta(a,b,c);
value a,b,c; integer a,b,c;
begin
delta:=sqrt(factorial((a+b-c)//2)*
  factorial((a-b+c)//2)*
  factorial((b+c-a)//2)/
  factorial((a+b+c+2)//2));
end delta;

if j1+j2<j3 or abs(j1-j2)>j3 or j1+j2+j3<>
  2*((j1+j2+j3)//2) or
  j1+l2<j3 or abs(j1-l2)>l3 or j1+l2+l3<>
  2*((j1+l2+l3)//2) or
l1+j2<l3 or abs(l1-j2)>l3 or l1+l2+l3<>
  ((l1+j2+l2)//2) or
  l1+l2<j3 or abs(l1-l2)>j3 or l1+l2+j3<>
  ((l1+l2+j3)//2) then sjs:=0 else
begin
  omega:=0;
  wmin:=j1+j2+j3;
  if wmin<j1+l2+l3 then wmin:=j1+l2+l3;
  if wmin<l1+j2+l3 then wmin:=l1+j2+l3;
  if wmin<l1+l2+j3 then wmin:=l1+l2+j3;
  wmax:=j1+j2+l1+l2;
  if wmax>j2+j3+l2+l3 then wmax:=j2+j3+l2+l3;
  if wmax>j3+l1+l3+l1 then wmax:=j3+j1+l3+l1;

  for w:=wmin step 2 until wmax do
    omega:=omega+(if w=4*(w//4) then 1 else -1)
    *factorial(w//2+1)/(factorial((w-j1-j2-j3)//2)
    *factorial((w-j1-l2-l3)//2)
    *factorial((w-l1-j2-l3)//2)
    *factorial((w-l1-l2-j3)//2)
    *factorial((j1+j2+l1+l2-w)//2)
    *factorial((j2+j3+l2+l3-w)//2)
    *factorial((j3+j1+l3+l1-w)//2));
  sjs:=delta(j1,j2,j3)*delta(j1,l2,l3)*
       delta(l1,j2,l3)*delta(l1,l2,j3)*omega;
  end
end sjs;

real procedure njs(j11,j12,j13,j21,j22,j23,j31,j32,j33);
value j11,j12,j13,j21,j22,j23,j31,j32,j33;
integer j11,j12,j13,j21,j22,j23,j31,j32,j33;
begin
integer k,kmin,kmax;
real nj;

if j11+j12<j31 or abs(j11-j21)>j31 or j11+j21+j31
    <>2*((j11+j21+j31)//2) or
  j21+j22<j23 or abs(j21-j22)>j23 or j21+j22+j23
    <>2*((j21+j22+j23)//2) or
  j31+j32<j33 or abs(j31-j32)>j33 or j31+j32+j33
    <>2*((j31+j32+j33)//2) or
j11+j12<j13 or abs(j11-j12)>j13 or j11+j12+j13
    <>2*((j11+j12+j13)//2) or
  j12+j22<j32 or abs(j12-j22)>j32 or j12+j22+j32
    <>2*((j12+j22+j32)//2) or
  j13+j23<j33 or abs(j13-j23)>j33 or j13+j23+j33
    <>2*(j13+j23+j33)//2 then njs:=0 else
begin
  nj:=0;
  kmin:=abs(j21-j32);
  if kmin<abs(j11-j33) then kmin:=abs(j11-j33);
  if kmin<abs(j12-j23) then kmin:=abs(j12-j23);
  kmax:=j21+j32;
  if kmax>j11+j33 then kmax:=j11+j33;
  if kmax>j12+j23 then kmax:=j12+j23;
  for k:=kmin step 2 until kmax do
    nj:=nj+(if k=2*(k//2) then 1 else -1)*(k+1)*
      sjs(j11,j21,j31,j32,j33,k)*
      sjs(j12,j22,j32,j21,  k,j23)*
      sjs(j13,j23,j33,  k,j11,j12);
  njs:=nj;
end;
end njs;

real procedure vcc(j1,j2,j,m1,m2,m);
value j1,j2,j,m1,m2,m;
integer j1,j2,j,m1,m2,m;
begin
integer z,zmin,zmax;
real cc;
if m1+m2<>m or abs(m1)>abs(j1) or abs(m2)>abs(j2) or
   abs(m)>abs(j) or j>j1+j2 or j<abs(j1-j2) or
   j1+j2+j<>2*((j1+j2+j)//2) then vcc:=0 else
begin
  zmin:=0;
  if j-j2+m1<0 then zmin:=-j+j1+m2;
  zmax:=j1+j2-j;
  if j2+m2-zmax<0 then zmax:=j2+m2;
  if j1-m1-zmax<0 then zmax:=j1-m1;
  cc:=0;
  for z:=zmin step 1 until zmax do
    cc:=cc+(if z=4*(z//4) then 1 else -1)/
      factorial(z//2)*
      factorial((j1+j2-j-z)//2)*
      factorial((j1-m1-z)//2)*
      factorial((j2+m2-z)//2)*
      factorial((j-j2+m1+z)//2)*
      factorial((j-j1-m2+z)//2);
  vcc:=sqrt((j+1)*factorial((j1+j2-j)//2)*
    factorial((j1-j2+j)//2)*
    factorial((j2+j-j1)//2)*factorial((j1+m2)//2)*
    factorial((j1-m1)//2)*factorial((j2+m2)//2)*
    factorial((j2-m2)//2)*factorial((j+m)//2)*
    factorial((j-m)//2)*factorial((j1+j2+j+2)//2))*cc;
  end;
end vcc;

real procedure fracpar_p(N,S,L,S1,L1);
value N,S,L,S1,L1; integer N,S,L,S1,L1;
begin
<*calculates fractional parantage coefficients for p
configurations.
N number of p electrons
S1L1 quantum numbers for N-1 electron core
SL total spin and angular momentum
Sobelman p. 104-106
*>
fracparp:=0;
if N>=1 and N<=6 and S>=0 and S<=3 and L<=6 and
S1<=3 and L1<=6 then
begin
case N of begin
  <*N=1*> if S=0 and L=0 and S1=0 and L1=0 then fracparp:=1;
  <*N=2 Table 18.*> if S1=1 and L1=2 then
  begin
  if S=0 and (L=0 or L=4) then fracparp:=1 else
  if (S=2 and L=2) then fracparp:=1;
  end;
  <*N=3 Table 19.*> 
  if S=1 and L=2 and S1=0 and L1=0 then fracparp:=sqrt(2)/3 else
  if S=3 and L=0 and L1=2 and S1=2 then fracparp:=1 else
  if S=1 and S1=2 and L>0 and L1>0 then
  begin
  fracparp:=case L//2+L1-2 of (-1/sqrt(2),
     1/sqrt(2),-sqrt(5/18),-sqrt(1/2));
  end;
  <*N=4 Table 20.*>
  if S1=3 and L1=0 and S=2 and L=2 then fracparp:=-sqrt(1/3) else
  if S=0 and L=0 and S1=1 and L1=2 then fracparp:=1 else
  if L>0 and L1>0 and S1=1 then
  begin
  if (S=2 and L=2) or (S=0 and L=4) then
    begin
    fracparpr:=case L//2+L1-2 of(-1/2,-1/2,sqrt(5/12),sqrt(3/4));
    end;
  end;
  <*N=5*> if S=1 and L=2 then
  begin
  if S1=0 and L1=0 then fracparp:=1/sqrt(15) else
  if S1=2 and L1=2 then fracparp:=sqrt(3/5) else
  if S1=0 and L1=4 then fracparp:=sqrt(1/3);
  end;
  <*N=6*>
  if S=0 and L=0 then fracparp:=1;
  end case N;
end;
end fracparp;

▶EOF◀