DataMuseum.dk

top - metrics - download

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

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦e6c2bcfa6⟧ »cryprog« 
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
    └─⟦4334b4c0b⟧ 
        └─⟦e6c2bcfa6⟧ »cryprog« 
            └─⟦this⟧ 

TextFile

<*                  rydiagpr   calculation of diagonal matrixelements      *>

procedure rydiag(bsname,parterm);
integer array parterm;
array bsname;
if chargesegdes(13)<0 and diag  then
begin
integer i,k,sq,ngr,em,ll1,lc,lastsegm,coresold;
real f0,f1,f2,f3,f4,df,del,del12,x,h,j,r2,eta,
      fak,ov,ovl,jint,iint,id,g,jd,Ep,Ek,Vc,Ve,grovl,ns2,
      cut;
algol list.off copy.statevar;
array pot(1:imax),pc,grf(1:if extype>0 then imax else 1),
    pe(1:if extype >1 then imax else 1),bsdianame(1:3);
integer array diagval(1:(chargestate(6)+1)*diarecsize//2);
array field dia;

prno:=6;
first:=true;
keystat:=survey:=false;
if diagkey then begin
  write(out,"nl",2,<:segment description record:>);
  writerec(out,segdes,0,0,segdesrecsize);
  end;
initrydfile(bsname(1) shift 24,real <:dia:>,segdes.diasize,false);
for i:=1,2 do bsdianame(i):=segdes.dianame(i);

if -,diagexact then begin
if extype=2 then fak:=1.5 else
if extype=3 or (extype>=10 and extype <=110) then fak:=1 else
if extype>=1000 then fak:=(extype-1000)/1000;
etam:=0;
if S=2 then fak:=-fak;
for k:=1,2,3,4 do
  case k of begin
    henta(ryp,6*psegm,pot,psegm);
    if extype>0 then henta(ryp,psegm,grf,psegm);
    if extype>0 then henta(ryp,3*psegm,pc,psegm);
    if extype>1 then henta(ryp,4*psegm,pe,psegm);
    end;
end;
algol list.off copy.stateloop;
if cores=0 then coresold:=0;
  ll1:=l*(l+1)//4;
ncur:=n; lcur:=l; jcur:=J;
setrydwhere(<:rydiag:>,n,l,0,0);
operator(false);
cpu:=systime(1,0,time);
if diagexact then nstar:=n else begin
if extype >=10 and extype <=110 then
  henta(ryexp,state*psegm,pe,psegm);
computed:=getcomputed(states,state,cutcase,cut);
end;
if computed then
  begin 
   array ryd(1:(if diagexact then 1 else (segm-1)*128));
  chargesegdes(14):=-1;
if -,diagexact then begin
  henta(ryf,(state+stateindex)*segm+1,ryd,segm-1);
    E:=(-Ip+Ecm)/Econv;
    r2:=h:=iint:=jint:=ovl:=Ep:=Ek:=Vc:=Ve:=grovl:=0;
    fak:=fak*(1.5/(arctan(1)**2))**(1/3);
    NE(-1):=0;
    icut:=cut/del0;
    f0:=f1:=f2:=ryd(NE(q-1));
    f3:=ryd(NE(q-1)-1);
    for sq:=q-1 step -1 until 0 do begin
    del:=DE(sq); del12:=1/12/del;
    x:=B(sq+1);
    ngr:=NE(sq-1)+1;
    if ngr<1 then ngr:=1;
    f0:=f1;
    f1:=(f2+f1)/2;
    for k:=NE(sq) step -1 until ngr do
    begin
      f4:=(if k>2 then ryd(k-2) else 0);
      g:=f2*del; ov:=f2*g; ovl:=ovl+ov;
      jd:=ov/x; jint:=jint+jd;
      if k>ngr and abs (k-icut)>3 then df:=(8*(f3-f1)+f0-f4)*del12 else
      if k>icut+2 and k=ngr and k>3 then df:=(8*(f3-f1)+f0-ryd(k-3))*del12
       else df:=0;
      id:=g*df; iint:=iint+id;
      h:=h+ov*x;
      r2:=r2+ov*x*x;
      Ek:=Ek+df*df*del+ll1*jd/x;
      if k<=imax then begin
       Ep:=Ep+pot(k)*ov;
    if extype>0 then begin
      Vc:=Vc+pc(k)*ov;
      grovl:=grovl+ov*grf(k);
      if extype>1 then
        begin 
         Ve:=Ve+ov*pe(k) end;
      end;
      end else Ep:=Ep-jd;
      f0:=f1; f1:=f2; f2:=f3; f3:=f4;
      x:=x-del;
    end;
    end;
end else begin
comment exact integrals;
ns2:=nstar**2;
Ep:=-1/ns2;
E:=-1/2/ns2;
Ek:=1/2/ns2;
h:=1/2*(3*ns2-ll1);
r2:=ns2/2*(5*ns2+1-3*ll1);
jint:=1/ns2;
iint:=2/ns2;
end exactdiag;
Ep:=Z*Ep;
E:=Z*Z*E;
Ek:=Z*Z*Ek;
Vc:=Z*Vc;
Ve:=Z*Ve;
jint:=Z*jint; iint:=Z*iint;
h:=h/Z; r2:=r2/Z/Z;
eta:=-Ep/Ek;
etam:=etam+eta;
if diag then begin
if first then begin
  page:=1;
  writepage(page,cores,parterm);
  head(r); 
  if diagexact then write(out,"nl",1,<:exact values:>) else
  writeparam;
  write(out,"nl",1,<:expectation values:>);
  end;
if first or cores<>coresold or
    (l<>lc and lin+series.ser(3)-series.ser(2)+2>maxlines) then
begin
  if -,first then writepage(page,cores,parterm);
  write(out,"nl",2,"sp",if finestruct then 9 else 5,
    <:<1/r>:>);
  write(out,<:<df/dr>:>);
  write(out,"sp",5,<:<r>:>,"sp",4,<:<r**2>:>,
    "sp",4,<:E:>,"sp",6,<:2*Ek:>,"sp",5,<:Ep:>,"sp",2);
  if overl then write(out,<: overlap  :>);
  if extype>0 then write(out,<:coulomb:>,
     "sp",3,if extype>1 then <:exchange  :> else <: :>,
     <:<fg!fnl>:>,"sp",4,if extype>1 then <:eta      :>
     else <:  :>);
  if cput then write(out,<:cpu-time:>);
  if nuitr=0 and -,exact then write(out,<:  cut:>);
  first:=false;
  lin:=lin+2;
  lc:=l;
  coresold:=cores;
end sideoverskrift;
write(out,"nl",if l<>lc then 2 else 1);
writestate(out,L,n,l,J);
   write(out,<<-d.ddd>,"sp",1,jint,"sp",1,iint-jint,
   string layr,h,
   string layr2,r2,<<-dd.ddd>,
   "sp",1,E,"sp",1,Ek,"sp",1,Ep);
if extype>0 then write(out,<< -d.dddddd>,Vc);
if extype>1 then write(out,<< -d.dddddd>,fak*Ve);
if overl then write(out,string laye,ovl);
if extype>0 then write(out,<< -d.dddddd>,grovl); 
if extype>1 then write(out,<< -d.dddddd>,eta);
if cput then write(out,<< -d.dd>,systime(1,time,time)-cpu);
if nuitr=0 and -,exact then write(out,<<dd.ddd>,cut);
lc:=l;
lin:=lin+1; 
end writediagonal values;
end computation ;
dia:=(state+stateindex)*diarecsize;
for i:=1 step 1 until 6 do
diagval.dia(i):=case i of (jint,iint-jint,h,r2,Ek,Ep);
end end end end stateloop;
ncur:=0;
comment if extype>=1000 then
  write(out,"nl",3,<:etamean = :>,<<-d.dddd>,etam/lexi(n,l));
lastsegm:=putstruct(bsdianame,diagval,
    0,chargestate(6),diarecsize,chargesegdes(9));
end rydiagonal;
▶EOF◀