;pro spectrademo

;+
; NAME:
;  spectrademo
; PURPOSE: (one line)
;   Demonstrate the use of the spectra library
; DESCRIPTION:
;   A heavily-commented proceedure to show how to use the spectral library
; CATEGORY:
;  Spectra
; CALLING SEQUENCE:
;  spectrademo
; INPUTS:
;  none
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  None.
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
;  Prints and plots
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2000 October, by Leslie Young, SwRI
;-

on_error, 2

; Load the physical constants
physconstants

; Load the spectralib variables.
spectralib

; Set the CH4 filename and open it
fn = !spectralib.by_molecdir+'06_HIT96.PAR'
openr, fp, fn, /get_lun

; read all the lines from 2 to 4.2 micron
numin = 1.e4/4.2 
numax = 1.e4/2.0
nmax = 10000
ln=new_specline()
lna=replicate(ln, nmax)
n=0L
rd_next, fp, ln, 6, 1, -1, -1,numin, numax, -1
while ((ln.mol gt 0) and (n lt nmax) ) do begin &$
  lna(n)=ln &$
  n=n+1 &$
  rd_next, fp, ln, 6, 1, -1, -1,numin, numax, -1 &$
end
close, fp

lnb = replicate(ln, n)
lnb(0:n-1) = lna(0:n-1)
delvar, lna

; Change the temperature and pressure to that of the
; Jovian lower thermosphere
t = 160
p = 1e4 ; 10 mbar
n = p/(!phys.k*t)
h = 24.45e5 ; scale height
ncol = sqrt(2.*!pi*71000.e5*h) * n

changept_specline_array, lnb, p, t

; make a spectrum
dnu = 0.0002  ; 1/5 pressure width
nu0 = 2384.9
nu1 = 4666.6
nnu = (nu1-nu0)/dnu
nu = nu0 + findgen(nnu)*dnu
tau = 0*nu
for i=3L, 6334-3-3 do begin &$
  j = (lnb(i).nu0 - nu0)/dnu &$
  tau(j-5000:j+5000) = tau(j-5000:j+5000) + $
      (lnb(i).s * ncol/(!pi*lnb(i).gamma_air)) * 1/( ((nu(j-5000:j+5000)-lnb(i).nu0)/lnb(i).gamma_air)^2 + 1) &$
end

;trans = exp(-tau)
;
;nusam = rebin( nu(0:11400000-1), 11400/2)
;transsam = rebin( trans(0:11400000-1), 11400/2)
;
;plot, 1e4/nusam, transsam

nnu = n_elements(nu)

; f = factor to use to adjust for 
f = 1.8e-3 * 1.5  ; CH4 fraction, and half-light's at 15 mbar

tau = tau * f

r = 73000.e5
h = 23.e5
v = 17e5
phi=[0.01,0.02,0.03,0.04,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,.99]
x = -alog( (1/phi-1) )
rho=r+x*h-h*exp(-x)
time = (rho-(r-h))/v
plot, time, phi, ps=-4

nx = n_elements(x)

ns = 1140L
tarr = fltarr(ns, nx)
for ix = 0,nx-1 do begin &$
   for is=0L,ns-1 do begin &$
     trans = exp(-exp(-x(ix))*tau(is*10000:is*10000+9999)) &$
     tarr(is,ix) = total(trans)/10000. &$
   endfor &$
endfor
lam=fltarr(ns)
for is=0L,ns-1 do lam(is) = total(nu(is*10000:is*10000+9999))/10000
lam = 1e4/lam

plot, lam, tarr(*,9), xr=[2,2.5], yr=[-.1,1.1]
plot, 1e4/smooth(lam,5), smooth(tarr(*,9),5), xr=[4000,5000]
oplot, 1e4/smooth(lam,5), smooth(tarr(*,4),5), /li
oplot, 1e4/smooth(lam,5), smooth(tarr(*,14),5), /li

plot, time, phi, xtit=lam(800), yr=[-.1,1.1], xr=[-50,10], /xs
oplot, time, tarr(800,*)*phi, /li
oplot, time, 0*time, li=2

save, file='ch4lines3.sav'

gd = where(lam gt 3.2 and lam lt 4)
plot, time, total(tarr(gd,*),1) * phi/312, xr=[-50,10],/xs, /li
oplot, time+3, phi
oplot, time, phi
s = sort(tarr(gd,9))

gd = where(lam gt 3.2 and lam lt 4)


; interpolate over time
ni = 1070
timei = findgen(107/0.1)*0.1-100.
rhoi = v*timei + (r-h)
xi = interpol(x, rho,rhoi)
phii = 1/(1+exp(-xi))
tarri = fltarr(ns, ni)
for si=0,ns-1 do tarri(si,*) = interpol(tarr(si,*), x,xi)
lci = tarri
for i=0,ni-1 do lci(*,i) = tarri(*,i) * phii(i)

plot, smooth(lam,5), smooth(total(tarri(*,1000:1009),2)/10,5), xr=[2.35,2.5]
for i=0,1050, 20 do $
  oplot, smooth(lam,5), smooth(total(tarri(*,i:i+10),2)/10,5), color=220-abs((i-1000))/2

; make a 2-sec lightcurve
t_2sec = fltarr(54)
for i=0,1060, 20 do t_2sec(i/20) = total(timei(i:i+9))/10
lc900 = reform(tarri(900,*))*phii
lc900_2sec = fltarr(54)
for i=0,1060, 20 do lc900_2sec(i/20) = total(lc900(i:i+9))/10
lc850 = reform(tarri(850,*))*phii
lc850_2sec = fltarr(54)
for i=0,1060, 20 do lc850_2sec(i/20) = total(lc850(i:i+9))/10

plot, t_2sec, lc900_2sec, ps=-4, xr=[-20,10]
oplot, [-10],[0.5],ps=4 
xyouts, -10,.5,lam(900)
oplot, t_2sec, lc850_2sec, ps=-1
oplot, [-10],[0.4],ps=1 & xyouts, -10,.4,lam(850)

plot, t_2sec, lc900_2sec, ps=-4, xr=[-20,10]
oplot, [-10],[0.5],ps=4 
xyouts, -10,.5,lam(900)
oplot, t_2sec, lc850_2sec, ps=-1
oplot, [-10],[0.4],ps=1 & xyouts, -10,.4,lam(850)

t_3sec = fltarr(36)
for i=0,1050, 30 do t_3sec(i/30) = total(timei(i:i+19))/20
lc900_3sec = fltarr(36)
for i=0,1050, 30 do lc900_3sec(i/30) = total(lc900(i:i+19))/20
lc850_3sec = fltarr(36)
for i=0,1050, 30 do lc850_3sec(i/30) = total(lc850(i:i+19))/20

plot, t_3sec, lc900_3sec, ps=-4, xr=[-20,10]
oplot, [-10],[0.5],ps=4 
xyouts, -10,.5,lam(900)
oplot, t_3sec, lc850_3sec, ps=-1
oplot, [-10],[0.4],ps=1 & xyouts, -10,.4,lam(850)

plot, smooth(lam,5), smooth(total(lci(*,990:990+19),2)/20,5), xr=[2.35,2.6], th=2, yr=[-.1,1.1], /ys
for i=900,1050, 30 do $
   oplot, smooth(lam,5), smooth(total(lci(*,i:i+19),2)/20,5)
i=0
oplot, smooth(lam,5), smooth(total(lci(*,i:i+19),2)/20,5)


;end