Update code for ELOSS fortran subroutine

TELoss/ELOSS.f90

@@ -4,320 +4,256 @@
 !== File: ELOSS.f90
 !== Date: 08/12/2000
 !======================================================================
-!  ELOSS.f90 
+!  ELOSS.f90
 !
 !  FUNCTIONS/SUBROUTINES exported from ELOSS.dll:
-!	ELOSS      - subroutine 
+!  ELOSS      - subroutine
+!======================================================================
+!== Date 2024-05-03
+!   Some reworking of code:
 !
-subroutine ELOSS(NEL,ZEL,AEL,WEL,DEN,ZP,AP,NE,ETAB,RE,ZW,AW)
+! * 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_')
 
-  ! Expose subroutine ELOSS to users of this DLL
-  !
-  !DEC$ ATTRIBUTES DLLEXPORT::ELOSS
 
-	IMPLICIT REAL*8(A-H,O-Z)
-	INTERFACE
-	SUBROUTINE rde(E, R, R1, R2, IRE)
-		REAL*8   E, R, R1, R2
-		INTEGER  IRE
-	END SUBROUTINE rde
-	END INTERFACE
+  IMPLICIT none
 
   ! Variables
 
   ! Input:
-	INTEGER NEL, NE
-	REAL*8  ZEL(NEL),AEL(NEL),WEL(NEL), DEN,ZP,AP, ETAB(NE)
+  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*8  RE(NE),ZW,AW
-
-  ! Local
-	COMMON adj(5),as(5),z(5),f(5),um(5),oz,ro,azm,aj,zion,aion,&
-	       zmed,amed,jc
+  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
-	
-	!== two parameters with question
-	ire  = 1
-	oz   = 1.
-	DO i = 1, NEL
-		as(i) = AEL(i)
-		z(i)  = ZEL(i)
-		um(i) = WEL(i)
-		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=0.
-	zam=0.
-	umed=0.
-	zme=0.
-	aj=0.
-	DO i = 1, jc
-		uma=uma+um(i)*as(i)
-		umed=umed+um(i)
-		zme=zme+um(i)*z(i)
-	END DO
-
-	amed=uma/umed
-	zmed=zme/umed
-	
-	DO i = 1, jc
-		f(i)=um(i)*as(i)/uma
-	END DO
-
-	DO i = 1, jc
-		zam=zam+f(i)*z(i)/as(i)
-	END DO
-
-	azm=1./zam
-	aj=0.
-	
-	DO i = 1, jc
-		aj=aj+f(i)*z(i)*log(adj(i))/as(i)
-	END DO
-	
-	aj=aj*azm
-	aj=exp(aj)
-	ZW = zmed
-	AW = amed
-
-	DO j = 1, NE
-		!== result [mg/cm^2]
-		CALL rde(ETAB(j), RE(j), rel1, rel, ire)
-	END DO
+  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)
-	IMPLICIT REAL*8(A-H,O-Z)
-
-	INTERFACE 
-		REAL*8 FUNCTION c(x)
-			REAL*8 x
-		END FUNCTION c
-
-		REAL*8 FUNCTION c1(x)
-			REAL*8 x
-		END FUNCTION c1
-	END INTERFACE
-
-
-
-! calculates range and de/dx for compounds
-	dimension ala(3),ala1(3)
-	dimension a(3,3),a1(4,4),a2(4,4),b(2,3),cc(5)
-	common	adj(5),as(5),z(5),f(5),um(5),oz,ro,azm,aj,zion,aion,&
-			zmed,amed,jc
-	dimension alaj1(4),altau1(4)
-
-	data a/	-.75265, .073736, .040556, 2.5398,-.312,.018664, &
-			-.24598, .11548, -.0099661/
-
-	data a1/-8.0155,     0.36916,    -1.4307e-02,  3.4718e-03,	&
-			 1.8371,    -1.452e-02,  -3.0142e-02,  2.3603e-03,	&
-			 0.045233,  -9.5873e-04,  7.1303e-03, -6.8538e-04,	&
-			-5.9898e-03,-5.2315e-04, -3.3802e-04,  3.9405e-05/
-
-	data a2/-8.725,		 0.8309,	 -0.13396,	    0.012625,	&
-			 1.8797,	 0.11139,	 -0.064808,		0.0054043,	&
-			 0.74192,	-0.528805,	  0.1264232,   -0.00934142,	&
-			 0.752005,  -0.5558937,   0.12843119,  -0.009306336/
-
-	data b/	0.422377e-06,   3.858019e-09, 0.0304043e-06, -0.1667989e-09, &
-           -0.00038106e-06, 0.00157955e-09/
-
-	if(e.le.0.002)then
-		range=0.
-		rel1=0.
-		rel=0.
-		return
-	endif
-
-! this a flag
-	ibis=-1
-
-	en	=	e/aion
-	tau	=	en*1.008
-
-	altau	=	log(tau)
-	alaj	=	log(aj)
-	alaj1(1)=	1.
-	altau1(1)=	1.
-	DO kk = 2,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)
-
-
-	s1=0.
-	DO i=1,4
-		DO j=1,4
-			s1	=	s1 + a2(j,i)*alaj1(j)*altau1(i)
-		END DO
-	END DO
-
-	ala(1) = azm*exp(s1)
-	
-	s2=0.
-	DO i=2,4
-		DO j=1,4
-			s2	=	s2 + (i-1)*a2(j,i)*alaj1(j)*altau1(i-1)
-		END DO
-	END DO
-	ala1(1)=ala(1)*s2/tau
-
-	s1=0.
-	DO i=1,4
-		DO j=1,4
-			s1	=	s1 + a1(j,i)*alaj1(j)*altau1(i)
-		END DO
-	END DO
-
-	ala(3) = azm*exp(s1)
-
-	s2=0.
-	DO i=2,4
-		DO j=1,4
-			s2	=	s2 + (i-1)*a1(j,i)*alaj1(j)*altau1(i-1)
-		END DO
-	END DO
-	ala1(3)=ala(3)*s2/tau
-
-	s1=0.
-	DO i=1,3
-		DO j=1,3
-			s1	=	s1	+	a(j,i)*alaj1(j)*altau1(i)
-		END DO
-	END DO
-	ala(2)=azm*exp(s1)/1000.
-
-	s2=0.
-	DO i=2,3
-		DO j=1,3
-			s2	=	s2 + (i-1)*a(j,i)*alaj1(j)*altau1(i-1)
-		END DO
-	END DO
-	ala1(2)=ala(2)*s2/tau
-
-25	continue
-	coefa=3.
-	coefb=1.
-
-	alaa=(ala(1)*(tanh(coefa*(.98 - en))+1.)/2.			&
-		+ ala(2)*(tanh(coefa*(en - .98))+1.)/2.)		&
-		*        (tanh(coefb*(8.0 - en))+1.)/2.			&
-		+ ala(3)*(tanh(coefb*(en  - 8.))+1.)/2.
-
-
-	alim1=0.
-	alim2=0.
-	if(coefa*(en-.98).lt.85)then
-		alim1=1.008/cosh(coefa*(.98-en))/cosh(coefa*(.98-en))
-	endif
-	if(coefb*(en-8.).lt.85)then
-		alim2=1.008/cosh(coefb*(8.-en))/cosh(coefb*(8.-en))
-	endif
-
-	alaa1=(ala1(1)*(tanh(coefa*(.98-en))+1.)/2.+		&
-		ala1(2)*(tanh(coefa*(en-.98))+1.)/2.)*			&
-		(tanh(coefb*(8.-en))+1.)/2.+					&
-		ala1(3)*(tanh(coefb*(en-8.))+1.)/2.+			&
-		coefa/2.*(tanh(coefb*(8.-en))+1.)/2.*			&
-		(ala(2)*alim1-ala(1)*alim1)+					&
-		coefb/2.*(ala(3)*alim2-							&
-		(ala(1)*(tanh(coefa*(.98-en))+1.)/2.+			&
-		ala(2)*(tanh(coefa*(en-.98))+1.)/2.)*alim2)
-
-
-	hi=137.*beta/zion
-	bz=(31.8+3.86*(aj**.625))*azm*.000001
-	bz=bz*(zion**2.666667)*c(hi)
-	bz1=(4.357+.5288*(aj**.625))*azm*.001
-	bz1=bz1*(zion**1.666667)*c1(hi)
-	bep=beta*beta
-	rel1=zion*zion/(alaa1+bz1*((1.-bep)**1.5)/931.141/beta)
-	rel1=rel1/1000.
-	range=(alaa+bz)*aion/1.008/zion/zion
-	range=range*1000.
-
-! Atention!! this version do not work correctly for ix.ne.1
-	if(ix.ne.1)return
-	ank=.153574*ro/azm
-	z23=zion**.666667
-	abet=beta*125./z23
-	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
-	del1=del1/ro
-	del2=2.*ank*zef2*(log(oz/aj)-cbet)/bep
-	del2=del2/ro
-	rel2=rel1-del1
-	rel3=rel1-rel2+del2
-	if(del1)8,8,9
-8	rel=rel1
-	goto 12
-9	if(del1+del2-rel2)10,10,11
-10	rel=rel3
-	goto 12
-11	if(del1.lt.rel1)then
-	 rel=rel2
-	 else
-	 rel=rel1
-	endif
-12	return
-	END
-
-
-	REAL*8 FUNCTION c(x) 
-	REAL*8 x
-	IF(x .LE. 0.2) THEN
-		c = -0.00006 + (0.05252 + 0.12857*x)*x
-	ELSE IF(x .LE. 2.) THEN
-		c = -0.00185 + (0.07355 + (0.07172 - 0.02723*x)*x)*x
-	ELSE IF(x .LE. 3.) THEN
-		c = -0.0793  + (0.3323  - (0.1234  - 0.0153*x)*x)*x
-	ELSE
-		c = 0.22
-	END IF
-	END
-
-	REAL*8 FUNCTION c1(x)
-	REAL*8 x
-
-	IF(x .LE. 0.2) THEN
-		c1 = 0.05252+.25714*x
-	ELSE IF(x .LE. 2.0) THEN
-		c1 = 0.07355 + (0.14344 - 0.08169*x)*x
-	ELSE IF(x .LE. 3.0) THEN
-		c1 = 0.3323  - (0.2468  - 0.0459*x)*x
-	ELSE
-		c1 = 0.
-	END IF
-	
-	END
-
-
+  ! 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