DataMuseum.dk

top - metrics - download

    Length: 8448 (0x2100)
    Types: TextFile
    Names: »simcalc«

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

TextFile

<*                                  simcalc                       *>
<*calculation of simulated decay curves*>

procedure simcalc(out,dip,states,dipindex,N0,dt,tmin,tmax,curstate,points);
value dt,tmin,tmax,curstate,points;
real dt,tmin,tmax;
integer curstate,points;
zone out;
integer array dip,states,dipindex;
array N0;
begin
array Ni,dNi(0:maxstates-1),pvalue(1:maxdiagram);
integer array Ewritten(0:maxstates-1);
algol copy.statevar;
integer i,j,state2,dti,tottrans,intv;
real dn,lambda,expdecayv,s0,tau,x,dx,xmin,xmax;
integer array field dipi;
array field dipf;
ewritten(Ewritten);
for i:=0 step 1 until maxstates-1 do Ni(i):=N0(i);

if proc and -,generate then
begin
  write(out,"nl",2);
  intv:=0;
  repeat intv:=intv+1;
    dx:=dxval(intv);
    xmin:=xval(intv);
    xmax:=xval(intv+1);
   j:=0;
    repeat
    x:=xmin+j*dx;
    j:=j+1;
    expdecayv:=simproc(x,maxdiagram,N0,Ewritten,pvalue);
    message simproc in preceeding line may be undeclared;
    if list then 
    write(out,"nl",1,<<dd.ddd>,x);
    for i:=1 step 1 until maxdiagram do
      write(out,<< -d.dddd'+zd>,pvalue(i));
  until x>=xmax-dx or x>=tmax;
until x>=tmax;
end simproc
else
for dti:=1 step 1 until points do
begin
  for i:=0 step 1 until maxstates-1 do dNi(i):=0;
  tottrans:=0;
  for state:=0 step 1 until maxstates-1 do
  begin
    tau:=gettau(states,state);
    s0:=if delt=0 then 1 else sinh(delt/tau)/delt*tau;
    for i:=gettransno(states,state) step -1 until 1 do
    begin
      dipi:=dipf:=diprecsize*tottrans;
      lambda:=dip.dipf(9);
      if lambda<-smallreal then
      begin
        dn:=Ni(state)*dip.dipf(7)/10*dt*s0;
        dNi(state):=dNi(state)-dn;
        if state<>dip.dipi(1) then alarm(<:state error :>,state,dip.dipi(1));
        state2:=dip.dipi(2);
        dNi(state2):=dNi(state2)+dn;
      end lambda<0;
      tottrans:=tottrans+1;
    end calculate dn;
  end state;
  for i:=0 step 1 until maxstates-1 do
    Ni(i):=Ni(i)+dNi(i);
  if list then write(out,"nl",if (dti-1) mod 3=0 then 1 else 0,
    "sp",if (dti-1) mod 3 <>0 then 3 else 0,
    << -d.dddddddd'+d>,dt*dti,Ni(curstate));
  if dti=dtt and test then
  begin
  write(out,"nl",3,<:snapshot of Ni at :>,dtt*dt);
  for i:=0 step 1 until maxstates-1 do
    write(out,"nl",if i mod 5=0 then 1 else 0,
      "sp",if i mod 5<>0 then 3 else 0,i,Ni(i));
  outendcur(10);
  end test;
  end time loop;
end simcalc;

procedure simmodel(bsname,curstate,parterm,N0,N00,list,model,beta);
value curstate,N00,list,model,beta;
integer curstate,model; real N00,beta; boolean list;
array bsname,N0;
integer array parterm;
if model>=0 then
begin
integer trans,tottrans;
real cut,tau;
algol list.off copy.statevar;

algol list.off copy.allstatel;
end for state;
getstate(states,curstate,n,J,seriesindex,nstar,Ecm,L,app);
ser:=seriesindex*seriessize;
l:=series.ser(1);
if model=0 then N00:=N00/N0(curstate);

if model>1 and model<6 then N00:=N00*nstar**3;
if model=3 or model=4 then N00:=N00*nstar**2;
if finestruct and model>1 then N00:=N00/(J+1);
if model=3 then N00:=N00/(l+1);
if model=5 and l mod 2=1 then N00:=N00/3;
if model=6 then N00:=N00*((beta*cau)**2*nstar+Z*Z/nstar)**3;
end stateloop;

if lookupentry(<:model6:>)=0 then
  begin
  write(out,"nl",3,<:test of model 6 :>,
     "nl",1,<:Z= :>,Z,
     "nl",1,<:beta=:>,beta,1/beta,
     "nl",1,<:N(0)= :>,N00);
  writechargerec(out);
  write(out,"nl",2);
  end testmodel6;

algol list.off copy.allstatel;
  case model+1 of begin
  begin
    comment model 0, renormalize to curstate;
    N0(stateindex+state):=N0(stateindex+state)*N00*
     (if stateindex+state<>curstate then cascadefactor else 1);
  end model 0;
  begin
   comment model 1, only curstate populated;
   if stateindex+state=curstate then N0(stateindex+state):=N00;
  end model1;
  begin
    comment model 2, (2*j+1)/n**3;
    N0(stateindex+state):=N00/(nstar**3)*
     (if finestruct then J+1 else 1);
  end model2;
  begin
    comment model 3, (2*J+1)*(l+1)/n**5;
    N0(state+stateindex):=N00/(nstar**5)*
    (if finestruct then J+1 else 1)*(l+1);
  end model3;
  begin
    comment model 4, (2*J+1)/n**5;
    N0(state+stateindex):=N00/(nstar**5)*
    (if finestruct then J+1 else 1);
  end model4;
  begin
    comment model 5, (2*J+1)/n**3*(l odd then 3 else 1);
    N0(state+stateindex):=N00/nstar**3*
    (if finestruct then J+1 else 1)*(if l mod 2=1 then 3 else 1);
  end model5;
  begin
  comment model 6, plasma model of D. G. Ellis,
    N:=(beta**2*cau**2*nstar+Z**2/nstar)**-3,
    lmax:=2*Z/cau/beta,
    beta:=13.9/300*sqrt(E in MeV/ Mass in amu);
    N0(state+stateindex):=N00/(beta*beta*cau*cau*nstar+
      Z*Z/nstar)**(-3)*
      (if finestruct then J+1 else 1);
  end model6;
begin
  comment model 7, (J+1)*(l+1)/n**2;
  N0(state+stateindex):=N00/(nstar**2)*
  (if finestruct then J+1 else 1)*(l+1);
end model7;
begin
  comment model 8, (J+1)*(l+1)/n**3;
  N0(state+stateindex):=N00/(nstar**3)*
  (if finestruct then J+1 else 1)*(l+1);
end model8;
begin
  comment model 9, (J+1)*(l+1)/n**4;
  N0(state+stateindex):=N00/(nstar**4)*
  (if finestruct then J+1 else 1)*(l+1);
end model9;
  end case;
  if list then
  begin
    write(out,"nl",1);
    writestate(out,L,n,l,J);
    write(out,"sp",4,<< -d.dddddd>,nstar,
     << -d.dddd'-dd>,N0(stateindex+state));
  end list;
end end stateloop;
outendcur(0);
end simmodel;

procedure beam_foil_decay_curve(bsname,parterm,curstate,tmax,dt);
value curstate,tmax,dt;
integer curstate; real tmax,dt;
array bsname;
integer array parterm;
if chargesegdes(17)>=0 then
  write(out,"nl",1,<:*dipole matrixelements not resorted in :>,
    string inc(bsname))
else
begin
integer i,trans,tottrans,sta;
real tau,lambda,N00,cut;
array bsdipname(1:2),N0(0:maxstates-1);
algol list.off copy.statevar;
algol list.off copy.statevar2;
integer array field dipi;
array field dipf;

initrydfile(bsname(1) shift 24, real <:dip:>,segdes.dipsize,true);
for i:=1,2 do bsdipname(i):=segdes.dipname(i);
keystat:=lookupentry(<:rysimkey:>)=0;

begin
integer array dipval(1:
     if chargestate(7)=0 or (proc and -,generate) then 1 else
     chargestate(7)*diprecsize//2),
  dipindex(0:maxstates);
if -,(proc and -,generate) then
begin
  getalldipval(bsdipname,dipval);
  calcdipindex(dipval,dipindex);
end;
N00:=e;

if keystat then
begin
  for i:=0 step 1 until maxstates-1 do write(out,dipindex(i));
  for i:=0 step 1 until chargestate(7)-1 do
    writediprec(dipval,i);
end keystat;

for i:=0 step 1 until maxstates-1 do N0(i):=0;
if readmodel then
begin
  model:=0;
  readstateaux(bsname,parterm,N0);
  N0(curstate):=N0(curstate)*tauf;
end readmodel;
simmodel(bsname,curstate,parterm,N0,N00,list,model,beta);

N00:=N0(curstate);

algol list.off copy.allstatel;
end for state;

  getstate(states,curstate,n,J,seriesindex,nstar,Ecm,L,app);
  ser:=seriesindex*seriessize;
  l:=series.ser(1);
  tau:=gettau(states,curstate);
  if list then
  begin
    write(out,"nl",1,<:model :>,model,
      <:N0 = :>,N00,
      <:tau= :>,tau,
      <:tf = :>,tauf);
  end list;
  outendcur(10);
  connectcuro(output);
  write(out,"nl",1);
  writeatsym(out,S,atno,Z);
  write(out,"nl",1);
  writestate(out,L,n,l,J);
  if -,generate then
  begin
    if proc then
    begin
      nsmax:=simproc(-2,0,N0,parterm,N0);
      brmin:=simproc(-3,0,N0,parterm,N0);
    end proc;
    write(out,"nl",1,<< -d.dddddd>,
         tau*tauf, points,tmax,"nl",1,
         model,delt,brmin,nsmax);
  write(out,if proc then maxdiagram else 0);
  end -,generate;
  if diagram and generate then
    writediagramproc(out,bsname,dipval,states,series,dipindex,N0,curstate,delt) else
  simcalc(out,dipval,states,dipindex,N0,dt,tmin,tmax,curstate,(tmax-tmin)/dt);
  closeout;
  end stateloop;
end dipval;
end beam foil decay curve;
▶EOF◀