ccccccccccccccccccccccccc     Program 2.10     ccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM SCATTERING
C
C This is the main program for the scattering problem.
C
      PARAMETER(M =  21)
      PARAMETER(N = 10001)
      DIMENSION FI(N),THETA(M),SIG(M),SIG1(M)
      COMMON /CCNST/ B,E,A
      B0 = 0.01
      DB = 0.5
      DX = 0.01
      E  = 1.0
      A  = 100.0
      DL = 1.E-06
      DO      100  I = 1, M
        B  = B0+(I-1)*DB
C
C Calculate the first term of theta
C
        DO     80  J = 1, N
          X = B+DX*J
          FI(J) = 1.0/(X*X*SQRT(FB(X)))
   80   CONTINUE
        CALL SIMP(N,DX,FI,G1)
C
C Find r_m from 1-b*b/(r*r)-U/E=0
C
        X0  = B
        DX0 = DX
        CALL SECANT(DL,X0,DX0,ISTEP)
C
C Calculate the second term of theta
C
        DO     90   J = 1, N
          X = X0 + DX*J
          FI(J) = 1.0/(X*X*SQRT(F(X)))
   90   CONTINUE
        CALL SIMP(N,DX,FI,G2)
        THETA(I) = 2.0*B*(G1-G2)
  100 CONTINUE
C
C Calculate d_theta/d_b
C
      CALL THREE(M,DB,THETA,SIG,SIG1)
C
C Put the cross section in log form with the exact result of
C the Coulomb scattering (RUTH)
C
      DO      120  I = M, 1, -1
        B      = B0 + (I-1)*DB
        SIG(I) = B/ABS(SIG(I))/SIN(THETA(I))
        RUTH   = 1.0/SIN(THETA(I)/2.0)**4/16.0
        SI     = ALOG(SIG(I))
        RU     = ALOG(RUTH)
        WRITE (6, 999) THETA(I),SI,RU
  120 CONTINUE
      STOP
  999 FORMAT (3F16.8)
      END
C
      FUNCTION F(X)
      COMMON /CCNST/ B,E,A
        F = 1.0-B*B/(X*X)-U(X)/E
      RETURN
      END
C
      FUNCTION FB(X)
      COMMON/CCNST/ B,E,A
        FB = 1.0-B*B/(X*X)
      RETURN
      END
C
      FUNCTION U(X)
      COMMON/CCNST/ B,E,A
        U = 1.0/X*EXP(-X/A)
      RETURN
      END
C
      SUBROUTINE SECANT (DL,X0,DX,ISTEP)
C
C Subroutine for the root of f(x)=0 with the secant method.
C
      ISTEP = 0
      X1 = X0 + DX
      DO    100  WHILE (ABS(DX).GT.DL)
        D  = F(X1) - F(X0)
        X2 = X1 - F(X1)*(X1-X0)/D
        X0 = X1
        X1 = X2
        DX = X1 - X0
        ISTEP = ISTEP + 1
  100 END DO
      RETURN
      END
C
      SUBROUTINE THREE(N,H,FI,F1,F2)
C
C Subroutine for 1st and 2nd order derivatives with the three-
C point formulas. Extrapolations are made at the boundaries.
C FI: input f(x); H: interval; F1: f'; and F2: f"
C
      DIMENSION FI(N),F1(N),F2(N)
C
C f' and f" from three-point formulas
C
      DO      100  I = 2, N-1
        F1(I) = (FI(I+1) - FI(I-1))/(2.*H)
        F2(I) = (FI(I+1) - 2.*FI(I) + FI(I-1))/(H*H)
  100 CONTINUE
C
C Linear extrapolation for the boundary points
C
      F1(1) = 2.*F1(2) - F1(3)
      F1(N) = 2.*F1(N-1) - F1(N-2)
      F2(1) = 2.*F2(2) - F2(3)
      F2(N) = 2.*F2(N-1) - F2(N-2)
      RETURN
      END
C
      SUBROUTINE SIMP (N,H,FI,S)
C
C Subroutine for integration over f(x) with the Simpson rule.
C FI: integrand f(x); H: interval; S: integral.
C
      DIMENSION FI(N)
C
      S  = 0.
      S0 = 0.
      S1 = 0.
      S2 = 0.
      DO      100 I = 2, N-1, 2
        S1 = S1 + FI(I-1)
        S0 = S0 + FI(I)
        S2 = S2 + FI(I+1)
  100 CONTINUE
      S = H*(S1 + 4.*S0 + S2)/3.
C
C If N is even, add the last slice separately
C
      IF(MOD(N,2).EQ.0) S = S
     *  + H*(5.*FI(N) + 8.*FI(N-1) - FI(N-2))/12.
      RETURN
      END