;+
; NAME:
;  bz1980_vp
; PURPOSE: (one line)
;  Return vapor pressure, optionally latent heat
; DESCRIPTION:
;  Return vapor pressure, optionally latent heat from Brown & Ziegler 1980
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;  vp = bz1980_vp(species, t, latentheat, energyratio, phase)
; INPUTS:
;  species: string 
;    Valid species are: p-H2, O2, N2, F2, CO, CO2, Ne, Ar, Kr, Xe, CH4, C2H4, C2H6
;  t : temperature, K
; OUTPUTS:
;  vp : vapor pressure in microbar
; OPTIONAL OUTPUT
;  latentheat : latent heat of sublimation or vaporization in erg/g
;  energyratio : ratio of latent to kinetic energy, unitless.
;    latent heat per molecule / (k*T)
;  phase : a string representing the phase
;  key: string giving designation in  bz1980_vp_table3.txt
;  ap: the polynomial coefficients A0..A6 in bz1980_vp_table3.txt
; RESTRICTIONS:
;  The files "bz1980_vp_table1.txt" and "bz1980_vp_table3.txt"
;   must be in the same directory as bz1980_vp.pro
; PROCEDURE:
;  This function implements the tables of 
;   Brown, G.N., Ziegler, W.T., 1980.
;   Vapor pressure and heats of vaporization and
;   sublimation of liquids and solids of interest in cryogenics below 1atm
;    pressure. Adv. Cryog. Eng. 25, 662–670.
;   (note there is confustion about the year of this publication, as
;   the conference was 1979 and the procedings were published in
;   1980).
;  Latent heat is calculated from the derivative of the vapor
;  pressure, as John Stansberry pointed out that the BZ latent heat
;  tables are incorrect.
; MODIFICATION HISTORY:
;  Written 2011 Apr 10, by Leslie Young, SwRI
;  2016 Mar 22 Cleaned up documentation, LAY 
;-

function bz1980_vp, species, t, latentheat, energyratio, phase, key, ap

    physconstants
    torr = 1.33322d3  ; microbar per torr
    cal = 4.184e7     ; erg per cal

    scalar = ( size(t,/n_dim) eq 0)

    ; get the directory that has this file, 'bz1980_vp'
    findpro, 'bz1980_vp', dirlist=dirlist, /noprint

    ; get the phase and Polynomial designation
    table1fn = filepath('bz1980_vp_table1.txt',root=dirlist[0])
    if not isfile(table1fn) then begin
       print, 'bz1980_vp: check bz1980_vp_table1.txt is in the same directory as bz1980_vp.pro'
       return, 0.
    endif
    readcol, table1fn, substance, state, t_high, t_low, p_high, p_low, designation, $
             for='A,A,F,F,D,D,A', skip=11, /silent
    p_high = p_high * torr
    p_low = p_low * torr
    if scalar then begin
       ok = where( (species eq substance) and (t_high gt t) and (t_low le t), nok)
       if nok eq 0 then begin
          print, 'bz1980_vp: no data for species ', species, ' at temperature ', t
stop
          return, 0.
       endif
       i = ok[0]
       phase = state[i]
       key = designation[i]
    endif else begin
       nt = n_elements(t)
       phase = strarr(nt)
       key = strarr(nt)
       p = dblarr(nt)
       latentheat = dblarr(nt)
       energyratio = dblarr(nt)
       for it= 0, nt-1 do begin
          ok = where( (species eq substance) and (t_high gt t[it]) and (t_low le t[it]), nok)
          if nok gt 0 then begin
             i = ok[0]
             phase[it] = state[i]
             key[it] = designation[i]         
          endif else begin
             print, 'bz1980_vp: no data for species ', species, ' at temperature ', t[it]
             stop
          endelse
       endfor
    endelse

    ; get the molecular weight
    substance = ['p-H2', 'O2',   'N2',     'F2',    'CO',    'CO2',   'Ne',    'Ar',    'Kr',    'Xe',    'CH4',   'C2H4', 'C2H6']
    molweight = [2.0159, 31.9989, 28.0134, 37.9968, 28.0102, 44.0096, 20.1798, 39.9481, 83.7982, 131.293, 16.0426, 28.0533, 30.0692]
    i = matchfirst(species, substance)
    if i eq -1 then begin
       print, 'bz1980_vp: no molecular weight data for species ', species
stop
       return, 0.
    endif
    mu = molweight[i]

    ; get the phase and Polynomial designation
    table3fn = filepath('bz1980_vp_table3.txt',root=dirlist[0])
    if not isfile(table3fn) then begin
       print, 'bz1980_vp: check bz1980_vp_table3.txt is in the same directory as bz1980_vp.pro'
       return, 0.
    endif
    readcol, table3fn, designation, degree, ap0, ap1, ap2, ap3, ap4, ap5, ap6, $
             for='A,I,D,D,D,D,D,D,D', skip=14, /silent
    if scalar then begin
       ok = where(key eq designation, nok)
       if nok eq 0 then begin
          print, 'bz1980_vp: no pressure coeficients for species ', species, ' with designation', key
stop
          return, 0.
       endif
       i = ok[0]
       ap = [ap0[i], ap1[i], ap2[i], ap3[i], ap4[i], ap5[i], ap6[i]]
       p = exp( total((t^(-findgen(7))) * ap))* torr       
       al = -ap[1:6] *  (findgen(6)+1)
       latentheat = total( t^(- (findgen(6)) ) * al) * !phys.k / (mu * !phys.m_u)       
    endif else begin
       key_uniq = key[uniq(key, sort(key))]
       for ik = 0, n_elements(key_uniq)-1 do begin
          ok = where(key_uniq[ik] eq designation, nok)
          if nok gt 0 then begin
             i = ok[0]
             ap = [ap0[i], ap1[i], ap2[i], ap3[i], ap4[i], ap5[i], ap6[i]]
             al = -ap[1:6] *  (findgen(6)+1)
             k = where(key eq key_uniq[ik], ntk)
             for itk = 0, ntk-1 do begin
                tk = t[k[itk]]
                p[k[itk]] = exp( total(( tk^(-dindgen(7))) * ap))* torr       
                latentheat[k[itk]] = total( tk^(- (dindgen(6)) ) * al) * !phys.k / (mu * !phys.m_u)       
             endfor
          endif else begin
             
          endelse
       endfor
    endelse

    energyratio = latentheat * mu * !phys.m_u / (!phys.k * t)

    return, p

end

pro bz1980_vp_test
; check we get the right pressures at the ends of the polynomial
; ranges

    physconstants
    torr = 1.33322e3            ; microbar per torr
    cal = 4.184e7               ; erg per cal
    
    ; get the directory that has this file, 'bz1980_vp'
    findpro, 'bz1980_vp', dirlist=dirlist, /noprint

    ; get the phase and Polynomial designation
    table1fn = filepath('bz1980_vp_table1.txt',root=dirlist[0])
    readcol, table1fn, substance, state, t_high, t_low, p_high, p_low, designation, $
             for='A,A,F,F,D,D,A', skip=11, /silent
    p_high = p_high * torr
    p_low = p_low * torr

    ;-------------------------------------------------------
    ; Print the calculated pressure ranges
print, $
" Substance    phase             temperature (K)  pressure (torr)      designation Pcalc/Ptab     "+$
"A0              A1               A2                A3              A4                A5               A6"

;    for i = 0, n_elements(designation) do begin
    for i = 0,-1 do begin
;       pp = bz1980_vp(substance[i], [t_high[i] - 1d-6, t_low[i] + 1d-6], l, a, phaselow, k,coef)
       p_highi = bz1980_vp(substance[i], t_high[i] - 1d-4, l, a, phasehigh, k,coef)
       p_lowi = bz1980_vp(substance[i], t_low[i] + 1d-4, l, a, phaselow, k,coef)
       pp = [p_highi, p_lowi]
       if phasehigh ne phaselow then stop
       print, substance[i], phasehigh, t_high[i], t_low[i], pp/torr, designation[i], pp/[p_high[i],p_low[i]], coef, $
              for='(" ", A-13, A-15, 2F9.3, 2E11.3, " ", A-8, 2F8.3, 7E17.8)'
;       print, substance[i], t_high[i], t_low[i], '  '+phasehigh, [ p_high[i], pp[0], p_highi, p_low[i], pp[1], p_lowi ] / torr
;print, coef, for='(7E17.8)'
    endfor

    ;-------------------------------------------------------
    ; Print the calculated latent heats
    table3fn = filepath('bz1980_vp_table3.txt',root=dirlist[0])
    readcol, table3fn, designation3, degree3, ap0, ap1, ap2, ap3, ap4, ap5, ap6, $
             for='A,I,D,D,D,D,D,D,D', skip=14, /silent
    table5fn = filepath('bz1980_vp_table5.txt',root=dirlist[0])
    readcol, table5fn, designation5, degree5, bl0, bl1, bl2, bl3, bl4, bl5, $
             for='A,I,D,D,D,D,D,D', skip=14, /silent
    subst_mw = ['p-H2', 'O2',   'N2',     'F2',    'CO',    'CO2',   'Ne',    'Ar',    'Kr',    'Xe',    'CH4',   'C2H4', 'C2H6']
    molweight = [2.0159, 31.9989, 28.0134, 37.9968, 28.0102, 44.0096, 20.1798, 39.9481, 83.7982, 131.293, 16.0426, 28.0533, 30.0692]
    for i = 0,30 do begin
       key = designation[i]
       im = matchfirst(substance[i], subst_mw)
       mu = molweight[im]
       thi = t_high[i]
       tlo = t_low[i]
       p_highi = bz1980_vp(substance[i], t_high[i] - 1d-4, lhi, a, phasehigh, k,coef)
       p_lowi = bz1980_vp(substance[i], t_low[i] + 1d-4, llo, a, phaselow, k,coef)
       ok = where(key eq designation3, nok) & i3 = ok[0]
       ap = [ap0[i3], ap1[i3], ap2[i3], ap3[i3], ap4[i3], ap5[i3], ap6[i3]]
       al = -ap[1:6] *  (findgen(6)+1)
       l3 = [total( thi^(- (findgen(6)) ) * al),total( tlo^(- (findgen(6)) ) * al)] * !phys.k / (mu* !phys.m_u)  
       ok = where(key eq designation5, nok)
       if nok gt 0 then begin
          i5 = ok[0]
          bl = [bl0[i5], bl1[i5], bl2[i5], bl3[i5], bl4[i5], bl5[i5]]
          l5 = [total( thi^((findgen(6)) ) * bl), total( thi^((findgen(6)) ) * bl)] * cal / mu
          print, substance[i], phasehigh, t_high[i], t_low[i], lhi, llo, k, l5/[lhi,llo], bl, $
                 for='(" ", A-13, A-15, 2F9.3, 2E11.3," ", A-8, 2F8.3, 6E17.8)'
       endif
    endfor

end