!!****if* source/physics/materialProperties/MagneticResistivity/MagneticResistivityMain/SpitzerHighZ/MagneticResistivity_fullState
!!
!! NAME
!!  MagneticResistivity_fullState
!!
!! SYNOPSIS
!!  call MagneticResistivity_fullState(real(in)    :: solnVec(NUNK_VARS),
!!                                     real(out)   :: resPar,
!!                            OPTIONAL,real(out)   :: resPerp)
!!
!! DESCRIPTION
!!
!! Computes the Spitzer electron Magnetic Resistivity for all materials,
!! including those with Z > 1. The specific equations used here all
!! come from "The Physics of Inertial Fusion" by Atzeni.
!!
!!  Returns Magnetic Resistivity, parallel and perpendicular components.
!!
!! ARGUMENTS
!!
!!   solnVec  :   solution state, a vector from UNK with all variables
!!   resPar   :   parallel component of Magnetic Resistivity
!!   resPerp :    perpendicular component of Magnetic Resistivity
!!
!!***

#include "Flash.h"
#include "constants.h"  


subroutine MagneticResistivity_fullState(solnVec,resPar, resPerp)
  use MagneticResistivity_interface, ONLY: MagneticResistivity
  use MagneticResistivity_data, ONLY: mag_useMagneticResistivity, &
       res_mele, res_qele, res_navo, res_speedlt, res_boltz, res_hbar
  use MagneticResistivity_data, ONLY: res_mUnit
  use MagneticResistivity_data, ONLY: res_coef
  use Eos_interface, ONLY: Eos_getAbarZbar

  implicit none

  real, intent(in)  :: solnVec(NUNK_VARS)
  real, intent(out) :: resPar
  real, OPTIONAL, intent(out) :: resPerp

  real :: dens
  real :: tele
  real :: nele
  real :: abar
  real :: zbar
  real :: ll
  real :: resPerpLoc


  call Eos_getAbarZbar(solnVec=solnVec,abar=abar,zbar=zbar)
  dens = solnVec(DENS_VAR)
  nele = zbar * dens * res_navo / abar

#ifdef FLASH_3T
  tele = solnVec(TELE_VAR)
#else
  tele = solnVec(TEMP_VAR)
#endif

  call logLambda(tele, nele, zbar, ll)
  resPerpLoc = (1.1*zbar +0.6) * sqrt(res_mele) * res_qele**2 * ll & 
      / (res_boltz * tele)**1.5

  if (present(resPerp)) resPerp = resPerpLoc
  ! This formula is only valid when the magnetic field is strong (and
  ! only for hydrogen):
  resPar = resPerpLoc/1.96

  !! Scale resistivity depending on units
  if (res_mUnit == "SI" .or. res_mUnit == "si" ) then
     resPerpLoc = resPerpLoc*1.e7/(4.*PI)
     resPar  = resPar *1.e7/(4.*PI)
  elseif (res_mUnit == "CGS" .or. res_mUnit == "cgs" ) then
     resPerpLoc = resPerpLoc*res_speedlt**2/(4.*PI)
     resPar  = resPar *res_speedlt**2/(4.*PI)
  else !no unit
     ! Do nothing
     resPerpLoc = resPerpLoc*res_speedlt**2/(4.*PI)
     resPar  = resPar *res_speedlt**2/(4.*PI)
  end if

  if (present(resPerp)) resPerp = resPerpLoc

contains

  subroutine loglambda(tele, nele, zbar, ll)
    implicit none

    ! This subroutine computes the Coulomb logarithm. The formula used
    ! comes from Atzeni.
    
    real, intent(in)  :: tele ! electron temperature [K]
    real, intent(in)  :: nele ! electron number density [cm^-3]
    real, intent(in)  :: zbar ! the average ionization [unitless]
    real, intent(out) :: ll   ! the coulomb logarithm [unitless]

    real :: bmax, bmin, bmin_classic, bmin_quantum
    real, parameter :: ll_floor = 1.0
    
    bmax = sqrt(res_boltz * tele / (4*PI * res_qele**2 * nele))
    
    bmin_classic = zbar * res_qele**2 / (3*res_boltz*tele)
    bmin_quantum = res_hbar / (2*sqrt(3*res_boltz*tele*res_mele))
    bmin = max(bmin_classic, bmin_quantum)

    ll = max(log(bmax/bmin), ll_floor)

  end subroutine loglambda

end subroutine MagneticResistivity_fullState