set_ceq_from_chease Subroutine

  pure subroutine set_ceq_from_chease(self)
    use constants, only: pi
    use read_chease, only: npsi_chease, nchi_chease, b0exp_chease
    use read_chease, only: psi_chease, f_chease, q_chease, p_chease
    use read_chease, only: r_chease, z_chease
    implicit none
    class(ceq_type), intent(in out) :: self
    real :: R_geo
    integer :: i, j

    !     nr == number of actual grid points in radial array
    self%nr = npsi_chease

    !     nt == number of theta grid points in theta eq grid
    self%nt = nchi_chease + 1

    call self%alloc_arrays()

    !     eqpsi(1:nr) == values of psi on the radial grid
    self%eqpsi = psi_chease

    !     fp(1:nr) == the function that satisfies
    !              B = fp grad zeta + grad zeta x grad psi
    self%fp = f_chease
    self%qsf = q_chease

    !     pressure(1:nr) == pressure profile on the radial grid,
    !     normalized to the value at the magnetic axis.
    self%pressure = p_chease

    !     psi_0 == value of psi at the magnetic axis
    !     psi_a == value of psi at the boundary of the plasma
    self%psi_0 = self%eqpsi(1) ; self%psi_a = self%eqpsi(self%nr)
    self%psi_bar = (self%eqpsi - self%psi_0) / (self%psi_a - self%psi_0)

    do j = 1, self%nt
      do i = 1, self%nr
        self%R_psi(i, j) = r_chease(i, chease_chi_index(nchi_chease, j))
      end do
    end do

    do j = 1, self%nt
      do i = 1, self%nr
        self%Z_psi(i,j) = z_chease(i, chease_chi_index(nchi_chease, j))
      end do
    end do

    !     aminor    == half diameter of last flux surface at elevation of Z_mag
    R_geo = (self%R_psi(self%nr, self%nt/2+1) + self%R_psi(self%nr, 1))/2.
    self%aminor = (self%R_psi(self%nr, self%nt/2+1) - self%R_psi(self%nr, 1))/2.
    self%B_T = self%fp(self%nr)/R_geo
    self%R_psi = self%R_psi / self%aminor
    self%Z_psi = self%Z_psi / self%aminor
    self%B_psi = 0.0

    !     R_mag == R position of magnetic axis
    !     Z_mag == Z position of magnetic axis
    self%R_mag = self%R_psi(1,1)
    self%Z_mag = self%Z_psi(1,1)
    self%beta = 8. * pi * 1.e-7 * self%pressure / b0exp_chease**2  !not correct definition yet
    self%beta_0 = self%beta(1)
    self%pressure = self%pressure / self%pressure(1)
  end subroutine set_ceq_from_chease