ADAS Subroutine nlthes
SUBROUTINE NLTHES(Z0,ZEFF,N,L,E0)
IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C PURPOSE: PROVIDES BINDING ENERGY OF TERM CENTRE FOR OUTER ELECTRON
C IN LITHIUM LIKE IONS
C
C FROM EDLEN (1979) PHYSICA SCRIPTA, 19, 255.
C
C FINE STRUCTURE FOR L>0 MUST BE ADDED EXTERNALLY. INFINITE MASS VALUE
C FOR RYDBERG CONSTANT IS USED. (NOT NOW SEE BELOW!!)
C ---
C 5 SEPT 1985 CHANGED R =109737.3 (I.E INFINITY MASS VALUE) TO
C ----------- THE Z DEPENDANT EQUATION SEE BELOW.......J.SPENCE
C
C ALL NEWLLPS3 RUNS FROM 5 SEPT 1985 ONWARDS HAVE THIS
C CHANGE ADDED IN, BUT %DIFF. IS SO VERY VERY SLIGHT THAT
C EVEN IF THE FINAL RESULTS DO CHANGE THEY WILL BE IN
C THE LAST FEW DEC. PLACES.
C E.G BEFORE 3.11321270 NOW GET 3.11321269 ===> APPROX. N/C
C
C NEWLLPS3.FORT(NLTHES2) ====> IS THE ORIGINAL ROUTINE WITH
C RZ=109737.318.
C (THIS ROUTINE EXISTS IN SOURCE FORM
C ONLY. I.E IT IS NOT COMPILED.)
C
C-----------------------------------------------------------------------
C VERSION : 1.1
C DATE : 18-03-1999
C MODIFIED : ???
C
C VERSION : 1.2
C DATE : 05-10-2000
C MODIFIED : ???
C - Removed junk from columns > 72
C
C VERSION : 1.3
C DATE : 16-05-2007
C MODIFIED : Allan Whiteford
C - Updated comments as part of subroutine documentation
C procedure.
C
C-----------------------------------------------------------------------
EL0(Z)=0.25D0*(Z*Z-3.18244D0*Z+2.0038D0+0.208015D0/(Z-1.3833D0))
&-6.3789D-6*RCORR*(Z-2.0D0)**2
SIG0(Z)=0.141441D0*(Z-1.7025D0-0.768371D0*(1.0D0-0.090333D0/(Z-
&2.184D0))/(Z-0.975D0))-1.9137D-6*RCORR*(Z-2.0D0)**2
P(XN,XL)=(3.0D0*XN*XN-XL*(XL+1.0D0))/(2.0D0*XN**5*(XL-0.5D0)*
&XL*(XL+0.5D0)*(XL+1.0D0)*(XL+1.5D0))
QP(XN,XL)=(35.0D0*XN**4-5.0D0*XN*XN*(6.0D0*XL*(XL+1.0D0)-5.0D0)+
&3.0D0*(XL-1.0D0)*XL*(XL+1.0D0)*(XL+2.0D0))/(8.0D0*XN**7*(XL-1.5D0)
&*(XL-1.0D0)*(XL-0.5D0)*XL*(XL+0.5D0)*(XL+1.0D0)*(XL+1.5D0)*(XL+
&2.0D0)*(XL+2.5D0))
DRS(Z)=0.456536D0*(Z-1.2808D0)**4+1.2763D-5*(Z-1.2808D0)**6+
&4.34D-10*(Z-1.2808D0)**8
DRP(Z)=0.21305D0*(Z-2.241D0)**4+0.466D-5*(Z-2.241D0)**6+
&1.48D-10*(Z-2.241D0)**8
DLS(Z)=4.5246D-3*(Z-1.6D0)**4*(-2.179D0-2.0D0*DLOG(7.29729D-3*
&(Z-1.6D0))+5.26427D-2*(Z-1.6D0)-5.32504D-5*(Z-1.6D0)**2*(3.0D0*
&(DLOG(7.29729D-3*(Z-1.6D0)))**2+8.695D0*DLOG(7.29729D-3*(Z-1.6D0))
&+19.081D0))
DLP(Z)=4.525D-3*(Z-2.0D0)**4*(3.0D-2-2.6412D-5*(Z-2.0D0)**2*DLOG(
&7.29729D-3*(Z-2.0D0)))
C
C ................. IMPROVEMENTS ....................
C
C R=109737.318D0 <<<<<<< OLD LINE BEFORE 5 SEPT 1985.
C REPLACEMENT FOR ABOVE LINE IS GIVEN BELOW.
R=109737.318D0-60.200D0/(2.00D0*Z0)
C
C EL0( ) AND SIG0( ) ABOVE HAD SOME CONSTANTS DIVIDED BY 109737.3
C THIS CORRECTION VALUE PUTS 109737.3 BACK IN AND INSTEAD DIVIDES
C THOSE CONSTANTS BY THE Z-DEPT R ABOVE.
RCORR=109737.3/R
C
C WRITE(6,1414)R,RCORR
1414 FORMAT(' R,RCORR = ',1P2E15.7)
EN=N
XL=L
T1=DRS(Z0)
T2=DRP(Z0)
T3=DLS(Z0)
T4=DLP(Z0)
C WRITE(6,101)T1,T2,T3,T4
IF(L-1)10,20,40
C S STATES
10 T2S=R*EL0(Z0)+DRS(Z0)-DLS(Z0)
C WRITE(6,101)T2S
T2S=T2S/(R*(Z0-2.0D0)**2)
C WRITE(6,101)T2S
D2S=2.0D0-1.0D0/DSQRT(T2S)
C=0.0828D0/Z0-0.2283D0/Z0**2
B=-5.5D-4*Z0+5.963D-3+0.19404D0/(Z0-0.36D0)-0.3368D0/(Z0-0.36D0)*
&*2
A=D2S-B*T2S-C*T2S*T2S
C WRITE(6,101)A,B,C
GO TO 30
C P STATES
20 T2P1=R*EL0(Z0)+DRS(Z0)-DLS(Z0)
T2P2=R*SIG0(Z0)+DRS(Z0)-DRP(Z0)-DLS(Z0)+DLP(Z0)
C WRITE(6,101)T2P1,T2P2
T2P=(T2P1-T2P2)/(R*(Z0-2.0D0)**2)
D2P=2.0D0-1.0D0/DSQRT(T2P)
C=-2.603D-2/(Z0-2.37D0)+1.326D-2/(Z0-2.37D0)**2
B=-1.2D-3*Z0+2.1237D-2-8.905D-2/(Z0-1.74D0)+4.803D-2/(Z0-1.74D0)*
&*2
A=D2P-B*T2P-C*T2P*T2P
C WRITE(6,101)A,B,C
30 U=0.0D0
31 U0=U
V=EN-U0
TNL=1.0D0/V**2
U=A+TNL*(B+TNL*C)
C WRITE(6,101)TNL,U,U0
101 FORMAT(1P4D15.7)
IF(DABS(U-U0).LE.1.0D-6)GO TO 35
GO TO 31
35 TNL=((Z0-2.0D0)/(EN-U))**2
GO TO 50
C L>2 CASES
40 S=0.3397D0+0.102D0/(Z0-0.4D0)
A=9.0D0*((Z0-2.0D0)/(Z0-S))**4
AK=0.2113D0*Z0+0.598D0-2.4D0/Z0
DELTAP=A*P(EN,XL)+A*DSQRT(A)*AK*QP(EN,XL)
C TRR=(Z0-2.0D0)**2*(1.0D0+5.32504D-5*(Z0-2.0D0)**2*(EN/(XL+0.5D0)-
C &0.75D0)/(EN*EN))/(EN*EN)
C WRITE(6,7845)TRR*R,DELTAP*R
7845 FORMAT(' TRR*R,DELTAP*R = ',1P2E15.7)
TNL=(Z0-2.0D0)**2*(1.0D0+5.32504D-5*(Z0-2.0D0)**2*(EN/(XL+0.5D0)-
&0.75D0)/(EN*EN))/(EN*EN)+DELTAP
50 E0=TNL
C WRITE(7,100)Z0,N,L,TNL
100 FORMAT(1H ,F5.1,2I5,1PD15.7)
C WRITE(6,7842)V
7842 FORMAT(' N STAR = ',1PE15.7)
C IF(L.LT.2)GOTO 7844
C WRITE(6,7843)EN,XL,R*P(EN,XL),QP(EN,XL)/P(EN,XL)
7843 FORMAT(' EN,XL,P,Q = ',1P4E15.7)
7844 RETURN
END
INTEGER L, N
REAL*8 E0, Z0, ZEFF