FALSE [IF]

Program to solve the nonrelativistic Schroedinger equation
for bound states ( E < 0, L = 0 ) in the Woods-Saxon plus
repulsive Coulomb potential

---------------------------------------------------
    (C) Copyright 1999  Julian V. Noble.
      Permission is granted by the author to
      use this software for any application pro-
      vided this copyright notice is preserved,
      as per GNU Public License agreement.
---------------------------------------------------

This is an ANS Forth compatible program with the following
environmental dependence:

      ANS FLOAT and FLOAT EXT wordsets
      ANS TOOLS EXT wordsets

      Assumes separate floating point stack

      Uses a FORmula TRANslator for ease of porting
      to other languages

[THEN]


MARKER  -sch

include ftran111.f      \ load FORmula TRANslator
include ansfalsi.f      \ load root finder
include flocals.f       \ load fp locals wordset

: fvariables    ( n --)  0 DO  FVARIABLE  LOOP  ;

9  fvariables kappa alpha r0 r dr psi psi0 psi1 drsq
3  fvariables chi chi0 chi1


3  fvariables  Amass   Z   a0
12e0 Amass F!
12e0  Z F!          \ fake charge (should be 6)
0.5e0  a0 F!

: Coul   ( f: radius -- Coul)   \ point-charge Coulomb potential
    FRAME| aa |             \ local name for radius is aa
    f" 0.0695*Z / aa "
    |FRAME    ;             \ finish up flocals

: U     ( f: radius -- U)   \ Woods-Saxon + Coulomb
    FRAME| aa |             \ local name for radius is aa
    aa F@  r0 F@  F<        \ inside charge dist?
    IF      f" Coul(r0)*(1.5 - 0.5*(aa/r0)^2) "  \ compute inside Coul. pot.
    ELSE    f" Coul(aa) "                        \ compute tail Coul. pot.
    THEN
    f" -2.414/(1+exp((aa-r0)/a0) ) "   F+        \ and add to W-S
    |FRAME                  \ finish up flocals
;

: startup               \ initialize parameters
    f" r0 = 1.2*Amass^(1/3)"
    f" dr = 0.1"
    f" r = dr"
    f" drsq = dr^2"
    f" psi0 = 0"  f" psi1 = dr"
    f" chi0 = 0"  f" chi1 = dr"
;

: psi_step              \ integrate one step in Numerov scheme
    f" psi = ( 2+drsq*( U(r)+kappa^2) )*psi1 - psi0"
    f" chi = ( 2+drsq*( Coul(r)+kappa^2) )*chi1 - chi0"
    f" psi0 = psi1"     f" psi1 = psi"
    f" chi0 = chi1"     f" chi1 = chi"
    f" r=r+dr"
;

: display       \ for debugging, as in  display psi_step
    r f@ f.  psi1 f@ f.  psi0 f@ f.   chi1 f@ fs.
    f" psi1 / abs(chi1) " fs.
;


2  fvariables OldB  NewB        \ some extra workspace--could use stack

: beta  ( f: kappa -- B[k])
    kappa F!
    startup
    psi_step
    f" NewB = psi1/abs(chi1)"
    BEGIN   f" OldB = NewB "
            psi_step
            f" NewB = psi1/abs(chi1)"
            f" abs(NewB-OldB)/(NewB+OldB) "  1e-7 F<
    UNTIL
    NewB F@
;

: energy   ( f: kappa -- energy)   f^2   20.71e0 F*   ;

: norm       ( f: -- norm)
    startup
    0e0                                 \ sum = 0.0
    begin   r F@ 10e0 F<                \ r < 10 fm ?
    while   psi_step  f" psi1^2"  F+    \ sum = sum + psi^2
    repeat  dr F@  F*  FSQRT   ;        \ return sqrt(sum*dr)


: .wf                \ print normalized wave function
    norm                                \ compute norm
    FDUP F.                             \ display it and leave on fstack
    startup                             \ initialize
    BEGIN   r F@ 10e0 F<                \ r < 10 fm ?
    WHILE   CR   r F@ F.                \ output r
                 psi1 F@  FOVER F/  F.  \ and psi(r)
                 psi_step               \ next step
    REPEAT  FDROP                       \ drop norm
;