National Bureau of Standards

Library, N.W. Bids

NBS I Il— i

1 K 1 J 1x1 Alt I I

410

Improved Fortran Program for

Single Particle Energy Levels

and Wave Functions in Nuclear

Structure Calculations

RANDALL S. CASWELL

\

f

Kt"1 <"

Q

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

UNITED STATES DEPARTMENT OF COMMERCE • John T. Connor, Secretary

NATIONAL BUREAU OF STANDARDS • A. V. Astin, Director

NBS TECHNICAL NOTE 410

ISSUED SEPTEMBER 30, 1966

Improved Fortran Program for

Single Particle Energy Levels

and Wave Functions in Nuclear

Structure Calculations

Randall S. Caswell

Institute for Basic Standards

National Bureau of Standards

Washington, D.C., 20234

NBS Technical Notes are designed to supplement the Bureau's regular publications program. They provide a means for making available scientific data that are of transient or limited interest. Technical Notes may be

listed or referred to in the open literature.

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C., 20402

Price 25 cents.
Contents

Page

1. Introduction 1

2. Description of the Calculation 2

3. Search Procedure 5

4. Sample Results 7

5. Accuracy and Checks of the Code 7

6. References 8

Appendix 1. Evaluation of f .(outside) for Scattering States -- 9

Appendix 2. Instructions for use of EIGENV subroutine and meaning of input variables 11

Appendix 3. Fortran II listing 13

Appendix 4. Fortran IV listing 21

ii
Improved Fortran Program for Single Particle Energy Levels and

Wave Functions in Nuclear Structure Calculations

Randall S. Caswell

An improved program has been developed for

numerical calculations of single particle energy

levels and wave functions for Woods- Saxon or other

real potential wells. The program is designed to

be used as a subroutine of a larger program for

nuclear structure calculation but may be used

separately if desired. Improvements over the

previous program include: (1) for bound states

the wave function is calculated by integrating

outward from zero radius and inward from a maximum

radius beyond the nuclear potential well with matching at an intermediate radius; (2) a first order

correction has been made to the starting conditions

at small radius for the integration of the wave function; (3) a large step size is used until the calculation has nearly converged on the eigenvalue; then a small step size is used to provide maximum accuracy;

(4) the correction of a small error in the KuttaRunge integration procedure has been made. As a re- sult of the changes the program is approximately five

times faster than the previous version. Checks with

harmonic oscillator potentials give an accuracy of

.01% in the energy level values. Listings of Fortran II

and Fortran IV versions are given.

KeyWords: Eigenvalues, eigenfunctions , energy levels , wave functions, single particle shell model.

1. Introduction

Nuclear shell model calculations usually have employed harmonic

oscillator wave functions for the evaluation of radial integrals for

both bound and continuum states. Deeply bound states can be well

represented by harmonic oscillator wave functions. For continuum

wave functions it appears desirable to use more realistic wave functions, for example, those of a real Woods- Saxon well. A previous

program (Caswell, 1962) was developed to provide these eigenvalues

and wave functions. The older program, however, used a method of

matching . inside and outside wave functions at successively larger radii

which is inherently slow. In some large nuclear structure programs
the calculation of the single particle wave functions occupies an undesirably large fraction of the computer time. The program described

here speeds up the calculation by a factor of roughly 5» provides

increased accuracy, and gives normalized wave functions at up to 201

radial points.

2. Description of the Calculation

The chief calculational steps are: (l) calculation of the WoodsSaxon real potential well including a spin-orbit term; (2) numerical

integration of the wave function, u(r) = r|(r), from zero to a matching radius, RMATCH and, for bound states, from a maximum radius beyond

the nuclear potential, RMAX, backwards in radius to the same matching

radius -- for scattering states the maximum radius is used as the

matching radius; and (3) calculation of the logarithmic derivatives,

f (inside) and f (outside). For scattering states f (outside) is

determined from the criterion of a 90° phase shift. It is also possible

to read in from the input data values for f .(outside) which are taken

JC as the value at the maximum radius, RMAX„ The approximate location of

the eigenvalues may be determined by a survey versus energy of f (inaj side), f (outside) , and the difference between them, called "signature

of the well". An automatic search procedure finds the energy for which

f (inside) = f (outside), and the wave function vs. radius is calculated

Xj AJ and stored in memory.

The potential is calculated by subroutine VR3 as a function of

radius at up to 399 points, which is the number required for the numerical integration routine for up to 201 radial points for the wave

functions. Spacing between points is RMAX/NINTVL where RMAX is the

maximum radius and NINTVL is the number of intervals (one less than

the number of points). The potential used is:

V(r) = V

c

p(r) - „V

c (£)* <*3)i g (1)

where p(r) = -— , the "Saxon" potential; a is the strength

1 + expf
of the spin-orbit interaction compared to the "Thomas" interaction

for a nucleon; M is the nucleon mass; V is the central real potential;

ft/Mc is the Compton wave length for a nucleon; Z is the orbital angular

—

>

momentum of the neutron; a is the Pauli spin operator of the incident

neutron; R is the nuclear radius; and a is the "diffuseness" parao

r meter of the potential „ If a negative value is given for a> it is

2 taken as V , the depth of the spin-orbit well in MeV where a,V (ttt" SO / a \2 C\cdrlC

is replaced r byJ V ( ) , m being the mass of the jt-meson. If

soVm c/ at

it j = Z + \t then a»i = I, and the values of the variable VPLUS (nuclear

potential) for different radii are calculated. If j = JL - \t then

a* Z = -(j6+1), and potential is calculated for the spin-orbit antiparallel case. Other potential shapes can be used by changes in this

subroutine. For example, the program has been tested with a harmonic

oscillator potential.

The calculation of f (inside) and f (outside), and the search

procedure for matching the f values is controlled by EIGENV, the main

subroutine of the program. First, the radial part of the non-relativistic Schrb'dinger wave equation is integrated out in radius for the

appropriate Z:

d

2 u.(r) 2m r

Z(Z+l)ft2 -1— +-

dr3 tf

E V(r)lu (r) = 0. (2)

2mr3 J l

The values of the function and the first derivative at the first

radial point (radius = H ) are determined from the relation u (r) =

JL+l

C r as r -> 0, where C is a constant. We make a first order correc- tion to the starting value of the wave function, u, and the derivadu , -. . , . -, tive j— at the first radial point a; rewritten (omitting subscript i,'s):

tive j— at the first radial point as follows: Equation (2) may be

£u = Ki+lju m 2h (e.v)u> (3)

dr3 r

3

fi2
As r -» 0, the last term on the right is negligible compared to the

first and our solution is just u = r as indicated above, taking

C equal to one. Taking into account the second term and using our

JO approximate solution for u, we may find the derivative at the first

point r (= HI) by integration of (3).

du , „ -v I 2m (E-V)r m

d7

= u+1)r -W (ilL) • w

Integrating again,

_ i+1 . 2m (E-V) 1+25 f5

x u_r W (l+2)U+3)

r

* t5;

Equations (1^) and (5) are used as the initial values of wave function

and derivative to be used by the wave function integration subroutine,

V is taken as the value of the potential at the point r = H . A similar problem arises in the backward integration from a radius EMAX for bound states. At this distance the centrifugal potential term is small and the differential equation may be written:

= f

1 [V(r)-E]u = 7

2 u (6)

where the expression in the square bracket is positive for a bound

state. If we consider V(r) as constant, then the descending ex- ponential solution to this equation is u = exp(-7r), du/dr =

-y exp(-7r), and (du/dr) /r = -7. Since normalization is initially

arbitrary, we set a = 1 and (du/dr) = -7 at RMAX for starting the

backward integration.

The second order Runge-Kutta method is used for the numerical

integration (see Zurmuhl, I96I) . An error in the previous code

(Caswell, I962) which led to errors of about 1 percent in the eigenvalues and eigenfunctions by causing a slightly wrong potential to be
ased at certain steps has been corrected. The integration is con- trolled by subroutine INTEG3 and carried out by subroutines KRR1

and FR.

Modes of the calculation . The calculation may be run in three alternative modes. When MODE = 1, if the energy is negative (bound state),

the wave function is integrated outward from zero to RMATCH and inward

from RMAX to RMATCH, using the starting conditions given above. For con- vergence on an eigenvalue it is required that the "inside" and "outside"

f 's match at RMATCH. If the energy is positive (continuum state), the

criterion of a 90 phase shift is applied to the wave function at the

radius RMAX where the inside and outside wave functions must join smoothly (as indicated by having the same value of f.)« For the 90 phase

Xj shift continuum case, the "outside" or asymptotic solutions are spherical Neumann functions, G. = -xn (x) , since Coulomb effects are not

Xj Xj considered. The argument x = kr. The evaluation of the logarithmic

derivatives, f , is given in Appendix 1.

Xj When MODE = 2, values of f (outside) at a radius beyond the nuclear

Xj

I potential radius, RMAX, may be given the EIGENV program by the calling

I program. The variable FLOUT is used to read in the value of f (outXj side). In the case of a bound state, after integration inward from

RMAX starting at the value FLOUT, the variable FLOUT is actually

i matched to the f (inside).

When MODE = 3, the program just calculates and normalizes the wave function at the given energy EN. No search is carried out. EN must

have been determined on a previous run. This mode is provided to minimize running time by avoiding unnecessary searches.

3. Search Procedure

To begin the search procedure, one needs a trial energy. In

principle this may be selected by guess. The guess will usually con- verge to some eigenvalue if enough search steps are allowed. A more systematic way to find approximate energies for starting the automatic

search procedure is to use the method called the "signature of the well."

This is carried out by setting the variable NOMAX = 1, and making a run with a large number of trial energies, for example, at 1 MeV steps

5
throughout the energy range of interest. When NOMAX = 1, the radius of

the well at the energy, EN, is used as the matching radius for bound

states, and RMAX is used for scattering states. Values are calculated

for f (inside), f (outside), and for the difference between the two

Xj Xj (see, for example, Figure 1). The integration in radius is carried out

only once for each value of EN. The energy values at which the difference, DIF, equals zero are close to the energy of the final converged

eigenvalue which is to be found, and they serve as suitable starting

values for the automatic search procedure.

Automatic search procedure . First an integration of the wave function is carried out at the trial eigenvalue energy (the value of EN

when EIGENV is called). A second integration is carried out at a nearby

energy (one percent difference in energy). Using the differences (between the inside and outside f * s found at the matching radius) found

Xj in these two trial calculations, a linear interpolation or extrapolation

is made to an energy for which the difference is predicted to be zero (corresponding to smooth joining of the inside and outside wave functions). Integration in radius is carried out at this predicted eigenvalue energy, and the difference is again found. We now know three

energies and the three corresponding differences. Using this information, the program predicts by quadratic interpolation a fourth energy

at which the difference should be closer to zero. The calculation

continues, energy prediction being made by quadratic interpolation on

the last three values, until either the difference between f.'s is

Xj closer to zero than the difference specified by the value of EPS given

to the subroutine, or until the maximum number of tries, NOMAX, has

been used up. The wave functions during the convergence procedure are shown in Figure 2 for a typical case. If the maximum number of tries

has been used up without converging on an eigenvalue, the value of EN

and values of the wave function at radial points, WAVEFN(I), are set

equal to zero to insure that the calling program does not use incorrect

information. When the converged eigenvalue has been found, the part of

the wave function between RMATCH and RMAX is renormalized so that it has

the proper absolute magnitude to join smoothly to the inner wave

6
function at RMATCH. A radial integration of the Ju

3 dr is then carried

out by subroutine SIMPSN, and the complete wave function is normalized

so that

subroutine ill

J

u

2 dr = 1. o To save computer time an integration step size five times the final

step size is used for the initial steps of the search (when DIFX5.1).

Integration step size is controlled by the variable KSTEP which is set

jequal to 5 in the early, rough part of the search and set equal to 1 for

the final determination of the eigenvalue and wave function.

4. Sample Results

Sample results for eigenvalues and wave functions are shown in

Figure 3. Eigenvalues are shown for the lsi, 2si , 3si, lp^» lPi» Ids

state. Wave functions are shown for the 15i, 25i, and Ids. states.

5. Accuracy and Checks of the Code

To test the accuracy of the program, subroutine VR3 was rewritten

to use an harmonic oscillator potential:

V(r) = -100 I 1 - (r/5)3 MeV

where r is in fermis. For this particular potential hco = 18.810 MeV

Ifor AMASS = 15.0. Based on this value we may compare eigenvalues given

by the code for various I values with the usual relation E = (n+j-Ohoo,

all energies being in MeV.

Harmonic Oscillator

Levels s states p states d states

-15.355 -15.352(2p)

-34.165 -34.162(2s) -34.l63(ld)

-52.975 -52.974(lp)

-71.785 -71.784(ls)

It can be seen that the code predicts the level separations to an

accuracy of about 0.01% which is about the accuracy of various constants

in the code.
6. References

Abramowitz, M. and Stegun, I. A., (1964). Handbook of Mathematical

Functions , National Bureau of Standards Applied Mathematics

Series No. 55.

Blatt, J. M. and Weisskopf, V. F. , (1952). Theoretical Nuclear Physics , 334, Wiley (New York).

Caswell, R. S. (1962). A Fortran code for calculation of eigenvalues

and eigenfunctions in real potential wells, NBS Technical Note 159.

The present report replaces this reference, which should no longer

be used.

Schiff, L. I. (1955). Quantum Mechanics , 77, McGraw-Hill (New York).

ZurmUhl, R. (1961). Praktische Mathematic fur Ingenieure and Physiker , 412, Springer-Verlag (Berlin, Gottingen, Heidelberg).

The author would like to thank Michael Danos for suggestions and

his continued interest in this problem.
Appendix 1. Evaluation of f (outside) for Scattering States

For continuum states (EX)) we may write for the logarithmic derivative:

(du /dr) [~(d/dr)(rn.)l

R f^(outside) = R £ ^ = R ±= Ui , (!)

(uA [rnA

where n, is the spherical Neumann function, which corresponds to a 90°

phase shift, as defined in Schiff (1955).

But — (xn^) = — (x*

+1x~^(x)) = (jH-l)[n

A

(x)] - **mW. (2)

/2m xi

If X = kR = (—^ E) R, we have for f , (outside) V J l

Xn 4^(X) , 6/ . f, (3) = jH-1 ^± (= x-^-1,

I Vx> G/

where the prime refers to differentiation with respect to X. The

quantities F and G are the "regular" and "irregular" wave functions

which for neutrons (no Coulomb interaction) are F = X j.(X) and

G

i

= - x vx) - For values of X>1, "forward" recursion formulae are used (see,

for example, Abramowitz and Stegun, 1964):

cosX

n

Q

(X) = , X

cosX sinX

n,(X) = 5 , and XXn (X) n

x+1

(X) = (2X+1)-— nM(X).
For I < 9 the values accurately reproduce tabular values from the same reference, there appearing to be no need for use of "backward" recursion

formulas (from large to small I).

For X < 1 a series expansion is used (see Abramowitz and Stegun,

1964)

1-3' 5'... (2j&-1) , h? (|x2

) a n/

X) = "1+1 \}~ 11(1-21)

+ 2l(l-2JL)(3-2Z) " " J

* X

10
Appendix 2. Instructions for use of EIGENV subroutine and meaning of

input variables.

To call the EIGENV subroutine values must be assigned to the following list of variables (with meanings given) and then the statement CALL

EIGENV calls the subroutine. The subroutine returns as the answer a value of EN which is the energy eigenvalue, and of a normalized WAVEFN(I)

at up to 201 equally spaced points. The wave function is rf.

TYPICAL VALUE MEANING

3.15 radius of the Saxon potential in fermis,

typically 1.25 A 1/3 .

.65 diffuseness parameter of the Saxon potential in fermis.

50.0 central real potential in MeV (is negative

for attractive potentials).

20.0 magnitude of the spin-orbit potential ex- pressed in number of times larger than for

the Thomas term of a nucleon. If a negative

number is given, it will be taken as V the

depth of the spin-orbit terra in MeV.

so

15.0 atomic mass in a.m.u. of the nucleus excluding the single neutron whose states are being considered.

8.0 maximum radius to which wave function will

be calculated in fermis (typically about

5 fra greater than the average radius of the

well).

,0001 the difference between the value of f . (inside) and f .(outside) must be less than EPS

for the eigenvalue to be considered determined.

normal search for bound or scattering states,

If a bound state, f .(outside) is determined

by integration inward from RMAX to the

matching radius RMATCH (which is automatically determined). For a scattering state

the matching is at RMAX and f .(outside) is

determined from the criterion for a 90°

phase shift.

f .(outside) is read in from the input data.

It is considered as the value of f .(outside)

at radius RMAX.

11
VARIABLE

MODE

NOMAX

TYPICAL VALUE MEANING

20

just calculate the wave function for the

given energy, EN. Do not search for a better eigenvalue.

maximum number of automatic search steps

seeking the eigenvalue. If N0MAX=1, one point is determined for the "signature of

the well."

NINTVL 100

JDBLD

JPRINT

EN

1

2

11.0

number of intervals in radius over which the wave function is to be determined. In the

example 101 radial points will be calculated, equally spaced between the first point

at R=0 and the last point at R=RMAX.

NINTVL must be a multiple of 10.

orbital angular momentum of the single

particle state sought.

twice the j -value of the state sought. In

the example j=-f so "j-doubled" is 3.

no printing of intermediate results.

print intermediate results for code checking.

value given to subroutine is the initial

energy. Value returned by subroutine is

the eigenvalue (if successful).

Value is in MeV.

FLOUT 4. read-in value of f .(outside) for the radius

RMAX.

12
Appendix 3. Fortran II listing.

The subroutine and their functions are:

MAINEV A main program for testing the EIGENV

program and a sample for how to call the

program.

EIGENV Controlling subroutine which performs matching and eigenvalue prediction.

VR3 Calculates V(r).

INTEG3 Controls integration in radius.

KRR1 Kutta-Runge integration subroutine.

FR Calculates f(r) for differential equation

d

3 u/dra

= f(r).

SIMPSON Simpson rule integration subroutine for

normalization of the wave function.

The DIMENSION and COMMON statements appearing in the subroutines

must also appear in the calling program.

13
* LABEL

CMAINEV

» LIST

C MAIN PROGRAM TO TEST EIGENV SINGLE PARTICLE EIGENVALUE AND

C WAVE FUNCTION PROGRAM

DIMENSION RPRINT(5) , WAVE (5) ,ENRGY( 10O),FL(lO0)

DIMENSION WAVEFN(201),VPLUS(3 99),XINTGD(201)

COMMON WAVEFN, RADIUS, A, VC, ALPHA, AMASS, RMAX, EPS, MODE, NOMA X,NI NT VL , 1L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,VPLUS,V,NMAX,

2F,RMATCH,I , XPLUS , XPLUSP ,FRR, INOUT, KSTEP , REFL , XINTGD , YINTG,XI NT 1 ,H,

3XINT2,X

5 WRITE OUTPUT TAPE 6,7

7 FORMAT(lHl)

READ INPUT TAPE 5,10

10 F0RMATI72H

1 ) WRITE OUTPUT TAPE 6,10

WRITE OUTPUT TAPE 6,4501

4501 F0RMAT(70H RADIUS A VC ALPHA(VSO) AMASS

1RMAX EPS ) READ INPUT TAPE 5 , 1 , RAD IUS, A, VC , ALPHA , AM ASS , RMAX, EPS

WRITE OUTPUT TAPE 6 , 1 ,RAD IUS, A, VC, ALPHA , AMASS, RMAX, EPS

1 F0RMAT(7F10.5)

WRITE OUTPUT TAPE 6,4503

4503 F0RMAT(70H ENERGIES MODE NOMAX NINTVL L J 1DBLD JPRINT ) READ INPUT TAPE 5 ,20 , KZ , MODE, NOMAX, N INT VL ,L , JDBLD , JPRINT

WRITEOUTPUTTAPE 6 ,2 0, KZ , MODE, NOMAX, M NT VL , L , JDBLD, JPRINT

20 F0RMAT(7I10)

NP0INT=NINTVL+1

GO TO (4081, 4082, 4081), MODE

4081 WRITE OUTPUT TAPE 6,4504

4504 F0RMATU5H TRIAL ENERGIES ) READ INPUT TAPE 5 , 1 , ( ENRGY ( KY ) , KY=1 ,KZ

)

WRITEOUTPUT TAPE6 , 1 ,( ENRGY < KY ), KY=1 ,KZ)

GO TO 4085

4082 WRITE OUTPUT TAPE 6,4083

4083 F0RMAT(25H ENERGY FLOUT, RMAX ) READ INPUT TAPE 5, 1 ,< ENRGY <KY) , FL ( KY ) ,KY=1 ,KZ

)

WRITE OUTPUT TAPE 6, 1, ( ENRGY ( KY ) ,FL ( KY ), KY=1 ,KZ

)

4085 DO 100 KY=1,KZ

EN=ENRGY(KY)

FLOUT=FL(KY)

WRITE OUTPUT TAPE 6,6180

6180 FORMAT( 68H ENERGY DIFFERENCE FL INSIDE FL OUT

1SIDE RMATCH ) CALL EIGENV

IF(ABSF(EN)-1.E-10) 100,4088,4088

4088 WRITE OUTPUTTAPE 6, 4089, EN

4089 FORMAT (34H WAVE FUNCTION VS RADIUS, ENERGY=, F10. 5 , 4H MEV)

R=-HI

14
INMAX=5

DO 200 N=1,NP0INT,5

DO 150 IN=1,INMAX

NN=N+IN-1

WAVE(IN)=WAVEFN(NN)

R=R+HI

RPRINT{ IN) =R

IF(NN-NPOINT) 150,150,155

150 CONTINUE

GO TO 159

155 INMAX=IN-1

159 WRITE OUTPUT TAPE 6 , 160, ( RPR I NT ( IN ) , WAVE ( IN )

,

IN = 1., INMAX)

160 F0RMAT(5(F9.3,E15.6))

200 CONTINUE

GO TO (100,170),JPRINT

170 WRITE OUTPUT TAPE 6,180,YINTG

180 F0RMAT(7H YINTG= , F10.5

)

100 CONTINUE

READ INPUT TAPE 5,20,NDUMMY

IF (NDUMMY) 1000, 1000,5

1000 CALL EXIT

END

• LABEL

CEIGENV

* LIST

SUBROUTINE EIGENV

DIMENSION Y(ll)

DIMENSION WAVEFN(201),VPLUS(399),XINTGD(201)

COMMON W AVEFN, RADIUS, A, VC, ALP HA, AMASS, RMAX, EPS, MODE, NOMAX,NI NTVL , 1L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,VPLUS,V,NMAX,

2F,RMATCH,I,XPLUS,XPLUSP,FRR,IN0UT,KSTEP,REFL,XINTGD,YINTG,XINT1,H,

3XINT2,X

K=L + 1 HI=RMAX/FLOATF(NINTVL)

NP0INT=NINTVL+1

IF(ALPHA) 1,2,2

1 ALPHA=ALPHA*181.1/VC

2 CALL VR3

KSTEP=1

H=HI

GO TO <15, 15, 6130), MODE

15 IF(EN) 180,180,170

170 RMATCH=RMAX

GO TO 190

180 RMATCH=(5.*FL0ATF(NINTVL/15+1))*HI

190 MPRINT=2

IF(NOMAX-l) 25,30,25

25 KSTEP=5

H=5«*HI

30 CALL INTEG3

N0 = 1 203 ENERGY=EN

IF<EN) 213,213,220

213 FLOMT=R«XPLUSP/XPLUS

DIF=REFL-FLOMT

GO TO 6001

220 PM0M=SQRTF(.04826*F*EN)

DIF=REFL-FL0UT

FL0MT=FL0UT

GO TO (2201, 6001, 6130), MODE

15
2201 XX=R*PMQM

KPLUS1=K+1

IF(XX-1.) 300,300,400

300 NPASS=1

301 SUM=1.

TERM=1.

DO 310 N=l,20

TERM=-TERM»0.5*XX**2/(FL0ATFtN)»FL0ATF(2*N-l-2*L)

)

SUM=SUM+TERM

IF(ABSF(TERM)-1.E-8)320,310,310

310 CONTINUE

320 IF(K-l) 325,335,325

325 00 330 IK=2,K

AK=IK

330 SUM=SUM*(2.*AK-3„)

335 GO Tl) (340,350),NPASS

340 YL=-SUM/(XX**K)

K=K + 1 L = L + 1 NPASS=2

GO TO 301

350 YLPLUS=-SUM/(XX**K)

K=K-1

L = L-1

GO TO 500

400 Y( 1)=-C0SF(XX)/XX

Y(2)=-C0SF(XX)/XX**2-SINF(XX)/XX

DO 410 IK=3,KPLUS1

AK = IK

410 Y(IK)=(2.*AK-3.)*Y( IK-1 ) /XX-Y ( IK-2

)

YL=Y(K)

YLPLUS=Y(KPLUS1)

500 FL0MT=FL0ATF(L)+1.-XX*YLPLUS/YL

DIF=REFL-FLOMT

6001 IF(NO-NOMAX) 6004,6100,6100

6C04 IF(N0-2) 6005,6080,6090

6005 EN=0 o 99*EN

GO TO 6100

6080 EN =(ENERG1*DIF-ENERGY*DIFM1)/(DIF-DIFM1 J GO TO 6100

6090 EN =(((ENERG2*DIFM1-ENERG1«DIFM2)»DIF/(DIFM1-DIFM2))-((ENERG1

1*DIF-ENERGY*DIFM1)*DIFM2/(DIF-DIFMD) ) / ( DIF-DIFM2 ) 6100 IF(ABSF((EN / ENERGY )- 1. )-. 2 ) 6110,6110,6101

6101 EN= ENERGY+ ABSF ( ENERGY )*S I GNF( 0,1, (EN- ENERGY)

)

6110 WRITE OUTPUT TAPE 6 , 6160, ENERGY ,DIF ,REFL , FLOMT ,R

6160 F0RMAT(F15.5,3E15.7,F8.3)

IF (NO-NOMAX) 6114,6140,6140

6114 IF(ABSF(DIF )-EPS ) 6115,6115,6120

6115 IF(KSTEP-l) 105,6116,105

6116 EN=ENERGY

GO TO 6130

6120 IF(KSTEP-l) 100,110,100

100 IF(ABSF(DIF )-.l) 105,110,110

105 KSTEP=1

H=HI

110 CALL INTEG3

DIFM2=DIFM1

DIFM1=DIF

ENERG2=ENERG1

ENERG1=ENERGY

16
N0=N0+1

GO TO 203

6130 MPRINT=1

CALL INTEG3

GO TO 6150

6140 EN = 0.

DO 6145 N*1,NP0INT

6145 WAVEFN(N)=0.

6150 CONTINUE

RETURN

END

• LABEL

CVR

* LIST

SUBROUTINE VR3

DIMENSION WAVEFN(201),VPLUS(399),XINTGD(201)

COMMON WAV EFN, RADIUS , A, VC r ALPHA, AMASS, RMAX, EPS, MODE, NOMAX ,N INTVL , 1L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,VPLUS,V,NMAX,

2F,RMATCH,I,XPLUS,XPLUSP,FRR,IN0UT,KSTEP,REFL,XINTGD,YINTG,XINT1,H,

3XINT2,X

R = HI

H=HI/2.

GO TO (2090, 2000), JPRINT

2000 WRITE OUTPUT TAPE 6,2091

2091 F0RMAT(25H V R ) 2090 C= ,0110270/A

NMAX=2*NINTVL-1

DO 2210 I=1,NMAX

REXP=EXPF( (R-RADIUS)/A)

VCC=VC/(1.0+REXP)

VS=ALPHA*VC*C*REXP/( ( ( 1 .O+REXP ) **2 ) *R

)

AK=K

IF(JD8LD-2*L) 22,22,21

21 VPLUSl I )=VCC+VS*(AK-1.Q)

GO TO 25

122

VPLUS ( I )=VCC-VS*AK

25 GO TO (40, 28), JPRINT

28 WRITE OUTPUT TAPE 6 , 30, VPLUSt I ) ,

R

30 FORMAT( E15.7*F7.3)

40 R=R+H

2210 CONTINUE

RETURN

END

'• LABEL

CINTEG

'• LIST

SUBROUTINE INTEG3

DIMENSION WAVEFN(201),VPLUS(399),XINTGD(201)

COMMON WAVEFN^RADIUS, A, VC^ALPHA, AMASS, RMAX, EPS, MODE, NOMAX, NINTVL,

1L, ELL, K,JDBLD,JPR INT, MPRINT.NPO INT, EN, FLOUT, H I, R,RR, VPLUS, V,NMAX,

2F,RMATCH,I,XPLUS,XPLUSP,FRR,IN0UT»KSTEP,REFL,XINTGD,YINTG,XINT1,H,

3XINT2,X

GO TO (1,15) ,MPRINT

1 GO TO (5,5,15),M0DE

5 IF(EN) 5215,4221,4221

15 WAVEFN(1)=0.

AK=K

F= AMASS/ (AMASS+1. 00898)

R=H

RR = R

17
V=VPLUS(2*KSTEP-1)

XPLUS=R**K

X=XPLUS

ELL=0.

CALL FR

ELL=K-1

XPLUS=XPLUS+FRR*R**2/( ( ELL+2. )*(ELL+3. ) ) WAVEFN(2)=XPLUS

XPLUSP =AK*(R **(K-1)) +FRR*R/ ( ELL+2. ) GO TO (4090, 4090, 4080)»M0DE

4080 N0MAX=2

4090 IN0UT=1

DC 4210 I=KSTEP,1000,KSTEP

1 = 1 CALL KRR1

R=R+H

4127 IF(NGMAX-l) 4130,4128,4130

4128 RMATCH=RMAX

IF IR-0.2*RADIUS ) 4130,4129,4129

4129 IF(EN-V) 10,10,4132

10 RMATCH=R

GO TO 4220

4130 WAVEFNC I+2)=XPLUS

4132 IFIR-RMATCH+.001) 4210,4220,422

4210 CONTINUE

4220 REFL=R*XPLUSP/XPLUS

WMATCH=XPLUS

IF(EN) 4211,4216,4216

4211 XPLUS=U

WAVEFN(NP0INT)=1.

GO TO (4212, 4213, 4213), MODE

4212 XPLUSP=-SQRTF(VPLUS(NMAX)-EN)

GO TO 4214

4213 XPLUSP=FLOUT/RMAX

4214 H=-H

R=RMAX

IN0UT=-1

DO 5210 J=1,10Q0,KSTEP

I=NINTVL-J

CALL KRR1

R = R + H WAVEFN( I+1)=XPLUS

IF(R-RMATCH-.OOl) 4215,4215,5210

5210 CONTINUE

4215 H=-H

GO TO (5215, 4230), MPRINT

5215 1 = 1 + 1 DO 5216 J=I,NP0INT

5216 WAVEFNt J)=WAVEFN( J ) *WMATCH/XPLUS

4216 GO TO (4221, 4230), MPRINT

4221 YINTG=0.

XINT1=0.

XINT2=RMAX

DO 4225 N=1,NP0INT

4225 XINTGD(N)=WAVEFN(N)*»2

CALL SIMPSN

DO 4228 N=1,NP0INT

4228 WAVEFN(N)=WAVEFN(N)/SQRTF(YINTG)

4230 RETURN

END

18
* LABEL

CKRR

• LIST

SUBROUTINE KRR1

DIMENSION WAVE FN ( 20 1),VPLUS( 399), XINTGD (201)

COMMON WAVEFN,RADIUS,A,VC,ALPHA,AMASS,RMAX,EPS,MODE,NOMAX,NINTVL,

1L,ELL,K,JDBLD, JPRINT, MPRINT.NPQINT, EN.FLOUT, HI,R,RR,VPLUS,V, NMAX,

2F, RMATCH, I, XPLUS, XPLUSP, FRR, INOUT , K STEP , REFL , XINTGD , YINTG, XI NT 1 ,H

,

3XINT2,X

GO TO (56,40) ,JPRINT

40 WRITE OUTPUT TAPE 6 ,45 ,H , R,XPLUS, XPLUSP , AI , EN , V, RMATCH

45 F0RMATL3H H=,F8.5,3H R=,F8.5,7H XPLUS=, E 12. 5 , 8H XPLUSP= ,E 12. 5 ,4H A 1I=,F8.5,4H EN=,F10.4,3H V=,F9.4,8H RMATCH=, F6. 3

)

56 RR=R

INDEX=2*I-IN0UT

V=VPLUS( INDEX)

X=XPLUS

CALL FR

AI=FRR*H**2/2.0

RR=R+H/2.0

X=XPLUS+ XPLUSP*H/2.0+AI/4.0

INDEX=2*I+INQUT*(KSTEP-1)

V=VPLUS(INDEX)

CALL FR

AH =FRR*(H*»2)/2oO

RR=R+H

X=XPLUS+XPLUSP*H+AI

I

INDEX=2*I+IN0UT*(2*KSTEP-1)

V=VPLUS( INDEX)

CALL FR

AI II=FRR*H**2/2.0

XPLUS=XPLUS+H*XPLUSP + (AI+2.0*AI I )/3.0

XPLUSP =XPLUSP +(AI +4.0*AII+AIII)/(3.0*H)

GO TO (4125, 190), JPRINT

190 IF(A8SF(RR-RMATCH)-.001) 195,4125,412 5 195 WRITE OUTPUT TAPE 6 ,45, H, RR ,X PLUS, XPLUSP , AI , EN , V, RMATCH

4125 CONTINUE

RETURN

END

'* LABEL

CFR

* LIST

SUBROUTINE FR

DIMENSION WAV EFN ( 201 ),VPLUS( 399), XINTGD (201)

COMMON WAVEFN,RADIUS,A,VC,ALPHA,AMASS,RMAX,EPS,MODE,NOMAX,NINTVL,

1L,ELL,K,JDBLD,-JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,VPLUS,V,NMAX,

2 F, RMATCH, I, XPLUS, XPLUSP, FRR, I NOUT,K STEP, REFL, XINTGD, YINTG, XI NT 1,H,

3XINT2,X

FRR=-.04826*F *(EN -20.72 1*ELL* ( ELL+1.0) /( F *RR**2)-V)*X

RETURN

END

* LABEL

CSIMPSN

* LIST

SUBROUTINE SIMPSN

DIMENSION W AVEFN ( 20 1),VPLUS( 399), XINTGD (201

)

COMMON WAVEFN,RADIUS,A,VC,.ALPHA,AMASS,RMAX,EPS,MODE,NOMAX,NINTVL,

lL,ELL,K,JDBLD,.JPRINT,MPRINT,NPOINT,EN,FLaUT,HI,R,RR,VPLUS,V,NMAX,

2 F, RMATCH, I , XPLUS , XPLUSP ,FRR, INOUT, K STEP , REFL , XINTGD , YINTG, XI NT 1 ,H,

3XINT2,X

19
DELINT=(XINT2-XINT1)/FL0ATF(NINTVL)

NPaiNT=NINTVL+l

YINTG3=XINTGD(1)

NPT=l

DO 20 NPT = 2,NINTVU2

YINTG3=yiNTG3+XINTGD{NPT)*4,0

20 CONTINUE

DO 30 NPT=3,NINTVL,2

YINTG3=YINTG3+XINTGD(NPT)*2.0

30 CONTINUE

YINTG3=YINTG3+XINTGQ(NP0INT)

YINTG=YINTG+YINTG3*DELINT/3.0

RETURN

END

20
Appendix 4. Fortran IV listing.

A block of COMMON called MEMEIG is common to MAINEV, EIGENV and to the subroutines called by EIGENV. It must appear in the calling program.

$IBFTC MAINEV LIST,M94,XR7

C MAIN PROGRAM TO TEST EIGENV SINGLE PARTICLE EIGENVALUE AND WAVE

C FUNCTION PROGRAM

DIMENSION RPRINT(5) , WAVE ( 5 ) , ENRGY

(

1C0),FL( IOC)

COMMON/MEMEIG/WAVEFN(201),RA0IUS,A,VC,ALPHA,AMASS,RMAX,EPS,MODE,

1N0MAX,NINTVL,L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,

2H,VPLUS(399) , V ,-NKAX , F„RMATCH, I , XPLUS , XPLUSP t FRR, INOUT, KSTEP,

3REFL,XINTGD(201) , YINTG, XINT1, XI NT2 ,

X

5 WRITE(6,7)

7 FORMAT! 1H1

)

READ(5,10)

10 FORMAT! 72H

1 ) WRITE(6,10)

WRITE(6,4501)

4501 F0RMAK70H RADIUS A VC ALPHA(VSO) AMASS

1RMAX EPS ) RE AD (5,1) RADIUS,A,VC,ALPHA,AMASS,RMAX,EPS

WRITE(6,l)RADIUS,A,VC,ALPHA,AMASS,RMAX,fcPS

1 FORMAT(7F10.5)

WRITE(6,4503)

4503 FORMAT! 70H ENERGIES MODE NOMAX NINTVL L J 1DBLD JPRINT ) READ (5,20) KZ, MODE, NOMAX ,-NINTVL , L , JD8LD, JPRINT

WRITE (6, 20) KZ,MODE,NOMAX,NINTVL,L, J D8LD, JPRINT

20 FORMAT(7I10)

NP0INT=NINTVL+1

GO TO (4081, 4082, 4081), MODE

4081 WRITE(6,4504)

4504 F0RMAT(15H TRIAL ENERGIES ) READ(5,1) (ENRGY(KY),KY=1,KZ)

WRITE(6,1) (ENRGY(KY) ,KY=1,KZ)

GO TU 4085

4082 WRIT£(6,4083)

4083 FORMAT! 25H ENERGY FLOUT, RMAX ) READ(5,1) ( ENRGY (KY),FL(KY),KY=1,KZ)

WRITE (6,1) (ENRGY(KY),FL(KY) ,KY=1,KZ

)

4085 DO 100 KY=1,KZ

EN=ENRGY(KY)

FLOUT=FL(KY)

WRITE (6, 6 180)

6180 FORMAT! 68H ENERGY DIFFERENCE FL INSIDE FL OUT

1SIDE RMATCH ) CALL EIGENV

IF(ABS (EN)-l.E-lO) 100,4088,4088

4088 WRITE(6,4089) EN

4089 FORMAT ( 34H WAVE FUNCTION VS RADIUS, ENERGY=, F 10. 5 , 4H MEV)

R=-HI

INMAX=5

DO 200 N=1,NPCINT,5

DO 150 IN=1,INMAX

NN=N+IN-1

21
WAVE(IN)=WAVEFN(NN)

R=R+HI

RPRINT( IN)=R

IF (NN-NPO I NT) 150,150,15 5 150 CONTINUE

GO TO 159

155 INMAX=IN-1

159 WRITc(6,160) ( RPRINT ( IN ) , WAVE ( I N) , I N=l, INMAX

)

160 F0RMAT(5(F9.3,E15.6) ) 200 CONTINUE

GO TO (100,170), JPRINT

170 WRITE(6,180) YINTG

180 F0RMAT(7H YINTG= , F 12, 5

)

100 CONTINUE

GO TO 5 END

SIBFTC EV LIST,M94,XR7

SUBROUTINE EIGENV

DIMENSION Y(ll)

C0MM0N/MEMEIG/WAVEFN(201),.RADIUS,A,VC,ALPHA,AMASS,RMAX,EPS,MODE,

1N0MAX,NINTVL,L,.ELL,K,JDBLD, JPRI NT, MPR INT ,NP0INT , EN, FLOUT,HI ,R,RR

,

2H,VPLUS(399) , V ,.NMAX , F*RMATCH, I, XPLUS , XPLUSP, FRR, INOUT,KSTEP,

3REFL,XINTGD(201) , Y I NTG, X I NT1 , X INT2,

X

K=L+1

Hl=RMAX/FLOAT(NINTVL)

NP0INT=NINTVL+1

IF(ALPHA) 1,2,2

1 ALPHA=ALPHA*181.1/VC

2 CALL VR3

KSTEP=1

H=HI

GO TO (15, 15,6130), MODE

15 IF(EN) 180,180,170

170 RMATCH=RMAX

GO TO 190

130 RMATCH=(5.»FL0AT(NINTVL/15+1) )*HI

190 MPRINT=2

IF(NOMAX-l) 25,30,25

25 KSTEP=5

H=5.»HI

30 CALL INTEG3

N0 = 1 203 ENERGY=EN

IF(EN) 213,213,220

213 FLOMT=R*XPLUSP/XPLUS

DIF=REFL-FL0MT

GO TO 6001

220 PM0M=SQRT(.04826*F»EN)

DIF=REFL-FLOUT

FLOMT=FLOUT

GO TO (2201, 6001, 6130), MODE

2201 XX=R*PMOM

KPLUS1=K+1

IFUX-1.) 300,300,400

300 NPASS=1

301 SUM=1.

TERM=1.

DO 310 N=l,20

TERM=-TERM*0.5#XX**2/( FLOAT (N)*FL0AT(2«N-1-2*L)

)

SUM=SUM+TERM

22
IF(ABS(TERM)-l.E-8) 320,310,310

310 CONTINUE

320 IF(K-l) 325,335,325

325 00 330 IK=2,K

AK=IK

330 SUM=SUM»(2.*AK-3.

)

335 GO TO (340,350),NPASS

340 YL=-SUM/(XX**K)

K=K+1

L = L+1

NPASS=2

GO TO 301

350 YLPLUS=-SUM/(XX**K)

K=K-1

L = L-1

GO TO 500

400 Y(1)=-C0S(XX)/XX

Y(2)=-COS(XX)/XX**2-SIN(XX)/XX

DO 410 IK=3,KPLUS1

AK = IK

, 410 Y( IK)=(2.*AK-3. )*Y( IK-1 ) /XX-Y ( IK-2

)

, YL=Y(K)

YLPLUS=Y(KPLUS1)

500 FL0MT=FL0AT(L)+1.-XX*YLPLUS/YL

DIF=REFL-FL0MT

5001 IF(NO-NOMAX) 6004,6100,6100

5004 IF(N0-2) 6005,6080,6090

5005 EN=0.99*EN

GO TO 6100

5080 EN =(ENERG1*DIF-ENERGY*DIFM1)/(DIF-DIFM1)

GO TO 6100

5090 EN =( <iENERG2*DIFMI-ENERGl*DIFM2)*DIF/(DIFMl-DIFM2))-((ENERGl

1*DIF-ENERGY*DIFM1)*DIFM2/(DIF-DIFM1)) ) / ( DIF-DI FM2

)

5100 IF(ABS((EN /ENERGY)-1. )-,2 ) 6110,6110,6101

5101 EN=ENERGY+ABS( ENERGY ) *S IGN< 0. 1, (EN-ENERGY)

)

5110 WRITE (6,6160) ENERGY, DIF ,REFL, FLOMT ,

R

6160 F0RMAT(F15.5,3E15.7,F8.3)

IF (NO-NOMAX) 6114,6140,6140

5114 IF(ABS(OIF )-EPS ) 6115^6115,6120

6115 IF(KSTEP-l) 105,6116,105

6116 EN=ENERGY

GO TO 6130

6120 IF(KSTEP-l) 100,110,100

100 IF(ABS(DIF )-.l) 105,110,110

105 KSTEP=1

H=HI

110 CALL INTEG3

DIFM2=DIFM1

DIFM1=0IF

ENERG2=ENERG1

ENERG1=ENERGY

N0=N0+1

GO TO 203

5130 MPRINT=1

CALL INTEG3

GO TO 6150

5140 EN=0.

DO 6145 N=1,NP0INT

6145 WAVEFN(N)=0.

6150 CONTINUE

23
RETURN

END

$IBFTC VR LIST>M94,XR7

SUBROUTINE VR3

COMMON/ MEME I G/WAVEFN ( 20 1),. RADIUS, A, VC, ALPHA, AMASS, RMAX, EPS, MODE,

1N0MAX,NINTVL,L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,

2H,VPLUS(399> , V ,.NKAX , FrRMATCH, I , XPLUS, XPLUSP, FRR , IN0UT,KSTEP,

3REFL,XINTGD(201) , YINTG, XINT1 ,XINT2,

X

R=HI

H=HI/2.

GO TO (2090,2000)»JPRINT

2000 WRITEI6,2091)

2091 FORMAT! 25H V R ) 2090 C= .0110270/A

NMAX=2*NINTVL-1

DO 2210 I=1,NMAX

REXP=EXP UR-RADIUS)/A)

VCC=VC/(1.0+REXP)

VS=ALPHA*VC*C*REXP/( ( ( 1 .Q+REXP) »»2 ) *R

)

AK=K

IF(JDBLD-2*L) 22,22,21

21 VPLUSU )=VCC+VS»(AK-1.0)

GO TO 25

22 VPLUS ( I)=VCC-VS*AK

25 GO TO (40,28),JPRINT

28 WRITE<6,30) VPLUSU)*R

30 FORMATl E15.7*F7.3)

40 R=R+H

2210 CONTINUE

RETURN

END

$IBFTC INTEG LIST,M94,XR7

SUBROUTINE INTEG3

COMMON/MEME IG/ WAVEFN( 201), RADIUS, A, VC , ALPHA, AMASS, RMAX, EPS, MODE

,

lNCMAX,NINTVL,L,ELL,K,JDBLD,JPRINT f MPRINT,NPOINT,EN,FLOUT,HI,R,RR,

2H,VPLUS(399) ,V„NMAX,F,.RMATCH, I, XPLUS, XPLUSP, FRR, IN0UT,KSTEP

,

3REFL,XINTGD(201) , Y INTG, X INT1 , XINT2,

X

GO TO (1,15),MPRINT

1 GO TO (5, 5,15), MODE

5 IF ( EN) 5215,4221,4221

15 WAVEFN(1)=0.

AK = K F=AMASS/( AMASS+ 1.00898)

R=H

RR = R V=VPLUS(2»KSTEP-1)

XPLUS=R«*K

X=XPLUS

ELL=0.

CALL FR

ELL=K-1

XPLUS=XPLUS+FRR*R«»2/((ELL+2. )*(ELL + 3.) ) WAVEFN(2)=XPLUS

XPLUSP =AK*(R **(K-D) +FRR*R/ ( ELL + 2. ) GO TO (4090, 4090,4080), MODE

4080 N0MAX=2

4090 IN0UT=1

DO 4210 I=KSTEP,1000,KSTEP

1=1

CALL KRR1

24
R = R+H

4127 IF(NOMAX-l) 4130,4128,4130

4128 RMATCH=RMAX

IF (R-0.2*RADIUS) 4130,4129,4129

4129 IF(EN-V) 10,10,4132

10 RMATCH=R

GO TO 4220

4130 WAVEFN(I+2)=XPLUS

4132 IF(R-RMATCH+.001) 4210,4220,4220

4210 CONTINUE

4220 REFL=R*XPLUSP/XPLUS

WMATCH=XPLUS

IF(EN) 4211,4216,4216

4211 XPLUS=1.

WAVEFN(NP0INT)=1.

GO TO (4212, 4213, 4213), MODE

4212 XPLUSP =-SQRT(VPLUS(NMAX)-EN)

GO TO 4214

4213 XPLUSP=FLOUT/RMAX

4214 H=-H

R=RMAX

IN0UT=-1

DO 5210 J=1,1000,KSTEP

I=NINTVL-J

CALL KRR1

R=R+H

WAVEFN( I+1)=XPLUS

IF(R-RMATCH-.OOl) 4215,4215,5210

5210 CONTINUE

4215 H=-H

GO TO (5215, 4230), MPRINT

5215 1=1+1

DO 5216 J=I,NPOINT

5216 WAVEFN(J)=WAVEFN(J)*WMATCH/XPLUS

4216 GO TO (4221, 4230), MPRINT

4221 YINTG=0.

XINT1=0.

XINT2=RMAX

DO 4225 N=1,NP0INT

4225 XINTGD(N)=WAVEFN(N)**2

CALL SIMPSN

DO 4228 N=1,NP0INT

4228 WAVEFN(N)=WAVEFN(N)/SQRT(YINTG)

4230 RETURN

END

iSIBFTC KRR LIST

SUBROUTINE KRR1 _ „ u4M cnc unr. c COMMON/MEMEIG/WAVEFNf 201), RADIUS, A, VC, ALPHA, AMASS, RMAX, EPS, MODE,

1N0MAX,NINTVL,L^ELL,K,JDBLD,JPRI NT, MPRINT, NPOINT, EN, FLOUT,HI,R,RR,

2H,VPLUS(399),V,NMAX,F rRMATCH,I, XPLUS , XPLUSP, FRR , INOUT,KSTEP,

3REFL,XINTGD(201) , Y INTG, XINT1, X I NT2,

X

GO TO (56,40) ,JPRINT

40 WRITE(6,45) H, R, XPLUS , XPLUSP, AI , EN, V, RMATCH

45 F0RMAT(3H H=,F8.5,3H R=,F9.5,7H XPLLS=, E12. 5, 8H XPLUSP= ,E12. 5,4H A

1I=,F9.5,4H EN=,F10.4,3H V=,F9.4,8H RMATCH=, F6.3

)

56 RR=R

INDEX=2*I-IN0UT

V=VPLUS( INDEX)

X=XPLUS

CALL FR

25
AI=FRR*H**2/2.0

RR=R+H/2,0

X=XPLUS+ XPLUSP»H/2.0+AI/4.0

INDEX=2*I+IN0UT*(KSTEP-1)

V=VPLUS(INDEX)

CALL FR

AH =FRR»(H**2)/2.0

RR=R+H

X = XPLUS+XPLUSP*H+AI I 1NDEX=2*I+IN0UT*(2*KSTEP-1)

V=VPLUS(INDEX)

CALL FR

AI II=FRR*H**2/2.0

XPLUS=XPLUS+H»XPLUSP+(AI+2.0*AII)/3a0

G0

L TO

P <:i2S

U

i q

P n;

(AI **-0-*II*Ai«I?lS.O«l ku TO (4125,190) ,JPRINT

190 IF(A8S(RR-RMATCH)-.00D 195,4125,4125

RETURN

END

$IBFTC FR LIST*M94,XR7

SUBROUTINE FR

C0MM0N/MEMEIG/WAVEFN(2Ol > RAnnic a ur *.„,.

3REFL,XINTGDI201),YINTG?«;n:« NT2 .x

' XPLUSP ' FRR ' '"OUT, KSTEP . RE?WN°

4826*F * <EN -^WJl.ELLMELL^.o./IF .RR» 2 ,-v,. X END

SIBFTC SIMPSN LIST,M94,XR7

SUBROUTINE SIMPSN

DEL INT=(XINT2-XINT1) /FLOAT (NINTVL) YlNrG3=XINTGD(l)

NPT=1

DO 20 NPT=2, NINTVL,

2

YINTG3=YINTG3+XINTG0(NPT)#4.0

20 CONTINUE

DO 30 NPT=3, NINTVL,

2

YINTG3=YINTG3+XINTGD(NPT)*2.0

30 CONTINUE

YINTG3=YINTG3+XINTGD(NP0INT)

YINTG=YINTG+YINTG3*DELINT/3.C

RETURN

END

-EPS, MODE,

26
Figure Captions

Results of the "signature of the well" run for the £=0 well

shown in Figure 3. Results are given for f (inside), f (outXj Hi side), and for the difference (DIF) between them. Note that

the energy at which the difference is zero is close to the

energy of the converged eigenvalue, indicated by the short

arrow. Improvement of the wave function match during the automatic

search procedure. Note that for proper normalization of the

bound state wave function that RMAX must be sufficiently large

2

that there is negligible contribution to m dr beyond RMAX.

Single particle states and the lsj, 2si, and ld-5 wave functions

8 2 2 in a Woods-Saxon potential well. The parameters of the well

are: V = -50 MeV, RADIUS =3.15 fermis, a = 0.650 fermis, and

c the spin-orbit term is 20 times the Thomas term for a nucleon

(ALPHA = 20). The shape of the well shown is for 4=0. During

the "signature of the well" preliminary run, the matching

radius is given by the heavy black line - the edge of the well

for bound states, and the maximum radius (8 fermis in this

example) for continuum states.

27
30N3d3dd
o

< gi- u-< ^UJ^ UJ_I

QUIh O;

X

o

I-<-2

HE

in

/fi J =n 'NOIlONnd 3AVM
FIGURE 3

3 4 5

RADIUS r, fm

»U.S. GOVERNMENT PRINTING OFFICE: 1966-251-019/27