Function filename,v
; Returns the appropriate data file name for the OH J and wavenumber
; data for the specified rotational level.
CASE 1 OF
v LE 5: Return,'c:\data\spectra\c_oh2p'+StrTrim(v,2)+'.dat'
v GT 5: Return,'c:\data\spectra\cf_oh2p'+StrTrim(v,2)+'.dat'
  ELSE: BEGIN
           Message,'Bad vibrational level: '+StrTrim(v,2),/Informational
           Return,'NUL'
           END
    ENDCASE
END

Function hl_r,jp,lp
; Returns an R-branch (Honl-London) line strength (delta lambda = 0)
Return,(jp+lp)*(jp-lp)/(Float(jp))
END

Function hl_q,jp,lp
; Returns an Q-branch (Honl-London) line strength (delta lambda = 0)
Return,(2.*jp+1)*lp*lp/(Float(jp)*(jp+1))
END

Function hl_p,jp,lp
; Returns an P-branch (Honl-London) line strength (delta lambda = 0)
Return,(jp+1+lp)*(jp+1-lp)/(Float(jp+1))
END

;+
; NAME:
;	MkOHFile
; PURPOSE:
;	Generate input data for use in an OH Meinel synthetic spectrum
;       generation program of the type provided by Abas Sivjee.
; CATEGORY:
;	Synthetic spectra.
; CALLING SEQUENCE:
;	MkOHFile, vp, vpp
; INPUTS:
;	vp:     Vibrational quantum number of the upper vibrational
;               level (v');
;       vpp:    Vibrational quantum number of the lower vibrational
;               level (v'').
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;	None.
; OUTPUTS:
;	None.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	Two files of OH rotational-vibrational energy levels are read
;       from C:\DATA\SPECTRA, to get the needed data for the specified
;       levels. A data file containing the generated data is produced
;       in C:\IDL\SYNSPEC.
; RESTRICTIONS:
;	The input data files for the specified vibrational levels must
;       exist. The data files so far produced cover only rotational
;       levels J = 0.5 - 9.5.
; PROCEDURE:
;	Relatively straightforward. The term values for the f1 and f2 
;       levels are used to generate energy differences (in wavenumbers).
;       The Honl-London formulae for delta lambda = 0 are used to
;       approximate the line strengths. The energy levels for use in
;       the Boltzmann distribution are generated by subtracting the
;       lowest term value in the given vibrational level (assumed to be
;       that of f1e(J=1.5)). The wavenumber differences, line strengths,
;       and energy level differences (in cm-1) are merged over the R-,
;       Q-, and P-branches, the vacuum wavelengths are calculated, and
;       the arrays are sorted in order of increasing wavelength. The
;       wavelengths are converted to their values in air. The
;       wavelengths, line strengths, and energy level differences are
;       then written out to a disk file.
; EXAMPLE:
;	mkohfile,8,3    ; should replicate C:\IDL\SYNSPEC\OH83WSF.CHC
; SEE ALSO:
;	oh83.pro
; MODIFICATION HISTORY:
; 	Written by:	DP Steele, March 10, 1997.
;-

PRO mkohfile,vp,vpp

; Check validity of inputs
IF N_Params() NE 2 THEN BEGIN
    Doc_Library,'mkohfile'
    Return
    ENDIF

IF vp LE vpp THEN BEGIN
    Message,"v' must be greater than v''",/Informational
    Return
    ENDIF

; Load the needed data
vpdata=DDRead(filename(vp))
jp=Reform(vpdata(0,*))
njp=N_Elements(jp)
f1ejp=Reform(vpdata(1,*))
f1fjp=Reform(vpdata(2,*))
f2ejp=Reform(vpdata(3,*))
f2fjp=Reform(vpdata(4,*))

vppdata=DDRead(filename(vpp))
jpp=Reform(vppdata(0,*))
njpp=N_Elements(jpp)
f1ejpp=Reform(vppdata(1,*))
f1fjpp=Reform(vppdata(2,*))
f2ejpp=Reform(vppdata(3,*))
f2fjpp=Reform(vppdata(4,*))

nj=njp<njpp

; R branch
delj=1
nu1r=DblArr(nj)
nu2r=nu1r
shl1r=nu1r
shl2r=nu1r
f1r=nu1r
f2r=nu1r
FOR ijp=1,nj-1 DO BEGIN
    IF ijp GT 1 THEN BEGIN
        nu1er=f1ejp(ijp)-f1ejpp(ijp-delj)
        nu1fr=f1fjp(ijp)-f1fjpp(ijp-delj)
        nu1r(ijp)=(nu1er+nu1fr)/2.      ; average two values
        shl1r(ijp)=hl_r(jp(ijp),1.5)    ; Honl-London line strength
        f1r(ijp)=(f1ejp(ijp)+f1fjp(ijp))/2.-f1ejp(1)    ; cm-1 above ref. level
        ENDIF
    IF ijp LT nj-1 THEN BEGIN
        nu2er=f2ejp(ijp)-f2ejpp(ijp-delj)
        nu2fr=f2fjp(ijp)-f2fjpp(ijp-delj)
        nu2r(ijp)=(nu2er+nu2fr)/2.
        shl2r(ijp)=hl_r(jp(ijp),0.5)
        f2r(ijp)=(f2ejp(ijp)+f2fjp(ijp))/2.-f1ejp(1)
        ENDIF
    ENDFOR

; Q branch
delj=0
nu1q=DblArr(nj)
nu2q=nu1q
shl1q=nu1q
shl2q=nu1q
f1q=nu1q
f2q=nu1q
FOR ijp=0,nj-1 DO BEGIN
    IF ijp GT 0 THEN BEGIN
        nu1eq=f1ejp(ijp)-f1ejpp(ijp-delj)
        nu1fq=f1fjp(ijp)-f1fjpp(ijp-delj)
        nu1q(ijp)=(nu1eq+nu1fq)/2.
        shl1q(ijp)=hl_q(jp(ijp),1.5)
        f1q(ijp)=(f1ejp(ijp)+f1fjp(ijp))/2.-f1ejp(1)
        ENDIF
    IF ijp LT nj-1 THEN BEGIN
        nu2eq=f2ejp(ijp)-f2ejpp(ijp-delj)
        nu2fq=f2fjp(ijp)-f2fjpp(ijp-delj)
        nu2q(ijp)=(nu2eq+nu2fq)/2.
        shl2q(ijp)=hl_q(jp(ijp),0.5)
        f2q(ijp)=(f2ejp(ijp)+f2fjp(ijp))/2.-f1ejp(1)
        ENDIF
    ENDFOR

; P branch
delj=-1
nu1p=DblArr(nj)
nu2p=nu1p
shl1p=nu1p
shl2p=nu1p
f1p=nu1p
f2p=nu1p
FOR ijp=0,nj-2 DO BEGIN
    IF ijp GT 0 THEN BEGIN
        nu1ep=f1ejp(ijp)-f1ejpp(ijp-delj)
        nu1fp=f1fjp(ijp)-f1fjpp(ijp-delj)
        nu1p(ijp)=(nu1ep+nu1fp)/2.
        shl1p(ijp)=hl_p(jp(ijp),1.5)
        f1p(ijp)=(f1ejp(ijp)+f1fjp(ijp))/2.-f1ejp(1)
        ENDIF
    IF ijp LT nj-2 THEN BEGIN
        nu2ep=f2ejp(ijp)-f2ejpp(ijp-delj)
        nu2fp=f2fjp(ijp)-f2fjpp(ijp-delj)
        nu2p(ijp)=(nu2ep+nu2fp)/2.
        shl2p(ijp)=hl_p(jp(ijp),0.5)
        f2p(ijp)=(f2ejp(ijp)+f2fjp(ijp))/2.-f1ejp(1)
        ENDIF
    ENDFOR

; Merge the results 
nu=[nu1r,nu2r,nu1q,nu2q,nu1p,nu2p]
shl=[shl1r,shl2r,shl1q,shl2q,shl1p,shl2p]
f=[f1r,f2r,f1q,f2q,f1p,f2p]

; Eliminate null wavenumber entries.
nz=Where(nu GT 0.)
nu=nu(nz)
shl=shl(nz)
f=f(nz)

; Sort them into wavelength order
w=1.0E8/nu
sw=Sort(w)
ws=w(sw)
nus=nu(sw)
shls=shl(sw)
fs=f(sw)

; Convert vacuum wavelengths to air.
VacToAir,ws
Stop

; Write the data to a file for inspection.
outname=String(vp,vpp,Format='("oh",2I1,"wsf.dat")')
OpenW,lun,'c:\idl\synspec\'+outname,/Get_LUN
ofmt='(F7.2,F7.3,F7.1)'
FOR i=0,N_Elements(w)-1 DO $
    PrintF,lun,ws(i),shls(i),fs(i),Format=ofmt
Free_LUN,lun

Return
END