pro geom_tests, j, k, ampl, fn, xL,nL,pL,tL,refrL,phaseL, angleL, dangleL, yL, fluxL

;+
; NAME:
;  geom_tests
; PURPOSE: (one line)
;   Make files for testing the wavelet code, geometric optics
; DESCRIPTION:
;   
;
; CATEGORY:
;  wavelet01 project
; CALLING SEQUENCE:
;  geom_tests, j, k, ampl, fn, zL,nL,pL,tL,refrL,phaseL,angleL, dangleL, yL, fluxL
; INPUTS:
;  j       = binary dilation parameter 
;  k       = binary translation parameter
;  ampl    = amplitude of wavelet
;  fn      = filename
; OPTIONAL INPUT PARAMETERS:
;  None.
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  xL      = array of planet radii (cm)
;  nL      = array of number densities (molecule cm^-3)
;  pL      = array of pressures (microbar)
;  tL      = array of temperatures (K)
;  refrL   = array of refractivities (unitless)
;  phaseL  = array of phase delays (unitless)
;  angleL  = array of bending angle (radians)
;  dangleL = array of derivative of bending angle (radians/cm)
;  yL      = array of shadow radii (cm)
;  fluxL   = array of fluxes (unitless)
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
;  Runs astrolib and physconstants, loading Astronomy Library & !phys system variables
;  Plots t vs altitude and flux vs shadow radius
;  Writes the file with specified filename
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2001 November, by Leslie Young, SwRI
;-

astrolib
physconstants  ; load !phys, to use !phys.k (Boltzmann constant, erg K^-1)

; Define some useful constants
km = 1e5  ; cm per km
dist = 19.935 * 1.496e8 * km; planet-Earth distance, cm
nuSTP =  0.000132 ;Refractivity at STP : 0.000132
mu = 2.2000 ;Mean molecular weight (gm/mole) : 2.2000
xh = 25938.7 * km ; planet radius of half-light for isothermal case
nlosch = 2.68719e19 ; Loschmidt's number (molecule/cm^3 at STP)
lambda = 2.2e-4 ; wavelength (in cm)
twopi = 2. * !pi
H = 60. * km ; scale height (cm)
d = H/32  ; spacing (cm).  The fft's nicer if (H/d) is a power of 2

; half-light values, for the isothermal case
dangleh = 1/dist  ; bending angle at half light, radians
angleh = -H * dangleh; bending angle at half light, radians/cm
phaseh = -(twopi*H/lambda) * angleH ; phase delay at half light
refrh = -angleh/sqrt(twopi * xh/H) ; refr. at half light
nh = (nlosch/nuSTP)*refrh; number density at half light, molecule/cm^2
Th = 142. ; temperature at half-light, K
ph = nh * !phys.k * Th 

; Calculate all the arrays
  psirefrL = refr_jk_xL(H, d, xh, j,k, xL)
  zL = (xL-xH)/H
  c = ampl / max(psirefrL)
  refrL = refrh * exp(-zL) * (1 + c * psirefrL)
  nL = (nh/refrh)*refrL
  
  psipL = p_jk_xL(H, d, xh, j,k, xL)
  pL = ph * exp(-zL) * (1 + c * psipL)

  tL = pL / (nL * !phys.k)
    
  psiphaseL = phase_jk_xL(H, d, xh, j,k, xL)
  phaseL = phaseh * exp(-zL) * (1 + c * psiphaseL)
  
  psiangleL = angle_jk_xL(H, d, xh, j,k, xL)
  angleL = angleH * exp(-zL) * (1 + c * psiangleL)
  yL = xL + dist*angleL
  
  psidangleL = dangle_jk_xL(H, d, xh, j,k, xL)
  dangleL = dangleh * exp(-zL) * (1 + c * psidangleL)
  fluxL = 1./abs((1. + dist*dangleL))

; Plot it
  window, 0, xs=1000,ys=500
  !p.multi=[0,2,0]
  plot, tL, zL, xr=[min(tL)-5, max(tL)+5], $
    title='temperature profile', $
    xtitl='Temperature (K)', $
    ytitle='Altitude in scale heights from half-light'
  plot, (yL-(xh-H))/H, fluxL, xr=[-20,10], $
    title='lightcurve', $
    xtitl='Planet radius in scale heights from half-light', $
    ytitle='Flux'

; Print it
  if fn eq 'stdout' then lun=-1 else openw, lun, fn, /get_lun
  n = n_elements(xL)
  i0 = min(where (yL gt 0) )
  i=i0

  printf, lun, dist/km, form='("dist =",E11.3," ; planet-Earth distance, km")'
  printf, lun, nuSTP, form='("nuSTP =",E11.3," ;Refractivity at STP")'
  printf, lun, xh/km, form='("xh =",F11.3," ; planet radius of half-light for isothermal case (km)")'
  printf, lun, lambda*1e4, form='("lambda =",E11.3,"; wavelength (micron")'
  printf, lun, H/km, form='("H =",F11.3,"; scale height (km)")'
  printf, lun, d/km, form='("dist =",F11.3,"; spacing (km).")'

  printf, lun, j, form='("j =",I3," ; binary dilation parameter")'
  printf, lun, k, form='("k =",I3," ; binary translation parameter")'
  printf, lun, ampl, form='("ampl =",I3," ; amplitude of wavelet")'
  
  printf, lun, 'x         = array of planet radii (km)'
  printf, lun, 'n         = array of number densities (molecule cm^-3)'
  printf, lun, 'p         = array of pressures (microbar)'
  printf, lun, 't         = array of temperatures (K)'
  printf, lun, 'refr      = array of refractivities (unitless)'
  printf, lun, 'phase     = array of phases (unitless)'
  printf, lun, 'angle     = array of bending angle (radians)'
  printf, lun, 'dangle/dx = array of derivative of bending angle (radians/cm)'
  printf, lun, 'y         = array of shadow radii (km)'
  printf, lun, 'flux      = array of fluxes (unitless)'
  formh='(A8,A9,A9,A7, A10,A10,A11,A11, A8,A6)'
  printf, lun,'x', 'n',    'p',   'T', 'refr', 'phase', 'angle', 'dangle/dx', 'y',  'flux', form=formh
  printf, lun,'km','cm^-3','ubar','K', '-',    '-',     'rad',   'rad/cm',    'km', '-', form=formh
  form='(F8.1,E9.2,E9.2,F7.1, E10.3,E10.3,E11.3,E11.3, F8.1,F6.3)'
 for i=i0,n-1 do $
    printf, lun,xL[i]/km, nL[i], pL[i], tL[i], $
              refrL[i], phaseL[i], angleL[i], dangleL[i], yL[i]/km, fluxL[i], form=form
  if fn ne 'stdout' then close, lun
end