program read_table
implicit none
integer   :: alpha           ! Number of bins
integer   :: Nvarchimie      ! Number of species
integer   :: nvar            ! Equal to Nvarchimie
integer   :: nion=9          ! Number of non-grain species

integer,parameter::dp=kind(1.0d0) ! default
!!!! Grains info
real(dp), allocatable, dimension(:) :: x_g        ! Abundance of grains per bins
real(dp), allocatable, dimension(:) :: r_g        ! Radius of grains per bins
real(dp), allocatable, dimension(:) :: m_g        ! Mass of grains per bins

!!!! Constants
real(dp), parameter :: pi=3.1415927410125732422_dp  
real(dp), parameter :: mp=1.6726d-24     ! Proton mass in g
real(dp), parameter :: me=9.1094d-28     ! Electron mass in g
real(dp), parameter :: e=4.803204d-10    ! Electron charge in cgs
real(dp), parameter :: Kb = 1.3807d-16   ! Boltzmann constant erg/K
real(dp), allocatable, dimension(:) :: q ! Electric charge 
real(dp), allocatable, dimension(:) :: m ! Mass

!!!! bins of grains
real(dp), parameter :: rho_s=2.3_dp        ! Bulk density g.cc
real(dp), parameter :: a_0=0.0375d-4       ! Reference radius cm
real(dp), parameter :: a_min=0.0181d-4     ! Minimum radius cm
real(dp), parameter :: a_max=0.9049d-4     ! Maximum radius cm
real(dp), parameter :: zeta=a_min/a_max    ! a_min/a_max
real(dp), parameter :: lambda_pow=-3.5d0   ! Coeff power law
real(dp)            :: rho_gtot            ! Grain density
real(dp) :: Lp1,Lp3,Lp4

! resistivites (cf Kunz & Mouschovias 2009, Marchand et al. 2016)
real(dp), allocatable, dimension(:)  :: sigma             
real(dp), allocatable, dimension(:)  :: zetas              
real(dp), allocatable, dimension(:)  :: phi
real(dp), allocatable, dimension(:)  :: tau_sn,tau_gp,tau_gm
real(dp), allocatable, dimension(:,:)  :: tau_inel
real(dp), allocatable, dimension(:)  :: gamma_zeta,gamma_omega
real(dp), allocatable, dimension(:)  :: omega,omega_bar
real(dp), parameter        :: c_l=299792458.d2     ! c in cm/s
real(dp)            :: sigv, muuu        
integer  :: iB,iH,iT,iX
real(dp) :: B,bmaxchimie,nH,T,xi,sigH,sigO,sigP

integer :: cmp_X

! For loops
integer :: tchimie   ! tchmie steps in temperature
real(dp) :: dtchimie   ! dtchmie step in temperature
real(dp) :: tminchimie   ! min tchmie 
integer :: nchimie   ! nchmie steps in density 
real(dp) :: dnchimie   ! dnchmie step in density 
real(dp) :: nminchimie   ! min nchmie  
integer :: bchimie   ! bchmie steps in B field
real(dp) :: dbchimie   ! dbchmie step in B field
real(dp) :: bminchimie   ! min bchmie 
integer :: xichimie ! Steps in ionisation rate
real(dp) :: dxichimie ! dxichimie step in ion rate
real(dp) :: ximinchimie ! min xichimie
integer :: nislin,tislin,xiislin ! Linear scale or not
real(dp),allocatable,dimension(:,:,:,:,:)::conductivities ! resistivities
real(dp),allocatable,dimension(:,:,:,:)::abundances ! abundances
integer :: i,j,k,ij


real(dp) :: sig_p, sig_perp, sig_h, eta_ohm, eta_hall, eta_ad


   open(42,file='Table_abundances.dat', status='old')
   read(42,*) nchimie, tchimie, xichimie, nvar       ! Number of  steps in density, temperature and ionisation rate
   read(42,*) nislin, tislin, xiislin                      ! 1: Constant step in logscale, 0: Not constant step in logscale
   read(42,*)
   read(42,*)
   read(42,*)
   allocate(abundances(-2:nvar,nchimie,tchimie,xichimie)) 
   ! abundances(-2,:,:,:)=density
   ! abundances(-1,:,:,:)=temperature
   ! abundances(0,:,:,:)=ionisation rate
   ! abundances(n>°,:,:,:)=relative abundance of species n


   do ij=1,xichimie
     do i=1,tchimie
        do j=1,nchimie  
           read(42,*)abundances(-2:nvar,j,i,ij)
        end do
        read(42,*)
        read(42,*)
     end do
   end do
   close(42)
   nminchimie=(abundances(-2,1,1,1))   !min density
   tminchimie=(abundances(-1,1,1,1))   !min temperature
   ximinchimie=(abundances(0,1,1,1))   !min ionisation rate
   ! Density, temperature and IR steps in log (if nislin, tislin, xiislin = 1)
   dnchimie=(log10(abundances(-2,nchimie,1,1))-log10(abundances(-2,1,1,1)))/(nchimie-1)
   dtchimie=(log10(abundances(-1,1,tchimie,1))-log10(abundances(-1,1,1,1)))/(tchimie-1)
   dxichimie=(log10(abundances(0,1,1,xichimie))-log10(abundances(0,1,1,1)))/(xichimie-1)


   !Allocating arrays
   alpha=(nvar-nion)/3

   allocate(x_g(alpha))
   allocate(r_g(alpha))
   allocate(m_g(alpha))
   allocate(q(nvar))
   allocate(m(nvar))
   allocate(sigma(nvar))
   allocate(zetas(nvar))
   allocate(phi(nvar))
   allocate(tau_sn(nvar))
   allocate(tau_gp(nvar))
   allocate(tau_gm(nvar))
   allocate(gamma_zeta(nvar))
   allocate(gamma_omega(nvar))
   allocate(omega(nvar))
   allocate(omega_bar(nvar))
   allocate(tau_inel(nvar,nvar))
   

   ! Computation of grains radii and species mass
   Lp1=dble(lambda_pow+1)
   Lp3=dble(lambda_pow+3)
   Lp4=dble(lambda_pow+4)

    ! Grain radius
   if  (alpha==1) then
     r_g(1)=a_0
   else
     do  i=1,alpha    ! cf Marchand et al. 2016
       r_g(i)=a_min*zeta**(-dble(i)/dble(alpha)) * &
            & dsqrt( Lp1/Lp3* (1d0-zeta**(Lp3/dble(alpha)))/(1d0-zeta**(Lp1/dble(alpha))))
     end do
   end if


   q(:)=1.d0*e    ! cations
   q(1)=-1.d0*e   ! electrons
   do  i=nion+1,nvar
      if (mod(i-nion,3)==1) q(i)=1.d0*e   ! g+
      if (mod(i-nion,3)==2) q(i)=-1.d0*e  ! g-
      if (mod(i-nion,3)==0) q(i)=0.d0     ! g0
   end do
   m(:)=0.d0
   m(1) = me           ! e-
   m(2) = 23.5d0*mp    ! ions metalliques
   m(3) = 29.d0*mp     ! ions moleculaires
   m(4) = 3*mp         ! H3+
   m(5) = mp           ! H+
   m(6) = 12.d0*mp     ! C+
   m(7) = 4.d0*mp      ! He+
   m(8) = 39.098*mp    ! K+
   m(9) = 22.99d0*mp   ! Na+
   do  i=1,alpha       ! masse des grains
      m_g(i)=4.d0/3.d0*pi*r_g(i)**3*rho_s
      m(nion+3*(i-1)+1:nion+3*i)=m_g(i)
   end do


   ! values for Btable
   bminchimie=1.d-10
   bmaxchimie=1.d10               ! ok for first core in nimhd. maybe not enough for second core.
   bchimie=150
   dbchimie=(log10(bmaxchimie)-log10(bminchimie))/real((bchimie-1),dp)

   allocate(conductivities(0:3,1:nchimie,1:tchimie,1:xichimie,1:bchimie))
   !conductivities(1,:,:,:)= sigma_//
   !conductivities(2,:,:,:)= sigma_perp
   !conductivities(3,:,:,:)= sigma_hall
   !conductivities(0,:,:,:)= sigma_hall sign


  ! Computation of resistivities

   tau_sn      = 0.0_dp
   omega       = 0.0_dp
   sigma       = 0.0_dp
   phi         = 0.0_dp
   zetas       = 0.0_dp
   gamma_zeta  = 0.0_dp
   gamma_omega = 0.0_dp
   omega_bar   = 0.0_dp


do  iB=1,bchimie
do  iX=1,xichimie
do  iT=1,tchimie
do  iH=1,nchimie

    nh=abundances(-2,iH,iT,iX)                ! density (.cc) of current point
    B =10.0d0**(log10(bminchimie)+dble(iB-1)*dbchimie)  ! Magnetic field
    T =abundances(-1,iH,iT,iX)                ! Temperature
    xi =abundances(0,iH,iT,iX)                ! Ionisation rate


    do i=1,nion
      if  (i==1) then  ! electron
        sigv=3.16d-11 * (dsqrt(8d0*kB*1d-7*T/(pi*me*1d-3))*1d-3)**1.3d0
        tau_sn(i) = 1.d0/1.16d0*(m(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nH/2.d0*sigv)
      else if (i>=2 .and. i<=nion) then ! ions
        muuu=m(i)*2d0*mp/(m(i)+2d0*mp)
        if (i==2 .or. i==3) then
          sigv=2.4d-9 *(dsqrt(8d0*kB*1d-7*T/(pi*muuu*1d-3))*1d-3)**0.6d0
        else if (i==4) then
          sigv=2d-9 * (dsqrt(8d0*kB*1d-7*T/(pi*muuu*1d-3))*1d-3)**0.15d0
        else if (i==5) then
          sigv=3.89d-9 * (dsqrt(8d0*kB*1d-7*T/(pi*muuu*1d-3))*1d-3)**(-0.02d0)
        else
          sigv=1.69d-9
        end if
        tau_sn(i) = 1.d0/1.14d0*(m(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nH/2.d0*sigv)
      end if
      omega(i) = q(i)*B/(m(i)*c_l)
      sigma(i) = abundances(i,iH,iT,iX)*nH*(q(i))**2*tau_sn(i)/m(i)
      phi(i) = 0.d0
      zetas(i) = 0.d0
      gamma_zeta(i) = 1.d0
      gamma_omega(i) = 1.d0
      omega_bar(i) = 0.d0
    end do
    
    do  i=1,alpha   ! grains
      tau_sn(nion+1+3*(i-1))=1.d0/1.28d0*(m_g(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nH/2.d0*(pi*r_g(i)**2*(8.d0*Kb*T/(pi*2.d0*mp))**0.5))
      omega(nion+1+3*(i-1)) = q(nion+1+3*(i-1))*B/(m_g(i)*c_l)
      sigma(nion+1+3*(i-1)) = abundances(nion+1+3*(i-1),iH,iT,iX)*nH*(q(nion+1+3*(i-1)))**2*tau_sn(nion+1+3*(i-1))/m_g(i)
    
      tau_sn(nion+2+3*(i-1))=tau_sn(nion+1+3*(i-1))
      omega(nion+2+3*(i-1)) = q(nion+2+3*(i-1))*B/(m_g(i)*c_l)
      sigma(nion+2+3*(i-1)) = abundances(nion+2+3*(i-1),iH,iT,iX)*nH*(q(nion+2+3*(i-1)))**2*tau_sn(nion+2+3*(i-1))/m_g(i)
    
    end do
  
    sigP=0.d0
    sigO=0.d0
    sigH=0.d0

    do i=1,nvar
       sigP=sigP+sigma(i)
       sigO=sigO+sigma(i)/(1.d0+(omega(i)*tau_sn(i))**2)
       sigH=sigH-sigma(i)*omega(i)*tau_sn(i)/(1.d0+(omega(i)*tau_sn(i))**2)
    end do

    conductivities(1,iH,iT,iX,iB)=log10(sigP)
    conductivities(2,iH,iT,iX,iB)=log10(sigO)
    conductivities(3,iH,iT,iX,iB)=log10(abs(sigH))
    conductivities(0,iH,iT,iX,iB)=sign(1.0_dp,sigH)

    !conductivities(:,:,:,:) contains the log10(conductivities)


    eta_ohm  = 1d0/sigP
    eta_hall = sigH/(sigO**2 + sigH**2)
    eta_ad   = sigO/(sigO**2 + sigH**2)-1d0/sigP

end do
end do
end do
end do



end program
