double FCOUL,ZTHIRD,Z137,Z137SQ,Z23,ALOGZ23,ALOGZ3;
double CTIP1,CTIP2,CTIP3,CTIP4,DEL,TOL;

#include "brem.h"

double photon(double ebeam, double egam, double trad, double zrad, double arad)
{  
  if(egam>ebeam) 
    {
      printf("!!!You sould not see this message\n");
      return 0;
    }
  
  extern double FCOUL,ZTHIRD,Z137,Z137SQ,Z23,ALOGZ23,ALOGZ3;
  extern double CTIP1,CTIP2,CTIP3,CTIP4,DEL,TOL;
  double PI,piz137,g,f,ce,cn;

  //generate a bunch of terms that are a function of Z only
  PI     = 3.14159265359;
  Z137   = zrad/137.;
  piz137 = PI*Z137;
  Z137SQ = Z137*Z137;
  ALOGZ3 = 1./3.*log(zrad);
  ALOGZ23= ALOGZ3*2;
  ZTHIRD = pow(zrad,1./3.);
  Z23    = pow(zrad,2./3.);

  FCOUL  = Z137SQ*(1./(1.+Z137SQ)+0.20206-0.0369*Z137SQ
		   +0.0083*pow(Z137,4)-0.002*pow(Z137,6))+ALOGZ3;

  CTIP1 = 4./15.*PI*Z137;				       // this block of
  CTIP2 = (23./9.-139./720.*PI*PI)*Z137SQ;		       // terms is used
  g     = sqrt(1.-Z137SQ);			               // in tip only
  f     = pow(gamf(g+1.)/(g*g*gamf(2.*g))*pow(2.*Z137,g-1.),2);// ref 1, eq 6
  CTIP3 = 4.*piz137*exp(-piz137)*f;                            // 
  CTIP4 = -2.*Z137*acos(Z137);				       // ref 1, eq 7
    
  deltol(zrad,DEL,TOL);
  
  ce = 3.4894e-4*zrad/arad;
  cn = ce*zrad;

  return trad*(cn*bethe(ebeam,egam)+ce*phie(ebeam,egam))/egam;
}

//
//************************************************************************
//                                                                       *
//     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.               *
//                                                                       *
//************************************************************************
//                                                                       *
//         note: the common /BETH/ must be loaded before using any of    *
//               these modules.                                          *
//                                                                       *
//************************************************************************
//
// here is only phin and tip
// deltol and bethe are in "brem.h"


double phin(double t0, double k, double kmatch, double phitip)
{
  if(t0 < k) return 0;
  if(k<kmatch) return bethe(t0,k);
  double b = bethe(t0,kmatch);
  return (phitip-b)/(t0-kmatch)*(k-kmatch)+b;
}


int tip(double t0, double& phitip, double& kmatch)
{
  extern double ZTHIRD,CTIP1,CTIP2,CTIP3,CTIP4,DEL,TOL;
  double E0,delt,beta,pbeta,rbeta,tang,slope;
  int    it;

  double mbeta,kbeta,khi,klo,k15;

  double ms_e = 0.5110034;

  E0    = t0+ms_e;
  delt  = E0/ms_e;
  beta  = sqrt(1.-pow(delt,-2));
  mbeta = 4./3.+(delt-2.)*(delt-1.)/(beta*beta*delt)*
    (1.+1./(2.*beta*delt*delt)*log((1.-beta)/(1.+beta)));
  kbeta = 1./pow(beta,3)*(
			  delt-17./2.+63./(4.*delt)
			  -25./(4.*delt*delt)-2./pow(delt,3)
			  -15./(8.*beta*pow(delt,3))*(delt-2.)*(delt-1.)
			  *log((1.-beta)/(1.+beta)));
  pbeta = pow(beta*delt,3)/pow(t0/ms_e,4)
    *exp(CTIP4*(1./beta-1.))*mbeta;
  rbeta = kbeta/mbeta;
  phitip= CTIP3*pbeta*(1.-CTIP1*rbeta+CTIP2);
  
  //	k15 is where the screening parameter gamma=15 at which point
  //	bh changes from screened to unscreened, possible discontinuity
  //	in bh at this point
  
  k15  = 15.*pow(t0+ms_e,2)*ZTHIRD/(100.*ms_e+15.*(t0+ms_e)*ZTHIRD);
  delt = DEL*t0/2.;
  khi  = t0;
  klo  = khi-10.;
  it   = 0;

  //*******end of initialization, now look for match

  while (it<100)
    {
      kmatch = (khi+klo)/2.;
      if( k15<(kmatch-delt)||k15>(kmatch+delt))
	{
	  tang  = 0.5*(bethe(t0,kmatch-delt)-bethe(t0,kmatch+delt))/delt;
	  slope = (bethe(t0,kmatch)-phitip)/(t0-kmatch);
	}
      else
	{  
	  tang  = 0.25*(bethe(t0,kmatch-3.*delt)-bethe(t0,kmatch-delt)+
			bethe(t0,kmatch+delt)-bethe(t0,kmatch+3.*delt))/delt;
	  slope = 0.5*((bethe(t0,kmatch-delt)-phitip)/(t0-(kmatch-delt))+
		       (bethe(t0,kmatch+delt)-phitip)/(t0-(kmatch+delt)));
	}
      
      if(fabs((tang-slope)/(0.5*(tang+slope)))<TOL) return 0;
      if(tang>slope) khi = kmatch;
      else klo = kmatch;
      it++;
    }

  printf(" **tip**   match not found after 100 iterations.\n");
  return 1;
}