!======================================================================
!== Energy loss in material
!==
!== File: ELOSS.f90
!== Date: 08/12/2000
!======================================================================
!  ELOSS.f90
!
!  FUNCTIONS/SUBROUTINES exported from ELOSS.dll:
!  ELOSS      - subroutine
!======================================================================
!== Date 2024-05-03
!   Some reworking of code:
!
! * remove obsolete fortran instructions;
! * replace common block with module;
! * change precision to standard C double;
! * remove restriction for 5 elements in material.

module m_eloss
  use iso_c_binding, only: c_int,c_double

  real(c_double),allocatable,dimension(:),private:: adj,as,z,f,um
  real(c_double),private:: oz,ro,azm,aj,zion,aion, zmed,amed
  integer,private :: jc

contains
  subroutine ELOSS(NEL,ZEL,AEL,WEL,DEN,ZP,AP,NE,ETAB,RE,ZW,AW)&
       bind(C,name='eloss_')


  IMPLICIT none

  ! Variables

  ! Input:
  INTEGER(c_int),intent(in) ::NEL ! number of elements in the material
  INTEGER(c_int) ::NE ! number of points in E -- Range table
  real(c_double),dimension(nel),intent(in):: zel
  real(c_double),dimension(nel),intent(in):: ael
  real(c_double),dimension(nel),intent(in):: wel
  real(c_double),intent(in) :: den
  real(c_double),intent(in) :: ZP,AP
  real(c_double),dimension(ne),intent(in):: Etab

  ! Output:
  real(c_double),dimension(ne),intent(out) :: RE
  REAL(c_double),intent(out) :: ZW,AW

!!$ local variables
  integer :: i,j,ire
  real(c_double) :: rel,rel1,uma,umed,zam,zme
  ! Body of ELOSS

  zion = ZP
  aion = AP
  jc   = NEL
  ro   = DEN
  allocate(adj(jc),as(jc),z(jc),f(jc),um(jc))

  !== two parameters with question
  ire  = 1
  oz   = 1.

  as=ael;z = ZEL; um = WEL
  DO concurrent (i = 1:jc)
     IF(z(i) .LE. 13.) THEN
        adj(i) = 12.*z(i) + 7.
     ELSE
        adj(i) = 9.76*z(i) + 58.8 / z(i)**0.19
     END IF
  END DO

  uma=sum(um*as)
  zme=sum(um*z)
  zam=zme/uma;  azm=1./zam
  umed=sum(um)
  aj=exp(sum(um*z*log(adj))/uma*azm)

  amed=uma/umed
  zmed=zme/umed

  f=um*as/uma

  ZW = zmed
  AW = amed

  DO j = 1, NE
     !== result [mg/cm^2]
     CALL rde(ETAB(j), RE(j), rel1, rel, ire)
  END DO
  deallocate(adj,as,z,f,um)
end subroutine ELOSS

SUBROUTINE rde(e,range,rel1,rel,ix)
  ! calculates range and de/dx for compounds
  IMPLICIT none !REAL(c_double) (A-H,O-Z)

  real(c_double),intent(in) ::E
  integer,intent(in) ::ix
  real(c_double),intent(out) ::range
  real(c_double),intent(out) ::rel1
  real(c_double),intent(out) ::rel

  integer :: i,j,k,kk

  real(c_double),parameter,dimension(3,3)::a=reshape([&
       -0.75265, .073736, .040556,&
       +2.53980,-.312000, .018664, &
       -0.24598, .115480,-.0099661],shape=[3,3])
  real(c_double),parameter,dimension(4,4)::a1=reshape([&
       -8.0155   , 0.36916  ,-1.4307e-2, 3.4718e-3, &
       +1.8371   ,-1.4520e-2,-3.0142e-2, 2.3603e-3, &
       -4.5233e-2,-9.5873e-4, 7.1303e-3,-6.8538e-4, &
       -5.9898e-3,-5.2315e-4,-3.3802e-4, 3.9405e-5],shape=[4,4])
  real(c_double),parameter,dimension(4,4)::a2=reshape([&
       -8.725000, 0.8309000,-0.13396000, 0.012625, &
       -1.879700, 0.1113900,-0.06480800, 0.0054043, &
       -0.741920,-0.5288050, 0.12642320,-9.341420e-3,&
       -0.752005,-0.5558937, 0.12843119,-9.306336e-3 ],shape=[4,4])
  real(c_double),parameter,dimension(2,3)::b=reshape([&
       +4.223770e-07, 3.858019e-09, 3.04043e-08,&
       -1.667989e-10,-3.810600e-10, 1.57955e-12],shape=[2,3])

  real(c_double),allocatable,dimension(:) :: cc
  real(c_double),dimension(3) :: ala, ala1
  real(c_double) :: alaa,alaa1
  real(c_double),dimension(4) :: alaj1, altau1
  real(c_double) :: abet, alaj, altau, ank, bep, beta, bz, bz1,cbet, t, tau, tt
  real(c_double),parameter ::  coefa=3.,  coefb=1.
  real(c_double) :: alim1, alim2, del1, del2, en, hi, om, rel2, rel3, zef, zef2

  if(e.le.0.002)then
     range=0.
     rel1=0.
     rel=0.
     return
  endif

  en = e/aion
  tau = en*1.008

  altau = log(tau)
  alaj = log(aj)
  DO concurrent (kk = 1:4)
     alaj1(kk)=alaj**(kk-1)
     altau1(kk)=altau**(kk-1)
  END DO
  t = tau/938.256
  tt = 2.*t + t*t
  beta= sqrt(tt)/(1.+t)


  ala(1) = azm*exp(dot_product(alaj1,matmul(a2,altau1)))
  ala(2) = azm*exp(dot_product(alaj1(1:3),matmul(a,altau1(1:3))))/1000.
  ala(3) = azm*exp(dot_product(alaj1,matmul(a1,altau1)))

  ala1(1)=ala(1)*dot_product(alaj1,matmul(a2(:,2:4),altau1(1:3)*[1,2,3]))/tau
  ala1(2)=ala(2)*dot_product(alaj1(1:3),matmul(a(:,2:3),altau1(1:2)*[1,2]))/tau
  ala1(3)=ala(3)*dot_product(alaj1,matmul(a1(:,2:4),altau1(1:3)*[1,2,3]))/tau

  associate(ca=>coefa*(.98 - en),cb=>coefb*(8.0 - en))
    associate(tca=>(1+tanh(ca))/2.,tcb=>(1+tanh(cb))/2.)
      alaa=(ala(1)*tca + ala(2)*(1.0-tca))*tcb + ala(3)*(1.0-tcb)

      alim1=0.; alim2=0.
      if(-ca.lt.85)then
         alim1=1.008*(cosh(ca)**(-2))
      endif
      if(-cb.lt.85)then
         alim2=1.008*(cosh(cb)**(-2))
      endif

      alaa1=(ala1(1)*tca+ ala1(2)*(1.0-tca))*tcb + ala1(3)*(1.0-tcb)+&
           coefa*tcb* (ala(2)-ala(1))*alim1/2.+   &
           coefb*(ala(3)-(ala(1)*tca+ala(2)*(1.-tca)))*alim2/2.
    end associate
  end associate

  hi=137.*beta/zion
  bz=(31.8+3.86*(aj**.625))*azm* 1e-6 *(zion**2.666667)*c(hi)
  bz1=(4.357+.5288*(aj**.625))*azm* .001 *(zion**1.666667)*c1(hi)
  bep=beta*beta
  rel1=zion*zion/(alaa1+bz1*((1.- beta**2)**1.5)/931.141/beta)/1000.
  range=(alaa+bz)*aion/(1.008*zion**2)*1000.

  ! Atention!! this version do not work correctly for ix.ne.1
  if(ix.ne.1) then
     return
  end if

  allocate(cc(size(as)))

  ank=.153574*ro/azm
!  z23=zion**.666667
  abet=beta*125.*(zion** (-2.0/3.0))
  zef=zion*(1.-exp(-abet))
  zef2=zef*zef
  om= 1022000. *bep/(1.-bep)
  cbet=0.
  DO k=1,jc
     cc(k)=0.
     DO i=1,2
        DO j=1,3
           cc(k)=cc(k)+b(i,j)*((1./bep-1.)**j)*((adj(k)**(i+1)))
        END DO
     END DO
     cbet=cbet+f(k)*cc(k)/as(k)
  END DO
  cbet=cbet*azm
  del1=ank*zef2*(log(om/oz)-bep)/bep/ro
  del2=2.*ank*zef2*(log(oz/aj)-cbet)/bep/ro
  rel2=rel1-del1
  rel3=rel1-rel2+del2
!!$ here is a replasement for an Arithmetic IF construct
  if (del1 <=0.0 ) then
     rel=rel1
  else if (del1+del2-rel2<=0) then
     rel=rel3
  else if(del1.lt.rel1)then
     rel=rel2
  else
     rel=rel1
  endif
  deallocate(cc)
  return
END SUBROUTINE rde


pure FUNCTION c(x) result(res)
   REAL(c_double),intent(in) :: x
   real(c_double) :: res
  IF(x .LE. 0.2) THEN
     res = -0.00006 + (0.05252 + 0.12857*x)*x
  ELSE IF(x .LE. 2.) THEN
     res = -0.00185 + (0.07355 + (0.07172 - 0.02723*x)*x)*x
  ELSE IF(x .LE. 3.) THEN
     res = -0.0793  + (0.3323  - (0.1234  - 0.0153*x)*x)*x
  ELSE
     res = 0.22
  END IF
END FUNCTION c

pure FUNCTION c1(x) result(res)
  REAL(c_double),intent(in) :: x
  real(c_double) :: res
  IF(x .LE. 0.2) THEN
     res = 0.05252+.25714*x
  ELSE IF(x .LE. 2.0) THEN
     res = 0.07355 + (0.14344 - 0.08169*x)*x
  ELSE IF(x .LE. 3.0) THEN
     res = 0.3323  - (0.2468  - 0.0459*x)*x
  ELSE
     res = 0.
  END IF

END FUNCTION c1

end module m_eloss