/****************************************************************/
/* Program to calculate the elastic and total cross-section     */
/* of neutron-nuclide reaction using square well Optical Model. */
/* See 22.113 Problem Set 4.1           Li Ju. Apr.10,1995      */
/****************************************************************/

#include <iostream.h>
#include <stdio.h>
#include <math.h>
#include "complex.h"

#define  V0    42.       /* Real depth of the well = 42 Mev       */
#define  W0    1.26      /* Imaginary part is 1.26 Mev            */
#define  zeta  0.03      /* Ratio of imaginary to real is 0.03    */
#define  A     63.54     /* Copper atom as nuclide                */
#define  rr    1.45E-13  /* radius of the well per amu in cm      */
#define  kk    2.198E12  /* wave vector unit in cm^-1 for neutron */
                         /* per Mev^(1/2)                         */

#define  LMAX 20         /* MAX partial wave number               */
#define  PI 3.141592653589793
#define IMAX(a,b) ((a) > (b) ? (a) : (b))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

/*
** These are copied from Numerical Recipes.
** Because direct linking to C++ didn't work. (maybe I'm stupid)
*/
float chebev(float a,float b,float c[],int m,float x);
void sphbes(int n,float x,float *sj,float *sy,float *sjp, float *syp);
void bessjy(float x,float xnu,float *rj,float *ry,float *rjp,float *ryp);
void beschb(double x,double *gam1,double *gam2,double *gampl,double *gammi);
void nrerror(char a[])
{ 
return;
}


void main()
{
  FILE * OUTPUT;
  double Delta[LMAX+1],s[LMAX+1];
  Complex qint[LMAX+1];
  double  tau[LMAX+1];
  Complex eita[LMAX+1];
  double  sigma_T,sigma_EL;     
  /* Total and Elastic cross-section in Barn = 1E-24 cm^2 */

  int L;
  double a,b,c;
  float  jL0,jL1,nL0,nL1;
  double E,k,xi;  /* Incoming energy in Mev, k=sqrt(E)*kk */
  Complex K,Z;

  OUTPUT = fopen ("cs.out","w");

  for (E=0.01;E<=30.;E+=(E<2.)?0.05:0.1)
    {
      K  = sqrt(Complex(E+V0,W0))*kk;  /* in cm^-1 now */
      Z  = K*cbrt(A)*rr;   /* cbrt(A)*rr is the width of the well in cm */ 
      k  = sqrt(E)*kk;     /* Real outgoing wave vector in cm^-1        */
      xi = k*cbrt(A)*rr;   /* k*r0                                      */
  
      Delta[0] = 0.;
      s[0]     = xi;

      for (L=1;L<=LMAX;L++)
	{
	  a = (L-Delta[L-1])*(L-Delta[L-1])+s[L-1]*s[L-1];
	  Delta[L] = xi*xi*(L-Delta[L-1])/a-L;
	  s[L] = xi*xi*s[L-1]/a;
	}
        
      qint[0] = Z*cos(Z)/sin(Z);
  
      for (L=1;L<=LMAX;L++)
	qint[L] = Z*Z/(L-qint[L-1])-L;
	
      sigma_T  = 0.;
      sigma_EL = 0.;
      for (L=0;L<=LMAX;L++)
	{
	  sphbes (L,xi,&jL0,&nL0,&jL1,&nL1);
	  /* regular and irregular spherical Bessel function */

	  tau [L] = atan(jL0/nL0);
	  eita[L] = exp(Complex(0,-2*tau[L]))*(qint[L]-Complex(Delta[L],-s[L]))
                    /(qint[L]-Complex(Delta[L],s[L]));
	  sigma_T  = sigma_T  + (2*L+1)*(1-eita[L].real);
	  sigma_EL = sigma_EL + (2*L+1)*(((1-eita[L])*(1-eita[L])).Abs());
	}

      sigma_T  *= 2*PI/k/k*1E24;   
      sigma_EL *= PI/k/k*1E24;     
      /* will be in barn */
      
      fprintf (OUTPUT,"%f  %f  %f\n",E,sigma_T,sigma_EL);
    }
}




/*************************************************/
/* Down below all copied from Numerical Recipes  */
/*************************************************/

#define RTPIO2 1.2533141
void sphbes(int n,float x,float *sj,float *sy,float *sjp, float *syp)
{
	float factor,order,rj,rjp,ry,ryp;

	if (n < 0 || x <= 0.0) nrerror("bad arguments in sphbes");
	order=n+0.5;
	bessjy(x,order,&rj,&ry,&rjp,&ryp);
	factor=RTPIO2/sqrt(x);
	*sj=factor*rj;
	*sy=factor*ry;
	*sjp=factor*rjp-(*sj)/(2.0*x);
	*syp=factor*ryp-(*sy)/(2.0*x);
}
#undef RTPIO2

#define EPS 1.0e-10
#define FPMIN 1.0e-30
#define MAXIT 10000
#define XMIN 2.0

void bessjy(float x,float xnu,float *rj,float *ry,float *rjp,float *ryp)
{
	int i,isign,l,nl;
	double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2,
		fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl,
		rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1,
		temp,w,x2,xi,xi2,xmu,xmu2;

	if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessjy");
	nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu-x+1.5)));
	xmu=xnu-nl;
	xmu2=xmu*xmu;
	xi=1.0/x;
	xi2=2.0*xi;
	w=xi2/PI;
	isign=1;
	h=xnu*xi;
	if (h < FPMIN) h=FPMIN;
	b=xi2*xnu;
	d=0.0;
	c=h;
	for (i=1;i<=MAXIT;i++) {
		b += xi2;
		d=b-d;
		if (fabs(d) < FPMIN) d=FPMIN;
		c=b-1.0/c;
		if (fabs(c) < FPMIN) c=FPMIN;
		d=1.0/d;
		del=c*d;
		h=del*h;
		if (d < 0.0) isign = -isign;
		if (fabs(del-1.0) < EPS) break;
	}
	if (i > MAXIT) nrerror("x too large in bessjy; try asymptotic expansion");
	rjl=isign*FPMIN;
	rjpl=h*rjl;
	rjl1=rjl;
	rjp1=rjpl;
	fact=xnu*xi;
	for (l=nl;l>=1;l--) {
		rjtemp=fact*rjl+rjpl;
		fact -= xi;
		rjpl=fact*rjtemp-rjl;
		rjl=rjtemp;
	}
	if (rjl == 0.0) rjl=EPS;
	f=rjpl/rjl;
	if (x < XMIN) {
		x2=0.5*x;
		pimu=PI*xmu;
		fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu));
		d = -log(x2);
		e=xmu*d;
		fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e);
		beschb(xmu,&gam1,&gam2,&gampl,&gammi);
		ff=2.0/PI*fact*(gam1*cosh(e)+gam2*fact2*d);
		e=exp(e);
		p=e/(gampl*PI);
		q=1.0/(e*PI*gammi);
		pimu2=0.5*pimu;
		fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2);
		r=PI*pimu2*fact3*fact3;
		c=1.0;
		d = -x2*x2;
		sum=ff+r*q;
		sum1=p;
		for (i=1;i<=MAXIT;i++) {
			ff=(i*ff+p+q)/(i*i-xmu2);
			c *= (d/i);
			p /= (i-xmu);
			q /= (i+xmu);
			del=c*(ff+r*q);
			sum += del;
			del1=c*p-i*del;
			sum1 += del1;
			if (fabs(del) < (1.0+fabs(sum))*EPS) break;
		}
		if (i > MAXIT) nrerror("bessy series failed to converge");
		rymu = -sum;
		ry1 = -sum1*xi2;
		rymup=xmu*xi*rymu-ry1;
		rjmu=w/(rymup-f*rymu);
	} else {
		a=0.25-xmu2;
		p = -0.5*xi;
		q=1.0;
		br=2.0*x;
		bi=2.0;
		fact=a*xi/(p*p+q*q);
		cr=br+q*fact;
		ci=bi+p*fact;
		den=br*br+bi*bi;
		dr=br/den;
		di = -bi/den;
		dlr=cr*dr-ci*di;
		dli=cr*di+ci*dr;
		temp=p*dlr-q*dli;
		q=p*dli+q*dlr;
		p=temp;
		for (i=2;i<=MAXIT;i++) {
			a += 2*(i-1);
			bi += 2.0;
			dr=a*dr+br;
			di=a*di+bi;
			if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN;
			fact=a/(cr*cr+ci*ci);
			cr=br+cr*fact;
			ci=bi-ci*fact;
			if (fabs(cr)+fabs(ci) < FPMIN) cr=FPMIN;
			den=dr*dr+di*di;
			dr /= den;
			di /= -den;
			dlr=cr*dr-ci*di;
			dli=cr*di+ci*dr;
			temp=p*dlr-q*dli;
			q=p*dli+q*dlr;
			p=temp;
			if (fabs(dlr-1.0)+fabs(dli) < EPS) break;
		}
		if (i > MAXIT) nrerror("cf2 failed in bessjy");
		gam=(p-f)/q;
		rjmu=sqrt(w/((p-f)*gam+q));
		rjmu=SIGN(rjmu,rjl);
		rymu=rjmu*gam;
		rymup=rymu*(p+q/gam);
		ry1=xmu*xi*rymu-rymup;
	}
	fact=rjmu/rjl;
	*rj=rjl1*fact;
	*rjp=rjp1*fact;
	for (i=1;i<=nl;i++) {
		rytemp=(xmu+i)*xi2*ry1-rymu;
		rymu=ry1;
		ry1=rytemp;
	}
	*ry=rymu;
	*ryp=xnu*xi*rymu-ry1;
}
#undef EPS
#undef FPMIN
#undef MAXIT
#undef XMIN
#undef PI

#define NUSE1 5
#define NUSE2 5

void beschb(double x,double *gam1,double *gam2,double *gampl,double *gammi)
{
	float xx;
	static float c1[] = {
		-1.142022680371172e0,6.516511267076e-3,
		3.08709017308e-4,-3.470626964e-6,6.943764e-9,
		3.6780e-11,-1.36e-13};
	static float c2[] = {
		1.843740587300906e0,-0.076852840844786e0,
		1.271927136655e-3,-4.971736704e-6,-3.3126120e-8,
		2.42310e-10,-1.70e-13,-1.0e-15};

	xx=8.0*x*x-1.0;
	*gam1=chebev(-1.0,1.0,c1,NUSE1,xx);
	*gam2=chebev(-1.0,1.0,c2,NUSE2,xx);
	*gampl= *gam2-x*(*gam1);
	*gammi= *gam2+x*(*gam1);
}
#undef NUSE1
#undef NUSE2

float chebev(float a,float b,float c[],int m,float x)
{
	float d=0.0,dd=0.0,sv,y,y2;
	int j;

	if ((x-a)*(x-b) > 0.0) nrerror("x not in range in routine chebev");
	y2=2.0*(y=(2.0*x-a-b)/(b-a));
	for (j=m-1;j>=1;j--) {
		sv=d;
		d=y2*d-dd+c[j];
		dd=sv;
	}
	return y*d-dd+0.5*c[0];
}