|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 18432 (0x4800)
Types: TextFile
Names: »rydiffint«
└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
└─⟦b2ec5d50f⟧
└─⟦this⟧ »rydiffint«
<* rydiffint solves the differential equation *>
<*81-01-13*>
real procedure fb2(fb1,l);
value fb1,l; real fb1; integer l;
fb2:=fb1*(2**(l+1));
real procedure ufb2(ufb1,energi,q,B,DE);
value ufb1,energi,q; real ufb1,energi;
integer q;
array B,DE;
ufb2:=ufb1*(if (1/B(q+1)+energi)>0 then 1 else
exp(DE(q)*sqrt(-2*(1/B(q+1)+energi))));
real procedure rval(i);
value i; integer i;
begin integer sq;
NE(-1):=0;
sq:=-1;
for sq:=sq+1 while i>NE(sq) do;
rval:=B(sq)+DE(sq)*(i-NE(sq-1));
end;
real procedure exactwf(n,l,B,DE,NE,ryd);
value n,l; integer n,l;
integer array NE;
array B,DE,ryd;
begin
real rho,drho;
integer i,imin,sq;
for sq:=q-1 step -1 until 0 do begin
imin:=NE(sq-1)+1;
if imin<1 then imin:=1;
drho:=DE(sq);
rho:=B(sq+1);
for i:=NE(sq) step -1 until imin do begin
comment ryd(i):=hydrogenwf(n,l,rho);
rho:=rho-drho;
end;
end;
exactwf:=-.5/n/n;
NE(-2):=n-l-1;
NE(-3):=8;
end;
real procedure renorm(f,DE,B,NE,q,r1m);
value q; integer q;
real r1m;
array f,DE,B;
integer array NE;
begin
integer sq,ngr,k;
real del,overlap,ovl,y,r,dr,y2,x,dr1m;
NE(-1):=0;
r1m:=r:=overlap:=0;
for sq:=q-1 step -1 until 0 do
begin
del:=DE(sq);
x:=B(sq+1);
ngr:=NE(sq-1)+1;
if ngr<1 then ngr:=1;
dr1m:=dr:=ovl:=0;
for k:=NE(sq) step -1 until ngr do
begin
y:=f(k);
y2:=y*y;
ovl:=ovl+y2;
dr:=dr+x*y2;
dr1m:=dr1m+y2/x;
x:=x-del;
end;
r1m:=r1m+dr1m*del;
overlap:=overlap+ovl*del;
r:=r+dr*del;
end;
r1m:=r1m/overlap;
renorm:=r/overlap;
overlap:=1/sqrt(overlap);
for k:=NE(q-1) step -1 until ngr do f(k):=f(k)*overlap;
end renorm;
\f
procedure potential(energy,nuccharge,q,n,l,lold,NE,DE,hj,pot);
value energy,nuccharge,q,n,l,lold;
real energy,nuccharge; integer q,n,l,lold;
array DE,hj,pot;
integer array NE;
begin
own integer lq;
if l<>lold or q<>lq then begin
integer i,m,k,l2,nmin;
real x,xf,nucfak,xmin,xmax;
henta(ryp,6*psegm,pot,psegm);
l2:=l*(l+1);
NE(-1):=0;
nucfak:=2*nuccharge;
for m:=q-1 step -1 until 0 do
begin
nmin:=NE(m-1)+1;
del:=DE(m);
xf:=B(m+1);
for k:=NE(m) step -1 until nmin do
begin
hj(k):=l2/xf/xf+ (if k>imax then -nucfak/xf
else nucfak*pot(k));
xf:=xf-del;
end;
end effective potential calculation;
gema(peff,0,hj,segm-1);
lq:=q;
end else henta(peff,0,hj,segm-1);
end potential;
\f
procedure lavpot(r,n,l);
value r,n,l; real r; integer n,l;
begin
if extype>=10 and extype<=110 then <* HFex(n,l)*>;
begin
real fak;
integer pst;
array pot(1:imax);
henta(ryp,2*psegm,pot,psegm);
if extype<>0 then begin
array pc,pe(1:imax);
if extype<5 or extype>=1000 then begin
henta(ryp,3*psegm,pc,psegm);
henta(ryp,4*psegm,pe,psegm);
if extype=1 then begin
for pst:=imax step -1 until 1 do
pot(pst):=2*pot(pst)+pc(pst);
end else begin
if extype =2 then fak:=1.5 else
if extype =3 then fak:=1 else
if extype>=1000 then fak:=(extype-1000)/1000;
if S= 2 then fak:=-fak;
for pst:=imax step -1 until 1 do
pot(pst):=2*pot(pst)+pc(pst)+fak*pe(pst);
end;
end Slater-type potentials;
if extype>=10 and extype<=110 then begin
if S=0 then fak:=1 else
if S=2 then fak:=-1;
henta(ryp,4*psegm,pc,psegm);
comment henta(ryexp,(lexi(n,l)-1)*psegm,pe,psegm);
if l=0 then begin
array off(1:imax);
henta(ryp,5*psegm,off,psegm);
for pst:=imax step -1 until 1 do
pot(pst):=2*pot(pst)+pc(pst)+fak*pe(pst)+off(1);
end l=0 loop else begin
for pst:=imax step -1 until 1 do begin
pot(pst):=2*pot(pst)+pc(pst)+fak*pe(pst)
end for loop;
end l<>0 loop;
end extype<>0;
end Hartree-Fock potentials;
gema(ryp,6*psegm,pot,psegm);
end main;
end lavpot;
\f
real procedure numug(q,n,l,lold,engu,zp,B,DE,NE,f,p,maxitr,Z);
value q,n,l,lold,engu,zp,maxitr,Z; integer q,n,l,lold,zp,maxitr,Z;
real engu;
array B,DE,f,p;
integer array NE;
begin
comment version 6-8-73;
message numug;
boolean b1,b2;
real energy,delen,f1,f2,f3,y1,y2,y3,del,del2,scai,scau,
int,iint,uint,norm,yfit;
integer i,itr,sq,elem,p0,nr,errorp,ilimit,
up_lim;
prno:=2;
scai:=scau:=1;
delen:=0;
nr:=nunr; elem:=NE(q-1); NE(-1):=itr:=1;
potential(engu,1,q,n,l,lold,NE,DE,p,f);
cleararray(f);
f(1):=DE(0); f(elem+nr):=exp(-B(q));
if f(elem+nr)<'-100 then f(elem+nr):='-100;
for energy:=-2*engu,energy+delen while
(itr<2 or abs(delen)>nuerror*energy) and itr<=maxitr do
begin
comment udefra integration;
setrydwhere(<:numug:>,ncur,lcur,0,0);
operator(false);
sq:=q;
b1:=b2:=true;
uint:=0;
y1:=f(elem+nr):=abs(scau*f(elem+nr));
y2:=f(elem+nr-1):=ufb2(y1,-.5*energy,sq-1,B,DE);
up_lim:=(NE(q-1)+3*NE(q-2))//4;
repudefra:
i:=NE(sq-1);
del:=DE(sq-1);
del2:=del*del;
f1:=del2*(energy+p(i));
f2:=del2*(energy+p(i-1));
int:=0;
ilimit:=NE(sq-2)-1;
if sq=1 then ilimit:=.5/DE(0);
for i:=i-2,i-1 while (b1 or b2) and i>=ilimit do
begin
if b1 then begin
if y2<y1 and i<up_lim then
begin
b1:=false;
yfit:=y2/2;
end;
end else
if y2<yfit then b2:=false;
f3:=del2*(energy+p(if i=ilimit then i-1 else i));
y3:=((24+10*f2)*y2+(f1-12)*y1)/(12-f3);
f1:=f2; f2:=f3; y1:=y2; y2:=f(i+nr):=y3;
errorp:=i;
int:=int+y3*y3;
end;
i:=ilimit+1;
uint:=uint+int*del;
if errorp>NE(0)-2 and errorp=ilimit then begin
comment next interval;
f(i+nr-1):=y2:=if nuipol then(f(i+nr-1)*105
+f(i+nr)*420
-f(i+nr+1)*210
+f(i+nr+2)*84
-f(i+nr+3)*15)/384
else (y1+y2)/2;
uint:=uint-(y3*y3-(y2*y2-y1*y1)/2)*del;
sq:=sq-1;
goto repudefra;
end;
sq:=p0:=0;
i:=1;
y1:=f(1):=abs(scai*f(1));
y2:=f(2):=fb2(y1,l);
iint:=(y1**2+y2**2)*DE(0);
repindefra:
ilimit:=NE(sq);
del:=DE(sq); del2:=del*del;
f1:=del2*(energy+p(if i=1 then 1 else i-1));
f2:=del2*(energy+p(if i=1 then 2 else i+1));
int:=0; up_lim:=(if ilimit<errorp+nr then ilimit
else errorp+nr-1);
for i:=i+2 step 1 until up_lim do
begin
f3:=del2*(energy+p(i));
y3:=((24+10*f2)*y2+(f1-12)*y1)/(12-f3);
f1:=f2; f2:=f3; y1:=y2; y2:=f(i):=y3;
if i<errorp then int:=int+y3*y3;
if y1>0 == y2<0 then begin
if i>.5/DE(0) then p0:=p0+1; end;
end;
i:=ilimit-1;
iint:=iint+int*del;
if errorp+nr-1>ilimit then
begin
comment next interval;
y1:=f(i-1);
sq:=sq+1;
goto repindefra;
end;
scau:=1/f(errorp+3*nr//2);
uint:=uint*scau**2;
scai:=1/f(errorp+nr//2);
iint:=iint*scai**2;
int:=iint+uint;
itr:=itr+1;
delen:=-(f(errorp+nr-1)*scai-f(errorp+2*nr-1)*scau
-f(errorp)*scai+f(errorp+nr)*scau)/int/del/(nr-1);
if abs delen>.1*energy then delen:=sign(delen)*energy*.1;
if p0<>zp then delen:=sign(p0-zp)*abs(delen);
if nukey or (maxitr<>0 and itr>=maxitr) then
begin
write(out,"nl",1,<:errorp:>,errorp,<: p0:>,p0,
<: delen:>,delen,
<: int:>,int,<: uint:>,uint,<: iint:>,iint,
"nl",1,<:sq:>,sq,<: scai:>,scai,<: scau:>,scau);
outendcur(10);
outendcur(10);
if -,nukey then goto ENDN; end;
end;
ENDN:
norm:=sign(scai)/sqrt(int);
scai:=norm*scai;
scau:=norm*scau;
for i:=1 step 1 until errorp-1 do f(i):=f(i)*scai;
for i:=errorp step 1 until elem do f(i):=f(i+nr)*scau;
NE(-3):=errorp;
NE(-2):=p0;
NE(-1):=itr-1;
numug:=-.5*energy;
end;
\f
real procedure nucut(q,n,l,lold,j,energy,B,DE,NE,f,p,nstar,zeroes,Z);
value q,n,l,lold,j,energy,nstar,Z; integer q,n,l,lold,j,Z;
real energy,nstar;
array B,DE,f,p;
integer array NE,zeroes;
begin
message nucut 81-01-29;
real f1,f2,f3,y0,y1,y2,y3,del,del2,r1m,r1me,
int,uint,norm,yfit,A,rns,delint,rint,x,overl,
error,newerror,ll2,dr,dnorm,dr1m,delta_r,d_ilim;
integer i,sq,elem,elem0,ilimit,imin,imax,
up_lim,zeroc,zp,cutcase;
boolean min,max,zero,nml2,nsint,nsgl1,rlist,absmax;
real procedure delr(x,del,i,dnorm,dr1m,list);
value x,del,i,list; integer i;
boolean list;
real x,del,dnorm,dr1m;
if rlfit and i>1 then
begin
real A,Ad;
A:=f(i)**2;
Ad:=A*del;
dnorm:=A*x;
dr:=delr:=(dnorm/(2*l+4)-Ad)*x;
dr1m:=(dnorm/(2*l+2)-Ad)/x;
dnorm:=dnorm/(2*l+3)-Ad;
if list and nukey then write(out,"nl",1,<:i,x,dr,dn,d1/r ,Ad :>,
i,x,dr,dnorm,dr1m,Ad);
end delr else delr:=dnorm:=dr1m:=0;
procedure removesol;
if i<up_lim and l>0 then
begin
ilimit:=i;
cutcase:=3;
if nukey then write(out,"nl",1,<:div y2,yfit :>,
<< -d.ddd'+ddd>,abs y2,yfit,<< dd.dd>,rval(i));
end;
prno:=3;
rlist:=lookupentry(<:nurlist:>)=0;
d_ilim:=0;
if lookupentry(<:nudilim:>)=0 then
begin
integer array t(1:10);
lookuptail(<:nudilim:>,t);
d_ilim:=t(6)-t(7);
end;
absmax:=false;
min:=zero:=false;
max:=true;
nsint:=abs(round nstar-nstar)<'-3;
ll2:=l*(l+1)/2;
nsgl1:=nstar>l+1;
if nsgl1 then d_ilim:=0;
elem0:=NE(0)-2;
zp:=NE(-2);
if nsint then zp:=nstar-l-1;
cutcase:=zeroc:=0;
if nukey then
begin
write(out,"nl",2,"*",67,"nl",1,"*",26,"sp",4,<:nucut:>,"sp",5,"*",27,
"nl",1,<:charge :>,Z,"nl",1,<:l,lold :>,l,lold,"nl",1);
writestate(out,-1,n,2*l,-1);
write(out,"nl",1);
end;
elem:=NE(q-1);
potential(energy,1,q,n,l,lold,NE,DE,p,f);
energy:=-2*energy;
NE(-1):=2;
r1me:=1/nstar/nstar;
rns:=1/2*(3*nstar**2-l*(l+1));
nml2:=-,randr1mcrit and ((l<=3 or
(l>3 and abs(round nstar - nstar)<.05)) and l>0);
imin:=if l<2 then 1.0 else
(if nsgl1 then abs(nstar-round nstar) else
if nstar-l>.75 then 4*DE(0) else
.5)/DE(0);
if l=2 then imin:=imin//2;
if l>2 and imin<.5/DE(0) then imin:=.5/DE(0);
if imin*DE(0)<l-4 then imin:=(l-4)/DE(0);
if l>3 then imin:=(l+0)/DE(0);
if -,autcut then imin:=1;
rint:=overl:=0;
setrydwhere(<:nucut:>,ncur,lcur,0,0);
operator(false);
yfit:=maxreal/100;
y1:=if nstar<l+1 then '-300 else exp(-B(q));
if y1<'-300 then y1:='-300;
f(elem):=y1;
y2:=f(elem-1):=ufb2(y1,-.5*energy,q-1,B,DE);
uint:=(y2*y2+y1*y1)*DE(q-1);
sq:=q+1;
if nukey then write(out,"nl",1,<:min ilimit :>,imin);
for sq:=sq-1 while sq>=1 and cutcase=0 do begin
i:=NE(sq-1);
up_lim:=NE(0)-2;
del:=DE(sq-1);
x:=B(sq)-del;
del2:=del*del;
f1:=del2*(energy+p(i));
f2:=del2*(energy+p(i-1));
int:=0;
ilimit:=NE(sq-2)-1;
if sq=1 then begin
ilimit:=if -,autcut then cut/DE(0) else imin;
if nukey then write(out,"nl",1,<:intlimit :>,ilimit,
<< dd.ddd>,rval(ilimit));
end;
if ilimit<1 then ilimit:=1;
for i:=i-2,i-1 while i>=ilimit and cutcase=0 do
begin
f3:=del2*(energy+
p(if i<>ilimit then i else if i>1 then i-1 else 1));
y3:=((24+10*f2)*y2+(f1-12)*y1)/(12-f3);
f1:=f2; f2:=f3; y0:=y1; y1:=y2; y2:=f(i):=y3;
delint:=y3*y3;
int:=int+delint;
if autcut then begin
if absmax and abs y2>yfit then removesol else
if max then begin
if y1>=y0 and y1>=y2 then
begin comment a true max found;
zero:=true; max:=false;
if -,absmax then yfit:=abs y2;
absmax:=true;
if nukey then write(out,"nl",1,<:max found :>,
"sp",6,<< ddd.d0>,rval(i),<< -d.ddd'+ddd>,yfit);
end;
end search for max else
if zero then begin
if (y2>0 and y0<0) or (y2<0 and y0>0) then
begin comment a zero found;
zero:=false;
zeroc:=zeroc+1;
zeroes(zeroc):=i;
if nukey then
write(out,"nl",1,<:zero found :>,"sp",6,<< ddd.d0>,rval(i),
<< d>,zeroc,<:.zero:>);
if y2>0 then max:=true else min:=true;
if zeroc>zp or (l>3 and zeroc>entier(nstar+.95)-l-1) then begin
ilimit:=i;
zeroc:=zeroc-1;
cutcase:=4;
end;
end zero found;
if i<elem0 and (nml2 or nsgl1) then begin
if y1<=y2 and y1<=y0 then begin
comment a false min found;
ilimit:=i;
cutcase:=1;
if nukey then
write(out,"nl",1,<:false min found :>,<< ddd.d0>,rval(i));
end false min found;
if y1>=y2 and y1>=y0 then begin
comment a false max found;
ilimit:=i;
cutcase:=2;
if nukey then
write(out,"nl",1,<:false max found :>,<< ddd.d0>,rval(i));
end false max found;
end test false min and max in innermost;
end search for zero else
if min then begin
if y2>=y1 and y1<=y0 then
begin comment a true min found;
zero:=true; min:=false;
if nukey then write(out,"nl",1,<:min found :>,
"sp",6,<< ddd.d0>,rval(i),<< -d.ddd'+ddd>,yfit);
end;
end search for min;
end search for a cut criteria;
end main integration loop;
i:=ilimit+1;
uint:=uint+(int-y3*y3)*del;
if sq>1 and cutcase=0 then begin
comment next interval;
f(i-1):=y2:=if nuipol then(f(i-1)*105
+f(i)*420
-f(i+1)*210
+f(i+2)*84
-f(i+3)*15)/384
else (y1+y2)/2;
uint:=uint+(y2*y2+y1*y1)*del;
end;
end for sq;
dnorm:=dr1m:=0;
for i:=1 step 1 until ilimit-1 do f(i):=0;
if cutcase=0 and autcut and l>0 and
(nstar<l+1 or abs y2>yfit) then removesol;
rint:=renorm(f,DE,B,NE,q,r1m);
x:=rval(ilimit);
dr:=delr(x,del,ilimit,dnorm,dr1m,true);
if nukey and l=0 then
begin
write(out,"nl",1,<:s-state after integration:>,
"nl",1,<:<r>,<r>e :>,rint,rns,
"nl",1,<:<1/r>,<1/r>e :>,r1me,r1m);
end;
if nukey then
begin
write(out,"nl",1,"*",2,<:after integration :>);
writestate(out,-1,n,2*l,j);
write(out,"nl",1,<:cut :>,rval(ilimit),
"nl",1,<:cutcase :>,cutcase,
"nl",1,<:zeroes :>,zeroc,
"nl",1,<:<r>e <r> :>,rns,(rint+dr)/(1+dnorm),<: <dr> <dnorm> :>,dr,dnorm);
end;
if rns<rhomin then rns:=rhomin;
if autcut and l>0 and -,nsint then
begin
if nukey then write(out,"nl",1,
<:************************** <r> adjust ***************************:>);
overl:=1;
x:=rval(ilimit);
dr:=delr(x,del,ilimit,dnorm,dr1m,true);
deltar:=f(ilimit)**2*del;
if nukey then write(out,"nl",1,<:<r>e,<r> :>,rns,(rint+dr)/(1+dnorm),
<: delta <r> :>,deltar*x/(1-deltar),
<:<10> cut :>,rval(ilimit));
i:=0;
cutcase:=6;
y1:=y2:=0;
sq:=0;
up_lim:=NE(0);
x:=del:=DE(0);
imax:=.9*elem;
error:=maxreal/2;
newerror:=maxreal/4;
<*adjustment will only take place when <r>e><r> *>
for i:=i+1 while newerror>2'-3 and abs newerror<= abs error and
-,((y2>0 and y1<0) or (y2<0 and y1>0)) and i<=imax do begin
error:=newerror;
if i>up_lim then begin
sq:=sq+1;
del:=DE(sq);
x:=B(sq+1);
up_lim:=NE(sq);
end;
y1:=y2;
y2:=f(i);
int:=y2*y2*del;
overl:=overl-int;
rint:=rint-x*int;
if i>=ilimit then
begin
dr:=delr(x,del,i,dnorm,dr1m,false);
newerror:=(1-(rint+dr)/(overl+dnorm)/rns);
if rlist then
begin
write(out,"nl",1,"*",1,i,
rval(i),(rint+dr)/(overl+dnorm),rint,dr,overl,
dnorm,<<-d.ddd'd>,newerror);
end;
if newerror>2'-3 and abs newerror<abs error then
begin
if i>1 then f(i-1):=0;
ilimit:=i;
end;
end i>=ilimit;
x:=x+del;
end;
if nstar<l+1 and zeroc=0 and cutcase<4 then cutcase:=5;
rint:=renorm(f,DE,B,NE,q,r1m);
dr:=delr(rval(ilimit),del,ilimit,dnorm,dr1m,true);
if nukey then write(out,"nl",1,<:after <r> adjustment <r>e,<r> :>,
rns,(rint+dr)/(1+dnorm),
<: cut :>,rval(ilimit),<: error :>,newerror);
end <r> adjust;
if nukey then write(out,"nl",1,
if r1mc then <:<1/r> criterion :> else <: :>,
<:<1/r> <1/r>e :>,(r1m+dr1m)/(1+dnorm),r1me,
<: cut :>,<<dd.dd>,rval(ilimit),
<: n* :>,<<d.ddd>,nstar,<: l :>,<<d>,l);
if autcut and l>0 and -,nsgl1 and r1mc then
begin
if nukey then write(out,"nl",1,
<:************************** <1/r> adjust ***************************:>);
cutcase:=5;
overl:=1; del:=DE(0); x:=del; i:=sq:=0;
up_lim:=NE(0);
for i:=i+1 while (r1m+dr1m)/(overl+dnorm)>=r1me and i<elem do begin
if i>up_lim then begin
sq:=sq+1; x:=B(sq+1); up_lim:=NE(sq);
del:=DE(sq);
end;
y2:=f(i);
int:=y2*y2*del;
overl:=overl-int;
r1m:=r1m-int/x;
if i>=ilimit then
begin
dr:=delr(x,del,i,dnorm,dr1m,false);
ilimit:=i;
end;
x:=x+del;
if i>1 then f(i-1):=0;
end adjust <1/r>;
rint:=renorm(f,DE,B,NE,q,r1m);
dr:=delr(rval(ilimit),del,ilimit,dnorm,dr1m,true);
if nukey then write(out,"nl",1,<:after adjustment <1/r>e <1/r> :>,r1me,
(r1m+dr1m)/(1+dnorm),
<: cut :>,rval(i));
end r1mcriterion;
if autcut and l>0 and randr1mcrit and -,nsgl1 then begin
if nukey then write(out,"nl",1,
<:************************** error adjust ***************************:>);
overl:=1; del:=DE(0); x:=del; i:=sq:=0;
up_lim:=NE(0);
newerror:=error:=abs((rint+dr)/(1+dnorm)-1.5/(r1m+dr1m)/(1+dnorm)+ll2);
cutcase:=6;
for i:=i+1 while newerror<=error and i<elem do
begin
if i>up_lim then begin
sq:=sq+1; x:=B(sq+1); up_lim:=NE(sq);
del:=DE(sq);
end;
y2:=f(i);
int:=y2*y2*del;
overl:=overl-int;
r1m:=r1m-int/x;
rint:=rint-int*x;
if i>=ilimit then
begin
dr:=delr(x,del,i,dnorm,dr1m,false);
ilimit:=i;
end;
x:=x+del;
if d_ilim=0 and i>1 then f(i-1):=0;
error:= newerror;
newerror:=abs((rint+dr)/(overl+dnorm)-1.5/(r1m+dr1m)*(overl+dnorm)+ll2);
end adjust;
if d_ilim<>0 then
begin
for i:=1 step 1 until ilimit+dilim-1 do f(i):=0;
write(out,"nl",1,<:cut (error) :>,rval(ilimit),<:d_ilimit :>,d_ilim);
ilimit:=ilimit+d_ilim
end;
rint:=renorm(f,DE,B,NE,q,r1m);
dr:=delr(rval(ilimit),del,ilimit,dnorm,dr1m,true);
if nukey then write(out,"nl",1,<:after error minimisation:>,"nl",1,
<:<1/r>e <1/r> :>,r1me,(r1m+dr1m)/(1+dnorm),"nl",1,
<:<r>e <r> :>,rns,(rint+dr)/(1+dnorm),"nl",1,
<:<37>error :>,newerror/(rint+dr)*(1+dnorm)*100,"nl",1,
<:cut :>,rval(ilimit));
end nml2 and weighted <r> <1/r> criterion;
if autcut then cut:=rval(ilimit);
if rlfit and l>0 then begin
A:=f(ilimit)/(rval(ilimit)**(l+1));
sq:=0;
del:=x:=DE(sq);
up_lim:=NE(sq);
for i:=1 step 1 until ilimit-1 do
begin
if i>up_lim then
begin
sq:=sq+1;
del:=DE(sq);
up_lim:=NE(sq);
x:=B(sq+1);
end change step;
f(i):=A*x**(l+1);
x:=x+del;
end;
rint:=renorm(f,DE,B,NE,q,r1m);
if nukey then write(out,"nl",1,<:after r**l+1 fit :>,"nl",1,
<:<1/r> <r> :>,r1m,rint);
if nukey then write(out,"nl",1,
<:<37>error :>,100*abs(rint-1.5/r1m+ll2));
end rlfit;
NE(-1):=0;
NE(-2):=zeroc;
NE(-3):=cutcase;
if nukey then write(out,"nl",1);
nucut:=-.5*energy;
end;
▶EOF◀