PROGRAM finale
      IMPLICIT NONE

C*********************************************
C IMPORTANT REFERENCES
C F.A.Cucinotta et al., AIP Conference Proceedings 362 (1996) 245;
C C.Grimani et al., Class. Quantum Grav. 32 (2015) 035001;
C M.Armano, et al., Ap.J. 854 (2018) 113;
C*********************************************

c Position of ions and electrons
      REAL*8 xi,yi,zi,xe,ye,ze
c Initial ion production
      REAL*8 x0,y0,z0
c ion and electron energy
      REAL eni, en, en0
      REAL enel, enelnew
c electrons and ions angles
      REAL*8 thi,phi,the,phe,thenew,phen
c ion mass
      REAL M
C number of mass and atomic number of the ion 
      INTEGER A,z
      real cs,ss,sc
c parameters
      REAL*8 eta, random, dx, step, pi,sl,stepel,uno
c conversion joule/eV
      REAL cnv
c energy deposited by electrons in TM
      REAL energy
c TM mass (kg)
      REAL mass
c electrons and ion counters
      INTEGER ne,ntot,nu,ntota
      INTEGER i,n,nioni
c ion flux parameters
      REAL num,c1,c2,c3,minimo,massimo,fg
c various fixed parameters
      parameter (pi=3.14159)
      parameter (sl=29979200)
      parameter (uno=1.)
      parameter (mass=2.)
      parameter (cnv=0.00000000000000000016)
c Creates output file containing energy of incoming particle 
c and number of electrons produced
      OPEN(12,FILE='numero_elettroni.txt')
      CLOSE(12,STATUS='DELETE')
      OPEN(2,FILE='numero_elettroni.txt',FORM='FORMATTED',
     &POSITION="append")
c creates 
      OPEN(1,FILE='elettroni.txt',FORM='FORMATTED',POSITION="append")
      
      OPEN(3,FILE='elettroni_propagati.txt',FORM='FORMATTED',
     &POSITION='append')

c ELIMINATES existing output file containing measureted dose and 
c total number of electrons produced in the simulation and the number
c of electrons per second
      OPEN(8,FILE='risultati_is.txt')
      CLOSE(8,STATUS='DELETE')
c geometrical factor
      fg=0.047214
c Integration time (s)
      t=1 
c parameters for the ion flux (M.Armano et al.)
      num=18000.
      c1=1.19
      c2=3.66
      c3=0.87
c minimum of the ion flux
      minimo=6.698*10**-3*fg*t
c maximum of the ion flux
      massimo=1488.58*fg*t
c step simulation ions 1 nm
      step=0.001
c step simulation electrons 1 Angstrom
      stepel=0.0001
c initialization of electron counter
      ntot=0
c set ions parameter
      z=1
      A=1
c set ion mass (eV)
      M=9.38*10**8
      
      DO nioni=1, 100
c extracting ion energy
       en0=(0.1+(200.-0.1)*rand())
       func=num*(1./(en0+c1)**c2)*en0**c3*fg*t
       alpha=minimo+(massimo-minimo)*rand()
       DO WHILE(alpha.gt.func)
        en0=(0.1+(200.-0.1)*rand())
        func=num*(1./(en0+c1)**c2)*en0**c3*fg*t
        alpha=minimo+(massimo-minimo)*rand()
       END DO
       eni=en0*10**9
c extracting ion angles
       thi=pi*rand()
       phi=2*pi*rand()
c extracting ion position on the electrodes
       j=int(rand()*6)

       IF (j==0) THEN
        xi=26000.9
        yi=(2.*rand()-1.)*26000.9
        zi=(2.*rand()-1.)*26000.9
       ELSE IF (j==1) THEN
        xi=-26000.9
        yi=(2.*rand()-1.)*26000.9
        zi=(2.*rand()-1.)*26000.9
        ELSE IF (j==2) THEN
         xi=(2.*rand()-1.)*26000.9
         yi=26000.9
         zi=(2.*rand()-1.)*26000.9
         ELSE IF (j==3) THEN
          xi=(2.*rand()-1.)*26000.9
          yi=-26000.9
          zi=(2.*rand()-1.)*26000.9
          ELSE IF (j==4) THEN
           xi=(2.*rand()-1.)*26000.9
           yi=(2.*rand()-1.)*26000.9
           zi=26000.9
           ELSE IF (j==5) THEN
            xi=(2.*rand()-1.)*26000.9
            yi=(2.*rand()-1.)*26000.9
            zi=-26000.9
       END IF
       x0=xi
       y0=yi
       z0=zi
c contatore elettroni
       ne=0

c propagation of ions in the electrodes
       DO WHILE(
     &(26000.le.xi.and.xi.le.26000.9.and.
     &-26000.9.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.zi.and.zi.le.26000.9).or.
c
     &(-26000.9.le.xi.and.xi.le.-26000.and.
     &-26000.9.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.zi.and.zi.le.26000.9).or.
c
     &(26000.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.xi.and.xi.le.26000.9.and.
     &-26000.9.le.zi.and.zi.le.26000.9).or.
c
     &(-26000.9.le.yi.and.yi.le.-26000.and.
     &-26000.9.le.xi.and.xi.le.26000.9.and.
     &-26000.9.le.zi.and.zi.le.26000.9).or.
c
     &(26000.le.zi.and.zi.le.26000.9.and.
     &-26000.9.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.xi.and.xi.le.26000.9).or.
c
     &(-26000.9.le.zi.and.zi.le.-26000.and.
     &-26000.9.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.xi.and.xi.le.26000.9))
        eta=-log(rand())
        DO WHILE(
     &(26000.le.xi.and.xi.le.26000.9.and.
     &-26000.9.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.zi.and.zi.le.26000.9
     &.and.(eni.ge.0).and.((eta-step/lambda(eni,z,A,M)).ge.0)).or.
c
     &(-26000.9.le.xi.and.xi.le.-26000.and.
     &-26000.9.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.zi.and.zi.le.26000.9
     &.and.(eni.ge.0).and.((eta-step/lambda(eni,z,A,M)).ge.0)).or.
c
     &(26000.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.xi.and.xi.le.26000.9.and.
     &-26000.9.le.zi.and.zi.le.26000.9
     &.and.(eni.ge.0).and.((eta-step/lambda(eni,z,A,M)).ge.0)).or.
c
     &(-26000.9.le.yi.and.yi.le.-26000.and.
     &-26000.9.le.xi.and.xi.le.26000.9.and.
     &-26000.9.le.zi.and.zi.le.26000.9
     &.and.(eni.ge.0).and.((eta-step/lambda(eni,z,A,M)).ge.0)).or.
c
     &(26000.le.zi.and.zi.le.26000.9.and.
     &-26000.9.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.xi.and.xi.le.26000.9
     &.and.(eni.ge.0).and.((eta-step/lambda(eni,z,A,M)).ge.0)).or.
c
     &(-26000.9.le.zi.and.zi.le.-26000.and.
     &-26000.9.le.yi.and.yi.le.26000.9.and.
     &-26000.9.le.xi.and.xi.le.26000.9
     &.and.(eni.ge.0).and.((eta-step/lambda(eni,z,A,M)).ge.0)))
         eta=eta-step/lambda(eni,z,A,M)
c updates ion position
         IF(x0.ge.0) xi=xi-step*cos(phi)*sin(thi)
         IF(x0.le.0) xi=xi+step*cos(phi)*sin(thi)
         IF(y0.ge.0) yi=yi-step*sin(thi)*sin(phi)
         IF(y0.le.0) yi=yi+step*sin(phi)*sin(thi)
         IF(z0.ge.0) zi=zi-step*cos(thi)
         IF(z0.le.0) zi=zi+step*cos(thi)

       
         eni=eni-step*SPioni(eni,z,A,M)
        END DO
c an electron has been generated
        ne=ne+1
c extracting ion energy
        en=10.
        random=rand()
        n=1
        dx=1000.
        enel=0.
        DO WHILE((dx>0.0001) .and. (n<100) )
         enel=en-(cumulativa(z,A,M,en,eni)-random)/sigma(z,A,M,en,eni)
         dx=abs(enel-en)
         en=enel
         n=n+1
        END DO
c values for the electrons angles
        the=acos(sqrt(enel/we(A,M,eni)))-thi
        phe=acos(2.*rand()-1.)
c update ion energy
        eni=eni-enel
c write parameter of created electron on the file
        WRITE(1,*) xi,yi,zi,the,phe,enel,x0,y0,z0   
c sono uscito dagli elettrodi
 100    END DO 
        ntot=ntot+ne
        WRITE(*,*) 'Ion number ',nioni
        WRITE(*,*) 'Electrons produced',ne
        WRITE(2,*) en0, ne
c fine dell'iterzaione sul numero degli ioni
      END DO
c       WRITE(*,*) 'Numero totale elettroni', achar(9),ntot
       WRITE(*,*) 'Propagation of electrons'
      ne=0
      
      CLOSE(1)
      OPEN(1,FILE='elettroni.txt',STATUS='old')
c propagation of electrons in the electrodes
      iline=0
      DO    
       READ (1,*, END=10) xe,ye,ze,the,phe,enel,x0,y0,z0
       iline=iline+1
        DO WHILE(
     &(26000.le.xe.and.xe.le.26000.9.and.
     &-26000.9.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.ze.and.ze.le.26000.9.and.
     & enel.ge.0).or.
c
     &(-26000.9.le.xe.and.xe.le.-26000.and.
     &-26000.9.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.ze.and.ze.le.26000.9.and.
     & enel.ge.0).or.
c
     &(26000.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.xe.and.xe.le.26000.9.and.
     &-26000.9.le.ze.and.ze.le.26000.9.and.
     & enel.ge.0).or.
c
     &(-26000.9.le.ye.and.ye.le.-26000.and.
     &-26000.9.le.xe.and.xe.le.26000.9.and.
     &-26000.9.le.ze.and.ze.le.26000.9.and.
     & enel.ge.0).or.
c
     &(26000.le.ze.and.ze.le.26000.9.and.
     &-26000.9.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.xe.and.xe.le.26000.9.and.
     & enel.ge.0).or.
c
     &(-26000.9.le.ze.and.ze.le.-26000.and.
     &-26000.9.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.xe.and.xe.le.26000.9.and.
     & enel.ge.0))  
        eta=-log(rand())
c        WRITE(*,*) eta-stepel/lambdae(enel,z,A,M)
        DO WHILE(
     &(26000.le.xe.and.xe.le.26000.9.and.
     &-26000.9.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.ze.and.ze.le.26000.9.and.
     & enel.ge.0.and.eta-stepel/lambdae(enel,z,A,M).ge.0).or.
c
     &(-26000.9.le.xe.and.xe.le.-26000.and.
     &-26000.9.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.ze.and.ze.le.26000.9.and.
     & enel.ge.0.and.eta-stepel/lambdae(enel,z,A,M).ge.0).or.
c
     &(26000.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.xe.and.xe.le.26000.9.and.
     &-26000.9.le.ze.and.ze.le.26000.9.and.
     & enel.ge.0.and.eta-stepel/lambdae(enel,z,A,M).ge.0).or.
c
     &(-26000.9.le.ye.and.ye.le.-26000.and.
     &-26000.9.le.xe.and.xe.le.26000.9.and.
     &-26000.9.le.ze.and.ze.le.26000.9.and.
     & enel.ge.0.and.eta-stepel/lambdae(enel,z,A,M).ge.0).or.
c
     &(26000.le.ze.and.ze.le.26000.9.and.
     &-26000.9.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.xe.and.xe.le.26000.9.and.
     & enel.ge.0.and.eta-stepel/lambdae(enel,z,A,M).ge.0).or.
c
     &(-26000.9.le.ze.and.ze.le.-26000.and.
     &-26000.9.le.ye.and.ye.le.26000.9.and.
     &-26000.9.le.xe.and.xe.le.26000.9.and.
     & enel.ge.0.and.eta-stepel/lambdae(enel,z,A,M).ge.0)) 
         eta=eta-stepel/lambdae(enel,z,A,M)
c update electron energy
         enel=enel-stepel*SPE(enel)
c updating electron position
          IF(x0.ge.0) xe=xe-p0/dp*stepel*cos(phe)*sin(the)
          IF(x0.le.0) xe=xe+p0/dp*stepel*cos(phe)*sin(the)
          IF(y0.ge.0) ye=ye-p0/dp*stepel*sin(phe)*sin(the)
          If(y0.le.0) ye=ye+p0/dp*stepel*sin(phe)*sin(the)
          IF(z0.ge.0) ze=ze-p0/dp*cos(the)
          IF(z0.le.0) ze=ze+p0/dp*cos(the)
        END DO
        IF(wee(enel).ge.10) THEN
c a secondary electron has been created
         ne=ne+1
c extracting secondary electron energy
         en=10.
         random=rand()
         n=1
         dx=1000.
         enelnew=0.
         DO WHILE((dx>0.0001) .and. (n<100) )
          enelnew=en-(cumulativae(z,A,M,en,enel)-random)/
     &sigmae(z,A,M,en,enel)
          dx=abs(enelnew-en)
          en=enelnew
          n=n+1
         END DO
c values of secondary electrons angles
          thenew=acos(sqrt(enelnew/wee(enel)))-the
          phen=acos(2.*rand()-1.)
c update primary electron energy
          enel=enel-enelnew
c scrivo i parametri dell'elettrone creato alla fine del file
          call append(1,iline,xe,ye,ze,thenew,phen,enelnew,x0,y0,z0)
        END IF
       END DO
c writing on the file position and energy of the electron at the 
c end of its life (either ecause energy=0 or because it is out 
c of the electrodes)
       WRITE(3,*) xe,ye,ze,enel
      END DO 

c close and deleted useless files
  10   CLOSE(1,STATUS='DELETE')
  20   CLOSE(2)
       CLOSE(3)
       OPEN(4,FILE='elettroni_propagati.txt',STATUS='old')
c counting electrons       
       nu=0
       energy=0
       DO
        READ(4,*,END=30) xe,ye,ze,enel
        IF((enel.ge.0.and.
     &20000.le.xe.and.xe.le.26000.and.
     &-26000.le.ye.and.ye.le.26000.and.
     &-26000.le.ze.and.ze.le.26000).or.
c
     &(enel.ge.0.and.
     &-26000.le.xe.and.xe.le.-20000.and.
     &-26000.le.ye.and.ye.le.26000.and.
     &-26000.le.ze.and.ze.le.26000).or.
c
     &(enel.ge.0.and.
     &20000.le.ye.and.ye.le.26000.and.
     &-26000.le.xe.and.xe.le.26000.and.
     &-26000.le.ze.and.ze.le.26000).or.
c
     &(enel.ge.0.and.
     &-26000.le.ye.and.ye.le.-20000.and.
     &-26000.le.xe.and.xe.le.26000.and.
     &-26000.le.ze.and.ze.le.26000).or.
c
     &(enel.ge.0.and.
     &20000.le.ze.and.ze.le.26000.and.
     &-26000.le.xe.and.xe.le.26000.and.
     &-26000.le.ye.and.ye.le.26000).or.
c
     &(enel.ge.0.and.
     &-26000.le.ze.and.ze.le.-20000.and.
     &-26000.le.xe.and.xe.le.26000.and.
     &-26000.le.ye.and.ye.le.26000)) THEN 
        nu=nu+1
        energy=energy+enel 
       END IF
       END DO
c close and delete useless file     
  30   CLOSE(4,STATUS='DELETE')
c Some notification 
       WRITE(*,*) 'Number of electrons that exited the electrodes',
     &achar(9),nu
       WRITE(*,*) 'Number of electron per second',achar(9),nu/t
       WRITE(*,*) 'DOSE', energy*cnv/(mass*t)
c write results on the output file
       OPEN(6,FILE='risultati_is.txt',FORM='FORMATTED'
     &,POSITION='append')      
       WRITE(6,*) 'Total number of electrons',achar(9),achar(32), ntot
       WRITE(6,*) 'Number of electrons that exited the electrodes', 
     &achar(9),achar(32),nu
       WRITE(6,*) 'Number of electrodes per second ',nu/t
       WRITE(6,*) 'DOSE',achar(9),achar(9),achar(9),achar(9),achar(32),
     &achar(32),energy*cnv/(mass*t)
       CLOSE(6)
C**************************************
C end of program
C**************************************
C Necessary functions
c**************************************
       CONTAINS
c electrons stopping power

       FUNCTION SPE(en)
       REAL:: SPE,en
      
       IF(0.le.en.and.en.le.27.14) SPE=0.00053525*en**2.47*10**4
       IF(27.14.le.en.and.en.le.40.6) SPE=0.0015853*en**2.82*10**4
       IF(40.6.le.en.and.en.le.98.4) SPE=1.128733*en**0.4044*10**4
       IF(98.4.le.en.and.en.le.386.73) SPE=8.90786*en**0.364610**4
       IF(386.73.le.en.and.en.le.4691.68) SPE=78.7338*en**-0.38643*
     &10**4
       IF(4691.68.le.en.and.en.le.10**5) SPE=719.83*en**-0.65125*10**4
       
       RETURN
       END FUNCTION
c total ion stopping power
       FUNCTION sioni(en,A,M)
        REAL:: en,M,me,i,h,wp,rho, sioni
        INTEGER::Zg,Ag,A
        me=511000.
        Zg=79
        Ag=197
        rho=19.32
        i=790.
c        h=6.582*10**(-16)
        h=0.0000000000000006582
        wp=28.816*sqrt(rho*Zg/Ag)

        sioni=100*0.307075*rho*Zg/Ag*z**2./beta(A,M,en)**2*(0.5*
     &    log((2.*me*beta(A,M,en)**2*g(A,M,en)**2*Wmax(en,A,M))/i**2)-
     &  beta(A,M,en)**2-(log(h*wp/i)+log(beta(A,M,en)*g(A,M,en))-0.5))
       RETURN
       END FUNCTION
c energy deposited by ionization
       FUNCTION ecross(en,z,A,M)
        REAL::ecross,en,rho,M
        INTEGER::z,A
        rho=19.32

        ecross=rho*Zs(z,A,M,en)**2*(beta(A,M,en)*(-0.893177/
     &sqrt(we(A,M,en))*(-10./3.+we(A,M,en))*Zs(z,A,M,en)**2+
     &61.5813*beta(A,M,en))+we(A,M,en)*(6.15813
     &*log(we(A,M,en)/10.)+0.188286*Zs(z,A,M,en)**2*beta(A,M,en)
     &-6.15813*beta(A,M,en)**2))/(we(A,M,en)*beta(A,M,en)**2)
       RETURN
       END FUNCTION
c stopping power electrons without energy deposited by ionization
       FUNCTION SPioni(en,z,A,M)
        REAL::en,M,SPioni
        INTEGER::A,z

        SPioni=sioni(en,A,M)-ecross(en,z,A,M)
       RETURN
       END FUNCTION

c cross section of electron production by ion
       FUNCTION crossec(en,z,A,M)
        REAL::en, crossec,rho,M
        INTEGER::z,A
        rho=19.32
        crossec=0.0893117*rho*Zs(z,A,M,en)**2*(-68.951+6.8951*
     &     we(A,M,en)+(-6.3456+sqrt(we(A,M,en))+10./sqrt(we(A,M,en)))*
     &     Zs(z,A,M,en)**2*beta(A,M,en)-68.951*log(we(A,M,en)/10.)*
     &     beta(A,M,en)**2)/(we(A,M,en)*beta(A,M,en)**2)
       RETURN
       END FUNCTION

c cross section of electron production by electron
       FUNCTION crossece(en,z,A,M)
        REAL:: crossece,rho,M
        REAL::en
        INTEGER::z,A
        rho=19.32
        crossece=2*0.0893117*rho*Zse(en)**2*(-68.951+6.8951*
     &  wee(en)+(-6.3456+sqrt(wee(en))+10./sqrt(wee(en)))*
     &  Zse(en)**2*betae(en)-68.951*log(wee(en)/10.)*
     &  betae(en)**2)/(wee(en)*betae(en)**2)
       RETURN
       END FUNCTION

c electron range
       FUNCTION rangel(en)
        REAL::en, rangel

        IF(0.lt.en.and.en.lt.27.14) rangel=0.000503525*en**2.47
        IF(27.14.lt.en.and.en.lt.40.6) rangel=0.00158583*en**2.82
        IF(40.6.lt.en.and.en.lt.98.4) rangel=1.128733*en**0.4044
        IF(98.4.lt.en.and.en.lt.386.73) 
     &rangel=7.796946695*en**-0.00002646
        IF(386.73.lt.en.and.en.lt.4691.68) rangel=76.7338*en**0.38643
        IF(4691.68.lt.en.and.en.lt.10**5) rangel=719.83*en**-0.6512254
       RETURN
       END FUNCTION

c mean free path of ions
       FUNCTION lambda(en,z,A,M)
        REAL::M,lambda,en
        INTEGER::z,A
        REAL::aa
        aa=4.1
        lambda=aa/crossec(en,z,A,M)
       
       RETURN
       END FUNCTION

c mean free path of electrons
       FUNCTION lambdae(en,z,A,M)
        REAL::M,lambdae
        INTEGER::z,A
        REAL::en
        lambdae=1./crossece(en,z,A,M)
       
       RETURN
       END FUNCTION

c velocity/c of ion
       FUNCTION beta(A,M,T)
        INTEGER::A
        REAL::M,T,beta
        beta=sqrt(1.-(1./((A*T/M)+1.))**2)
       RETURN
	END

c velocity/c electrons
       FUNCTION betae(T)
        REAL::M,T
        REAL*8::betae
        M=511000.
        betae=sqrt(1.-(1./((1.*T/M)+uno))**2)
       RETURN
	END

c screened ion charge (see Cucinotta et al.)
       FUNCTION Zs(z,A,M,T)
        INTEGER::z,A
        REAL::M,T,Zs

        Zs=z*(1.-exp(-125.*beta(A,M,T)/z**(2./3.)))
       RETURN
       END

c screened electron charge (see Cucinotta et al.)
       FUNCTION Zse(T)
        INTEGER::z,A
        REAL::M,Zse
        REAL::T

        Zse=-1*(1.-exp(-125.*betae(T)/z**(2./3.)))
       RETURN
       END

c maximum energy transfer for primary electrons
c (Cucinotta et al.)
       FUNCTION we(A,M,T)
        INTEGER::A
        REAL::M,T,me,we
        me=511000.	
	
        we=(2*me*beta(A,M,T)**2)/(1-beta(A,M,T)**2)
        RETURN
	END

c maximum energy transfer for secondary electrons
c (Cucinotta et al.)
       FUNCTION wee(en)
        REAL::me,wee,en
        me=511000.
       
        wee=0.5*(2*(A*en/Me)*((A*en/Me)+2.))*me/
     &(1+2*((A*en/Me)+1.)*(me/Me)+(me/Me)**2)
        RETURN
        END FUNCTION

c gamma factor      
       FUNCTION g(A,M,T)
        REAL::M,T,g
        INTEGER::A
        g=1./sqrt(1.-beta(A,M,T)**2)
        
       END FUNCTION

c maximum energy transfer       
       FUNCTION Wmax(en,A,M)
        REAL::en,M,me,Wmax
        INTEGER::A
        me=511000.
       
        Wmax=(2*(A*en/M)*((A*en/M)+2.))*me/
     &(1+2*((A*en/M)+1.)*(me/M)+(me/M)**2)
        
       RETURN
       END FUNCTION

c cumulative of cross section for primary electrons production
       FUNCTION cumulativa(z,A,M,T,en)
        INTEGER::z,A
        REAL::M,T,en,cumulativa

        cumulativa=(6.8951*(T*beta(A,M,en)*((0.14503*
     &      sqrt(we(A,M,en))+1.4503/sqrt(we(A,M,en))-0.458627*
     &      sqrt(T/we(A,M,en)))*Zs(z,A,M,en)**2
     &      -10.*log(T/10.)*beta(A,M,en))+we(A,M,en)*(-10+T-
     &      0.458627*sqrt(T/we(A,M,en))*Zs(z,A,M,en)**2*beta(A,M,en))
     &      ))/T/(-68.951+6.8951*we(A,M,en)+(-6.32456+sqrt(we(A,M,en))
     &      +10./sqrt(we(A,M,en)))*Zs(z,A,M,T)**2*beta(A,M,en)-
     &      68.951*log(we(A,M,en)/10.)*beta(A,M,en)**2)
       RETURN
       END 

c cumulative of cross section for secondary electrons production
       FUNCTION cumulativae(z,A,M,T,en)
        INTEGER::z,A
        REAL::M,cumulativae
        REAL::en,T

        cumulativae=(6.8951*(T*betae(en)*((0.14503*
     &    sqrt(wee(en))+1.4503/sqrt(wee(en))-0.458627*
     &    sqrt(T/wee(en)))*Zse(en)**2
     &    -10.*log(T/10.)*betae(en))+wee(en)*(-10+T-
     &    0.458627*sqrt(T/wee(en))*Zse(en)**2*betae(en))
     &    ))/T/(-68.951+6.8951*wee(en)+(-6.32456+sqrt(wee(en))
     &    +10./sqrt(wee(en)))*Zse(T)**2*betae(en)-
     &    68.951*log(wee(en)/10.)*betae(en)**2)
       RETURN
       END 

c cross section for primary electrons production
       FUNCTION sigma(z,A,M,T,en)
        REAL::M,T,en,sigma
        INTEGER::A,z

		sigma=(68.951*we(A,M,en)*(1.+3.14*(T/we(A,M,en))**1.5*
     &      (-T+we(A,M,en))*Zs(z,A,M,en)**2*beta(A,M,en)/(137*T))-
     &      T*beta(A,M,en)**2/we(A,M,en))/T**2/(-68.951+6.8951*
     &      we(A,M,en)+(-6.32456+sqrt(we(A,M,en))+10./
     &      sqrt(we(A,M,en)))*Zs(z,A,M,en)**2*
     &      beta(A,M,en)-68.951*log(we(A,M,en)/10.)*beta(A,M,en)**2)
       RETURN
       END

c cross section for secondary electrons production
       FUNCTION sigmae(z,A,M,T,en)
        REAL::M,sigmae
        INTEGER::A,z
        REAL::T,en

		sigmae=(68.951*wee(en)*(1.+3.14*(T/wee(en))**1.5*
     &      (-T+wee(en))*Zse(en)**2*betae(en)/(137*T))-
     &      T*betae(en)**2/wee(en))/T**2/(-68.951+6.8951*
     &      wee(en)+(-6.32456+sqrt(wee(en))+10./
     &      sqrt(wee(en)))*Zse(en)**2*
     &      betae(en)-68.951*log(wee(en)/10.)*betae(en)**2)
       RETURN
       END
      
c definition of ion flux (Grimani et al. 2015)
       FUNCTION flusso(num,c1,c2,c3,en)
        REAL::num,c1,c2,c3,flusso,en
        flusso=num*(1./(en+c1)**c2)*en**c3
        RETURN
       END

c function needed to write on the files

      SUBROUTINE append(iunit,iline,xe,ye,ze,the,phe,enel,x0,y0,z0)
      IMPLICIT NONE
      INTEGER, intent(in)::iunit, iline
      REAL, intent(in)::enel
      REAL*8::xe,ye,ze,the,phe,x0,y0,z0
      INTEGER::i
c salta alla fine del file
      DO 
       READ(iunit,*,END=20)
      END DO
 20   CONTINUE
c salta indietro di una riga per evitare l'EOF
      BACKSPACE(iunit)
c appendere al file
      WRITE(iunit,*) xe,ye,ze,the,phe,enel,x0,y0,z0
c ritorno all riga in cui ero
      REWIND(iunit)
      i=0
      DO
       i=i+1
       READ(iunit,*)
       IF (i==iline) exit
      END DO
      END 

      END