'This program calculates the differential scattering cross section for free
'electrons incident on a central potential.  The potential is a model
'potential for an atomic hydrogen target.  The electron energy e is 100 eV.
'The electron energy can easily be changed.
'This program outputs the cross section for angles between 0 and 180
'degrees in units of angstrom^2.

DECLARE SUB bessel ()
COMMON SHARED offset%, lmax%, x
'These variables are shared with the subroutine bessel.
CLS

'Here constants are declared.
pi = 3.14159:    hbarm = 7.62003:    a0 = .529
'hbarm =(hbar^2)/m in units of eV*angstrom^2.
rmax = 3:         rext = 3.2:                 dr = .000002
'The constants rmax, rext, and dr may be changed to test the program.
'The constants rmax and rext should be chosen bigger for bigger atoms.
nr = rmax / dr:  nrext = rext / dr:           nang% = 36
e = 1000
'The program can handle electron energies up to 1000 eV.
k2 = 2 * e / hbarm: k = SQR(k2)
lmax% = k * rmax / 2 + 4
'For l%>lmax% the phase shift delta(l%) will be zero.
lstart% = 0

'Here we declare arrays and double precision variables.
DIM pl(lmax%, nang%), dsigma(nang%)
DIM SHARED jl(2 * lmax% + 1), nl(2 * lmax% + 1), delta(lmax%)
DIM psi1 AS DOUBLE, psi2 AS DOUBLE, psi3 AS DOUBLE, ps AS DOUBLE
DIM psinr AS DOUBLE


'This is a routine to calculate and store the Legendre polynomials for all
'l% from zero to lmax% and for all angles fom 0 to 180 degrees in 5 degree
'intervals.  The routine uses a recursion formula.
FOR i% = 0 TO nang%
        theta = i% * pi / nang%
        x = COS(theta)
        pl(0, i%) = 1
        pl(1, i%) = x
        FOR l% = 1 TO lmax% - 1
                pli = (2 * l% + 1) * x * pl(l%, i%)
                pl(l% + 1, i%) = (pli - l% * pl(l% - 1, i%)) / (l% + 1)
        NEXT l%
NEXT i%

'This routine calculates the energy-dependent variables. Here you
'can change the energy.
e = 100
k2 = 2 * e / hbarm: k = SQR(k2)
lmax% = k * rmax / 2 + 4
'Smaller energies get a smaller lmax%.
PRINT "lmax = ", lmax%

'Here we calculate the Bessel fuctions of rmax and rext for all l%
'from zero to lmax% by calling a subroutine. The subroutine uses
'a recursion formula.
x = k * rmax
offset% = 0
CALL bessel
x = k * rext
offset% = lmax% + 1
CALL bessel

'Here we open a file to store our cross section.
a$ = "e" + RIGHT$(STR$(e), 3) + ".dat"
OPEN a$ FOR OUTPUT AS #1
PRINT "l", "delta"
PRINT #1, "l", "delta"

'routine to integrate the Schroedinger equation using the Numerov method
cons = dr * dr / (hbarm * 6)
FOR l% = lstart% TO lmax%
PRINT l%,
PRINT #1, l%,
'We have to do a seperate integration for each l%.
'First we plant the seeds to start the Numerov method.
        ll = l% * (l% + 1) * hbarm / 2
        psi1 = 0: psi2 = 1D-25
        r = 4 * dr
        v = 14.3998 * (1 / r + 1 / a0) * EXP(-2 * r / a0)
        v = v + ll / (r * r)
        kim1 = cons * (e - v)
        r = 5 * dr
        v = 14.3998 * (1 / r + 1 / a0) * EXP(-2 * r / a0)
        v = v + ll / (r * r)
        ki = cons * (e - v)
'Now we integrate using the Numerov method.
        FOR i = 5 TO nrext - 1
                r = (i + 1) * dr
                IF i + 1 > nr THEN v = 0: GOTO outer
                v = 14.3998 * (1 / r + 1 / a0) * EXP(-2 * r / a0)
outer:          v = v + ll / (r * r)
                kip1 = cons * (e - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
                IF i + 1 = nr THEN psinr = psi3
        NEXT i


'Here we calculate the phase shifts. This is where we need the previously
'computed Bessel functions.
        g = rmax * psi3 / (rext * psinr)
        numer = g * jl(l%) - jl(l% + lmax% + 1)
        denom = g * nl(l%) - nl(l% + lmax% + 1)
        delta(l%) = ATN(numer / denom)
        PRINT delta(l%)
        PRINT #1, delta(l%)
NEXT l%


'Now that we have all the phase shifts we can calculate the differential
'cross section for a given energy as a function of angle.
'This routine calculates dsigma/domega.
FOR i% = 0 TO nang%
FREAL = 0: FIMAG = 0
FOR l% = lstart% TO lmax%
        SD = SIN(delta(l%))
        CD = COS(delta(l%))
        FREAL = FREAL + (2 * l% + 1) * SD * CD * pl(l%, i%)
        FIMAG = FIMAG + (2 * l% + 1) * SD * SD * pl(l%, i%)
NEXT l%
FREAL = FREAL / k
FIMAG = FIMAG / k
dsigma(i%) = FREAL * FREAL + FIMAG * FIMAG
NEXT i%


'Here we write the scattering cross sections to a file.
PRINT "angle", "dsigma"
PRINT #1, "angle", "dsigma"
FOR i% = 0 TO nang%
PRINT i% / .2, dsigma(i%)
PRINT #1, i% / .2, dsigma(i%)
NEXT i%
CLOSE


SUB bessel
'subroutine to calculate the spherical Bessel functions
nl(offset%) = -COS(x) / x
nl(1 + offset%) = -(COS(x) / (x * x)) - SIN(x) / x
FOR l% = (1 + offset%) TO (lmax% - 1 + offset%)      'forward recursion
        ltrue% = l% - offset%
        nl(l% + 1) = (2 * ltrue% + 1) * nl(l%) / x - nl(l% - 1)
NEXT l%
lupper% = x + 10
jp1 = 0
j = 9.999999E-21
FOR l% = (lupper% + offset%) TO (1 + offset%) STEP -1   'backward recursion
        ltrue% = l% - offset%
        jm1 = (2 * ltrue% + 1) * j / x - jp1
        IF (l% - 1) <= (lmax% + offset%) THEN jl(l% - 1) = jm1
        jp1 = j: j = jm1
NEXT l%
norm = SIN(x) / (x * jl(offset%))
FOR l% = (offset%) TO (lmax% + offset%)
        jl(l%) = norm * jl(l%)
NEXT l%
END SUB