module CB
   real, unit(b) :: B
   real, unit(e) :: E
   real, unit(a) :: A
end module CB
program SCATTERING
   use CB
   implicit none
   integer, parameter , unit(m) :: M = 21
   integer, parameter  :: N = 10001
   integer :: I,J,ISTEP
   real :: DL,B0,DB,DX,X0,X,DX0,F,FX,FB,FBX,G1,G2,RU,RUTH,SI
   real, dimension (N) :: FI
   real, dimension (M) :: THETA,SIG,SIG1
   DL = 1.E-06
   B0 = 0.01
   DB = 0.5
   DX = 0.01
   E = 1.0
   A = 100.0
   do I = 1, M, 1
      B = B0+(I-1)*DB
      do J = 1, N, 1
         X = B+DX*J
         FI(J) = 1.0/((X*X)*sqrt(FBX(X)))
      end do
      call SIMP(N,DX,FI,G1)
      X0 = B
      DX0 = DX
      call SECANT(DL,X0,DX0,ISTEP)
      do J = 1, N, 1
         X = X0+DX*J
         FI(J) = 1.0/((X*X)*sqrt(FX(X)))
      end do
      call SIMP(N,DX,FI,G2)
      THETA(I) = (2.0*B)*(G1-G2)
   end do
   call THREE(M,DB,THETA,SIG,SIG1)
   do 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 (number = 6,'"(3F16.8)"') THETA(I),SI,RU
   end do
end program SCATTERING
subroutine SIMP(N,H,FI,S)
   implicit none
   integer, intent(in)  :: N
   integer :: I
   real, intent(in)  :: H
   real :: S0,S1,S2
   real, intent(out)  :: S
   real, intent(in) , dimension (N) :: FI
   S = 0.0
   S0 = 0.0
   S1 = 0.0
   S2 = 0.0
   do I = 2, N-1, 2
      S1 = S1+FI(I-1)
      S0 = S0+FI(I)
      S2 = S2+FI(I+1)
   end do
   S = (H*((S1+4.0*S0)+S2))/3.0
   if (MOD(N,2)==0) then