C
C  Function PHOTON - Calculate bremsstrahlung yield per el,MeV
C                    for a given beam
C                    energy and photon energy
C

       DOUBLE PRECISION function photon(ebeam,egam,trad,zrad,arad)
       implicit none
C
       DOUBLE PRECISION trad,zrad,arad
       DOUBLE PRECISION egam,ebeam,CE,CN
C
C  Functions
C
       DOUBLE PRECISION bethe,phie
C
       if(egam.gt.ebeam) then
             write(6,*) 'You sould not see this message'
             photon = 0.0
             return
       endif
C       
        call phinit(zrad)

        CE=3.4894E-4*zrad/arad
	CN=CE*zrad
C
	photon = trad*(CN*BETHE(ebeam,egam)+CE*PHIE(ebeam,egam))/egam
C
        return
        end
C
C************************************************************************
C                                                                       *
C     This is PHIN, the nuclear bremsstrahlung package.  It contains    *
C     four modules:                                                     *
C                       1) DELTOL                                       *
C                       2) TIP                                          *
C                       3) BETHE                                        *
C                       4) PHIN                                         *
C                                                                       *
C       DELTOL(Z,DEL,TOL):                                              *
C               Subroutine. Provides values for DEL and TOL for         *
C               use by TIP.                                             *
C                                                                       *
C       TIP(T0,DEL,TOL,PHITIP,KMATCH)                                   *
C               Subroutine. Calculates the endpoint yeild (PHITIP) and  *
C               the photon energy (KMATCH) at which a straight line is  *
C               joined to the BETHE distribution. T0 is the incident    *
C               electron kinetic energy. DEL and TOL are matching       *
C               criteria.                                               *
C                                                                       *
C       BETHE(T0,K):                                                    *
C               Subroutine.  Calculates the BETHE distribution.         *
C                                                                       *
C       PHIN(T0,K,KMATCH,PHITIP)                                        *
C               Function routine. Calculates the nuclear bremsstrahlung *
C               yield at photon energy=K from electrons of kinetic      *
C               energy=T0.  If K .GT. KMATCH, it interpolates linearly  *
C               between (PHITIP,T0) and (BETHE(T0,KMATCH),KMATCH).      *
C               For K .LT. KMATCH, it simply calls BETHE.               *
C                                                                       *
C************************************************************************
C                                                                       *
C         note: the common /BETH/ must be loaded before using any of    *
C               these modules.                                          *
C                                                                       *
C************************************************************************
C
        SUBROUTINE PHINIT(Z)
        IMPLICIT NONE

        DOUBLE PRECISION Z
        DOUBLE PRECISION FCOUL,ZTHIRD,Z137,Z137SQ,Z23,ALOGZ23,ALOGZ3
        DOUBLE PRECISION CTIP1,CTIP2,CTIP3,CTIP4,DEL,TOL
        DOUBLE PRECISION X,PI,PIZ137,G,F
        DOUBLE PRECISION GAMF

        COMMON/BETH/    FCOUL,ZTHIRD,Z137,Z137SQ,
     *                  Z23,ALOGZ23
        COMMON/TIPCB/CTIP1,CTIP2,CTIP3,CTIP4,DEL,TOL

        GAMF(X)=1.-0.5748646*(X-1)+0.9512363*(X-1)**2   ! GAMMA FUNCTION, 1<X<2
     *          -0.6998588*(X-1)**3+0.4245549*(X-1)**4  ! USING POLYNOMIAL 
     *		-0.1010678*(X-1)**5 			! NBS MATHEMATICAL HANDBOOK) 

	PI = 3.14159 					! GENERATE A BUNCH OF
	Z137=Z/137.					! TERMS THAT ARE A
	PIZ137=PI*Z137					! FUNCTION OF Z ONLY
	Z137SQ=Z137*Z137
	ALOGZ3=1./3.*LOG(Z)
	ALOGZ23=ALOGZ3*2.
	ZTHIRD=Z**(1./3.)
	Z23=Z**(2./3.)

	FCOUL=Z137SQ*(1./(1.+Z137SQ)+0.20206-0.0369*Z137SQ	! REF 1, EQ 4
     *		+0.0083*Z137**4-0.002*Z137**6)+ALOGZ3		! USED IN BETHE

	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 =(GAMF(G+1.)/(G*G*GAMF(2.*G))*(2.*Z137)**(G-1.))**2	! REF 1, EQ 6
	CTIP3=  4.*PIZ137*EXP(-PIZ137)*F
	CTIP4=  -2.*Z137*ACOS(Z137)				! REF 1, EQ 7

	CALL DELTOL(Z,DEL,TOL)
	RETURN
	END


	FUNCTION PHIN(T0,K,KMATCH,PHITIP)
        IMPLICIT NONE
        DOUBLE PRECISION T0,KMATCH,PHITIP
	DOUBLE PRECISION K,B,PHIN,BETHE
	IF(T0.LT.K)GOTO 3
	IF(K.LT.KMATCH)GOTO 1
	B=BETHE(T0,KMATCH)
	PHIN=(PHITIP-B)/(T0-KMATCH)*(K-KMATCH)+B
	RETURN
   1	PHIN=BETHE(T0,K)
	RETURN
   3	PHIN=0.
	RETURN
	END


	SUBROUTINE TIP(T0,PHITIP,KMATCH)
C	COMPUTE TIP REGION OF BREMSSTRAHLUNG SPECTRUM USING DMH ENDPOINT WITH
C	JABBUR-PRATT-GAVRILA ELECTRON ENERGY DEPENDENCE
C	FIND POINT ON BETHE-HEITLER SPECTRUM WHERE TANGENT EQUALS SLOPE OF
C	STRAIGHT LINE TO ENDPOINT
C	INPUTS:Z,T0,DELL,TOL     OUTPUTS:PHITIP,KMATCH
        IMPLICIT NONE
        DOUBLE PRECISION T0,PHITIP
        DOUBLE PRECISION FCOUL,ZTHIRD,Z137,Z137SQ,Z23,ALOGZ23
        DOUBLE PRECISION CTIP1,CTIP2,CTIP3,CTIP4,DELL,TOL
        DOUBLE PRECISION X,Y,BH,TANG1,TANG2,SLOP1,SLOP2
        DOUBLE PRECISION BETHE
        DOUBLE PRECISION E0,DELT,BETA,PBETA,RBETA,TANG,SLOPE  
        INTEGER IT

	COMMON/BETH/	FCOUL,ZTHIRD,Z137,Z137SQ,
     *			Z23,ALOGZ23
	COMMON/TIPCB/CTIP1,CTIP2,CTIP3,CTIP4,DELL,TOL
	DOUBLE PRECISION MBETA,KBETA,KHI,KLO,KMATCH,K15
C*************************************************************
	BH(X)=BETHE(T0,X)
	TANG1(X,Y)=0.25*(BH(X-3.*Y)-BH(X-   Y)+
     *		         BH(X+   Y)-BH(X+3.*Y))/Y
	TANG2(X,Y)=0.5* (BH(X-Y)-BH(X+Y))/Y
	SLOP1(X,Y)=0.5*((BH(X-Y)-PHITIP)/(T0-(X-Y))+
     *			(BH(X+Y)-PHITIP)/(T0-(X+Y)))
	SLOP2(X  )=(BH(X)-PHITIP)/(T0-X)
C*************************************************************
	E0=T0+0.511

	DELT=E0/0.511
	BETA=SQRT(1.-DELT**(-2))

	MBETA=4./3.+(DELT-2.)*(DELT-1.)/(BETA**2*DELT)*		! REF 1, EQ 9 
     *		(1.+1./(2.*BETA*DELT**2)*LOG((1.-BETA)/(1.+BETA)))

	KBETA=1./BETA**3*(					! REF 1, EQ 10
     *			 DELT-17./2.+63./(4.*DELT)
     *			-25./(4.*DELT**2)-2./DELT**3
     *			-15./(8.*BETA*DELT**3)*(DELT-2.)*(DELT-1.)
     *			*LOG((1.-BETA)/(1.+BETA)))

	PBETA=(BETA*DELT)**3/((T0/0.511)**4)
     *		*EXP(CTIP4*(1./BETA-1.))*MBETA			! REF 1, EQ 7

	RBETA=KBETA/MBETA					! REF 1, EQ 8

	PHITIP=CTIP3*PBETA*(1.-CTIP1*RBETA+CTIP2)		! REF 1, EQ 6

C	K15 IS WHERE THE SCREENING PARAMETER GAMMA=15 AT WHICH POINT
C	BH CHANGES FROM SCREENED TO UNSCREENED, POSSIBLE DISCONTINUITY
C	IN BH AT THIS POINT

	K15=15.*(T0+0.511)**2*ZTHIRD/(51.1+15.*(T0+0.511)*ZTHIRD)
	DELT=DELL*T0/2.
	KHI=T0
	KLO=KHI-10.
	IT=0
C*******END OF INITIALIZATION, NOW LOOK FOR MATCH
   5	IT=IT+1
	KMATCH=(KHI+KLO)/2.
	IF(IT.GT.100)GOTO 20
	IF((K15.LT. (KMATCH-DELT)).OR.(K15.GT. (KMATCH+DELT)))THEN
		TANG =TANG2(KMATCH,DELT)
		SLOPE=SLOP2(KMATCH     )
	ELSE
		TANG =TANG1(KMATCH,DELT)
		SLOPE=SLOP1(KMATCH,DELT)
	END IF
	IF(ABS((TANG-SLOPE)/(0.5*(TANG+SLOPE))).LT.TOL)RETURN
	IF(TANG.GT.SLOPE)THEN
	KHI=KMATCH
	ELSE
	KLO=KMATCH
	END IF
	GOTO 5
  20	STOP ' **TIP**   MATCH NOT FOUND AFTER 100 ITERATIONS'
	END





	DOUBLE PRECISION FUNCTION BETHE(T0,K)
        IMPLICIT NONE
        DOUBLE PRECISION FCOUL,ZTHIRD,Z137,Z137SQ,Z23,ALOGZ23
        DOUBLE PRECISION E0,T0,E,GAMMA,PHI1,PHI2,R,P0,P,EPS0,EPS,EE,PP
	COMMON/BETH/	FCOUL,ZTHIRD,Z137,Z137SQ,
     *			Z23,ALOGZ23
	DOUBLE PRECISION K,L,temp,temp0
	IF((T0-.00002)-K)1,1,2 !.00002 is a fudge to keep routine from crashing
   1	BETHE=0.
	RETURN
   2	E0=T0+0.511
	E=E0-K
	GAMMA=51.1*K/(E0*E*ZTHIRD)				! REF 1, AFTER
	IF(GAMMA.GT.15.)GOTO 4					! EQ 2

C*SCREENED
	PHI1=19.24-4.*LOG(GAMMA+2./(GAMMA+3.))-		! REF 1, EQ 3
     *		0.12*GAMMA*EXP(-GAMMA/3.) 
	PHI2 = PHI1
	IF (GAMMA.LE.0.80)PHI2=PHI2-0.027-(0.8-GAMMA)**2 

	R=E/E0							! REF 1, EQ 2
	BETHE=4.*(	 (1.+R*R)*(1./4.*PHI1-FCOUL)
     *			-2./3.*R* (1./4.*PHI2-FCOUL))
	RETURN

C*NON-SCREENED							! REF 1, EQ 5
   4	E0=E0/0.511						! UNITS OF ELECTRON MASS
	K =K /0.511
	E= E0-K
	P0=SQRT(E0*E0-1.)
	P=SQRT(E*E-1.)
C
        temp0 = E0 - P0
        if(temp0.eq.0.0) then
          temp0 = 1.0/(2.0*E0) + 1.0/(8.0*E0*E0*E0)
        endif
        temp = E - P
        if(temp.eq.0.0) then
          temp = 1.0/(2.0*E) + 1.0/(8.0*E*E*E)
        endif
C
	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
	BETHE=P/P0*(	 4./3.-2.*EE*(P**2+P0**2)/PP**2
     *			+EPS0*E/P0**3+EPS*E0/P**3-EPS*EPS0/PP
     *			+L*(	 K/(2.*PP)*(
     *						 (EE+P0*P0)/P0**3*EPS0
     *						-(EE+P *P )/P **3*EPS
     *						+2.*K*EE/PP**2	)
     *				+8.*EE/(3.*PP)+	K*K*(EE**2+PP**2)/PP**3	))
	K=K*.511
	RETURN
	END


      SUBROUTINE DELTOL(Z, DEL, TOL)
      IMPLICIT NONE
      DOUBLE PRECISION Z, DEL, TOL,IZ
C      INTEGER IZ
c      IZ = IFIX(Z)
      IZ=Z
      TOL = 0.01
      IF (IZ.LE.4) GO TO 1
      IF (IZ.LE.6) GO TO 2
      IF (IZ.LE.13) GO TO 3
      IF (IZ.LE.29) GO TO 4
      IF (IZ.LT.73) GO TO 5
      TOL = 0.0075
      DEL = 0.003
      RETURN
 1    DEL = 0.00001
      TOL=.02
      RETURN
 2    DEL = 0.0002
      TOL=.02
      RETURN
 3    DEL = 0.0005
      TOL=.02
      RETURN
! 4    DEL = 0.00075 	! USED FOR 365A25.R1
 4    DEL = 0.0001	! Used for NE8 (ERK 19DEC87)
! 4    DEL = 0.001
      RETURN
! 5    DEL = 0.00075	! USED FOR 365A42.MA
! 5    DEL = 0.001	! USED FOR 365A42.R1,
			!  315A42.MA, AND 265A42.MA
 5    DEL = 0.002
      RETURN
      END
C************************************************************************
C									*
C	This is PHIE, the electron bremsstrahlung package.  It contains	*
C	one modules:							*
C			1) PHIE						*
C									*
C	PHIN(T0,K)							*
C		Function routine. Calculates the electron bremsstrahlung*
C		yeild at photon energy=K from electrons of kinetic	*
C		energy=T0.						*
C									*
C************************************************************************
C									*
C	note:	the common /BETH/ must be loaded before using PHIE	*
C									*
C************************************************************************
C
	DOUBLE PRECISION FUNCTION PHIE(T0,K)
        IMPLICIT NONE
CD	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DOUBLE PRECISION FCOUL,ZTHIRD,Z137,Z137SQ,Z23,ALOGZ23
        DOUBLE PRECISION E0,KX,E,P0,KMAX,RHO,SQETA,ETA,L,X
        DOUBLE PRECISION T0,EPSI,X2,X3,X4,X5,PSI
	COMMON/BETH/	FCOUL,ZTHIRD,Z137,Z137SQ,
     *			Z23,ALOGZ23
	DOUBLE PRECISION K
CD	REAL*8*8 K,L,KMAX,KX
	E0=(T0+0.51100340)/0.51100340		! UNITS OF ELECTRON MASS
	KX= K/0.51100340
	E0=(T0+0.5110)/0.5110		! UNITS OF ELECTRON MASS
	KX= K/0.5110
	E = E0-KX
	P0=SQRT(E0**2-1.00)
	KMAX=(E0-1.00)/(E0+1.00-P0)		! REF 1, AFTER EQ 11
	IF(KX.GE.KMAX)GOTO 15
	RHO=E0-KX*(1.00+E0-P0)-1.00
	ETA=RHO/(RHO+2.00)
	SQETA=SQRT(ETA)
	L=2.00*LOG(	(SQRT(E0-1.00)+SQRT(ETA*(E0+1.00)))/
     *			(SQRT(E0-1.00)-SQRT(ETA*(E0+1.00))))
	X=E/E0
	PHIE=	2.00*(1.00-2.00/3.00*X+X**2)*(L-SQETA)+	! REF 1, EQ 17
     *		SQETA*(1.00-L**2/(2.00*RHO)-1.00/RHO**2*
     *		(L/2.00-SQRT(RHO*(RHO+2.00)*(E0+1.00)
     *          /(E0-1.00)))**2)
C*SCREENED
	EPSI = 100.00*KX/(E0*E*Z23)
	IF (EPSI.GE.0.880)GOTO 10
	X = 0.880 - EPSI
	X2 = X*X
	X3 = X2*X
	X4 = X3*X
	X5 = X4*X
	PSI =19.700+4.11700*X-3.80560*X2+31.8350*X3 ! REF 1, EQ 15
     *		-58.6300*X4+40.7740*X5
	PHIE=PHIE*(0.250*PSI-1.0-ALOGZ23)/(3.7980-LOG(EPSI)-ALOGZ23) 
  10	RETURN
  15	PHIE=0.00
	RETURN
	END