⟦d22babf4a⟧

TextFile


<*      diagrampr generate and ALGOL procedure for a specific decay curve  *>
<*      redesigned 1981-02-17*>

procedure quicksort(lo,r,ns,idx,n,l,J,L,tr,tau);
value lo,r; integer lo,r;
integer array idx,n,l,J,L,tr;
array ns,tau;
begin
integer i,j,k,wi;
real x,w;
i:=lo; j:=r;
x:=ns((lo+r)//2);
while i<=j do begin
  while ns(i)<x do i:=i+1;
  while ns(j)>x do j:=j-1;
  if i<=j then begin
    w:=ns(i); ns(i):=ns(j); ns(j):=w;
    wi:=idx(i); idx(i):=idx(j); idx(j):=wi;
    wi:=n(i); n(i):=n(j); n(j):=wi;
    wi:=l(i); l(i):=l(j); l(j):=wi;
    wi:=J(i); J(i):=J(j); J(j):=wi;
    wi:=L(i); L(i):=L(j); L(j):=wi;
    wi:=tr(i); tr(i):=tr(j); tr(j):=wi;
    w:=tau(i); tau(i):=tau(j); tau(j):=w;
    i:=i+1; j:=j-1;
    end i<=j;
  end loop;
  if lo<j then quicksort(lo,j,ns,idx,n,l,J,L,tr,tau);
  if i<r then quicksort(i,r,ns,idx,n,l,J,L,tr,tau);
end quicksort;


procedure writediagramproc(out,
   bsname,dipval,states,series,dipindex,N0,curstate,delt);
value curstate,delt; integer curstate;
real delt;
comment delt is the time vindow in ns;
zone out;
integer array dipval,states,series,dipindex;
array bsname,N0;
begin
integer diag,diagram,lower,res,transcalc,transdrop,
  termstot,termsBcrit,statesincalc,i,j,id1,id2,
tottrans,trans;
real b,nstarlower,an,s0,tau,taumin,taulow;
integer array field dip;
array field dipr;
integer array nt,maxd,idx,nv,lv,Jv,Lv,trv(0:(maxstates-1)),
   tail(1:10);
array tauv,nsv(0:maxstates-1);
long array diagname(1:2);
array field raf;
algol list.off copy.statevar;

procedure generatediagramsum(curstate,diagram,max,allow,Av);
value curstate,diagram,max; integer curstate,diagram,max;
integer array allow;
array Av;
begin
integer array ntc,s(0:diagram);
integer i,j,k,p,high,last,nidx,qmin,
  trans,tottrans;
real cut,tau;
real A,B;
boolean more,nzero,first;
algol list.off copy.statevar;
algol list.off copy.statevar2;

procedure writediagram(m,low,diagram,curidx,Av);
value m,low,diagram; integer m,low,diagram;
integer array curidx;
array Av;
begin
integer i,j,k,p,q,chars;
real tau,tau2;
boolean nzero,Ball;
integer array Scoeff(1:diagram);
array Bcoeff(1:diagram);
nzero:=true;
if proc then
begin
  chars:=0;
  if first then
  begin
chars:=chars+write(out,<:
<*:>);
timing(out,<:diagrams:>,false);
if diagram>1 then
chars:=chars+write(out,"nl",1,<:written :>,termstot);
chars:=chars+write(out,"nl",1,<:dropped :>,termsBcrit);
chars:=chars+write(out,<:
*>
<10>:>,<<d>,<:if maxdiagram>=:>,diagram,
   <: then
begin
simproc:=r:=r:>);
outendcur(0);
if lookuptail(diagname,tail)=0 then
begin
  tail(6):=systime(7,0,0.0);
  tail(8):=termstot;
  tail(9):=termsBcrit;
  tail(10):=diagram;
  changetail(diagname,tail);
end;
end first;
  first:=false;
end else
if printdiagram then
begin
  n:=nv(m); l:=lv(m); J:=Jv(m); L:=Lv(m);
  chars:=write(out,<<d>,<:<10>+N0(:>);
  chars:=chars+writestate(out,L,n,l,J);
  chars:=chars+write(out,<:)*:>);
end;
A:=1;
for q:=1 step 1 until diagram-1 do begin
  i:=curidx(q);
  j:=curidx(q+1);
  A:=A*Av(i*(i+1)//2+j);
  if list then
  begin
    n:=nv(i); l:=lv(i); J:=Jv(i); L:=Lv(i);
    write(out,"nl",1,"*",1);
    writestate(out,L,n,l,J); write(out,"sp",1);
    n2:=nv(j); l2:=lv(j); J2:=Jv(j); L2:=Lv(j);
    writestate(out,L2,n2,l2,J2);
    write(out,A);
  end list;
  if printdiagram then begin
    n:=nv(i); l:=lv(i); J:=Jv(i); L:=Lv(i);
    chars:=chars+write(out,<<d>,<:A(:>);
    chars:=chars+writestate(out,L,n,l,J);
    chars:=chars+write(out,<:,:>);
    n2:=nv(j); l2:=lv(j); J2:=Jv(j); L2:=lv(j);
    chars:=chars+writestate(out,L2,n2,l2,J2);
    chars:=chars+write(out,<:):>);
    if chars>72-14 then chars:=write(out,"nl",1,"sp",5);
    chars:=chars+write(out,"*",1);
  end printdiagram;
  end q loop;
if printdiagram then chars:=chars+write(out,<<d>,<:(:>);
for j:=low step 1 until m do begin
for q:=1 step 1 until diagram do begin
  if curidx(q)=j then begin
  tau:=tauv(j);
  if j=lower then tau:=tau*tauf;
  if printdiagram then begin
    if chars>72-15 then chars:=write(out,"nl",1);
    chars:=chars+write(out,<<d>,<:+:>,if delt>0 then <:S*:> else <::>,<:exp(-a(:>);
    chars:=chars+writestate(out,L,n,l,J);
    chars:=chars+write(out,<:)*t):>);
  end printdiagram;
  B:=A;
  for k:=low step 1 until m do begin
  for p:=1 step 1 until diagram do begin
    if curidx(p)=k and k<>j then begin
      tau2:=tauv(k);
      if k=lower then tau2:=tau2*tauf;
      B:=B/(1/tau2-1/tau);
      if printdiagram then begin
        n:=nv(j); l:=lv(j); J:=Jv(j); L:=Lv(j);
        n2:=nv(k); l2:=lv(k); J2:=Jv(k); L2:=Lv(k);
        n:=nv(m);  l :=lv(m); J :=Jv(m); L :=Lv(m);
        if chars>72-15 then chars:=write(out,"nl",1,"sp",5);
        chars:=chars+write(out,<<d>,<:/(a(:>);
        chars:=chars+writestate(out,L2,n2,l2,J2);
        chars:=chars+write(out,<:)-a(:>);
        chars:=chars+writestate(out,L,n,l,J);
        chars:=chars+write(out,<:)):>);
      end printdiagram;
      p:=diagram;
      end curstate=k;
      end p;
      end k;
      Bcoeff(q):=B;
      Scoeff(q):=j;
    q:=diagram;
    end curstate=q;
    end q;
    end j;
if proc then
begin
  Ball:=false;
  for i:=1 step 1 until diagram do
  begin
    Ball:=Ball or Bcoeff(i)>Bmin;
  end diagram;
  if -,Ball then termsBcrit:=termsBcrit+1 else
    termstot:=termstot+1;
  if Ball  then
  begin
    chars:=write(out,"nl",1,<:+N:>,<<d>,m,<:*(:>);
    for i:=1 step 1 until diagram do
    begin
      B:=Bcoeff(i); j:=Scoeff(i);
      if chars>72-14 then chars:=write(out,"nl",1,"sp",3);
      chars:=chars+write(out,if B<0 then <:-:> else
        if i=1 then <::> else <:+:>,
        <:E:>,<<d>,j);
      if delt>0 and false then
      begin
      <*old version now this correction is multiplied directly on the exp*>
      an:=delt/tauv(j); <*tau in ns and delt in ns*>;
      B:=B*sinh(an)/an;
      end;
      if B<>1.0 then chars:=chars+write(out,<:*:>,<<d.dddddddd'd>,abs B);
    end diagram;
  chars:=chars+write(out,<:):>);
  end Ball ;
end proc;
if printdiagram then
  chars:=chars+write(out,<:):>);
end writediagram;

message start of diagramsum;
s(0):=lower-1;
s(1):=lower;
i:=1;
more:=true;
if diagram>1 then
begin
repeat i:=i+1;
  j:=allow(s(i-1),1);
  s(i):=j;
until i=diagram or j<0;
for i:=i+1 step 1 until diagram do s(i):=-1;
end diagram>1;
for i:=1 step 1 until diagram do ntc(i):=if s(i)>0 then 1 else -1;
nt(lower-1):=1;
first:=true;

if dsurvey and more then write(out,"nl",1,<:diagram sum:>,
   diagram,
   <: lower:>,lower,<: upper :>,s(diagram));

while  more do
begin
  high:=if s(diagram-1)>=0 then nt(s(diagram-1)) else -1;
  for nidx:=1 step 1 until high do 
  begin
     writediagram(s(diagram),curstate,diagram,s,Av);
     if nidx<high then s(diagram):=allow(s(diagram-1),nidx+1);
   end nidx;
  if dsurvey then write(out,"nl",1,<:survey:  :>,<: high:>,high);
  last:=diagram+1;
  while s(last-1)<0  and last>0 do last:=last-1;
  ntc(diagram):=high;
DIAGRAMERROR:
  repeat
    last:=last-1;
    ntc(last):=ntc(last)+1;
  until ntc(last)<=nt(s(last-1)) or last=1;
  if dsurvey then write(out,"nl",1,<:survey: last after adjust:>,last);
  if last>1 then
  begin
   if dsurvey then write(out,"nl",1,<:s(last):>,s(last));
   s(last):=allow(s(last-1),ntc(last));
  if dsurvey then write(out,"nl",1,<:new s(last):>,s(last));
   j:=last;
   repeat j:=j+1;
     i:=allow(s(j-1),1);
     ntc(j):=1;
     if i>=0 then s(j):=i;
  until i<0 or j=diagram;
  if i<0 then
  begin
    if dsurvey then write(out,"nl",1,<:diagram error:>,last,s(last-1));
    last:=j;
    goto DIAGRAMERROR;
  end;
  end last>=1;
  more:=last>1 and s(1)=1;
  if dsurvey then begin
    write(out,"nl",1,<:search last :>,last,
      <: more :>,if more then <:true:> else <:false:>,s(1));
    for i:=1 step 1 until diagram do
    begin
    k:=s(i);
    if k>0 then
    begin
    n:=nv(k); l:=lv(k); J:=Jv(k); L:=Lv(k);
    writestate(out,L,n,l,J);
    end k>0;
    end diagram;
    end test;
end more;
if proc   then write(out,<:;:>);
end generate sum;


termstot:=termsBcrit:=0;

raf:=0;
diagname(1):=long <:diagr:> add 'a';
diagname(2):=long <:m:>;
if lookuptail(diagname,tail)=0 then
begin
  tail(10):=0;
  changetail(diagname,tail);
end;
statesincalc:=-1;
timing(out,<:diagram generation start:>,true);
taulow:=gettau(states,curstate);
for i:=0 step 1 until maxstates-1 do
begin
  getstate(states,i,n,J,seriesindex,nstar,Ecm,L,0);
  computed:=getcomputed(states,i,0,0.0);
  tau:=gettau(states,i);
  if computed and nstar<=nsmax  and tau<2*tmax*taulow then
  begin
   statesincalc:=statesincalc+1;
   idx(statesincalc):=i;
   ser:=seriesindex*seriessize;
   l:=series.ser(1);
   nv(statesincalc):=n; lv(statesincalc):=l;
   Jv(statesincalc):=J; Lv(statesincalc):=L;
   nsv(statesincalc):=nstar; trv(statesincalc):=gettransno(states,i);
   tauv(statesincalc):=gettau(states,i);
  end;
end;
quicksort(1,statesincalc,nsv,idx,nv,lv,Jv,Lv,trv,tauv);

if dsurvey then
begin
write(out,"nl",3);
for i:=0 step 1 until statesincalc do
  write(out,"nl",1,
    <<dd>,<:state :>,i,<<dd.ddd>,<: n* :>,nsv(i));
end;

lower:=i:=-1;
for i:=i+1 while lower<0 and i<=statesincalc do
begin
  if idx(i)=curstate then lower:=i;
end;
if lower<0 then alarm(<:curstate not found:>);
if dsurvey or test then write(out,"nl",1,<:lower :>,lower);

begin <*local block for allow and Av*>
integer array allow(lower:statesincalc,1:2*nsmax);
array Av(1:(statesincalc+1)*statesincalc//2);

for id1:=lower step 1 until statesincalc do
begin
  nt(id1):=0;
  for id2:=2*nsmax step -1 until 1 do
  begin
    allow(id1,id2):=-1;
  end;
end;

for id1:=lower step 1 until statesincalc-1 do
for id2:=id1+1 step 1 until statesincalc do Av(id1*(id1+1)//2+id2):=0;

transcalc:=transdrop:=0;
for id1:=lower step 1 until statesincalc-1 do
begin
  i:=0;
  for id2:=id1+1 step 1 until statesincalc do
  begin
   if abs(lv(id1)-lv(id2))=2 then
   begin
     if abs(Jv(id1)-Jv(id2))<=2 then
     begin
     if abs(nsv(id1)-nsv(id2))>'-4 then
     begin
       if abs(Lv(id1)-Lv(id2))<=2 then
       begin
         tottrans:=dipindex(idx(id2));
         trans:=trv(id2);
         repeat
           dipr:=dip:=(tottrans-trans)*diprecsize;
           trans:=trans-1;
         until trans<=0 or idx(id1)=dipval.dip(2);

         if idx(id1)=dipval.dip(2) then
         begin
           if dipval.dipr(7)>brmin then
           begin
             transcalc:=transcalc+1;
             Av(id1*(id1+1)//2+id2):=dipval.dipr(5)/10;
             i:=i+1;
             allow(id1,i):=id2;
             nt(id1):=nt(id1)+1;
           end brmin else transdrop:=transdrop+1;
        end found;
       end allowed;
     end n*;
     end J;
   end l;
  end id2;
end id;

for id1:=lower step 1 until statesincalc do
begin
  maxd(id1):=0;
  j:=id1;
  repeat
    j:=allow(j,1);
    if j>0 then maxd(id1):=maxd(id1)+1;
  until j<0 or j>=statesincalc;
  if dsurvey then write (out,"nl",1,<:maxdiag :>,id1,maxd(id1));
end for id1;

if dsurvey then
begin
write(out,"nl",1,<:states connected:>);
for id1:=lower step 1 until statesincalc-1 do
begin
  write(out,"nl",1);
  writestate(out,Lv(id1),nv(id1),lv(id1),Jv(id1));
  write(out,<< dddd>,nt(id1),allow(id1,1));
end;
end test;

if proc and mindiagram>1 then
begin
  res:=connectcuro(<:simprocext:>);
  if res<>0 then unstackcuro;
end else
if proc then
begin
  res:=connectcuro(<:simproctxt:>);
  if res<>0 then unstackcuro;
  write(out,<:
;ali time 5 0
mode list.yes
lookup simproctxt
clear user simproc
simproc=set 25 disc3
scope user simproc
lookup rydlist
if ok.yes
mode 15.yes
simproc=algol details.8.8
:>);
writedate(out,0.0);
write(out,"nl",1);
writeatsym(out,S,atno,Z); write(out,"sp",4);
n:=nv(lower); l:=lv(lower); J:=Jv(lower); L:=Lv(lower);
writestate(out,L,n,l,J);
write(out,<:
external real procedure simproc(t,maxdiagram,N);
value t,maxdiagram; real t; integer maxdiagram;
array N;
if t<-.5 then simproc:=:>,curstate,<: else
begin
real :>);
for i:=lower step 1 until statesincalc do
begin
  write(out,"nl",if i mod 10 =7 then 1 else 0,
   <:E:>,<<d>,i,if i<10 then "sp" else false,1,
  <:,N:>,<<d>,i,if i<10 then "sp" else false,1,<:,:>);
end for i;
write(out,<:r;
r:=0;
comment maxdiagram= :>,maxdiagram,"nl",1);
writeatsym(out,S,atno,Z); write(out,"sp",4);
writestate(out,L,n,l,J);
if delt>0 then write(out,
   "nl",1,<:delta t(slit)           :>,delt,<: ns:>,
   "nl",1,<:n* max                  :>,nsmax,
   "nl",1,<:states in calculation   :>,statesincalc,
   "nl",1,<:transitions             :>,transcalc,
   "nl",1,<:transitions dropped     :>,transdrop,
   "nl",1,<:minimum branching ratio :>,brmin,
   "nl",1,<:minimum coefficient     :>,Bmin);
for i:=0 step 1 until statesincalc do
  begin
  write(out,"nl",if i mod 4=0 then 1 else 0,
            "sp",if i mod 4<>0 then 2 else 0);
  n:=nv(i); l:=lv(i); J:=Jv(i); L:=Lv(i); nstar:=nsv(i);
  writestate(out,L,n,l,J);
  write(out,"sp",1,<<dd.ddd>,nstar,"sp",1,
    <:=:>,<<ddd>,i,if i<statesincalc then <:,:> else <::>);
  end comment;
write(out,<:;:>);
for i:=lower step 1 until statesincalc do
begin
  tau:=tauv(i);
  if i=lower then tau:=tau*tauf;
  write(out,"nl",if i mod 1=0 then 1 else 0,
    "sp",if i mod 1<>0 then 1 else 0,
    <:E:>,<<d>,i,if i<10 then "sp" else false,1,
      <::=exp(-t/:>,<<ddd.dddddddd>,tau,
    <:):>);
    if delt>0 then write(out,<:*:>,<<d.ddddddd>,sinh(delt/tau)/delt*tau);
    write(out,<:; N:>,<<d>,i,if i<10 then "sp" else false,1,
    <::=N(:>,<<dd>,i,<:);:>);
end higher;
end proc;

for diag:=statesincalc-lower step -1 until lower do begin
diagram:=statesincalc-lower-diag+1;
if diagram<=maxdiagram and diagram>=mindiagram then begin
  if dsurvey or test then write(out,<<d>,<:<10><10>diagram :>,diagram);
  generatediagramsum(lower,diagram,statesincalc,allow,Av);
  if -,proc then outendcur(10);
  end;
end;
if proc then
begin
for i:=1,1 step 1 until maxdiagram do write(out,"nl",1,<:end:>);
timing(out,<:diagram:>,false);
write(out,<<  ddddd>,
  "nl",1,<:terms total        :>,termstot,
  "nl",1,<:terms dropped B    :>,termsBcrit,  <: B min   :>,<<d.d'-dd>,Bmin,
     "nl",1,<<d.dd'-dd>,<:brmin :>,brmin,
    <: tau max :>,<<  ddddd>,2*tmax,<: tau0:>,
  "nl",1,<:blocksread       :>,blocksread,
"nl",1,<:simproc;
end;
:>);
end;
if res=0 and proc then closeout;
if dsurvey or test then
write(out,
          "nl",1,<:blocksread          :>,blocksread);
timing(out,<:diagram end:>,false);
if proc then
write(out,"nl",2,<:maxdiagram           :>,maxdiagram,
          "nl",1,<:terms total          :>,termstot,
          "nl",1,<:terms dropped B      :>,termsBcrit,
          "nl",1,<:blocksread           :>,blocksread);
outendcur(10);
end local block for allow and Av;
end writediagramproc;
▶EOF◀
DataMuseum.dk

RC4000/8000/9000

⟦d22babf4a⟧ TextFile

Derivation

TextFile
