/*************************    Program 2.10   ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define f(x) (1-b*b/(x*x)-exp(-x/a)/(x*e))
#define fb(x) (1-b*b/(x*x))
#define NMAX 10001
#define MMAX    21
const double a=100,e=1;
double b;

main()
/* This is the main program for the scattering problem.
   Copyright (c) Tao Pang 1997. */
{
int i,j,n,m,istep;
double fi[NMAX],theta[MMAX],sig[MMAX],sig1[MMAX];
double x,x0,g1,g2,b0,db,dx,df;
double si,ruth,ru;
double dl=1e-6;
extern double b;
void simp();
void secant();
void three();

n = NMAX;
m = MMAX;
b0 = 0.01;
db = 0.5;
dx = 0.01;
for (i = 0; i < m; ++i)
  {
  b = b0+i*db;
   
/* Calculate the first term of theta */
 
  for (j = 0; j < n; ++j)
    {
    x = b+dx*(j+1);
    fi[j] = 1/(x*x*sqrt(fb(x)));
    }
  simp (n,dx,fi,&g1);

/* Find r_m from 1-b*b/(r*r)-U/E=0 */
 
   secant (dl,b,dx,&x0,&df,&istep);

/* Calculate the second term of theta */

  for (j = 0; j < n; ++j)
    {
    x = x0+dx*(j+1);
    fi[j] = 1/(x*x*sqrt(f(x)));
    }
  simp (n,dx,fi,&g2);
  theta[i] = 2*b*(g1-g2);
  }

/* Calculate d_theta/d_b */

three (m,db,theta,sig,sig1);

/* Put the cross section in log form with the exact result of
   the Coulomb scattering (RUTH) */

for (i = m-1; i >= 0; i = i-1)
  {
  b = b0+i*db;
  sig[i] = b/(fabs(sig[i])*sin(theta[i]));
  ruth = 1/(pow(sin(theta[i]/2),4)*16);
  si = log(sig[i]);
  ru = log(ruth);
  printf("%16.8lf %16.8lf %16.8lf\n", theta[i],si,ru);
  }
}

void simp (n,h,f,s)
/* Subroutine for integration over f(x) with the Simpson rule.
   f: integrand f(x); h: interval; s: integral.
   Copyright (c) Tao Pang 1997. */
int n;
double h,*s;
double f[];
{
int i;
double s0,s1,s2;

s0 = 0;
s1 = 0;
s2 = 0;
for (i = 1; i < n-1; i = i+2)
  {
  s0 = s0+f[i];
  s1 = s1+f[i-1];
  s2 = s2+f[i+1];
  }
*s = h*(s1+4*s0+s2)/3;

/* If n is even, add the last slice separately */

if (n%2 == 0)
  {
  *s = *s+h*(5*f[n-1]+8*f[n-2]-f[n-3])/12;
  }
}

void secant (dl,xi,di,xf,df,istep)
/* Subroutine for the root of f(x)=0 with the secant method.
   Copyright (c) Tao Pang 1997. */
int *istep;
double dl,xi,di;
double *xf,*df;
{
double x0,x1,x2,d;

x0 = xi;
x1 = xi+di;
*df = di;
*istep = 0;
while (fabs(*df) > dl)
  {
  d = f(x1)-f(x0);
  x2 = x1-f(x1)*(x1-x0)/d;
  x0 = x1;
  x1 = x2;
  *df = x1-x0;
  *istep = *istep+1;
  }
  *xf = x0;
}

void three(n,h,f,f1,f2)
/* Subroutine for 1st and 2nd order derivatives with the three-
   point formulas. Extrapolations are made at the boundaries.
   FI: input f(x); H: interval; F1: f'; and F2: f". 
   Copyright (c) Tao Pang 1997. */
int n;
double h;
double f[],f1[],f2[];
{
int i;
/* f' and f" from three-point formulas */

for (i = 1; i < n-1; ++i)
  {
  f1[i] = (f[i+1]-f[i-1])/(2*h);
  f2[i] = (f[i+1]-2*f[i]+f[i-1])/(h*h);
  }
/* Linear extrapolation for the boundary points */

f1[0]   = 2*f1[1]-f1[2];
f1[n-1] = 2*f1[n-2]-f1[n-3];
f2[0]   = 2*f2[1]-f2[2];
f2[n-1] = 2*f2[n-2]-f2[n-3];
}