pro recon_atm

; Paul Withers, 2001.07.20
; Open University, Britain

; Called by:				N/A
; Calls:				get_n, get_openchute, get_dt, 
;					get_a, get_m, get_cd, get_uncert,
;					get_relsigmacd, get_convergence
; Declares common blocks:		N/A
; Reads from:				aeroacc_res.dat, vrel_res.dat, 
;					height_res.dat
; Writes to:				rho_res.dat, press_res.dat, 
;					temp_res.dat, test_model_atm.dat
; Options:				range_0, height_0, gm, 
; 					mean_molecular_mass,
; 					universal_gas_constant
; Purpose:				Reconstruct atmospheric structure

common marsblk
; Properties of Mars and physical constants
; Defined in marsprops

common sctblk
; Properties of spacecraft and datasets
; Defined in sctprops

forward_function get_n, get_openchute, get_dt, get_a, get_m, get_cd, $
 get_relsigmacd, get_convergence

print, 'Starting recon_atm.pro'

convergence = sctpropsstruct.convergence
range_0 = sctpropsstruct.range_0
height_0 = sctpropsstruct.height_0
omega = marspropsstruct.omega
gm = marspropsstruct.gm
rref = marspropsstruct.rref
c20 = marspropsstruct.c20
mean_molecular_mass = marspropsstruct.mean_molecular_mass
universal_gas_constant = marspropsstruct.universal_gas_constant
avogadro = marspropsstruct.avogadro
moldiameter = marspropsstruct.moldiameter
ratioheats = marspropsstruct.ratioheats
n = sctpropsstruct.n
tchute = sctpropsstruct.tchute
tdetect = sctpropsstruct.tdetect
nchute = sctpropsstruct.nchute
ndetect = sctpropsstruct.ndetect
npadtop = sctpropsstruct.npadtop
r0 = sctpropsstruct.r0
sct = sctpropsstruct.sct
a = sctpropsstruct.a
m = sctpropsstruct.m
alimit = sctpropsstruct.alimit
nalpha = sctpropsstruct.nalpha
uncert = sctpropsstruct.uncert
nmonte_traj = sctpropsstruct.nmonte_traj
nmonte_atm = sctpropsstruct.nmonte_atm
sca = sctpropsstruct.aerodb.sca
alphaerr = sctpropsstruct.aerodb.alphaerr
nomt0 = sctpropsstruct.nomt0
snomt0 = sctpropsstruct.snomt0
outdir = sctpropsstruct.outdir
itlimit = sctpropsstruct.itlimit

; Declare arrays
; acc_array is the magnitude of the aerodynamic acceleration, ms-2
; vrel_array is the magnitude of the relative velocity, ms-1
; height_array is radial distance, m
aeroacc_array = dblarr(n)
lat_array = aeroacc_array
lon_array = aeroacc_array
t_array = aeroacc_array
vrel_array = aeroacc_array
height_array = aeroacc_array
axial_array = aeroacc_array
normal_array = aeroacc_array

if uncert eq 0 then begin 
; Read in data from results of recon_traj.pro
openr, lun1, outdir + 'traj_' + sct + '_aeroacc_nominal.dat', /get_lun
readf, lun1, aeroacc_array
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_axial_nominal.dat', /get_lun
readf, lun1, axial_array
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_normal_nominal.dat', /get_lun
readf, lun1, normal_array
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_vrel_nominal.dat', /get_lun
readf, lun1, vrel_array
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_height_nominal.dat', /get_lun
readf, lun1, height_array
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_t_nominal.dat', /get_lun
readf, lun1, t_array
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_lat_nominal.dat', /get_lun
readf, lun1, lat_array
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_lon_nominal.dat', /get_lun
readf, lun1, lon_array
free_lun, lun1
endif

if uncert ne 0. then begin
traj_monte_t = dblarr(nmonte_traj,n)
traj_monte_height = traj_monte_t
traj_monte_lat = traj_monte_t
traj_monte_lon = traj_monte_t
traj_monte_vrel = traj_monte_t
traj_monte_aeroacc = traj_monte_t
traj_monte_axial = traj_monte_t
traj_monte_normal = traj_monte_t

openr, lun1, outdir + 'traj_' + sct + '_monte_t.dat', /get_lun
readf, lun1, traj_monte_t
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_monte_height.dat', /get_lun
readf, lun1, traj_monte_height
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_monte_lat.dat', /get_lun
readf, lun1, traj_monte_lat
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_monte_lon.dat', /get_lun
readf, lun1, traj_monte_lon
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_monte_vrel.dat', /get_lun
readf, lun1, traj_monte_vrel
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_monte_aeroacc.dat', /get_lun
readf, lun1, traj_monte_aeroacc
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_monte_axial.dat', /get_lun
readf, lun1, traj_monte_axial
free_lun, lun1
openr, lun1, outdir + 'traj_' + sct + '_monte_normal.dat', /get_lun
readf, lun1, traj_monte_normal
free_lun, lun1
endif

; nchute here or later????
; How many data points in these arrays occur prior to opening of parachute?

;aa = where(t_array eq tchute)
;nchute = aa(0)
print, 'nchute = ', nchute

;aa = where(t_array eq tdetect)
;ndetect = aa(0)
print, 'ndetect = ', ndetect
 
n = nchute - ndetect
; Put this here so I can cut trajectory from full traj to short traj

if uncert eq 0. then begin
t_array = t_array(ndetect:nchute-1)
aeroacc_array = aeroacc_array(ndetect:nchute-1)
axial_array = axial_array(ndetect:nchute-1)
normal_array = normal_array(ndetect:nchute-1)
vrel_array = vrel_array(ndetect:nchute-1)
height_array = height_array(ndetect:nchute-1)
lat_array = lat_array(ndetect:nchute-1)
lon_array = lon_array(ndetect:nchute-1)
endif


if uncert ne 0. then begin
atm_monte_height = dblarr(nmonte_atm,n)
atm_monte_t = atm_monte_height
atm_monte_lat = atm_monte_height
atm_monte_lon = atm_monte_height
atm_monte_vrel = atm_monte_height
atm_monte_aeroacc = atm_monte_height
atm_monte_axial = atm_monte_height
atm_monte_normal = atm_monte_height
atm_monte_ca = atm_monte_height
atm_monte_cn = atm_monte_height
atm_monte_alpha = atm_monte_height
atm_monte_mmw = atm_monte_height
atm_monte_area = atm_monte_height
atm_monte_mass = atm_monte_height
atm_monte_rho = atm_monte_height
atm_monte_press = atm_monte_height
atm_monte_temp = atm_monte_height
atm_monte_kn = atm_monte_height
atm_monte_ma = atm_monte_height
atm_monte_grav = atm_monte_height
atm_monte_err = atm_monte_height
endif

a_array = a + dblarr(n)
m_array = m + dblarr(n)

j=0L
while j lt nmonte_atm do begin
print, 'j = ', j, ' nmonte_atm = ', nmonte_atm

itfailed = 0. ; change to 1 if iteration failed

err_array = dblarr(n) ; any errors?

if uncert ne 0. then begin
if nmonte_traj ne nmonte_atm then k = floor(randomu(seed) * nmonte_traj)
; Randomly pick one of my monte trajectories
if nmonte_traj eq nmonte_atm then k = j
; Or use 1-1 trajectories
aeroacc_array = traj_monte_aeroacc(k,ndetect:nchute-1)
axial_array = traj_monte_axial(k,ndetect:nchute-1)
normal_array = traj_monte_normal(k,ndetect:nchute-1)
vrel_array = traj_monte_vrel(k,ndetect:nchute-1)
t_array = traj_monte_t(k,ndetect:nchute-1)
height_array = traj_monte_height(k,ndetect:nchute-1)
lat_array = traj_monte_lat(k,ndetect:nchute-1)
lon_array = traj_monte_lon(k,ndetect:nchute-1)
endif

aratio_array = normal_array/axial_array
;aa = where(normal_array lt alimit)
;if aa(0) ne -1 then aratio_array(aa) = 0.
aratio_array(0:nalpha-ndetect) = 0.

g_r = double(-1.)*gm/(height_array^2) - double(3.)*gm*(rref^2)/(height_array^4)*double(0.5)*(double(3.)*(cos(!pi/180.*(90.-lat_array))^2)-double(1.))*c20

c_r = omega^2 * height_array * sin(!pi/180. * (90.-lat_array))^2

grav_array = gm / (height_array)^2
grav_array = abs(g_r + c_r)
; See 2005.07.06 notes for details on signs


; Assuming spherical symmetry, calculate acceleration due to gravity
; at each data point, ms-2
; More elegant procedures would link this in with get_grav

mmw_array = dblarr(n) + mean_molecular_mass
; Allowing option of mmw that varies with altitude

cd_array = dblarr(n) + 2.
; Dimensionless drag coefficient

; Read in reference area and mass of spacecraft
a_array = a + height_array*0.
; use diameter from Schoenenbeger et al. AIAA-2005-0056
; Put in get_a() soon

aoverm = a/m
print, '(m2) Area = ', a
print, '(kg) Mass = ', m

print, 'Calculating first-cut atmospheric density profile with CD=2'
rho_array = double(2.) *axial_array /( (vrel_array^2) * cd_array * aoverm)
; Solve basic force balance equation for atmospheric density, kg m-3
aa = where(axial_array lt 0.)
if aa(0) ne -1 then rho_array(aa) = 0.

newrho_array =dblarr(n)
; Another density array, useful for testing convergence of iteration
; to stable solution for density
compare_rho = 0
i=0
; Flag to mark whether convergence has been achieved
while compare_rho eq 0 do begin
print, 'Iterating with aerodynamic database...', i

press_array = dblarr(n) ; Pa
k=1L
while k lt n do begin
press_array(k) = press_array(k-1) - $
 0.5*(rho_array(k-1)+rho_array(k))*0.5*(grav_array(k-1)+grav_array(k)) * $
 (height_array(k)-height_array(k-1))
; Try this integration routine instead, faster
k=k+1
endwhile

aa = where(rho_array ne 0.)
grav_0 = interpol(grav_array(aa), height_array(aa), [height_0]) ; ms-2
rho_0 = interpol(rho_array(aa), height_array(aa), [height_0]) ; kg m-3
mmw_0 = interpol(mmw_array(aa), height_array(aa), [height_0]) ; kg
; gravity, mmw, and density at height_0

if uncert eq 0 then begin
aa= where(height_array lt height_0 and height_array gt (height_0-range_0) and $
 rho_array ne 0.)

result = poly_fit(height_array(aa), alog(rho_array(aa)), 1)
; Do a linear fit to log(rho) in the relevant altitude range
; to calculate the density scale height

density_scale_height_0 = double(-1.)/result(1) ; m
; Obtain density scale height at height_0
;print, '(m) Density Scale Height at height_0 = ', density_scale_height_0
temp_0 = mmw_0 * grav_0 * density_scale_height_0 / universal_gas_constant
endif ; if uncert eq 0

if uncert ne 0 then begin
temp_0 = nomt0 + snomt0 * randomn(seed)
density_scale_height_0 = temp_0 * universal_gas_constant / (mmw_0 * grav_0)
endif

press_0 = rho_0 * grav_0 * density_scale_height_0
; Obtain pressure at height_0, Pa

press_offset = interpol(press_array, height_array, [height_0])
; What did we calculate as the pressure at height_0
; without using a constant of integration, Pa?
press_array = press_array + press_0(0) - press_offset(0)
; Shift all pressures so that pressure at height_0 is what we told it to be

aa = where(press_array lt 0.)
if aa(0) ne -1 then press_array(aa) = 0.
; deal with nasty noisy data

; Solve equation of state for temperature, K
temp_array = press_array *  mmw_array/ (rho_array * universal_gas_constant)
aa = where(rho_array eq 0.)
if aa(0) ne -1 then temp_array(aa) = 0.

kn_array = mmw_array / $
 ( sqrt(2.)*!pi * moldiameter^2 * rho_array * avogadro * sqrt(4.*a_array/!pi) )
; Equation from my 1A Physics notes
aa = where(rho_array eq 0.)
if aa(0) ne -1 then kn_array(aa) = 0.

ma_array = vrel_array / sqrt(ratioheats * universal_gas_constant * temp_array/mmw_array)
aa = where(rho_array eq 0.)
if aa(0) ne -1 then ma_array(aa) = 0.

ca_array = dblarr(n)
cn_array = dblarr(n)
alpha_array = dblarr(n)
k=0
while k lt n do begin
if rho_array(k) ne 0. then begin
junk = get_cd(kn_array(k), ma_array(k), vrel_array(k), aratio_array(k), i)
;junk = get_cdx(kn_array(k), ma_array(k), vrel_array(k), aratio_array(k), i)
ca_array(k) = junk(0)
cn_array(k) = junk(1)
alpha_array(k) = junk(2)
endif
k=k+1
endwhile

newrho_array = double(2.) * axial_array / ( vrel_array^2 * ca_array * aoverm)
; Derive updated density profile - not cos(alpha) dependent!!!
aa = where(axial_array lt 0.)
if aa(0) ne -1 then newrho_array(aa) = 0.

aa = where(rho_array ne 0. and newrho_array ne 0.)

if min(newrho_array(aa)/rho_array(aa)) ge 1-convergence and $
 max(newrho_array(aa)/rho_array(aa)) le 1+convergence then compare_rho = 1
; Compare old and new density profiles

if i gt itlimit then begin
; escape with a dummy result
rho_array = exp(-height_array) 
compare_rho = 1. 
itfailed = 1.
endif

oldrho_array = rho_array
rho_array = newrho_array
; Relabel new as old in preparation for next steps
i=i+1
; Helpful to know how many times we iterated
endwhile
print, 'Atmospheric density profile converged'

; Solve equation of hydrostatic equilibrium for pressure at each datapoint
; This integration routine is inelegant
; Improve it if possible

print, 'Calculating atmospheric pressure profile'
; Integrate pressure with no constant of integration
; It will be added later

press_array = dblarr(n) ; Pa
k=1L
while k lt n do begin
press_array(k) = press_array(k-1) - $
 0.5*(rho_array(k-1)+rho_array(k))*0.5*(grav_array(k-1)+grav_array(k)) * $
 (height_array(k)-height_array(k-1))
; Try this integration routine instead, faster
k=k+1
endwhile

;;; PW 2004.01.06
aa= where((height_0 - height_array) lt range_0 and height_array lt height_0 and $
 rho_array ne 0.)

result = poly_fit(height_array(aa), alog(rho_array(aa)), 1)
; Do a linear fit to log(rho) in the relevant altitude range
; to calculate the density scale height

density_scale_height_0 = double(-1.)/result(1) ; m
grav_0 = interpol(grav_array, height_array, [height_0]) ; ms-2
rho_0 = interpol(rho_array, height_array, [height_0]) ; kg m-3
mmw_0 = interpol(mmw_array, height_array, [height_0]) ; kg
temp_0 = mmw_0 * grav_0 * density_scale_height_0 / universal_gas_constant
; Obtain density scale height, gravity, and density at height_0
print, '(m) Density Scale Height at height_0 = ', density_scale_height_0
print, '(K) Temperature at height_0 = ', temp_0

press_0 = rho_0 * grav_0 * density_scale_height_0
; Obtain pressure at height_0, Pa

press_offset = interpol(press_array, height_array, [height_0])
; What did we calculate as the pressure at height_0
; without using a constant of integration, Pa?
press_array = press_array + press_0(0) - press_offset(0)
; Shift all pressures so that pressure at height_0 is what we told it to be

aa = where(press_array lt 0.)
if aa(0) ne -1 then press_array(aa) = 0.
; deal with nasty noisy data

print, 'Calculating atmospheric temperature profile'
; Solve equation of state for temperature, K
temp_array = press_array *  mmw_array/ (rho_array * universal_gas_constant)
aa = where(rho_array eq 0.)
if aa(0) ne -1 then temp_array(aa) = 0.

kn_array = mmw_array / $
 ( sqrt(2.)*!pi * moldiameter^2 * rho_array * avogadro * sqrt(4.*a_array/!pi) )
; Equation from my 1A Physics notes
aa = where(rho_array eq 0.)
if aa(0) ne -1 then kn_array(aa) = 0.

ma_array = vrel_array / sqrt(ratioheats * universal_gas_constant * temp_array/mmw_array)
aa = where(rho_array eq 0.)
if aa(0) ne -1 then ma_array(aa) = 0.

aa = where(press_array eq 0. or alpha_array eq alphaerr or rho_array eq 0.)
if aa(0) ne -1 then err_array(aa) = 1. ; ignore az<0 only where it occurs
if aa(0) ne -1 then err_array(*) = 1. ; ignore az<0 in entire traj
if itfailed eq 1 then err_array(*) = 1.

; Introduce some uncertainty due to aerodb...
rho_array = rho_array * (1. + randomn(seed, n) * sca * uncert)
press_array = press_array * (1. + randomn(seed, n) * sca * uncert)
temp_array = temp_array
kn_array = kn_array * (1. + randomn(seed, n) * sca * uncert)
ma_array = ma_array
ca_array = ca_array * (1. + randomn(seed, n) * sca * uncert)
cn_array = cn_array
alpha_array = alpha_array * (1. + randomn(seed, n) * sca * uncert)

xx = (height_array-r0)/1e3
aa = where( abs(xx-30.) eq min(abs(xx-30.)))

print, '30 km state is'
print, '(s) t = ', t_array(aa)
print, '(m) r = ', height_array(aa)
print, '(km) z = ', (height_array(aa)-r0)/1e3
print, '(degN) lat = ', lat_array(aa)
print, '(degE) lon = ', lon_array(aa)
print, '(m s-1) vrel = ', vrel_array(aa)
print, '(kg m-3) rho = ', rho_array(aa)
print, '(Pa) press = ', press_array(aa)
print, '(K) temp = ', temp_array(aa)
print, ''

print, 'Chute state is'
print, '(s) t = ', t_array(n-1)
print, '(m) r = ', height_array(n-1)
print, '(km) z = ', (height_array(n-1)-r0)/1e3
print, '(degN) lat = ', lat_array(n-1)
print, '(degE) lon = ', lon_array(n-1)
print, '(m s-1) vrel = ', vrel_array(n-1)
print, '(kg m-3) rho = ', rho_array(n-1)
print, '(Pa) press = ', press_array(n-1)
print, '(K) temp = ', temp_array(n-1)
print, ''



if uncert eq 0 then begin
; Store results for later use
openw, lun1, outdir + 'atm_' + sct + '_height_nominal.dat', /get_lun
printf, lun1, height_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_t_nominal.dat', /get_lun
printf, lun1, t_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_lat_nominal.dat', /get_lun
printf, lun1, lat_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_lon_nominal.dat', /get_lun
printf, lun1, lon_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_aeroacc_nominal.dat', /get_lun
printf, lun1, aeroacc_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_axial_nominal.dat', /get_lun
printf, lun1, axial_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_normal_nominal.dat', /get_lun
printf, lun1, normal_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_rho_nominal.dat', /get_lun
printf, lun1, rho_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_press_nominal.dat', /get_lun
printf, lun1, press_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_temp_nominal.dat', /get_lun
printf, lun1, temp_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_ca_nominal.dat', /get_lun
printf, lun1, ca_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_cn_nominal.dat', /get_lun
printf, lun1, cn_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_alpha_nominal.dat', /get_lun
printf, lun1, alpha_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_kn_nominal.dat', /get_lun
printf, lun1, kn_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_ma_nominal.dat', /get_lun
printf, lun1, ma_array
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_err_nominal.dat', /get_lun
printf, lun1, err_array
free_lun, lun1

;openw, lun1, 'test_model_atm.dat', /get_lun
;printf, lun1, rho_array, height_array
;free_lun, lun1

endif

if uncert ne 0. then begin
atm_monte_t(j,*) = t_array(*)
atm_monte_height(j,*) = height_array(*)
atm_monte_lat(j,*) = lat_array(*)
atm_monte_lon(j,*) = lon_array(*)
atm_monte_vrel(j,*) = vrel_array(*) 
atm_monte_aeroacc(j,*) = aeroacc_array(*)
atm_monte_axial(j,*) = axial_array(*)
atm_monte_normal(j,*) = normal_array(*)
atm_monte_ca(j,*) = ca_array(*)
atm_monte_cn(j,*) = cn_array(*)
atm_monte_alpha(j,*) = alpha_array(*)
atm_monte_mmw(j,*) = mmw_array(*)
atm_monte_area(j,*) = a_array(*)
atm_monte_mass(j,*) = m_array(*)
atm_monte_rho(j,*) = rho_array(*)
atm_monte_press(j,*) = press_array(*)
atm_monte_temp(j,*) = temp_array(*)
atm_monte_kn(j,*) = kn_array(*)
atm_monte_ma(j,*) = ma_array(*)
atm_monte_err(j,*) = err_array(*)
atm_monte_grav(j,*) = grav_array(*)
endif ; if get_uncert() ne 0. 

j=j+1
endwhile

if uncert ne 0. then begin


; Pad top of arrays!

nnew = n + npadtop

xxx = atm_monte_height
atm_monte_height = dblarr(nmonte_atm,nnew)
atm_monte_height(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_t
atm_monte_t = dblarr(nmonte_atm,nnew)
atm_monte_t(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_lat
atm_monte_lat = dblarr(nmonte_atm,nnew)
atm_monte_lat(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_lon
atm_monte_lon = dblarr(nmonte_atm,nnew)
atm_monte_lon(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_vrel
atm_monte_vrel = dblarr(nmonte_atm,nnew)
atm_monte_vrel(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_aeroacc
atm_monte_aeroacc = dblarr(nmonte_atm,nnew)
atm_monte_aeroacc(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_axial
atm_monte_axial = dblarr(nmonte_atm,nnew)
atm_monte_axial(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_normal
atm_monte_normal = dblarr(nmonte_atm,nnew)
atm_monte_normal(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_ca
atm_monte_ca = dblarr(nmonte_atm,nnew)
atm_monte_ca(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_cn
atm_monte_cn = dblarr(nmonte_atm,nnew)
atm_monte_cn(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_alpha
atm_monte_alpha = dblarr(nmonte_atm,nnew)
atm_monte_alpha(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_mmw
atm_monte_mmw = dblarr(nmonte_atm,nnew)
atm_monte_mmw(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_area
atm_monte_area = dblarr(nmonte_atm,nnew)
atm_monte_area(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_mass
atm_monte_mass = dblarr(nmonte_atm,nnew)
atm_monte_mass(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_rho
atm_monte_rho = dblarr(nmonte_atm,nnew)
atm_monte_rho(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_press
atm_monte_press = dblarr(nmonte_atm,nnew)
atm_monte_press(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_temp
atm_monte_temp = dblarr(nmonte_atm,nnew)
atm_monte_temp(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_kn
atm_monte_kn = dblarr(nmonte_atm,nnew)
atm_monte_kn(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_ma
atm_monte_ma = dblarr(nmonte_atm,nnew)
atm_monte_ma(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_grav
atm_monte_grav = dblarr(nmonte_atm,nnew)
atm_monte_grav(*,npadtop:*) = xxx(*,*)
xxx = atm_monte_err
atm_monte_err = dblarr(nmonte_atm,nnew)
atm_monte_err(*,npadtop:*) = xxx(*,*)

openw, lun1, outdir + 'atm_' + sct + '_monte_t.dat', /get_lun
printf, lun1, atm_monte_t
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_height.dat', /get_lun
printf, lun1, atm_monte_height
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_lat.dat', /get_lun
printf, lun1, atm_monte_lat
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_lon.dat', /get_lun
printf, lun1, atm_monte_lon
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_vrel.dat', /get_lun
printf, lun1, atm_monte_vrel
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_aeroacc.dat', /get_lun
printf, lun1, atm_monte_aeroacc
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_axial.dat', /get_lun
printf, lun1, atm_monte_axial
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_normal.dat', /get_lun
printf, lun1, atm_monte_normal
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_ca.dat', /get_lun
printf, lun1, atm_monte_ca
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_cn.dat', /get_lun
printf, lun1, atm_monte_cn
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_alpha.dat', /get_lun
printf, lun1, atm_monte_alpha
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_mmw.dat', /get_lun
printf, lun1, atm_monte_mmw
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_area.dat', /get_lun
printf, lun1, atm_monte_area
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_mass.dat', /get_lun
printf, lun1, atm_monte_mass
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_rho.dat', /get_lun
printf, lun1, atm_monte_rho
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_press.dat', /get_lun
printf, lun1, atm_monte_press
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_temp.dat', /get_lun
printf, lun1, atm_monte_temp
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_kn.dat', /get_lun
printf, lun1, atm_monte_kn
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_ma.dat', /get_lun
printf, lun1, atm_monte_ma
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_err.dat', /get_lun
printf, lun1, atm_monte_err
free_lun, lun1
openw, lun1, outdir + 'atm_' + sct + '_monte_grav.dat', /get_lun
printf, lun1, atm_monte_grav
free_lun, lun1

n = nnew

best_height = dblarr(n)
best_t = best_height
best_lat = best_height
best_lon = best_height
best_vrel = best_height
best_aeroacc = best_height
best_axial = best_height
best_normal = best_height
best_ca = best_height
best_cn = best_height
best_alpha = best_height
best_mmw = best_height
best_area = best_height
best_mass = best_height
best_rho = best_height
best_press = best_height
best_temp = best_height
best_kn = best_height
best_ma = best_height
best_grav = best_height

best_averr = best_height

sigma_height = best_height
sigma_t = best_height
sigma_lat = best_height
sigma_lon = best_height
sigma_vrel = best_height
sigma_aeroacc = best_height
sigma_axial = best_height
sigma_normal = best_height
sigma_ca = best_height
sigma_cn = best_height
sigma_alpha = best_height
sigma_mmw = best_height
sigma_area = best_height
sigma_mass = best_height
sigma_rho = best_height
sigma_press = best_height
sigma_temp = best_height
sigma_kn = best_height
sigma_ma = best_height
sigma_grav = best_height

print, 'Calculating mean and uncertainty for results'
print, 'nmonte_atm = ', nmonte_atm
i=0L
while i lt n do begin
best_averr(i) = mean(atm_monte_err(*,i))
aa = where(atm_monte_err(*,i) eq 0.) ; where error flag=0

if n_elements(aa) gt 1 then begin
best_height(i) = mean(atm_monte_height(aa,i))
best_t(i) = mean(atm_monte_t(aa,i))
best_lat(i) = mean(atm_monte_lat(aa,i))
best_lon(i) = mean(atm_monte_lon(aa,i))
best_vrel(i) = mean(atm_monte_vrel(aa,i))
best_aeroacc(i) = mean(atm_monte_aeroacc(aa,i))
best_axial(i) = mean(atm_monte_axial(aa,i))
best_normal(i) = mean(atm_monte_normal(aa,i))
best_ca(i) = mean(atm_monte_ca(aa,i))
best_cn(i) = mean(atm_monte_cn(aa,i))
best_alpha(i) = mean(atm_monte_alpha(aa,i))
best_mmw(i) = mean(atm_monte_mmw(aa,i))
best_area(i) = mean(atm_monte_area(aa,i))
best_mass(i) = mean(atm_monte_mass(aa,i))
best_rho(i) = mean(atm_monte_rho(aa,i))
best_press(i) = mean(atm_monte_press(aa,i))
best_temp(i) = mean(atm_monte_temp(aa,i))
best_kn(i) = mean(atm_monte_kn(aa,i))
best_ma(i) = mean(atm_monte_ma(aa,i))
best_grav(i) = mean(atm_monte_grav(aa,i))

junk = moment(atm_monte_height(aa,i))
sigma_height(i) = sqrt(junk(1))
junk = moment(atm_monte_t(aa,i))
sigma_t(i) = sqrt(junk(1))
junk = moment(atm_monte_lat(aa,i))
sigma_lat(i) = sqrt(junk(1))
junk = moment(atm_monte_lon(aa,i))
sigma_lon(i) = sqrt(junk(1))
junk = moment(atm_monte_vrel(aa,i))
sigma_vrel(i) = sqrt(junk(1))
junk = moment(atm_monte_aeroacc(aa,i))
sigma_aeroacc(i) = sqrt(junk(1))
junk = moment(atm_monte_axial(aa,i))
sigma_axial(i) = sqrt(junk(1))
junk = moment(atm_monte_normal(aa,i))
sigma_normal(i) = sqrt(junk(1))
junk = moment(atm_monte_ca(aa,i))
sigma_ca(i) = sqrt(junk(1))
junk = moment(atm_monte_cn(aa,i))
sigma_cn(i) = sqrt(junk(1))
junk = moment(atm_monte_alpha(aa,i))
sigma_alpha(i) = sqrt(junk(1))
junk = moment(atm_monte_mmw(aa,i))
sigma_mmw(i) = sqrt(junk(1))
junk = moment(atm_monte_area(aa,i))
sigma_area(i) = sqrt(junk(1))
junk = moment(atm_monte_mass(aa,i))
sigma_mass(i) = sqrt(junk(1))
junk = moment(atm_monte_rho(aa,i))
sigma_rho(i) = sqrt(junk(1))
junk = moment(atm_monte_press(aa,i))
sigma_press(i) = sqrt(junk(1))
junk = moment(atm_monte_temp(aa,i))
sigma_temp(i) = sqrt(junk(1))
junk = moment(atm_monte_kn(aa,i))
sigma_kn(i) = sqrt(junk(1))
junk = moment(atm_monte_ma(aa,i))
sigma_ma(i) = sqrt(junk(1))
junk = moment(atm_monte_grav(aa,i))
sigma_grav(i) = sqrt(junk(1))
endif

i=i+1
endwhile

openw, lun1, outdir + 'atm_' + sct + '_err_mean.dat', /get_lun
printf, lun1, best_averr
free_lun, lun1


openw, lun1, outdir + 'atm_' + sct + '_height_mean.dat', /get_lun
printf, lun1, best_height
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_t_mean.dat', /get_lun
printf, lun1, best_t
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_lat_mean.dat', /get_lun
printf, lun1, best_lat
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_lon_mean.dat', /get_lun
printf, lun1, best_lon
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_vrel_mean.dat', /get_lun
printf, lun1, best_vrel
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_aeroacc_mean.dat', /get_lun
printf, lun1, best_aeroacc
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_axial_mean.dat', /get_lun
printf, lun1, best_axial
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_normal_mean.dat', /get_lun
printf, lun1, best_normal
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_ca_mean.dat', /get_lun
printf, lun1, best_ca
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_cn_mean.dat', /get_lun
printf, lun1, best_cn
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_alpha_mean.dat', /get_lun
printf, lun1, best_alpha
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_mmw_mean.dat', /get_lun
printf, lun1, best_mmw
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_area_mean.dat', /get_lun
printf, lun1, best_area
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_mass_mean.dat', /get_lun
printf, lun1, best_mass
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_rho_mean.dat', /get_lun
printf, lun1, best_rho
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_press_mean.dat', /get_lun
printf, lun1, best_press
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_temp_mean.dat', /get_lun
printf, lun1, best_temp
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_kn_mean.dat', /get_lun
printf, lun1, best_kn
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_ma_mean.dat', /get_lun
printf, lun1, best_ma
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_grav_mean.dat', /get_lun
printf, lun1, best_grav
free_lun, lun1



openw, lun1, outdir + 'atm_' + sct + '_height_sigma.dat', /get_lun
printf, lun1, sigma_height
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_t_sigma.dat', /get_lun
printf, lun1, sigma_t
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_lat_sigma.dat', /get_lun
printf, lun1, sigma_lat
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_lon_sigma.dat', /get_lun
printf, lun1, sigma_lon
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_vrel_sigma.dat', /get_lun
printf, lun1, sigma_vrel
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_aeroacc_sigma.dat', /get_lun
printf, lun1, sigma_aeroacc
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_axial_sigma.dat', /get_lun
printf, lun1, sigma_axial
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_normal_sigma.dat', /get_lun
printf, lun1, sigma_normal
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_ca_sigma.dat', /get_lun
printf, lun1, sigma_ca
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_cn_sigma.dat', /get_lun
printf, lun1, sigma_cn
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_alpha_sigma.dat', /get_lun
printf, lun1, sigma_alpha
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_mmw_sigma.dat', /get_lun
printf, lun1, sigma_mmw
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_area_sigma.dat', /get_lun
printf, lun1, sigma_area
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_mass_sigma.dat', /get_lun
printf, lun1, sigma_mass
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_rho_sigma.dat', /get_lun
printf, lun1, sigma_rho
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_press_sigma.dat', /get_lun
printf, lun1, sigma_press
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_temp_sigma.dat', /get_lun
printf, lun1, sigma_temp
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_kn_sigma.dat', /get_lun
printf, lun1, sigma_kn
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_ma_sigma.dat', /get_lun
printf, lun1, sigma_ma
free_lun, lun1

openw, lun1, outdir + 'atm_' + sct + '_grav_sigma.dat', /get_lun
printf, lun1, sigma_grav
free_lun, lun1

endif ; if uncert ne 0.

if uncert eq 0 then zkm = (height_array-r0)/1e3
if uncert ne 0 then zkm = (best_height-r0)/1e3

if uncert ne 0 then print, 'Fraction of profiles ignored ', best_averr(n-1)

print, 'End of recon_atm.pro'
stop
end