C An interface between Atima and python
C For details related to Atima contact: H.Weick@gsi.de tel. 2158
C For questions to the wrapper: P.Boutachkov@gsi.de  tel. 2346
c
c Compile, run: make clean; make
c
c Usage:
c in python
c > import atima
c > atima.eloss(...)
c
c example 238U into 1 mg/cm2 27Al at 100 MeV/u 
c > atima.eloss([238, 92], [[27,13,1]], 2.7, 0, 1, 100)
c
c example 197Au into 10 mg/cm2 27Al at 1000 MeV/u 
c > atima.eloss([ 196.96655, 79], [[27,13,1]], 2.7, 0, 10, 1000)
c
c Read the comments below for comments to the the in and out variables

        SUBROUTINE eloss(proj, pn, ebeam,
     &                targ, nnuca, tm, rhoa, igasa, thick,  
     &                ein, eout, range, dedxi,dedxo,remra,rstr,estr,
     &                astr, tof, intrthick )
       IMPLICIT NONE
c An array defining [A, Z] of the projectile, here pn should be 2 
       INTEGER pn
       REAL*8  proj(2)
Cf2py intent(in) proj
Cf2py integer intent(hide),depend(proj) :: pn=shape(proj,0)      

c Beam Energy. If positive this is the energy IN the target. If negative, this denotes the energy after the target [MeV/u]      
       REAL*8 ebeam
Cf2py intent(ebeam) ebeam

c An array defining [[A0, Z0, SR0], ... ] the target composition, here SRi
c is the stochastic ratio of Ai, Zi 
c here nnuca is number of target compounds and tm = 3
       INTEGER nnuca, tm
       REAL*8  targ(nnuca, tm)
Cf2py intent(in) targ
Cf2py integer intent(hide),depend(targ) :: nnuca=shape(targ,0), tm=shape(targ,1)       

c Material density in g/cm3
       REAL*8  rhoa
Cf2py intent(in) rhoa

c 1 for gas 0 for solid
       INTEGER  igasa      
Cf2py intent(in) igasa

c target thickness in [mg/cm2]
       REAL*8  thick      
Cf2py intent(in) thick

c Energy IN [MeV/u]
       REAL*8 ein
Cf2py intent(out) ein

c Energy OUT [MeV/u]
       REAL*8 eout
Cf2py intent(out) eout

c Range in mg/cm2
       REAL*8  range
Cf2py intent(out) range

c dE/dX in the first 1 mg/cm2 of material [MeV/u/mg/cm2]
       REAL*8 dedxi 
Cf2py intent(out) dedxi

c dE/dX in the last 1 mg/cm2 of material [MeV/u/mg/cm2]       
       REAL*8 dedxo
Cf2py intent(out) dedxo

c Remaining Range in [mg/cm2]              
       REAL*8 remra
Cf2py intent(out) remra

c Range straggling [mg/cm2]
       REAL*8 rstr
Cf2py intent(out) rstr

c Energy straggling [mg/cm2]
c Multiply by dedxo to get the straggling in MeV/u
       REAL*8 estr
Cf2py intent(out) estr

c Angular straggling in [rad]      
       REAL*8 astr
Cf2py intent(out) astr

c Time Of Flight through the target in [ns]       
       REAL*8 tof
Cf2py intent(out) tof

c Returns the interpolated target thickness in [mg/cm2]
c If the ions stops in the material thick = range       
       REAL*8 intrthick
Cf2py intent(out) intrthick       
       
ccccccccc
c These are the corresponding variables in the world of ATIMA

c See the assignment bellow for the relation to the input       
       INTEGER  zpa
       REAL*8  ma(nnuca), IValue(nnuca)
       INTEGER  za(nnuca), stnr(nnuca)
       REAL*8  fntpa, GIValue, mpa
       
c This maximum dimension comes from  atima.1.32/src/dedx/calcint.f 
       INTEGER fspc
       PARAMETER (fspc=60000)

c These are splines for, range, range straggling and angular straggling
c here t is the knot sequence and corresponding b-coefficient       
c the corresponding n is the length of the  b-coefficient array and
c k is the spline order

       INTEGER  n1,k1,n2,k2,n3,k3,n4,k4

       REAL*8  t1(fspc)
       REAL*8  bc1(fspc)
       REAL*8  t2(fspc)
       REAL*8  bc2(fspc)
       REAL*8  t3(fspc)
       REAL*8  bc3(fspc)
       REAL*8  t4(fspc)
       REAL*8  bc4(fspc)


c Work variables, FG is used to pass the spline dimension info
       INTEGER  FG(8)       
       
       INTEGER j
       
c File name for the spline and the system time
c used to make a unique file name. Here chzp will the be zp as a character
c The function uniquefile checks if the name is indeed unique id not it adds 0
c to the name
       CHARACTER*35 filename, uniquefile
       CHARACTER*10 time
       
c A function which returns the Ionization potential      
       REAL*8 ipot
       
c Check if the input is ok
       IF (pn .NE. 2) THEN
         PRINT *,
     &   'Wrong Projectile list format  it must be [Ap, Zp, SR0]'
         RETURN     
       ENDIF

       IF (tm .NE. 3) THEN
         PRINT *,
     &  'Wrong Target list format  it must be [[A0, Z0, SR0], ... ] '     
         RETURN 
       ENDIF


c Transfer the parameters to ATIMA notation
       fntpa = 1.0D0
      
c This a global ionization potential
       GIValue = 0.0d0
      
       zpa = NINT(proj(2))
       mpa = proj(1)
      
       DO 101 j = 1, nnuca
         za(j) = NINT(targ(j, 2))
         ma(j) = targ(j, 1)
         IValue(j) = ipot(za(j))
         stnr(j) = NINT(targ(j, 3))

 101   CONTINUE 

c First generate the splines, ranges, angular strangling , ... 
c for the target projectile combination.
c The splines will be written in filename

c Get the current time, to be used to make a unique file name for the spline
       CALL date_and_time(TIME=time)
c Open a file in the shared memory
       filename = '/dev/shm/atima.'//time
       
c Check if the file name is unique, if not add 0 to it
       filename = uniquefile(filename)

       CALL atima( filename, 
     &                    nnuca, igasa, fntpa, GIValue, rhoa, zpa, mpa, 
     &                    za, ma, IValue, stnr, FG )

c Now read the spline data
       n1 = FG(1)
       k1 = FG(2)
       
       n2 = FG(3)
       k2 = FG(4)
       
       n3 = FG(5)
       k3 = FG(6)
 
       n4 = FG(7)
       k4 = FG(8)
 
       CALL SplinesLesen( filename,
     &                           n1, k1, t1, bc1, n2, k2, t2, bc2, 
     &                           n3, k3, t3, bc3, n4, k4, t4, bc4 )

c Clean the spline data as we have it in memory
       CALL delfile(filename)

c If the Ebeam is equal smaller than zero assume this is the exit energy from the target.
c First calculate Ein, if unknown
      IF (ebeam .LE. 0.) THEN
        eout = -ebeam
        CALL Sp0GroessenNeu( n1, k1, t1, bc1, eout, thick, ein)
      ELSE
        ein = ebeam
      ENDIF

c Now new can calculate what ever we want
       eout = -1.0D0
       CALL SpGroessenNeu( n1, k1, t1, bc1, n2, k2, t2, bc2, 
     &                            n3, k3, t3, bc3, n4, k4, t4, bc4,
     &                            ein, thick, eout,
     &                            range,dedxi,dedxo,remra,rstr,estr,
     &                            astr,tof )
      
       intrthick =  thick
              
       END SUBROUTINE

c Function returns the iIonization potentials for Z atomic number 
       REAL*8 FUNCTION ipot( nz )
       IMPLICIT NONE
       
       INTEGER nz
c Ionization potentials, from NIST page (ICRU Report#49, 1994) + extrapol.
       REAL*8  i_bar(120)
      
       IF (( nz .LT. 1) .OR. ( nz .GT. 120)) THEN
         PRINT *, 'Invalid target Z'
         ipot = 0.0
         RETURN
       ENDIF 
      
       i_bar(1)=19.2
       i_bar(2)=41.8
       i_bar(3)=40.0
       i_bar(4)=63.7 
       i_bar(5)=76.0
       i_bar(6)=78.0
       i_bar(7)=82.0
       i_bar(8)=95.0
       i_bar(9)=115.0
       i_bar(10)=137.0
       i_bar(11)=149.0
       i_bar(12)=156.0
       i_bar(13)=166.0
       i_bar(14)=173.0
       i_bar(15)=173.0
       i_bar(16)=180.0
       i_bar(17)=174.0
       i_bar(18)=188.0
       i_bar(19)=190.0
       i_bar(20)=191.0
       i_bar(21)=216.0
       i_bar(22)=233.0
       i_bar(23)=245.0
       i_bar(24)=257.0
       i_bar(25)=272.0
       i_bar(26)=286.0
       i_bar(27)=297.0
       i_bar(28)=311.0
       i_bar(29)=322.0
       i_bar(30)=330.0
       i_bar(31)=334.0
       i_bar(32)=350.0
       i_bar(33)=347.0
       i_bar(34)=348.0
       i_bar(35)=343.0
       i_bar(36)=352.0
       i_bar(37)=363.0
       i_bar(38)=366.0
       i_bar(39)=379.0
       i_bar(40)=393.0
       i_bar(41)=417.0
       i_bar(42)=424.0
       i_bar(43)=428.0
       i_bar(44)=441.0
       i_bar(45)=449.0
       i_bar(46)=470.0
       i_bar(47)=470.0
       i_bar(48)=469.0
       i_bar(49)=488.0
       i_bar(50)=488.0
       i_bar(51)=487.0
       i_bar(52)=485.0
       i_bar(53)=491.0
       i_bar(54)=482.0
       i_bar(55)=488.0
       i_bar(56)=491.0
       i_bar(57)=501.0
       i_bar(58)=523.0
       i_bar(59)=535.0
       i_bar(60)=546.0
       i_bar(61)=560.0
       i_bar(62)=574.0
       i_bar(63)=580.0
       i_bar(64)=591.0
       i_bar(65)=614.0
       i_bar(66)=628.0
       i_bar(67)=650.0
       i_bar(68)=658.0
       i_bar(69)=674.0
       i_bar(70)=684.0
       i_bar(71)=694.0
       i_bar(72)=705.0
       i_bar(73)=718.0
       i_bar(74)=727.0
       i_bar(75)=736.0
       i_bar(76)=746.0
       i_bar(77)=757.0
       i_bar(78)=790.0
       i_bar(79)=790.0
       i_bar(80)=800.0
       i_bar(81)=810.0
       i_bar(82)=823.0
       i_bar(83)=823.0
       i_bar(84)=830.0
       i_bar(85)=825.0
       i_bar(86)=794.0
       i_bar(87)=827.0
       i_bar(88)=826.0
       i_bar(89)=841.0
       i_bar(90)=847.0
       i_bar(91)=878.0
       i_bar(92)=890.0
c extension by calc. + extrapolation, H.Weick 2012 for V1.3
       i_bar(93)=  900.0
       i_bar(94)=  910.0
       i_bar(95)=  920.0
       i_bar(96)=  930.0
       i_bar(97)=  940.0
       i_bar(98)=  950.0
       i_bar(99)=  960.0
       i_bar(100)= 970.0
       i_bar(101)= 980.0
       i_bar(102)= 990.0
       i_bar(103)=1000.0
       i_bar(104)=1010.0
       i_bar(105)=1020.0
       i_bar(106)=1030.0
       i_bar(107)=1040.0
       i_bar(108)=1050.0
       i_bar(109)=1060.0
       i_bar(110)=1070.0
       i_bar(111)=1080.0
       i_bar(112)=1090.0
       i_bar(113)=1100.0
       i_bar(114)=1110.0
       i_bar(115)=1120.0
       i_bar(116)=1130.0
       i_bar(117)=1140.0
       i_bar(118)=1150.0
       i_bar(119)=1160.0
       i_bar(120)=1170.0
       
       ipot = i_bar(nz)
      
       END FUNCTION
       
c This subroutine erases the file with name file name
c recursive was introduced in f90
       RECURSIVE FUNCTION uniquefile(filename) RESULT(newfile)
       IMPLICIT NONE
       
       CHARACTER*35 filename, newfile
       CHARACTER*10 time
       INTEGER stat
       
       newfile = filename
c Try to open the file       
       OPEN(UNIT=1234, IOSTAT=stat, FILE=newfile, STATUS='old')
       CLOSE(1234)
       
       IF (stat .EQ. 0) THEN
c If it is exists try again to generate a unique file name       
           CALL date_and_time(TIME=time)
           newfile = uniquefile(filename(1:25)//time)
       ENDIF

       END FUNCTION

c This subroutine erases the file with name file name
       SUBROUTINE delfile(filename)
       IMPLICIT NONE
       
       CHARACTER*35 filename
       INTEGER stat
       
       OPEN(UNIT=1234, IOSTAT=stat, FILE=filename, STATUS='old')
       IF (stat .EQ. 0) CLOSE(1234, STATUS='delete')
       
       END SUBROUTINE

c This subroutine is kept for compatibility with AtimaPy1.2
       SUBROUTINE calculate(proj, pn, ein,
     &                targ, nnuca, tm, rhoa, igasa, thick,  
     &                eout, range, dedxi,dedxo,remra,rstr,estr,
     &                astr, tof, intrthick )
       IMPLICIT NONE
c An array defining [A, Z] of the projectile, here pn should be 2 
       INTEGER pn
       REAL*8  proj(2)
Cf2py intent(in) proj
Cf2py integer intent(hide),depend(proj) :: pn=shape(proj,0)      

c Energy IN [MeV/u]      
       REAL*8 ein
Cf2py intent(in) ein

c An array defining [[A0, Z0, SR0], ... ] the target composition, here SRi
c is the stochastic ratio of Ai, Zi 
c here nnuca is number of target compounds and tm = 3
       INTEGER nnuca, tm
       REAL*8  targ(nnuca, tm)
Cf2py intent(in) targ
Cf2py integer intent(hide),depend(targ) :: nnuca=shape(targ,0), tm=shape(targ,1)       

c Material density in g/cm3
       REAL*8  rhoa
Cf2py intent(in) rhoa

c 1 for gas 0 for solid
       INTEGER  igasa      
Cf2py intent(in) igasa

c target thickness in [mg/cm2]
       REAL*8  thick      
Cf2py intent(in) thick

c Energy OUT [MeV/u]
       REAL*8 eout
Cf2py intent(out) eout

c Range in mg/cm2
       REAL*8  range
Cf2py intent(out) range

c dE/dX in the first 1 mg/cm2 of material [MeV/u/mg/cm2]
       REAL*8 dedxi 
Cf2py intent(out) dedxi

c dE/dX in the last 1 mg/cm2 of material [MeV/u/mg/cm2]       
       REAL*8 dedxo
Cf2py intent(out) dedxo

c Remaining Range in [mg/cm2]              
       REAL*8 remra
Cf2py intent(out) remra

c Range straggling [mg/cm2]
       REAL*8 rstr
Cf2py intent(out) rstr

c Energy straggling [mg/cm2]
c Multiply by dedxo to get the straggling in MeV/u
       REAL*8 estr
Cf2py intent(out) estr

c Angular straggling in [rad]      
       REAL*8 astr
Cf2py intent(out) astr

c Time Of Flight through the target in [ns]       
       REAL*8 tof
Cf2py intent(out) tof

c Returns the interpolated target thickness in [mg/cm2]
c If the ions stops in the material thick = range       
       REAL*8 intrthick
Cf2py intent(out) intrthick

c Beam Energy. If positive this is the energy IN the target. If negative, this denotes the energy after the target [MeV/u]      
       REAL*8 ebeam
       
c Make sure that we ebeam is the energy of the beam going into the target
       ebeam = ABS(ein) 
       CALL eloss(proj, pn, ebeam,
     &                targ, nnuca, tm, rhoa, igasa, thick,  
     &                ein, eout, range, dedxi,dedxo,remra,rstr,estr,
     &                astr, tof, intrthick )
     
       END SUBROUTINE
