PRO neutdens,hkm,ndens,nlayers=nlayers,del_h=del_h,debug=debug

; This little routine returns neutral density (MKS units) vs. geometric
; height over the height range 0 - 50 km.  The data are taken from the
; CIRA-86 grand mean (which is nominally for 30 deg N).
;
; DPSteele, Nov. 6, 1996

IF ((nlayers-0.5)*del_h) GT 51. THEN BEGIN
    Message,'Density data available only up to 51 km at present' $
           ,/Informational
    Return
    ENDIF

IF Keyword_Set(debug) THEN Stop

k=1.380658D-23  ; Boltzmann's constant

hgeom=[57.,2166.,4197.,6148.,8011.,9771.,11450.,13055. $        ; meters
      ,14604.,16116.,17625.,19151.,20697.,22267.,23860.,25478. $
      ,27120.,28788.,30479.,32198.,33948.,35737.,37563.,39425. $
      ,41325.,43265.,45231.,47213.,49198.,51170.]

p_mb=[1013.,788.93,614.42,478.51,372.66,290.23,226.03,176.03 $  ; mbar
     ,137.09,106.77,83.15,64.76,50.43,39.28,30.59,23.82 $
     ,18.55,14.45,11.25,8.7642,6.8255,5.3157,4.1399,3.2242 $
     ,2.5110,1.9556,1.5230,1.1861,0.9237,0.7194]

temp=[295.1,283.1,272.1,259.9,247.0,234.6,223.6,214.8 $         ; Kelvin
     ,208.5,205.8,207.1,210.0,212.8,215.6,218.3,221.0 $
     ,223.9,226.9,230.3,234.4,238.9,243.7,248.7,253.8 $
     ,258.8,263.0,266.0,267.1,266.5,263.5]

p_Pa=100.*p_mb          ; Pascals
ngeom=p_Pa/(k*temp)     ; m**-3

IF N_Elements(nlayers) EQ 0 THEN nlayers=25
IF N_Elements(del_h) EQ 0 THEN del_h=2.

hkm=Findgen(nlayers)*del_h+del_h/2.
ndens=Interpol(ngeom,hgeom,hkm*1000.)

Return
END

PRO ozpaths,debug=debug,ps=ps
; PATHS.PRO - Compute path lengths for NITEOWL s/c viewing
; RLG Oct/96
; Cosmetic changes added by D P Steele, October 28/96
; Major changes made by D P Steele, November 6, 1996

nh=25   ;  Assume 25 layers, each 2 km thick
ozden=DblArr(nh)
tpath=DblArr(nh)
dpath=DblArr(nh)
ozpath=DblArr(nh)
ozpath(*)=0.0
ozpaths=ozpath
neutpath=ozpath
re=6371.2D0  ; Earth radius, km
; centrht=Dindgen(nh)*2 +1
yn=' '
ozden(*)=[1.0,1.2,1.3,1.0,.5,.3,.3,.5,1.3,2.3,3.8,4.4,4.8, $
          4.5,3.7,3.0,2.2,1.7,1.2,.8,.4,.3,.2,.15,.1]*1.0D12    ; cm**-3

neutdens,centrht,neutden,nlayers=nh,del_h=2.,debug=debug
neutden=neutden/1.0E6   ; convert from m**-3 to cm**-3

CASE 1 OF
N_Elements(ps): psopen,'c:\usask\893\ozpaths.ps'
          ELSE: set_plot,'win'         ; Line added 961028 by D P Steele
    ENDCASE

; openw,1,'ozpath.txt'

; printf,1, ' Ozone Densities [10(12)] at 1,3,5, .... 49 km '
; printf,1,ozden
ozdens=ozden
ozdens(0:5)=0.0  ; Assume no tropospheric O3 to do comparison.

for th=0,24 do begin   ;  Compute [O3] weighted path for 25 layers
 curpath=0.0   ; Start with zero length at tangent point
 dpath(*)=0.0   &   tpath(*)=0.0  ; Clean up
 for h=th,24 do begin    ; Scan from current tangent upwards to top
  tpath(h)=sqrt((re+(h+1.D0)*2.D0)^2 - (re+th*2.D0)^2)*1.0D5    ; cm
  dpath(h)=tpath(h)-curpath  ; Store layer geometric length
  curpath=tpath(h)
  ozpath(th)=ozpath(th)+ozden(h)*dpath(h)     ; Strato and tropo O3
  ozpaths(th)=ozpaths(th)+ozdens(h)*dpath(h)  ; Strato O3 only
  neutpath(th)=neutpath(th)+neutden(h)*dpath(h) ; Neutral atmosphere
 endfor
; if th eq 0 then begin   ; A check print
;  printf,1,' Path Lengths (Km) and Deltas for Tangent = ',th*2+1,' Km'
;  printf,1,tpath
;  printf,1,dpath
; endif
endfor

; printf,1,' Column O3 [cm-2] vs Tangent Height With and Without Tropo [O3]'
; printf,1,ozpath
; printf,1,ozpaths

plot,ozpath,centrht,linestyle=0, $
title=' Column O3 Concentration (cm!e-2!n) vs Tangent Height', $
xtitle='Column Ozone Amount (cm!e-2!n)', $
ytitle='Tangent Height (km)'
oplot,ozpaths,centrht,linestyle=2

read,'Press <Return> to see the assumed ozone profiles.',yn

plot,ozden,centrht,linestyle=0,title=' Assumed Ozone Profiles', $
xtitle='Ozone Density (cm!e-3!n)', $
ytitle='Height (km)'
oplot,ozdens,centrht,linestyle=2

; close,1

Read,'Show absorption/scattering effects on a stellar spectrum? ([Y]/N) ' $
     ,yn
doabs=(StrUpCase(yn) NE 'N')
IF doabs THEN BEGIN
    Read,'Magnitude of stellar spectrum? ',starmag
    factor=10.D0^((2.486593-starmag)*0.4)
;    Read,'Altitude of tangent ray? (1..49 km) ',alt
    FOR alt=10,50,10 DO BEGIN
    ozcol=(Interpol(ozpath,centrht,alt))(0)
    neutcol=(Interpol(neutpath,centrht,alt))(0)

    IF alt EQ 10 THEN BEGIN
    Read,'Ozone absorption? ([Y]/N) ',yn
    ozabs=(StrUpCase(yn) NE 'N')
    Read,'Rayleigh scattering? ([Y]/N) ',yn
    raysc=(StrUpCase(yn) NE 'N')

    IF Keyword_Set(debug) THEN Stop

    ; Load star spectrum for G-type star
    OpenR,lun,'\usask\893\crossect\wavelgth.txt',/Get_LUN
    starwave=DblArr(1084)
    ReadF,lun,starwave
    Free_LUN,lun
    dwl=Median(starwave-Shift(starwave,1))
;    starflux=ddread('\usask\893\crossect\g1--2v.txt)
    OpenR,lun,'\usask\893\crossect\g1--2v.txt',/Get_LUN
    starflux=DblArr(1084)
    ReadF,lun,starflux
    Free_LUN,lun
    starflux=starflux*factor

    ; Load ozone absorption cross-section from Shettle et al., GRL
    IF ozabs THEN rdcols,'\usask\893\crossect\ozonecro.dat',no3,o3wl,o3xs
    o3wl=[Findgen(70)+4000.,o3wl]		; Extrapolate xsection to 4000 A.
    o3xs=[Replicate(2.0E-24,70),o3xs]
    no3=no3+70

    ; Load Rayleigh scattering cross-section from Jayson Eppler (MODTRAN)
    IF raysc THEN BEGIN
        rdcols,'\usask\893\crossect\ray2xsec.txt',nrs,rswl,rsxs
        rswl=10.*rswl   ; convert to Angstroms
        ENDIF

    ; Stay within wavelength range covered by all data
    ; wllow=Min(starwave)>Min(o3wl)>Min(rswl)
    wllow=4000.
    ; wlhigh=Max(starwave)<Max(o3wl)<Max(rswl)
    wlhigh=8000.
    ok=Where((wllow LE starwave) and (starwave LE wlhigh),nflux)
    starwave=starwave(ok)
    starflux=starflux(ok)

    ; Interpolate O3 and Rayleigh cross-sections onto star wavelengths
    IF ozabs THEN ozcross=Interpol(o3xs,o3wl,starwave)
    IF raysc THEN rscross=Interpol(rsxs,rswl,starwave)

    ENDIF

    ; Determine appropriate optical depths
    tau=Replicate(0.D0,nflux)
    IF Keyword_Set(debug) THEN Stop
    IF ozabs THEN tau=tau+ozcol*ozcross
    IF raysc THEN tau=tau+neutcol*rscross

    ; Determine effect on star spectrum
    transflux=starflux*Exp(-tau)

    ; Integrated the spectra over the required bin size
    nsamples=Round(50./dwl)     ; defaults to 5 nm per pixel
    sfluxperpix=Total(Reform(starflux,nsamples,nflux/nsamples),1)
    tfluxperpix=Total(Reform(transflux,nsamples,nflux/nsamples),1)
    binwave=Total(Reform(starwave,nsamples,nflux/nsamples),1)/nsamples

    ; Display results
    IF (alt EQ 10) THEN Plot,binwave/10,sfluxperpix,/nodata $
        ,XTitle='!17Wavelength (nm)',XRange=[400,800] $
        ,YTITLE='!17Spectral Flux (ph s!e-1!n pixel!e-1!n)' $
        ,PSym=10,Title='!17Horizontal Stellar Irradiance !18versus!17 Altitude'
    OPlot,binwave/10,tfluxperpix,linestyle=5-(alt/10);,PSym=10
    ENDFOR
    ENDIF

    legend,String((5-Indgen(5))*10,Format='("!17",I2," km")') $
          ,LineStyle=Indgen(5),Position=[0.72,0.91],Box=0
    
IF N_Elements(ps) EQ 1 THEN psclose         ; Line added 961028 by D P Steele

Return
END