1 gaskelld 1.1 	subroutine enerloss_new(len,dens,zeff,aeff,epart,mpart,typeflag,Eloss)
 2              
 3              	implicit none
 4              
 5              	real*8 thick,len,dens,zeff,aeff,epart,mpart,Eloss
 6              	real*8 x,chsi,lambda,gauss1,Eloss_mp,gamma
 7              	real*8 denscorr,CO,hnup,log10bg,I,beta,Eloss_mp_new
 8              
 9              	integer typeflag          !1=normal eloss (picked from distribution)
10                                                !2=min eloss
11              	                          !3=max eloss
12                                                !4=most probable eloss
13              	real*8 me
14                      parameter(me=0.51099906)
15              
16              	thick = len*dens
17              	gamma=epart/mpart
18                      beta = sqrt(1.-1./gamma**2)
19              
20              	if(zeff.eq.1) then	!Ionization potential in MeV
21              	  I = 21.8d-06
22 gaskelld 1.1 	else
23              	  I = (16.*zeff**0.9)*1.0d-06
24              	endif
25              
26              	hnup = 28.816d-06*sqrt(dens*zeff/aeff) !plasma frequency
27              	log10bg = log(beta*gamma)/log(10.)
28              	CO=log(hnup)-log(I)+0.5
29              
30              C DJG Get density effect correction (I got this from JV).
31              
32              	if(log10bg.lt.0.) then
33              	  denscorr=0.
34              	elseif(log10bg.lt.3.) then
35              	  denscorr=CO+log(10.)*log10bg+abs(CO/27.)*(3.-log10bg)**3
36              	elseif(log10bg.lt.4.7) then
37              	  denscorr=CO+log(10.)*log10bg
38              	else
39              	  denscorr=CO+log(10.)*4.7
40              	endif
41              	   
42              	if (thick.le.0.) then
43 gaskelld 1.1 	  Eloss = 0.
44              	else
45              	  Eloss_mp = 0.1536d-03 * zeff/aeff * thick * ( 19.26 +
46                   &          log(thick/dens) )
47              	  Eloss_mp_new = 0.1536d-03 * zeff/aeff *thick/beta**2* (
48                   &          log(me/I**2) + 1.063 + 2.*log(gamma*beta) + 
49                   &		log(0.1536*zeff/aeff*thick/beta**2)-beta**2-denscorr)
50              c	  write(6,*) 'ELOSS',Eloss_mp,Eloss_mp_new 
51              ! ........ convert to MeV, the unit of choice in THIS program
52              ! ........ (cf. EVCOIN where GeV prevail)
53              	  Eloss_mp = Eloss_mp_new*1000.
54              	  chsi = 0.307075/2.*zeff/aeff*thick/beta**2
55              	  if(typeflag.eq.1)then
56              	    x=abs(gauss1(10.0d0))
57              	  elseif(typeflag.eq.2)then
58              	    x=3
59              	  elseif(typeflag.eq.3)then
60              	    x=0.0067
61              	  elseif(typeflag.eq.4)then
62              	    x=1
63                        endif
64 gaskelld 1.1 	  if(x.gt.0.0) then
65              	    lambda = -2.0*log(x)
66              	  else
67              	    lambda = 100000.
68              	  endif
69              	  Eloss = lambda*chsi+eloss_mp
70              	endif
71              
72              	return
73              	end