//
//************************************************************************
//                                                                       *
//     This is PHIN, the nuclear bremsstrahlung package.  It contains    *
//     four modules:                                                     *
//                       1) DELTOL                                       *
//                       2) TIP                                          *
//                       3) BETHE                                        *
//                       4) PHIN                                         *
//                                                                       *
//       DELTOL(Z,DEL,TOL):                                              *
//               Subroutine. Provides values for DEL and TOL for         *
//               use by TIP.                                             *
//                                                                       *
//       TIP(T0,DEL,TOL,PHITIP,KMATCH)                                   *
//               Subroutine. Calculates the endpoint yeild (PHITIP) and  *
//               the photon energy (KMATCH) at which a straight line is  *
//               joined to the BETHE distribution. T0 is the incident    *
//               electron kinetic energy. DEL and TOL are matching       *
//               criteria.                                               *
//                                                                       *
//       BETHE(T0,K):                                                    *
//               Subroutine.  Calculates the BETHE distribution.         *
//                                                                       *
//       PHIN(T0,K,KMATCH,PHITIP)                                        *
//               Function routine. Calculates the nuclear bremsstrahlung *
//               yield at photon energy=K from electrons of kinetic      *
//               energy=T0.  If K .GT. KMATCH, it interpolates linearly  *
//               between (PHITIP,T0) and (BETHE(T0,KMATCH),KMATCH).      *
//               For K .LT. KMATCH, it simply calls BETHE.               *
//                                                                       *
//************************************************************************

// brem.h include 
//    1) bethe
//    2) phie
//    3) deltol
//    4) gamf


double bethe(double t0, double k)
{
  double E0,E,gamma,phi1,phi2,r,p0,p,eps0,eps,EE,pp;
  extern double FCOUL,ZTHIRD;
  double l,temp,temp0,tmp_bethe;
  
  double ms_e = 0.5110034;

  if( t0-0.00002 <= k ) return 0;

  E0 = t0 + ms_e;
  E  = E0-k;
  gamma = 100.*ms_e*k/(E0*E*ZTHIRD);

  // NON-SCREENED
  if(gamma>15.)
    {
      E0 = E0/ms_e;
      k  = k/ms_e;
      E  = E0-k;
      p0 = sqrt(E0*E*-1.);
      p  = sqrt(E*E-1.);
      
      temp0 = E0-p0;
      if(temp0==0) temp0 = 1./(2.*E0) + 1./(8.*pow(E0,3));
      temp = E-p;
      if(temp==0) temp = 1./(2.*E) + 1./(8.*pow(E,3));
      
      l    = 2.*log((E0*E+p0*p-1.)/k);
      eps0 = log((E0+p0)/temp0);
      eps  = log((E+p)/temp);
      EE   = E0*E;
      pp   = p0*p;
      
      tmp_bethe = p/p0*(4./3.-2.*EE*(p*p+p0*p0)/(pp*pp)
			+eps0*E/pow(p0,3)+eps*E0/pow(p,3)-eps*eps0/pp
     			+l*(k/(2.*pp)*(
				       (EE+p0*p0)/pow(p0,3)*eps0
				       -(EE+p*p)/pow(p,3)*eps
				       +2.*k*EE/pow(pp,2))
			    +8.*EE/(3.*pp)+k*k*(EE*EE+pp*pp)/pow(pp,3)));
      return tmp_bethe;
    }

  // SCREENED
  phi1 = 19.24-4.*log(gamma+2./(gamma+3.))-0.12*gamma*exp(-gamma/3.);
  phi2 = phi1;
  if (gamma <= 0.8) phi2=phi2-0.027-pow(0.8-gamma,2);
  r    = E/E0;
  tmp_bethe = 4.*((1.+r*r)*(1./4.*phi1-FCOUL)
		  -2./3.*r*(1./4.*phi2-FCOUL));
  return tmp_bethe;
}


//***********************************************************************
//									*
//	This is PHIE, the electron bremsstrahlung package.  It contains	*
//	one modules:							*
//			1) PHIE						*
//									*
//	PHIN(T0,K)							*
//		Function routine. Calculates the electron bremsstrahlung*
//		yeild at photon energy=K from electrons of kinetic	*
//		energy=T0.						*
//									*
//***********************************************************************

double phie(double t0, double k)
{
  double E0,kx,E,p0,kmax,rho,sqeta,eta,l,x;
  double epsi,psi;
  extern double Z23,ALOGZ23;
  double tmp_phie;

  double ms_e = 0.5110034;   // electron mass in MeV
  E0   = (t0+ms_e)/ms_e;     // units of electron mass
  kx   = k/ms_e;
  E    = E0-kx;
  p0   = sqrt(E0*E0-1.);
  kmax = (E0-1.)/(E0+1.-p0);

  if(kx>=kmax)  return 0;

  rho   = E0-kx*(1.+E0-p0)-1.;
  eta   = rho/(rho+2.);
  sqeta = sqrt(eta);
  l     = 2.*log((sqrt(E0-1.)+sqrt(eta*(E0+1.)))/
		 (sqrt(E0-1.)-sqrt(eta*(E0+1.))));
  x     = E/E0;
  tmp_phie = 2.*(1.-2./3.*x+x*x)*(l-sqeta)
    +sqeta*(1.-l*l/(2.*rho)-1./(rho*rho)
	    *pow(l/2.-sqrt(rho*(rho+2.)*(E0+1.)/(E0-1.)),2));

  // SCREENED
  epsi = 100.*kx/(E0*E*Z23);
  if(epsi>=0.88) return tmp_phie;

  x   = 0.88 - epsi;
  psi = 19.7+4.117*x-3.8056*x*x+31.835*pow(x,3)
    -58.63*pow(x,4)+40.774*pow(x,5);
  tmp_phie = tmp_phie*(0.25*psi-1.-ALOGZ23)/(3.798-log(epsi)-ALOGZ23);
  return tmp_phie;
}


void deltol(double z, double& del, double& tol)
{
  double iz;
  iz  = z;
  tol = 0.01;
  if(iz<=4)
    {
      del = 0.00001;
      tol = 0.02;
    }
  else if(iz<=6)
    {
      del = 0.0002;
      tol = 0.02;
    }
  else if(iz<=13)
    {
      del = 0.0005;
      tol = 0.02;
    }
  else if(iz<=29)
    {
      del = 0.0001;
    }
  else if(iz<73)
    {
      del = 0.002;
    }
  else
    {
      del = 0.003;
      tol = 0.0075;
    }
  return;
}

// GAMMA FUNCTION, 1<X<2 USING POLYNOMIAL NBS MATHEMATICAL HANDBOOK 
double gamf(double x)
{
  double tmpx = 1.-0.5748646*(x-1.)+0.9512363*pow(x-1.,2)
    -0.6998588*pow(x-1.,3)+0.4245549*pow(x-1.,4)
    -0.1010678*pow(x-1.,5);
  return tmpx;
}