; Set up plot

SET_PLOT, 'ps'

DEVICE, XSIZE=16, YSIZE=16, /ENCAPSULATED

; Constants

@constants.idl

; Solve the Saha equation over a grid of T and n

T = 1D4 + DINDGEN(1000)/999*(2e5-1e4)

rho = 1D-4
n = rho/(4.*U_ATOM)
helium_saha, n, T, z_1, z_2, z_0, z_e

!P.MULTI = [0,1,2]

DEVICE, FILENAME='2011-03-23-fig-saha-4.eps'

PLOT, T, z_0, $
      XSTYLE=1, YSTYLE=1, $
      YRANGE=[0,1.05], $
      XTITLE='T (K)', YTITLE='z!Di!N!X', TITLE='rho = 1E-4 g cm!E-3!N!X'

OPLOT, T, z_1, THICK=3, LINESTYLE=1
OPLOT, T, z_2, THICK=3, LINESTYLE=2

PLOT, T, z_e, $
      XSTYLE=1, YSTYLE=1, $
      YRANGE=[0,2.05], $
      XTITLE='T (K)', YTITLE='z!De!N!X', TITLE='rho = 1E-4 g cm!E-3!N!X'

DEVICE, /CLOSE

PRINT, 'rho :',rho
PRINT, 'Half ionization (i) :', T[MIN(WHERE(z_1 gt 0.5))]
PRINT, 'Half ionization (ii):', T[MIN(WHERE(z_2 gt 0.5))]

;

DEVICE, FILENAME='2011-03-23-fig-saha-6.eps'

rho = 1D-6
n = rho/(4.*U_ATOM)
helium_saha, n, T, z_1, z_2, z_0, z_e

!P.MULTI = [0,1,2]

PLOT, T, z_0, $
      XSTYLE=1, YSTYLE=1, $
      YRANGE=[0,1.05], $
      XTITLE='T (K)', YTITLE='z!Di!N!X', TITLE='rho = 1E-6 g cm!E-3!N!X'

OPLOT, T, z_1, THICK=3, LINESTYLE=1
OPLOT, T, z_2, THICK=3, LINESTYLE=2

PLOT, T, z_e, $
      XSTYLE=1, YSTYLE=1, $
      YRANGE=[0,2.05], $
      XTITLE='T (K)', YTITLE='z!De!N!X', TITLE='rho = 1E-6 g cm!E-3!N!X'

DEVICE, /CLOSE

PRINT, 'rho :',rho
PRINT, 'Half ionization (i) :', T[MIN(WHERE(z_1 gt 0.5))]
PRINT, 'Half ionization (ii):', T[MIN(WHERE(z_2 gt 0.5))]

;

DEVICE, FILENAME='2011-03-23-fig-saha-8.eps'

rho = 1D-8
n = rho/(4.*U_ATOM)
helium_saha, n, T, z_1, z_2, z_0, z_e

!P.MULTI = [0,1,2]

PLOT, T, z_0, $
      XSTYLE=1, YSTYLE=1, $
      YRANGE=[0,1.05], $
      XTITLE='T (K)', YTITLE='z!Di!N!X', TITLE='rho = 1E-8 g cm!E-3!N!X'

OPLOT, T, z_1, THICK=3, LINESTYLE=1
OPLOT, T, z_2, THICK=3, LINESTYLE=2

PLOT, T, z_e, $
      XSTYLE=1, YSTYLE=1, $
      YRANGE=[0,2.05], $
      XTITLE='T (K)', YTITLE='z!De!N!X', TITLE='rho = 1E-8 g cm!E-3!N!X'

DEVICE, /CLOSE

PRINT, 'rho :',rho
PRINT, 'Half ionization (i) :', T[MIN(WHERE(z_1 gt 0.5))]
PRINT, 'Half ionization (ii):', T[MIN(WHERE(z_2 gt 0.5))]

;

DEVICE, FILENAME='2011-03-23-fig-saha-10.eps'

rho = 1D-10
n = rho/(4.*U_ATOM)
helium_saha, n, T, z_1, z_2, z_0, z_e

!P.MULTI = [0,1,2]

PLOT, T, z_0, $
      XSTYLE=1, YSTYLE=1, $
      YRANGE=[0,1.05], $
      XTITLE='T (K)', YTITLE='z!Di!N!X', TITLE='rho = 1E-10 g cm!E-3!N!X'

OPLOT, T, z_1, THICK=3, LINESTYLE=1
OPLOT, T, z_2, THICK=3, LINESTYLE=2

PLOT, T, z_e, $
      XSTYLE=1, YSTYLE=1, $
      YRANGE=[0,2.05], $
      XTITLE='T (K)', YTITLE='z!De!N!X', TITLE='rho = 1E-10 g cm!E-3!N!X'

DEVICE, /CLOSE

PRINT, 'rho :',rho
PRINT, 'Half ionization (i) :', T[MIN(WHERE(z_1 gt 0.5))]
PRINT, 'Half ionization (ii):', T[MIN(WHERE(z_2 gt 0.5))]