! Copyright (c) 2022 Le-Cnam/INRiM.

! MIT-like license with commercial restrictions.

! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! “Software”), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute and sublicense copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! 1. The above copyright notice and this permission notice shall be
!    included in all copies or substantial portions of the Software.

! 2. The Software, or part of it, is not included or used in commercial
!    applications without prior written agreement from Le-Cnam/INRiM.

! 3. Except as contained in this notice, the name of Le-Cnam/INRiM shall
!    not be used in advertising or otherwise to promote the sale, use or
!    other dealings in this Software without prior written authorization
!    from Le-Cnam/INRiM.

! THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND,
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
! NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
! LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
! OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
! WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

! COMPILE the LIBRARY:
!--------------------------------------------------------------------------------
! On   LINUX:
! gfortran -o gfer.so -shared -static-libgfortran -static-libgcc -fPIC gfer.f08
! on WINDOWS:
! gfortran -o gfer.dll -shared -static-libgfortran -static-libgcc -fPIC gfer.f08

! ONLY FOR DEBUG:
! gfortran -o gfer.so -static-libgfortran -static-libgcc -fPIC gfer.f08 ../Helium_EoS/He4_LabView16/He4.f08


MODULE gfer

  ! COSTANTS
  Real(8)		::	pi 				= 4.0*atan(1.0)
  Real(8)	  ::	zAC(13)		= (/ 0.0, 4.4934094579,7.7252518369,10.904121659,14.066193912, &
       17.220755271,20.371302959,23.519452498,26.666054258, &
       29.811598790,32.956389039,36.100622244, 39.244432361 /)  
  Real(8)		::  zTE(12)		= (/ 4.4934094579,7.7252518369,10.904121659,14.066193912, &
       17.220755271,20.371302959,23.519452498,26.666054258, &
       29.811598790,32.956389039,36.100622244, 39.244432361 /)
  Real(8)		::  zTM(13)		= (/2.7437072699,6.1167642642,9.3166156285,12.485937368, &
       15.643866106,18.796253353,21.945518067,25.092847019, &
       28.238892111,31.384041635,34.528541442,37.672554848,40.816194611/)
  
	type constants
	  Real(8)		::	 R	= 8.3144598d0;					! J / (mol K)
	endtype constants
	type(constants) :: cnst
!end Module gfsr
contains

   ! Function for predicting the frequencies of the acoustic and electromagnetic modes triaxial-ellipsoids

   ! Return the radius of the equivalent sphere using nominal ellipsoid semi-axes
   function gfsr_Req(Rx, Ry, Rz) bind (C, name='Ell_Req')
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)	::	Rx, Ry, Rz, gfsr_Req
     gfsr_Req = (Rx*Ry*Rz)**(1.d0/3.0d0)
   end function gfsr_Req

   !------------------------------------------------------------------------------------------------------------------------------
   ! ACOUSTIC MODES

   function f_id(nn, w, RReq) bind(C, name='f_id')															 ! Frequency in an ideal resonator
     use, intrinsic :: iso_c_binding, only : c_double, c_int
     real(c_double)	:: w, RReq, Req, f_id
     integer(c_int) :: nn
     Req = RReq * 1d-3
     f_id = zAC(nn)/(2.0d0*pi*Req)*w
   end function f_id
   
   function d_t				(K, rho, cp, f)	bind(C,name='d_t')												! Thermal penetration length [m]
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::	K, rho, cp, f, d_t								
     ! d_t				=		sqrt(2*Dt/(2*pi*f))
     d_t				=	sqrt ( K/(pi*rho*cp*f) )
   end function d_t

   function d_v				(Nu, rho, f) bind(C,name='d_v')														! Viscous penetration length [m]
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::		Nu, f, d_v, rho  																				! kinetic viscosity eta/rho m^2/s, Hz
     ! d_v				=			sqrt(2*Dv/(2*pi*f))
     d_v				= sqrt ( Nu/(pi*rho*f) )
   end function d_v

   function I_t				(T, pp, cv, K, m0, h) bind(C,name='I_t')												! Thermal accomodation length [m]
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::		K, T, pp, p, cv, h, m0, I_t
     p = pp*1d3
     I_t				=			K/p * sqrt( (pi*m0*T)/(2d0*cnst%R) ) * (2d0-h)/h *&
          						(1 / (cv*m0/cnst%R + 0.5d0))
   end function I_t

   function df_th			(g, dt, dt_Cu, I_th, K_gas, K_Cu, aa, f)	bind(C,name='df_th')									! Thermal boundary layer correction [Hz]
     ! ToDO: Split contribution terms
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::		g, dt, dt_Cu, I_th, K_gas, K_Cu, a, aa, f, df_th
     Integer		::		n
     a = aa * 1d-3
     !df_th			=			-( (g-1d0)/(2d0*a)*dt*(1d0-2d0*I_th/dt))*f 							! Gavioso et al.
     df_th = -(g-1)/(2*a)*dt + (g-1)/a*I_th + (g-1)/(2*a)*dt_Cu*K_gas/K_Cu			! Relative value
     df_th = df_th * f																													! Absolute value in hertz
   end function df_th

   function df_dth	(g, aa, dth) bind(C, name='df_dth')
     use, intrinsic :: iso_c_binding, only: c_double
     Real(c_double)	::g, aa, a, dth, df_dth
     a	=	aa * 1d-3
     df_dth = -(g-1)/(2*a)*dth
   end function df_dth
   
	 function df_Ith	(g, aa, I_th) bind(C, name='df_Ith')
     use, intrinsic :: iso_c_binding, only: c_double
     Real(c_double)	::g, aa, a, I_th, df_dth
     a	=	aa * 1d-3
     df_Ith = -(g-1)/(2*a)*I_th
   end function df_Ith

 	 function df_CuTh	(g, aa, dt_Cu, K_g, K_Cu) bind(C, name='df_CuTh')
     use, intrinsic :: iso_c_binding, only: c_double
     Real(c_double)	::g, aa, a, dt_Cu, K_g, K_Cu
     a	=	aa * 1d-3
     df_CuTh = -(g-1)/(2*a)*dt_Cu*K_g/K_Cu
   end function df_Cuth

   function g_t 			(T,pp,g,d_th, dt_Cu, K_Gas, K_Cu, aa, f)	bind(C,name='g_t')											! Contribution of thermal boundary layer to resonance half-width [Hz]
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::		T,pp,p,g,d_th,dt_Cu, K_Gas, K_Cu, aa,f, g_t
     p = pp*1d3; 	a = aa * 1d-3;
     ! g_t 				=			(g-1)/(2*a)*d_th*f																			! Gavioso et al.
     g_t		= (g-1)/(2*a)*d_th + (g-1)/(2*a)*dt_Cu*K_gas/K_Cu - &
              0.5d0*(g-1)*(2*g-1)*(d_th/a)**2																	  ! Relative value
     g_t		=	g_t*f																															! Absolute value in hertz
   end function g_t

   function g_th 			(g, Req, d_th)	bind(C,name='g_th')												! Separated contribution of thermal boundary layer to resonance half-width [Hz]
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::		g, Req, d_th, a, g_th
     a = Req * 1d-3;
     ! g_t 				=			(g-1)/(2*a)*d_th*f																			! Gavioso et al.
     g_th		= (g-1)/(2*a)*d_th 						  																	  ! Relative value
   end function g_th

   function g_thO2 			(g, Req, d_th)	bind(C,name='g_thO2')										! Separated contribution of thermal boundary layer to resonance half-width [Hz]
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::		g, Req, d_th, a, g_thO2
     a = Req * 1d-3;
     g_thO2		=  0.5d0*(g-1)*(2*g-1)*(d_th/a)**2																  ! Relative value
   end function g_thO2
   
	 function g_Cu 		(g, Req, dt_Cu, K, K_Cu)	bind(C,name='g_Cu')								! Separated contribution of thermal boundary layer to resonance half-width [Hz]
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::		g, Req, d_th, a, dt_Cu, K, K_Cu, g_Cu
     a = Req * 1d-3
     g_Cu		=  (g-1)/(2*a)*dt_Cu*K/K_Cu																				  ! Relative value
   end function g_Cu
   
   function g_b				(T,pp,w,gamma, d_vs, d_th, f) bind(C,name='g_b')					! Contribution of bulk attenuation sound to resonance half-width [Hz]
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double)		::		T,pp,p,w,gamma, d_vs, d_th, f, g_b
     p = pp * 1d3
     g_b				=			pi**2/w**2 * (4d0/3d0*d_vs**2 + (gamma-1)*d_th**2) * f**3d0
   end function g_b
   
    subroutine df_ducts		(rr, LL, N, dns, cp, w, gamma, k_g, nu_g, RReq, f, nn, df_f, dg_f) bind(C, name='df_ducts')
			use, intrinsic :: iso_c_binding, only: c_double, c_int										! Ducts effect on frequency and half-width [df/f, dg/f]
			! Zt is assumed to be negligible
			Integer(c_int)		::		N, nn
			Real(c_double)		::		rr(N), r(N), LL(N), L(N), dns, cp, w, gamma, k_g, nu_g, Req, RReq, f
			Real(c_double)		::		df_f, dg_f
			Real(c_double)		::		A(N)
			Complex(8)				::		kv, kt, Fv(N), Ft(N), G(N), Z0(N), Za(N), Zb(N)
			Complex(8)				::		Ze													! Acoustic impedence

   		r = rr * 1d-3;		L = LL * 1d-3; Req = RReq * 1d-3
			! Tubes sections
			A			= pi*r**2
			! Physical effects
			kv	= (1d0, -1d0)/d_v(nu_g, dns, f);															kt = (1d0,-1d0)/d_t(k_g,dns, cp,f);
			Fv	= 2d0 * BesselJ(1, kv*r)/(kv*r*BesselJ(0,kv*r));    				  Ft = 2d0 * BesselJ(1, kt*r)/(kt*r*BesselJ(0,kt*r));
			G		=	(0d0, 1d0)*(2*pi*f/w)*sqrt( (1+(gamma-1)*Ft)/(1-Fv) )
			! Acoustic impedences
			Z0	= (dns*w/A) / sqrt( (1d0+(gamma-1)*Ft)*(1-Fv) )
			Za	= Z0 * tanh(G*L/2d0);																					Zb = Z0 / sinh(G*L)
		
			! Acoustic impedence calculated iteratively
			!Ze		= Za(N) + Zb(N)		 ! Zt -> infinity
			Ze	 = (0d0,0d0)
			do i=N,1,-1
				 Ze = Za(i) + (Za(i) + Ze)*Zb(i)/(Za(i) + Zb(i) + Ze)
			end do

			 Ze		= ((0d0,1d0)*dns*w/Ze)/(4d0*pi*Req**2*zAC(nn))																! Ze is used as buffer
			 df_f = realpart(Ze);		dg_f = imagpart(Ze)
		end subroutine df_ducts

   function df_mic(rr, aa, rho, w) bind(C,name='df_mic')														! Effects of microphone [df/f] 
     use, intrinsic :: iso_c_binding, only : c_double
     Real(c_double) :: chi = 1d-10
     Real(c_double)	:: rr, r, aa, a, rho, w, df_mic
     ! df_mic = - (r/(2d0*a))**2 * g*p*chi/a																		! Mehl et al.
     r = rr * 1d-3;			a = aa * 1d-3;
     df_mic = -rho*w**2/2d0 * r**2/a**3 * chi																		! NPL report
   end function df_mic
   
   function O2dZ2AC		(Eps,n)	bind(C,name='O2dZ2AC')														! Second order triaxial-ellipsoid perturbation on  z^2 for (0,x) modes [dz^2/z^2]
     use, intrinsic :: iso_c_binding, only : c_double, c_int
     !use gfsr
     Real(c_double)		::  	Eps(2), xi, O2dZ2AC
     Integer(c_int)		::  	n
     xi 					=		zAC(n)
     O2dZ22AC			=		8*xi**2/135d0 * (Eps(1)**2 - Eps(1)*Eps(2) + Eps(2)**2)
   end function O2dZ2AC

   function Cp_Cu(T) bind(C, name='cp_cu')
     use, intrinsic 	:: iso_c_binding, only: c_double
     real(c_double)		:: T, Cp_Cu, logCp

		 logCp = 1.131d0 - 9.454d0*log(T)/log(10.0) + 12.99d0*log(T)**2/log(10.0)**2 &
        - 5.501d0*log(T)**3/log(10.0)**3 + 0.7637d0*log(T)**4/log(10.0)**4
   		Cp_Cu = 10.0**logCp
   end function Cp_Cu

   function K_Cu(T, RRR) bind (C, name='k_cu')
     use, intrinsic 	:: iso_c_binding, only: c_double
     real(c_double)		:: T, RRR, K_Cu
     real(c_double)		:: Wo, Wi, Wio, Wc = 0.0
     real(c_double)		:: beta, beta_r
     real(c_double)		:: P1=1.754d-8, P2=2.763, P3=1102, P4=-0.165, &
          P5=70, P6=1.756, P7

     beta=0.634/RRR
     beta_r=beta/0.0003
     P7=0.838/(beta_r**0.1661)
     
     Wo		=	 beta/T
     Wi		=  P1*T**P2 / (1+ P1*P3*T**(P2+P4)*exp(-(P5/T)**P6)) + Wc
     Wio  =  P7*Wi*Wo / (Wi +Wo)

     K_Cu = 1d0/(Wo + Wi + Wio)
   end function K_Cu
           
   !------------------------------------------------------------------------------------------------------------------------------
   ! ELECTROMAGNETIC MODES
   function d_sk (mu_r, s, f) bind(C, name='d_sk')
     use, intrinsic		:: iso_c_binding, only: c_double
     real(c_double)		:: mu_r, s, f, d_sk
     real(c_double)		:: mu0 = 1.25663706212d-6											  ! CODATA
     d_sk	=	1d0/sqrt( pi*mu_r*mu0*s*f)																! May 2004
   end function d_sk

   subroutine df_skTE (d_sk, aa, df_f, dg_f) bind(C, name='df_skTE')
     use, intrinsic		:: iso_c_binding, only: c_double
     real(c_double)		:: d_sk, a, aa, df_f, dg_f
     complex(8)				:: dZ
     a    = 	aa * 1d-3
     dZ		=		d_sk/a*(-1d0, 1d0)/2d0																	! May 2004
     df_f	=		realpart(dZ);			dg_f = imagpart(dZ)
   end subroutine df_skTE

   subroutine df_skTM1n (d_sk, aa, nn, df_f, dg_f) bind(C, name='df_skTM1n')
     use, intrinsic		:: iso_c_binding, only: c_double, c_int
     real(c_double)		:: d_sk, aa, a, df_f, dg_f
     integer(c_int)		:: nn
     complex(8)				:: dZ
     a    = aa * 1d-3
     dZ		=		d_sk/a*(-1d0, 1d0)/2d0 * zTM(nn)**2/(zTM(nn)**2 - 2)				! May 2004
     df_f	=		realpart(dZ);			dg_f = imagpart(dZ)
   end subroutine df_skTM1n

   function df_dcTM1n (rr, aa, nn) bind(C, name='df_dcTM1n')
     use, intrinsic		:: iso_c_binding, only: c_double, c_int
     real(c_double)		:: r,rr,aa, a, df_dcTM1n
     integer(c_int)		:: nn
     ! ToDO: rescale r and check using Underwood, correction on each component
     r 		= rr * 1d-3
     a    = aa * 1d-3
     df_dcTM1n	=	-r**3/(4d0*pi*a**3) * (0.950*zTM(nn)**2 - 1.152)/(zTM(nn)**2-2)  ! Underwood 2010
   end function df_dcTM1n

   function df_dcTE (rr, aa) bind(C, name='df_dcTE')
     use, intrinsic		:: iso_c_binding, only: c_double
     real(c_double)		:: rr,r,aa, a, df_dcTE
     r    = rr * 1d-3
     a    = aa * 1d-3
     df_dcTE	=	-r**3/(4d0*pi*a**3) * 0.950						! Underwood 2010
   end function df_dcTE

   function hTM (xi) 
     Real(8)	:: xi, hTM
     hTM = (2*xi**2 + 4) / (15*xi**2 - 2)
   end function hTM

   function hTE (xi)
     Real(8)	:: xi, hTE
     hTE = (xi**2 + 4) / (xi**2 - 2) 
   end function hTE

   subroutine EpsTM (f1, f2, f3, n, ep1, ep2)	bind(C,name='Eps_TM1n')	 ! Epsilon from TM1x modes
     use, intrinsic :: iso_c_binding, only: c_double, c_int
     Real(c_double) :: f1,f2,f3,f0, xi
     Real(c_double) :: f(3), ep1, ep2
     Integer(c_int) :: m, n
   	 xi	= zTM(n)
     f 	= SortFreq(f1,f2,f3)
     f0	= sum(f)/3.0d0					
     ep1 = (f(2)**2 - f(1)**2)/( 3*f0**2*hTM(xi))
     ep2 = (f(2)**2 - f(3)**2)/( 3*f0**2*hTM(xi)) 
   end subroutine EpsTM

   subroutine EpsTE (f1, f2, f3, n, ep1, ep2) bind(C, name='Eps_TE1n') ! Epsilon from TE1x modes
		 use, intrinsic :: iso_c_binding, only: c_double, c_int
     Real(c_double) :: f1,f2,f3,f0, xi
     Real(c_double) :: f(3), ep1, ep2
		 Integer(c_int) :: n
   	 xi	= zTE(n)
     f 	= SortFreq(f1,f2,f3)
     f0	= sum(f)/3.0d0					
     ep1= (f(2)**2 - f(1)**2)/( 3*f0**2*hTE(xi))
     ep2= (f(2)**2 - f(3)**2)/( 3*f0**2*hTE(xi)) 
   end subroutine EpsTE

   function dzm2_zTM1n(ep1, ep2,n) bind(C, name='dzm2_zTM1n')								! Second order relative shift of z^2 for TM1x modes
     use, intrinsic :: iso_c_binding, only: c_double, c_int
     Real(c_double)		::  ep1, ep2,xi, O2dZm2TM
     Integer(c_int)		::  n
     xi 				=		zTM(n)
     O2dfm2TM		=		2d0*(33d0*xi**8 - 245d0*xi**6 + 714d0*xi**4 - 1152d0*xi**2 + 160d0) / (1125d0*(xi**2-2d0)**3) * &
     								(ep1**2 - ep1*ep2 + ep2**2)
   end function dzm2_zTM1n
   
   function dzm2_zTE1n(ep1, ep2,n) bind(C, name='dzm2_zTE1n')								! Second order relative shift of z^2 for TE1x modes
     use, intrinsic:: iso_c_binding, only: c_double, c_int
     Real(c_double)		::  ep1, ep2, xi, O2dZm2E
     Integer(c_int)		::  n
     xi 				=		zTE(n)
     O2dfm2TE		=		(22d0*xi**2/375d0 - 2d0/225d0) * (ep1**2 - ep1*ep2 + ep2**2)
   end function dzm2_zTE1n

    
   !------------------------------------------------------------------------------------------------------------------------------
   ! Supporting functions
   
   function SortFreq (f1, f2, f3) result(f)
     Real(8)	:: f1, f2, f3
     Real(8)  :: f(3), buf(7)
     integer	:: i,j
     
     ! Sort frequency values
     buf = -1 			! Buffer initialization
     buf(4) = f1
     if (f2 .lt. f1) then
        buf(2) = f2
     else
        buf(6)=f2
     end if

     do i=1,7
        if (f3 .lt. buf(i)) then
           buf(i-1)=f3
           exit
        end if
        if (i .eq. 7) then
           buf(7)=f3
        end if
     end do
     ! Copy sorted values
     j = 1
     do i=1,7
        if (buf(i) .ne. -1) then
           f(j) = buf(i)
           j = j + 1
        end if
     end do
   end function SortFreq
   
	 elemental function BesselJ (n,z)
			! Implementation of
			! David J. Sookne, "Bessel functions I and J of Complex Argument and integer order",
			!									J. Nat. Res. Bureau of Standards - B. Mathematical Science, Vol 77B, Nos 3&4, July-December 1973
			complex(8),intent(IN)		::  z
			integer, intent(IN)    	::  n
			complex(8)		::	mu, BesselJ
			integer				::	nN
			integer				::	M,i,j, LL
			complex(8)		::	p(3)
			real(8)				::	k, l
			integer,parameter ::	deg = 10
			real(8)				::	mL, rho
			real(8)				::	Q								! sqrt(2*10**deg*abs( p_L * p_(L+1) ))		eq.(4)
  	  !complex(8), allocatable :: y(:)			! eq.(1)
		  !complex(8)		:: y(500)			! eq.(1)
  	  complex(8), pointer :: y(:)			! eq.(1)

			M				=	floor(abs(z));					LL = max(M, n);

			p(1)			= 0d0
			p(2)			= 1d0
			p(3)			= 2*(M+1)/z

			do i=1,200					 
				p(1) = p(2);			p(2)	= p(3);				p(3)	= 2*(M+1)/z*p(2) - p(1)
				if (M+i-1 .eq. LL) then
					 Q = sqrt(2*10d0**deg*abs( p(2) * p(1)))
				end if
			
				k			=	abs(p(3)/p(2));	l			= (M+i)/abs(z) + sqrt( (M+i+1)/abs(z)**2 - 1)
				rho		= min(k, l)
			
				if ( abs(p(2)) .ge. (sqrt(rho/(rho**2-1))*Q) ) then
					 exit
				end if							
			end do
			nN = M+i

  	  BesselJ = NormBesselJ(n, z, nN, i, p(3))
   					
 end function BesselJ

pure function NormBesselJ(n, z, nN, i, p3)
   integer, intent(IN) 		:: n, nN, i
   complex(8),intent(IN)  :: z, p3
   complex(8)   :: mu, NormBesselJ
   integer			:: j
   complex(8)		:: y(nN+1)
 			! Backward iteration
			!allocate(y(nN+1))
			y(nN+1) = 0;				y(nN)=1d0/p3;
			do j=nN+1,3,-1											
				y(j-2) = 2*(j-2)/z*y(j-1) - y(j)
			end do
	 
			mu = 0
			! Normalization of Jn with Im(z)>0
			if (imag(z) .gt. 0d0) then
			 do j=1,nN
					mu = mu + 2*(0,-1)**j*y(j+1)
			 end do
			 mu = (mu + y(1))* exp((0,1)*z)
			end if
			! Normalization of Jn with Im(z)=0
			if (imag(z) .eq. 0d0) then
			 do j=1,floor(i/2d0)
					mu = mu + 2*y(j*2)
			 end do
			 mu = mu + y(1)
			end if
			! Normalization of Jn with Im(z)<0
			if (imag(z) .lt. 0d0) then
			 do j=1,i
					mu = mu + 2*(0,1)**j*y(j)
			 end do
			 mu = (mu + y(1))* exp((0,-1)*z)
			end if

   NormBesselJ = y(n+1)/mu
 end function NormBesselJ

 end module gfer
 
