!.....Composition module for the Hempel and Schaffner-Bielich EOS v1.03
!.....RMF parameterization: ASB
!.....Nuclear mass table: FRDM

!.....Note about usage:
!.....The main subroutine which calculates and interpolates the nuclear 
!.....distributions for arbitrary temprature t, baryon number density nb
!.....and electron fraction ye is 'sub_dist_interpol', whereas a 
!.....trilinear interpolation scheme is used for the variables 
!.....log(t), log(nb), ye. If the input variables are out of the range
!.....of the table, the success flag 'sflag' is set to zero, and the
!.....abundances of all nuclei and nucleons, too.
!.....There is also 'sub_dist_grid' if the composition shall be 
!.....calculated for a specific grid point of the EOS without 
!.....interpolation. 
!.....The other two subroutines in this module 'sub_coulomb'
!.....and 'sub_degen' are called automatically from 'sub_dist_***', 
!.....and they are needed only internally.
!.....Before 'sub_dist_***' can be used, you need to read in the data
!.....file by calling the subroutine 'compdata_readin' once.
!.....Apart from the output of 'sub_dist_***' you have to use the
!.....array 'az', because it specifies the mass and charge number of
!.....every nucleus. The variables are explained further where they
!.....appear for the first time.

!.....Note about units:
!.....This module uses natural units internally. This means that the
!.....units of all quantities are powers of MeV. If a number density
!.....with unit [MeV^3] has to be converted to [fm^-3] one has to devide
!.....by hqc**3. This is done only for the number densities which are
!.....the output of the main subroutine 'sub_dist'.
      module asb_frdm_composition_module

      implicit none

!.....hqc: reduced Planck`s constant times c [MeV fm]
!.....pi: have a guess
      double precision, parameter, private :: hqc = 197.3269631d0
      double precision, parameter, private :: pi = 4.0d0*datan(1.0d0)
      
!.....t_entries: number of temperature entries
!.....y_entries: number of electron fraction entries, not counting ye=0
!.....nb_entries: number of nb entries
!.....c_entries: number of entries per grid-point in the data file
      integer, parameter, private :: t_entries = 81
      integer, parameter, private :: y_entries = 60
      integer, parameter, private :: nb_entries = 326
      integer, parameter, private :: c_entries = 8          
      
!.....define the constants of the t, ye, and nb array
      double precision, parameter, private :: nbin = 1.0d-12
      double precision, parameter, private :: tin = 0.1d0
      double precision, parameter, private :: yin = 0.0d0
      double precision, parameter, private :: dnb = 0.04d0
      double precision, parameter, private :: dt = 0.04d0
      double precision, parameter, private :: dy = 0.01d0
        
!.....parameters, which depend on the used EOS table:
!.....kmax: number of nuclei
!.....nb0: saturation density of the RMF parameterization [MeV^3]
      integer, parameter :: kmax = 8140
      double precision, parameter, private :: 
     & nb0 = 1.5824150805d-01*hqc**3

!.....comp(i,j,k,l): array with the input data
!.....i=1, ..., t_entries: temperature index
!.....j=0, ..., y_entries: electron fraction index
!.....k=1, ..., nb_entries: baryon number density index
!.....l=1, ..., c_entries: composition data entry index
      double precision,
     & dimension(t_entries,0:y_entries,nb_entries,c_entries),
     & private :: comp

!.....various arrays with the properties of the nuclei
!.....mass: mass [MeV]
!.....bind: binding energy [MeV]
!.....degen: degeneracy factor []
!.....coul: Coulomb energy [MeV]
!.....az(inuc,1): mass number of nucleus  'inuc' []
!.....az(inuc,2): charge number of nucleus 'inuc' []
      double precision, dimension(kmax), private :: mass, bind, degen,
     & coul
      integer, dimension(kmax,2) :: az

      contains

!.....Subroutine which reads in the binary data file (generated with
!.....ifort on a 64 bit machine) for the calculation. Has to be called
!.....once.
        subroutine compdata_readin

        implicit none

        open(10,file='asb_frdm_comp_v1.03.bin',form='unformatted',
     & status='old')
        read(10)az,mass,bind,comp
        close(10)

        return
        end subroutine compdata_readin

!.....Main subroutine which calculates the nuclear distributions,
!.....i.e. the nuclear and nucleon densities and mass fractions for
!.....given temperature, electron fraction and baryon density
        subroutine sub_dist_interpol(t,ye,nb,xaz,xn,xp,naz,nn,np,sflag)

        implicit none

!.....input:   
!.....t: temperature [MeV]
!.....ye: electron fraction []
!.....nb: baryon number density [fm^-3]
        double precision, intent(in) :: t,ye,nb

!.....output:
!.....naz(inuc): number density of nucleus 'inuc' [fm^-3]
!.....nn: number density of neutrons [fm^-3]
!.....nn: number density of protons [fm^-3]
!.....xaz(inuc): mass fraction of nucleus 'inuc' []
!.....xn: mass fraction of neutrons []
!.....xp: mass fraction of protons []
!.....sflag: success flag
!.....sflag=0: input is not inside the table and the composition cannot
!.....         be calculated
!.....sflag=1: succesful call of the routine

        double precision, dimension(kmax), intent(out) :: naz, xaz
        double precision, intent(out) :: xn, xp, nn, np
        integer, intent(out) :: sflag

!.....internal variables:
!.....i: temperature index
!.....j: electron fraction index
!.....k: baryon number density index
        integer :: i, j, k, l, indx(1:8,1:3)      
!.....fnb: total baryon number density [MeV^3]
!.....fnp: total proton number density (nuclei and free protons) [MeV^3]
        double precision :: fnb, fnp       
!.....x, y, z: interpolation variables
        double precision :: x, y, z       
!.....various variables used for the interpolation        
        double precision  :: tlocal
        double precision, dimension(kmax) :: degen1, degen2
        double precision, dimension(1:8,kmax) :: nazlocal, alphalin
        double precision, dimension(1:8,1:2) :: nucleonslocal, alphalin2
        double precision, dimension(1:2) :: nucleons
!.....trilinear interpolation matrix
        integer, parameter :: trilm(1:8,1:8) =
     & reshape( (/
     & 1,-1,-1, 1,-1, 1, 1,-1,
     & 0, 1, 0,-1, 0,-1, 0, 1,
     & 0, 0, 1,-1, 0, 0,-1, 1,
     & 0, 0, 0, 1, 0, 0, 0,-1,
     & 0, 0, 0, 0, 1,-1,-1, 1,
     & 0, 0, 0, 0, 0, 1, 0,-1,
     & 0, 0, 0, 0, 0, 0, 1,-1,
     & 0, 0, 0, 0, 0, 0, 0, 1
     & /),(/8,8/))  
      
        sflag=0
        naz=0.0d0
        xaz=0.0d0
        xn=0.0d0
        xp=0.0d0
        nn=0.0d0
        np=0.0d0
        
!.....determine the indices of the grid corresponding to the cube around
!.....t, nb, ye
        if (t.lt.tin) then
          write(*,*)"Error: t too low"
          write(*,*)t,tin
          return
        end if
        if (t.gt.10.0d0**(dble(t_entries-1)*dt)*tin) then
          write(*,*)"Error: t too large"
          return
        end if
        i=int(dlog10(t/tin)/dt)+1
        if (i.ge.t_entries) then
          i=t_entries-1
        end if
        if (ye.lt.yin) then
          write(*,*)"Error: Y_e too low"
          write(*,*)ye,yin
          return
        end if
        if (ye.gt.dble(y_entries)*dy+yin) then
          write(*,*)"Error: Y_e too large"
          return
        end if
        j=int((ye-yin)/dy)
        if (j.eq.y_entries) then
          j=y_entries-1
        end if
        if (nb.lt.nbin) then
          write(*,*)"Error: n_B too low"
          write(*,*)nb,nbin
          return
        end if
        if (nb.gt.10.0d0**(dble(nb_entries-1)*dnb)*nbin) then
          write(*,*)"Error: n_B too large"
          return
        end if
        k=int(dlog10(nb/nbin)/dnb)+1
        if (k.eq.nb_entries) then
          k=nb_entries-1
        end if           
            
        indx(1,1)=i
        indx(1,2)=j
        indx(1,3)=k

        indx(2,1)=i+1
        indx(2,2)=j
        indx(2,3)=k

        indx(3,1)=i
        indx(3,2)=j+1
        indx(3,3)=k

        indx(4,1)=i+1
        indx(4,2)=j+1
        indx(4,3)=k

        indx(5,1)=i
        indx(5,2)=j
        indx(5,3)=k+1

        indx(6,1)=i+1
        indx(6,2)=j
        indx(6,3)=k+1

        indx(7,1)=i
        indx(7,2)=j+1
        indx(7,3)=k+1

        indx(8,1)=i+1
        indx(8,2)=j+1
        indx(8,3)=k+1

!.....calculate the internal degeneracy for the two different
!.....temperatures involved in the interpolation
        call sub_degen(tin*10.0d0**((i-1)*dt))
        degen1=degen
        call sub_degen(tin*10.0d0**((i)*dt))
        degen2=degen  
            
        nazlocal=0.0d0
        nucleonslocal=0.0d0             
        do l=1,8
!.....calculate the distributions for each of the eight corners if 
!.....nuclei are present    
          if (comp(indx(l,1),indx(l,2),indx(l,3),5).ne.1.0d0) then
!.....set the local temperature
            tlocal=tin*10.0d0**((indx(l,1)-1)*dt)
!.....set the internal degeneracy            
            if (indx(l,1).eq.i) then
              degen=degen1
            else
              degen=degen2
            end if
!.....calculate the coulomb energies
            call sub_coulomb(comp(indx(l,1),indx(l,2),indx(l,3),6))           
!.....calculate the distribution
            nazlocal(l,:)=(mass*tlocal/2.0d0/pi)**1.5d0
     & *dexp((dble(az(:,1))*(comp(indx(l,1),indx(l,2),indx(l,3),1)
     & -comp(indx(l,1),indx(l,2),indx(l,3),4)/nb0)+dble(az(:,2))
     & *comp(indx(l,1),indx(l,2),indx(l,3),2)+bind-coul)/tlocal)
     & *degen*comp(indx(l,1),indx(l,2),indx(l,3),3)
     & *(1.0d0-comp(indx(l,1),indx(l,2),indx(l,3),5))
          end if 
          
!.....set the nucleon densities          
          nucleonslocal(l,1)=comp(indx(l,1),indx(l,2),indx(l,3),7)
          nucleonslocal(l,2)=comp(indx(l,1),indx(l,2),indx(l,3),8)
        end do

!.....set the interpolation variables
        x=dlog10(t/tin)/dt-dble(i-1)
        y=(ye-yin)/dy-dble(j)        
        z=dlog10(nb/nbin)/dnb-dble(k-1)

!....trilinear interpolation        
        alphalin=matmul(trilm,nazlocal)
        alphalin2=matmul(trilm,nucleonslocal)
        nucleons=0.0d0
        do i=0,1
          do j=0,1
            do k=0,1
              naz(:)=naz(:)+alphalin(1+i+2*j+4*k,:)*x**i*y**j*z**k
              nucleons(:)=nucleons(:)
     & +alphalin2(1+i+2*j+4*k,:)*x**i*y**j*z**k
            end do
          end do
        end do      

!.....calculate the the total densities        
        fnb=dot_product(naz,dble(az(:,1)))+nucleons(1)+nucleons(2)
        fnp=dot_product(naz,dble(az(:,2)))+nucleons(2)
        
!.....calculate the mass fractions, sum is normalized to 1        
        xaz=naz*dble(az(:,1))/fnb
        xn=nucleons(1)/fnb
        xp=nucleons(2)/fnb

!.....convert the number densities which are output to [fm^-3]
        naz=naz/hqc**3
        nn=nucleons(1)/hqc**3
        np=nucleons(2)/hqc**3
        
        sflag=1 
        
        return
        end subroutine sub_dist_interpol


!.....Subroutine which calculates the nuclear distributions,
!.....i.e. the nuclear and nucleon densities and mass fractions for
!.....given temperature index i, electron fraction index j,
!.....and baryon number density index k
        subroutine sub_dist_grid(i,j,k,xaz,xn,xp,naz,nn,np)

        implicit none

!.....input:
!.....i,j,k: see above
        integer, intent(in) :: i, j, k

!.....output:
!.....naz(inuc): number density of nucleus 'inuc' [fm^-3]
!.....nn: number density of neutrons [fm^-3]
!.....nn: number density of protons [fm^-3]
!.....xaz(inuc): mass fraction of nucleus 'inuc' []
!.....xn: mass fraction of neutrons []
!.....xp: mass fraction of protons []
        double precision, dimension(kmax), intent(out) :: naz, xaz
        double precision, intent(out) :: xn, xp, nn, np

!.....internal variables:
!.....t: temperature [MeV]
!.....fnb: total baryon number density [MeV^3]
!.....fnp: total proton number density (nuclei and free protons) [MeV^3]
        double precision :: t, fnb, fnp

!.....calculate the temperature
        t=tin*10.0d0**((i-1)*dt)

        naz=0.0d0
        xaz=0.0d0
        xn=0.0d0
        xp=0.0d0
        fnb=0.0d0
        fnp=0.0d0

!.....calculate the distribution only if nuclei are present
        if (comp(i,j,k,5).ne.1.0d0) then
!.....calculate the temperature dependent degeneracy
          call sub_degen(t)
!.....calculate the Coulomb energies
          call sub_coulomb(comp(i,j,k,6))
!.....calculate the distribution
          naz=(mass*t/2.0d0/pi)**1.5d0
     & *dexp((dble(az(:,1))*(comp(i,j,k,1)-comp(i,j,k,4)/nb0)
     & +dble(az(:,2))*comp(i,j,k,2)+bind-coul)/t)
     & *degen*comp(i,j,k,3)*(1.0d0-comp(i,j,k,5))
!.....calculate the contribution of nuclei to the total densities
          fnb=dot_product(naz,dble(az(:,1)))
          fnp=dot_product(naz,dble(az(:,2)))
        end if

!.....add neutrons and protons to the total densities
        fnb=fnb+comp(i,j,k,7)+comp(i,j,k,8)
        fnp=fnp+comp(i,j,k,8)

!.....calculate the mass fractions, sum is normalized to 1
        xaz=naz*dble(az(:,1))/fnb
        xn=comp(i,j,k,7)/fnb
        xp=comp(i,j,k,8)/fnb

!.....convert the number densities which are output to [fm^-3]
        naz=naz/hqc**3
        nn=comp(i,j,k,7)/hqc**3
        np=comp(i,j,k,8)/hqc**3

        return
        end subroutine sub_dist_grid


!......Subroutine which calculates the Coulomb energies for given
!......ne: electron number density [MeV^3]
        subroutine sub_coulomb(ne)

        implicit none

!.....alpha: fine structure constant
        double precision, parameter :: alpha = 7.2973525376d-3
        
        double precision, intent(in) :: ne
        double precision, dimension(1:kmax) :: fkern, rkern

        coul=0.0d0
        if (ne.ne.0.0d0) then
          fkern=(ne*dble(az(:,1))/nb0/dble(az(:,2)))**(1.0d0/3.0d0)
          rkern=(0.75d0/pi*dble(az(:,1))/nb0)**(1.0d0/3.0d0)
          coul=0.5d0*dble(az(:,2))**2*alpha/rkern*(0.6d0*fkern**3
     & -1.8d0*fkern)
        end if

        return
        end subroutine sub_coulomb


!......Subroutine which calculates the internal degeneracy for given t
        subroutine sub_degen(t)

        implicit none

        double precision, intent(in) :: t
        double precision, dimension(1:kmax) :: spin
        double precision, dimension(1:kmax) :: u1, u2
        double precision, dimension(1:kmax) :: xm, xr, dx, x1, x2
        double precision :: w(5),x(5)
        SAVE w,x
        DATA w/.2955242247d0,.2692667193d0,.2190863625d0,.1494513491d0,
     &.0666713443d0/
        DATA x/.1488743389d0,.4333953941d0,.6794095682d0,.8650633666d0,
     &.9739065285d0/
        integer i

        degen=0.0d0
        if (t.ne.0.0d0) then
          u1=-dsqrt(dble(az(:,1))*0.125d0*
     & (1.0d0-0.8d0*dble(az(:,1))**(-1.0d0/3.0d0))*t*0.5d0)
          u2=u1+dsqrt(bind/t)
          xm=0.5d0*(u1+u2)
          xr=0.5d0*(u2-u1)
          do i=1,5
            dx=xr*x(i)
            x1=xm+dx
            x2=xm-dx
            degen=degen
     & +w(i)*(dexp(u1**2-x1**2)*(x1-u1)+dexp(u1**2-x2**2)*(x2-u1))
          end do
          degen=0.4d0/dble(az(:,1))**(5.0d0/3.0d0)*t*xr*degen
          if (ANY(degen.gt.1.0d300)) then
            write(*,*)'Error: Degen out of range!'
            stop
          end if
        end if
        do i=1,kmax
          if ((az(i,1)+1)/2*2-az(i,1).eq.0) then
            spin(i)=0.0d0
          else
            spin(i)=0.5d0
          end if
          if ((az(i,1).eq.2).and.az(i,2).eq.1) then
            spin(i)=1.0d0
          end if
        end do

        degen=degen+2.0d0*spin+1.0d0

        return
        end subroutine sub_degen
        
      end module asb_frdm_composition_module