pro makeodyaccresultsformro

; Paul Withers, 2006.01.06
; Center for Space Physics, Boston University

; Generate Rapid Data Products from ODY ACC data in support of MRO aerobraking


; Nav_recontr_thru_P338.xls has
; orb, date, time, 
; z, lat, lon, r-3397.2
; LST, Ls, vrel at peri, 
; r, a, e, i
; mass, period

area = 11. ; Reference area of ODY
; m2 from Smith and Bell (2005)
aerocoeff = 2. ; aerodynamic force coefficient
gm = 4.2828e13 ; GM of Mars
; m3 s-2 from Lodders and Fegley (1998)

datadir = 'rawdata/'
outdir = 'outdir/'
navfile = 'navstuff_v01.csv'
; Stripped-down version of file from JPL server MMONT1, 
; M01-AB/NAV/Reconstruction_Data/Nav_reconstr_thru_P0338.xls

;;;
; Find number of lines in Nav file
;;;
junk=''
j=0
openr, lun1, navfile, /get_lun
while not eof(lun1) do begin
readf, lun1, junk
j=j+1
endwhile
free_lun, lun1
norb = j

;;;
; Define variables useful for constant altitude data
;;;
nrep = 8
; number of columns per altitude level
z0list = [100, 110, 120, 130, 140]
; altitude levels, km
strz0list = strcompress(string(round(z0list)), /remove_all)
naltentries = 16+nrep*n_elements(z0list)*2
; 16 preamble entries, plus the actual results,
; gives the total number of entries used
; to print out the constant altitude results
bigaltarray = dblarr(norb,naltentries-1)
; put all the alt results in here for printing
dz = 5. ; fit to data between z0+dz and z0-dz, dz in km
altlimit = 3. ; only perform fit if data exist 
; below z0-dz and above z0+dz
maxi=2 ; 2 => inbound and outbound
inout = ['in', 'out']

;;;
; Create blank arrays for data in Nav file
;;;
orb = dblarr(norb) 		; orbit number (-)
utcdatetime = strarr(norb) 	; periapsis UTC date/time (-)
palt = orb 			; periapsis altitude above unknown reference surface (km)
prefalt = orb 			; periapsis altitude above 3397.2 km sphere (km)
pradius = orb 			; periapsis radius (km)
plat = orb 			; periapsis latitude (degN)
plon = orb 			; periapsis longitude (degE)
plst = orb 			; periapsis LST (hrs)
pls = orb 			; periapsis Ls (deg)
pvrel = orb 			; atmosphere-relative velocity at periapsis (km s-1)
psma = orb 			; periapsis semi-major axis (km)
pecc = orb			; periapsis eccentricity (-)
pinc = orb 			; periapsis inclination (deg)
pmass = orb 			; spacecraft mass (kg)
period = orb 			; period (hrs)
pnoise = orb 			; noise in raw acc measurements (m s-2)
prho1 = orb 			; density derived from raw acceleration (kg m-3)
psrho1 = orb			; 1-sigma uncertainty in density derived from raw acceleration (kg m-3)
prho7 = orb 			; density derived from 7-point mean acceleration (kg m-3)
psrho7 = orb 			; 1-sigma uncertainty in density derived from 7-point mean acceleration (kg m-3)
prho39 = orb 			; density derived from 39-point mean acceleration (kg m-3)
psrho39 = orb 			; 1-sigma uncertainty in density derived from 39-point mean acceleration (kg m-3)
ptime = orb 			; time after periapsis of datapoint that is closest to periapsis (s)
pratptime = orb 		; radius at this closest datapoint (km)

strorb = utcdatetime
; 3 character string with orbit number, eg 007

;;;
; Read in data from Nav file
;;;
junk=''
j=0
openr, lun1, navfile, /get_lun
while not eof(lun1) do begin
readf, lun1, junk
xx = str_sep(junk, ',')
orb(j) = double(xx(0))
yy = str_sep(xx(3), '"')
utcdatetime(j) = yy(1)
palt(j) = double(xx(7))
prefalt(j) = double(xx(8))
pradius(j) = double(xx(24))
plat(j) = double(xx(10))
plon(j) = double(xx(11))
plst(j) = double(xx(19))
pls(j) = double(xx(21))
pvrel(j) = double(xx(22))
psma(j) = double(xx(25))
pecc(j) = double(xx(26))
pinc(j) = double(xx(27))
j=j+1
endwhile
free_lun, lun1

; Make 3 character string with orbit number, eg 007
i=0
while i lt norb do begin
junk = strcompress(string(round(orb(i))), /remove_all)
while strlen(junk) lt 3 do junk = '0'+junk
strorb(i) = junk
i=i+1
endwhile

;;;
; Get ODY mass and period data
; Data comes from same Nav_reconstr_thru_P0338.xls
; just moved into a different stripped-down file
;;;
xx = dblarr(norb)
yy = xx
zz = xx
openr, lun1, 'odyorb.txt', /get_lun
readf, lun1, xx
free_lun, lun1
openr, lun1, 'odymass.txt', /get_lun
readf, lun1, yy
free_lun, lun1
openr, lun1, 'odyperiod.txt', /get_lun
readf, lun1, zz
free_lun, lun1
i=0
while i lt norb do begin
aa = where(orb eq xx(i))
pmass(aa(0)) = yy(i)
period(aa(0)) = zz(i)
i=i+1
endwhile

;;;
; Write out table of ancilliary information
;;;

outfileanc = 'ancinfo.txt'

formatanc = '(F9.1,",",A30,13(",",E12.5))'
openw, lun1, outdir+outfileanc, /get_lun
printf, lun1, ';;; start of header ;;;'
printf, lun1, '21 lines of header'
printf, lun1, strcompress(string(n_elements(orb)), /remove_all) + $
 ' lines of data'
printf, lun1, '14 columns of data separated by commas'
printf, lun1, 'Source of data is JPL server MMONT1, file ' + $
 'M01-AB/NAV/Reconstruction_Data/Nav_reconstr_thru_P0338.xls'
printf, lun1, '1: orbit number (-)'
printf, lun1, '2: periapsis UTC date/time (-)'
printf, lun1, '3: periapsis altitude above unknown reference surface (km)'
printf, lun1, '4: periapsis altitude above 3397.2 km sphere (km)'
printf, lun1, '5: periapsis radius (km)'
printf, lun1, '6: periapsis latitude (degN)'
printf, lun1, '7: periapsis longitude (degE)'
printf, lun1, '8: periapsis LST (hrs)'
printf, lun1, '9: periapsis Ls (deg)'
printf, lun1, '10: periapsis vrel (km s-1)'
printf, lun1, '11: periapsis semi-major axis (km)'
printf, lun1, '12: periapsis eccentricity (-)'
printf, lun1, '13: periapsis inclination (deg)'
printf, lun1, '14: spacecraft mass (kg)'
printf, lun1, '15: period (hrs)'
printf, lun1, ';;; end of header ;;;'
i=0
while i lt norb do begin
printf, lun1, format=formatanc, $
orb(i), utcdatetime(i), palt(i), prefalt(i), pradius(i), $
 plat(i), plon(i), plst(i), pls(i), $
 pvrel(i), psma(i), pecc(i), pinc(i), pmass(i), period(i)
i=i+1
endwhile
free_lun, lun1

;;;
; Prepare output file for alt results
;;;

outfilealt = 'allalt.txt'
junkfile = 'junk.txt'
spawn, 'ls -1 ' + datadir + ' > ' + junkfile
i=0
junk=''
openr, lun1, junkfile, /get_lun
while not eof(lun1) do begin
readf, lun1, junk
i=i+1
endwhile
nrawdatafiles = i
free_lun, lun1
spawn, '/bin/rm ' + junkfile

openw, lun2, outdir+outfilealt, /get_lun
printf, lun1, ';;; start of header ;;;'
printf, lun1, '42 lines of header'
printf, lun1, strcompress(string(nrawdatafiles), /remove_all) + $
 ' lines of data'
printf, lun1, strcompress(string(naltentries), /remove_all) + $
 ' columns of data separated by commas'
printf, lun1, 'Sources of data are JPL server MMONT1, file ' + $
 'M01-AB/NAV/Reconstruction_Data/Nav_reconstr_thru_P0338.xls'
printf, lun1, 'and pXXXacc.dat from Jim Murphy/PDS'
printf, lun1, '1: orbit number (-)'
printf, lun1, '2: periapsis UTC date/time (-)'
printf, lun1, '3: periapsis altitude above unknown reference surface (km)'
printf, lun1, '4: periapsis altitude above 3397.2 km sphere (km)'
printf, lun1, '5: periapsis radius (km)'
printf, lun1, '6: periapsis latitude (degN)'
printf, lun1, '7: periapsis longitude (degE)'
printf, lun1, '8: periapsis LST (hrs)'
printf, lun1, '9: periapsis Ls (deg)'
printf, lun1, '10: periapsis vrel (km s-1)'
printf, lun1, '11: periapsis semi-major axis (km)'
printf, lun1, '12: periapsis eccentricity (-)'
printf, lun1, '13: periapsis inclination (deg)'
printf, lun1, '14: spacecraft mass (kg)'
printf, lun1, '15: period (hrs)'
printf, lun1, '16: noise in raw acc measurements (m s-2)'
printf, lun1, '17-24: 100 km inbound results'
printf, lun1, '25-32: 100 km outbound results'
printf, lun1, '33-40: 110 km inbound results'
printf, lun1, '41-48: 110 km outbound results'
printf, lun1, '49-56: 120 km inbound results'
printf, lun1, '57-64: 120 km outbound results'
printf, lun1, '65-72: 130 km inbound results'
printf, lun1, '73-80: 130 km outbound results'
printf, lun1, '81-88: 140 km inbound results'
printf, lun1, '89-96: 140 km outbound results'
printf, lun1, 'Eight columns (A-H) for each altitude level'
printf, lun1, 'A: fitted density (kg m-3)'
printf, lun1, 'B: 1-sigma uncertainty in fitted density (kg m-3)'
printf, lun1, 'C: fitted density scale height (km)'
printf, lun1, 'D: 1-sigma uncertainty in fitted density scale height (km)'
printf, lun1, 'E: Reduced chi-squared value for fit (-)'
printf, lun1, 'F: Number of data points used in fit (-)'
printf, lun1, 'G: latitude at appropriate altitude (degN)'
printf, lun1, 'H: longitude at appropriate altitude (degE)'
printf, lun1, ';;; end of header ;;;'

formatalt = '(F9.1,",",A30,14(",",E12.5),10(",",E12.5,",",E12.5,",",E12.5,",",E12.5,",",E12.5,",",E12.5,",",E12.5,",",E12.5))'

;;;
; Read in an equipotential surface based
; on some emails from Bob Tolson in 2002
; This may or may not be used later
;;;

equipot_radius = dblarr(181,361) ; radius of equipotential surface, km
equipot_lat = equipot_radius ; latitude of equipotential surface, deg N
equipot_lon = equipot_radius ; longitude of equipotential surface, deg E
openr, lun1, 'equipot_radius.dat', /get_lun
readf, lun1, equipot_radius, $ 
 equipot_lat, equipot_lon
free_lun, lun1
equipot_lat=equipot_lat*180./!pi ; -90 to 90
equipot_lon=equipot_lon*180./!pi ; -180 to 180





;;;;;;;;;
;;;;;;;;;
; start loop for all orbits
;;;;;;;;;
;;;;;;;;;

ii = 0 
while ii lt norb do begin

;ii = 74
;while ii lt 75 do begin







;;;
; Start reading in ACC data for one aerobraking pass
;;;

junk=''
nheaderlines=4
i=0
infile = 'p' + strorb(ii) + 'acc.txt'
if file_test(datadir+infile) eq 1 then begin
; Skip this section if raw acc datafile does not exist

openr, lun1, datadir+infile, /get_lun
while not eof(lun1) do begin
readf, lun1, junk
i=i+1
endwhile ; while not eof(lun1) do begin
free_lun, lun1
n = i - nheaderlines
print, n
; number of lines of data in raw acc datafile

junkarr=strarr(n)
i=0
openr, lun1, datadir+infile, /get_lun
readf, lun1, junk
readf, lun1, junk
readf, lun1, junk
readf, lun1, junk
while i lt n do begin
readf, lun1, junk
junkarr(i) = junk
i=i+1
endwhile ; while i lt n do begin
free_lun, lun1

acc = dblarr(n) ; raw ay accelerometer data, ms-2
; Times are given in file as UTC(?) YY/DDD-HH:MM:SS.SSS
; eg 01/339-04:48:49.454
; split into hours, minutes, integer seconds, and fractions of a second
; ignore date
hr = acc
mint = acc
fullsec = acc
fracsec = acc
i=0
while i lt n do begin
acc(i) = double( strmid(junkarr(i), 32, 14))
junk1 = strmid(junkarr(i), 7, 12)
junk2 = str_sep(junk1, ':')
hr(i) = junk2(0)
mint(i) = junk2(1)
xx = str_sep(junk2(2), '.')
fullsec(i) = xx(0)
fracsec(i) = double( '0.' + xx(1) ) 
i=i+1
endwhile ; while i lt n do begin

;;;
; Finish reading in ACC data for one aerobraking pass
;;;

;;;
; Start time-tag datapoints by "seconds after periapsis"
;;;
; Split time of periapsis for this orbit
; into hour, min, sec, frac sec as well
xx = utcdatetime(ii)
yy = str_sep(xx, ' ')
zz = str_sep(yy(1), ':')
pp = str_sep(zz(2), '.')

; Convert periapsis time into number of seconds after midnight
; keep sec and frac sec separate
psecaftermidnightfull = double(pp(0)) + 60. * (double(zz(1) + 60. * (double(zz(0)) ) ) )
psecaftermidnightfrac = double( '0.' + pp(1) ) 

; Convert times of all acc datapoints into same scheme
orbsecaftermidnightfull = fullsec + 60.*(mint + 60.*hr)
orbsecaftermidnightfrac = fracsec

; Ensure times of acc datapoints are continuous
; through midnight
if min(orbsecaftermidnightfull) lt 1.*60.*60. and $
 max(orbsecaftermidnightfull) ge 23.*60.*60. then begin
; orbit spans midnight...
aa = where(orbsecaftermidnightfull ge 12.*60.*60.)
bb = where(orbsecaftermidnightfull lt 12.*60.*60.)
if psecaftermidnightfull ge 23.*60.*60. then $
 orbsecaftermidnightfull(bb) = orbsecaftermidnightfull(bb) + 24.*60.*60.
if psecaftermidnightfull lt 1.*60.*60. then $
 orbsecaftermidnightfull(aa) = orbsecaftermidnightfull(aa) - 24.*60.*60.
endif ; if min(orbsecaftermidnightfull) lt 1.*60.*60. and $
; max(orbsecaftermidnightfull) ge 23.*60.*60. then begin


; Convert times of acc datapoints into
; seconds after midnight, then include
; fractions of a second as well.
timeafterperiapsis = orbsecaftermidnightfull - $
 psecaftermidnightfull
timeafterperiapsis = timeafterperiapsis + $
 orbsecaftermidnightfrac - psecaftermidnightfrac

;;;
; Finish time-tag datapoints by "seconds after periapsis"
;;;

;;;
; Start calculating orbit using a, e, i, time
; Determine r, v, lat, and lon
; as functions of time after periapsis
; also determine a nominal altitude
;;;

dt = 0.1 ; timestep in seconds
nsteps = 1e4 ; we don't need the entire ellipse
tmax = dt * nsteps

rell = dblarr(nsteps) ; radial distance, m
thetaell = rell ; angle, radians
lell = rell ; curve length, m
tell = rell ; time from periapsis, sec
vell = rell ; speed, m s-1

a = psma(ii)*1e3 ; semi-major axis, m
e = pecc(ii) ; eccentricity
rell(0) = a * (1.-e) 
; calculated periapsis radius
; this agrees well with Nav periapsis radius, pradius
; <5 m difference due to rounding
tell(0) = 0.
thetaell(0) = 0.
lell(0) = 0.
vell(0) = sqrt( gm/a * (1.+e)/(1.-e) ) 
; orbital speed at periapsis,
; calculated from v2/2 -GM/r = -GM/2a
; energy equation

i=1
while i lt nsteps do begin
dtheta = rell(0) * vell(0) * dt / rell(i-1)^2
; r-squared dtheta/dt = constant
; Equal area law, conservation of angular momentum
thetaell(i) = thetaell(i-1) + dtheta
rell(i) = a * (1.-e*e) / (1. + e*cos(thetaell(i)))
; New r from ellipse equation
dl = sqrt( (rell(i)-rell(i-1))^2 + rell(i-1)^2*dtheta^2 )
; dl from usual dr and dtheta relation
lell(i) = lell(i-1) + dl
vell(i) = dl/dt
tell(i) = tell(i-1) + dt
i=i+1
endwhile ; while i lt nsteps do begin

; Convert some variables from m and ms-1 to km and kms-1
r_km = rell / 1e3
v_kmpersec = vell / 1e3
l_km = lell /1e3

;;;
; Start finding latitude and longitude
; around orbit - note that equations 
; may need to be changed if applied to any SH data
;;;

; Will need to use periapsis conditions
; to find conditions at ascending node
; that are needed in the standard formulae
perilat = plat(ii) ; degN
perilon = plon(ii) ; degE
periinc = pinc(ii) ; deg

; Whether ODY is heading north
; or south at periapsis is necessary
; information for finding the angle
; turned through between the ascending 
; node and periapsis
; i=124, 125, perilat = inc exactly
; so it's neither n-bound or s-bound
; in this sense at periapsis
; Angle turned through must be 90deg

iturnlo = 124
iturnhi = 125
if orb(ii) lt iturnlo then dirn='s'
if orb(ii) gt iturnhi then dirn='n'

;;;
; Orbit equations from my undergrad lecture notes
; Rees, Remote Sensing
;;;

;;;
; Find angle turned through to go from
; ascending node to periapsis and 
; longitude at ascending node
;;;
if orb(ii) ne iturnlo and orb(ii) ne iturnhi then $
periphi = 180./!pi * asin( sin(perilat*!pi/180.)/sin(periinc*!pi/180.) )
; asin returns -90 to 90
; periphi is angle turned through to go from
; ascending node to periapsis
; since we're in NH, there's no ambiguity in asin
if orb(ii) eq iturnlo or orb(ii) eq iturnhi then periphi = 90.
if dirn eq 's' then periphi = 180.-periphi
if dirn eq 'n' then periphi = periphi
print, periphi

if orb(ii) ne iturnlo and orb(ii) ne iturnhi then $
lml0 = 180./!pi * acos( cos(periphi*!pi/180.) / cos(perilat*!pi/180.) )
; acos returns 0 to 180
; lml0 is longitude at periapsis minus longitude at ascending node, degE
if periphi lt 90. and lml0 lt 90. then lml0 = 360. - lml0
if periphi gt 90. and lml0 gt 90. then lml0 = 360. - lml0
; periphi and lml0 must be in opposite quadrants
; thus resolving acos ambiguity
if orb(ii) eq iturnlo or orb(ii) eq iturnhi then lml0 = 270.
l0 = perilon - lml0
; longitude at ascending node, degE
print, l0

; r is symmetric about periapsis
; lat, lon are not
; so need to perform this bit twice
; inbound and outbound
; use subscripts m and p to indicate minus or plus
; in allphi
jmax=2
jsign =[-1,1]
j=0
while j lt jmax do begin

allphi = periphi + jsign(j) * thetaell*180./!pi
; phi = angle turned through to go from
; ascending node to specific position on ellipse
; allphi = phi for all positions on elliptical arc, deg
sinlatell = sin(allphi*!pi/180.) * sin(periinc*!pi/180.)
latell = 180./!pi * asin(sinlatell)
; latell = latitude for all positions on elliptical arc, degN
coslellml0 = cos(allphi*!pi/180.) / cos(latell*!pi/180.) 
lellml0 = 180./!pi * acos(coslellml0)
; acos returns 0 to 180
; lellml0 = longitude for all positions on elliptical arc
; minus longitude of ascending node, degE

aa = where(allphi lt 90. and lellml0 lt 90.)
if aa(0) ne -1 then lellml0(aa) = 360.-lellml0(aa)
aa = where(allphi gt 90. and lellml0 gt 90.)
if aa(0) ne -1 then lellml0(aa) = 360.-lellml0(aa)
; resolve acos ambiguity again
lell = lellml0 + l0
; lell= longitude for all positions on elliptical arc, degE
aa = where(lell lt 0.)
if aa(0) ne -1 then while lell(aa(0)) lt 0. do lell(aa) = lell(aa) + 360.
aa = where(lell gt 360.)
if aa(0) ne -1 then while lell(aa(0)) gt 360. do lell(aa) = lell(aa) - 360.
; Make sure lell is in 0-360 range

if j eq 0 then begin
latellm = latell
lellm = lell
endif
if j eq 1 then begin
latellp = latell
lellp = lell
endif

j=j+1
endwhile ; while j lt jmax do begin

; m is always inbound, p is always outbound
latelli = latellm
latello = latellp
lelli = lellm
lello = lellp

; Join inbound and outbound together
aa = where(timeafterperiapsis gt 0.)
to = timeafterperiapsis(aa)
lato = interpol(latello, tell, to)
lono = interpol(lello, tell, to)
aa = where(timeafterperiapsis lt 0.)
ti = (-1.)*timeafterperiapsis(aa)
lati = interpol(latelli, tell, ti)
loni = interpol(lelli, tell, ti)

latonpass = [lati, lato] 
; latitude of all positions on elliptical arc, degN
lononpass = [loni, lono]
; longitude of all positions on elliptical arc, degN

; Use times of acc datapoints and calculated orbit
; to get position and speed at each datapoint
abst = abs(timeafterperiapsis)
ronpass = interpol(r_km, tell, abst) ; km
vonpass = interpol(v_kmpersec, tell, abst) ; km s-1


; Use Tolson areoid to find
; radius of areoid beneath orbital path
xx = lononpass
aa = where(lononpass gt 180.)
if aa(0) ne -1 then xx(aa) = xx(aa) - 360.
yy = xx + 180.
zz = latonpass + 90.
requipotonpass = interpolate(equipot_radius, zz, yy)

; Nav_recontr_thru_P338.xls gives a simple areoid
; as well, if you manipulate it
; zonally symmetric
rnavareoid = pradius-palt
aa = where(plat eq max(plat))
latx = plat(aa(n_elements(aa)-1):*)
rnavareoidx = rnavareoid(aa(n_elements(aa)-1):*)
rnavareoidonpass = interpol(rnavareoidx, latx, latonpass)

;;; 
; Now calculate altitude along pass
; from radius along pass
; Several areoid choices, none are good
; If you know of a better areoid, use it!
;;;

;zonpass = ronpass - requipotonpass ; km
; Tolson areoid

;zonpass = ronpass - rnavareoidonpass ; km
; NAV crude areoid

zonpass = ronpass - pradius(ii) + palt(ii) ; km
; Best consistency with ODYA_0001 results on PDS
; if I define altitude as radius - periapsis radius + periapsis altitude
; based on Nav_recontr_thru_P338.xls information
; This is locally, but not globally accurate

;;;
; Finish calculating orbit using a, e, time
;;;



;;;
; Convert acc into rho
;;;

; Use unsmoothed data, 7 point mean, and 39 point mean
acc1 = acc
acc7 = smooth(acc, 7, /edge_truncate)
acc39 = smooth(acc, 39, /edge_truncate)
; several types of accelerations (m s-2)

; Ratio of density to acceleration from drag equation
; note that vonpass is in km
; units are kg m-3 (m s-2)-1
rhooveracc = 2. * pmass(ii) / ( area * aerocoeff * (vonpass*1e3)^2 )

; Densities, kg m-3, based on unsmoothed acc, 
; 7 point mean, and 39 point mean
rho1 = acc1 * rhooveracc
rho7 = acc7 * rhooveracc
rho39 = acc39 * rhooveracc

;;;
; Calculate uncertainties in acc and rho
;;;

; I would like to use pre-pass and post-pass
; acc data to determine noise levels, drift, and offset
; but some orbits have thruster activity after periapsis
; This ruins the post-pass data, so I just use pre-pass data
; Can't determine drift without post-pass data
; Don't determine offset, because I haven't worked out a way
; to deal with pathological cases that have acc measurements
; beginning when acc is already detecting atmosphere,
; plus I don't want to start correcting acc data yet
; Define noise as the standard deviation within
; a subset of pre-pass data

; Use arbitrary limits for noise subset
; Don't want to be close to start of time series,
; due to wake-up effects
; Don't want to be within atmosphere either -
; which is tricky for later, longer passes
istartnoise = 100
iendnoise = 200
junk = moment(acc1(istartnoise:iendnoise))
; Use this value as acc measurement uncertainty for
; all data points on this orbit

if ii eq 27 then junk = moment(acc1(20:120))
if ii eq 279 then junk = moment(acc1(500:600))
; acc data starts very close to atm detection
; timing must have been wrong

noise = sqrt(junk(1)) 
; noise in acc data for this profile, m s-2

; This is the zero offset, if trusted and needed
; zerolevel = junk(0) ; m s-2

pnoise(ii) = noise

; 1-sigma measurement uncertainty in
; unsmoothed, 7 point mean, and 39 point mean
; accelerations. Reduce noise by square root
; of n for the running means.
sacc1 = noise + acc1*0. ; m s-2
sacc7 = noise / sqrt(7.) + acc7*0. ; m s-2
sacc39 = noise / sqrt(39.) + acc39*0. ; m s-2

; 1-sigma uncertainty in densities derived from
; unsmoothed, 7 point mean, and 39 point mean
; accelerations (kg m-3)
srho1 = sacc1 * rhooveracc
srho7 = sacc7 * rhooveracc
srho39 = sacc39 * rhooveracc


;;;
; Reject measurements that are less than uncertainties
; More accurately, accept only those measurements
; that are greater than the uncertainties AND
; form a continuous series all the way to periapsis
; with no bad data points within the series
;;;

t = timeafterperiapsis ; seconds
; rename using short name

numaveragedaccs = 3
; unsmoothed, 7pt, and 39 pt = 3 types
j = 0
while j lt numaveragedaccs do begin

if j eq 0 then begin
xa = acc1
xsa = sacc1
xr = rho1
xsr = srho1
endif ; if j eq 0 then begin

if j eq 1 then begin
xa = acc7
xsa = sacc7
xr = rho7
xsr = srho7
endif ; if j eq 1 then begin

if j eq 2 then begin
xa = acc39
xsa = sacc39
xr = rho39
xsr = srho39
endif ; if j eq 0 then begin

; xa = accelerations, m s-2
; xsa = uncertainties in accelerations, m s-2
; xr = densities, kg m-3
; xsr = uncertainties in densities, kg m-3
; This lets me use the same variable names for
; unsmoothed, 7pt, and 39 pt data

fullacc = dblarr(n_elements(xa))
fullsacc = fullacc
fullrho = fullacc
fullsrho = fullacc
; these arrays will be clones of xa, xsa, etc
; with bad values set to zero

aa = where(xr ge xsr)
if aa(0) eq -1 then begin
; If all measurements are less than uncertainty
; then life is easy, do nothing at all
; all arrays remain full of zeroes

endif else begin ; if aa(0) eq -1 then begin

; split in and outbound
aa = where(t ge 0.)

; I don't think there's any use for ti, ri, and zi
; and their outbound counterparts
; but I'm not going to delete them in case
; I need them later

; i indicates inbound
acci = xa(0:aa(0)) ; acceleration, m s-2
sacci = xsa(0:aa(0)) ; 1 sigma uncertainty in acceleration, m s-2
rhoi = xr(0:aa(0)) ; density, kg m-3
srhoi = xsr(0:aa(0)) ; 1 sigma uncertainty in density, kg m-3
ti = t(0:aa(0)) ; time after periapsis, s
ri = ronpass(0:aa(0)) ; radial distance, km
zi = zonpass(0:aa(0)) ; altitude, km

; o indicates outbound
acco = xa(aa(0)+1:*) ; acceleration, m s-2
sacco = xsa(aa(0)+1:*) ; 1 sigma uncertainty in acceleration, m s-2
rhoo = xr(aa(0)+1:*) ; density, kg m-3
srhoo = xsr(aa(0)+1:*) ; 1 sigma uncertainty in density, kg m-3
to = t(aa(0)+1:*) ; time after periapsis, s
ro = ronpass(aa(0)+1:*) ; radial distance, km
zo = zonpass(aa(0)+1:*) ; altitude, km

flagi=0 ; this changes from 0 to 1 
; if there are no good inbound data points

bb = where(srhoi ge rhoi)
if bb(n_elements(bb)-1)+1 eq n_elements(rhoi) then begin
; there are no good inbound datapoints
; set everything to zero and flip the flag
goodacci=[0.] ; acceleration, m s-2
goodsacci = goodacci ; 1 sigma uncertainty in acceleration, m s-2
goodrhoi = goodacci ; density, kg m-3
goodsrhoi = goodacci ; 1 sigma uncertainty in density, kg m-3
goodti = goodacci ; time after periapsis, s
goodri = goodacci ; radial distance, km
goodzi = goodacci ; altitude, km
flagi=1
endif else begin ; if bb(n_elements(bb)-1)+1 eq n_elements(rhoi) then begin
; extract the good inbound datapoints, leave flag as is
goodacci = acci(bb(n_elements(bb)-1)+1:*) ; acceleration, m s-2
goodsacci = sacci(bb(n_elements(bb)-1)+1:*) ; 1 sigma uncertainty in acceleration, m s-2
goodrhoi = rhoi(bb(n_elements(bb)-1)+1:*) ; density, kg m-3
goodsrhoi = srhoi(bb(n_elements(bb)-1)+1:*) ; 1 sigma uncertainty in density, kg m-3
goodti = ti(bb(n_elements(bb)-1)+1:*) ; time after periapsis, s
goodri = ri(bb(n_elements(bb)-1)+1:*) ; radial distance, km
goodzi = zi(bb(n_elements(bb)-1)+1:*) ; altitude, km
endelse ; if bb(n_elements(bb)-1)+1 eq n_elements(rhoi) then begin

; Repeat for outbound
flago=0
bb = where(srhoo ge rhoo)
if bb(0) eq 0 then begin
goodacco=[0.]
goodsacco = goodacco
goodrhoo = goodacco
goodsrhoo = goodacco
goodto = goodacco
goodro = goodacco
goodzo = goodacco
flago=1
endif else begin ; if bb(0) eq 0 then begin
goodacco = acco(0:bb(0)-1)
goodrhoo = rhoo(0:bb(0)-1)
goodsacco = sacco(0:bb(0)-1)
goodsrhoo = srhoo(0:bb(0)-1)
goodto = to(0:bb(0)-1)
goodro = ro(0:bb(0)-1)
goodzo = zo(0:bb(0)-1)
endelse ; if bb(0) eq 0 then begin

;;;
; Now combine inbound and outbound data back together
;;;

; Good inbound data, good outbound data
if flagi eq 0 and flago eq 0 then begin
cleanacc = [goodacci, goodacco] ; acceleration, m s-2
cleansacc = [goodsacci, goodsacco] ; 1 sigma uncertainty in acceleration, m s-2
cleanrho = [goodrhoi, goodrhoo] ; density, kg m-3
cleansrho = [goodsrhoi, goodsrhoo] ; 1 sigma uncertainty in density, kg m-3
cleant = [goodti, goodto] ; time after periapsis, s
cleanr = [goodri, goodro] ; radial distance, km
cleanz = [goodzi, goodzo] ; altitude, km
endif ; if flagi eq 0 and flago eq 0 then begin

; Bad inbound data, good outbound data
if flagi eq 1 and flago eq 0 then begin
cleanacc = [goodacco]
cleansacc = [goodsacco]
cleanrho = [goodrhoo]
cleansrho = [goodsrhoo]
cleant = [goodto]
cleanr = [goodro]
cleanz = [goodzo]
endif ; if flagi eq 1 and flago eq 0 then begin

; Good inbound data, bad outbound data
if flagi eq 0 and flago eq 1 then begin
cleanacc = [goodacci]
cleansacc = [goodsacci]
cleanrho = [goodrhoi]
cleansrho = [goodsrhoi]
cleant = [goodti]
cleanr = [goodri]
cleanz = [goodzi]
endif ; if flagi eq 0 and flago eq 1 then begin

; Bad inbound data, bad outbound data
if flagi eq 1 and flago eq 1 then begin
cleanacc = [0.]
cleansacc = cleanacc
cleanrho = cleanacc
cleansrho = cleanacc
cleant = cleanacc
cleanr = cleanacc
cleanz = cleanacc
endif ; if flagi eq 1 and flago eq 1 then begin

;;;
; I now have only the good datapoints,
; plus maybe a zero or two for bad ones
; I need to order them against the original timeseries
; with zeroes for the bad datapoints
;;;

if max(abs(cleant)) ne 0. then begin
; Nothing happens if cleant = 0, corresponding
; to bad inbound and bad outbound data
; so arrays remain full of zeroes
i=0
while i lt n_elements(t) do begin
; for each datapoint in original time series,
; do we have a good datapoint at the same time?
; [Clearly this approach is inefficient...]
aa = where(cleant eq t(i))
if aa(0) ne -1 then begin
fullacc(i) = cleanacc(aa(0)) ; acceleration, m s-2
fullsacc(i) = cleansacc(aa(0)) ; 1 sigma uncertainty in acceleration, m s-2
fullrho(i) = cleanrho(aa(0)) ; density, kg m-3
fullsrho(i) = cleansrho(aa(0)) ; 1 sigma uncertainty in density, kg m-3
; What about r, z, and t?
; They keep their values regardless of whether datapoint is good or bad
endif ; if aa(0) ne -1 then begin
i=i+1
endwhile ; while i lt n_elements(t) do begin

endif ; if max(abs(cleant)) ne 0. then begin

endelse ; if aa(0) eq -1 then begin

; File the processed acc/rho values away before
; going onto next loop around
if j eq 0 then begin
tidyacc1 = fullacc ; acceleration, m s-2
tidysacc1 = fullsacc ; 1 sigma uncertainty in acceleration, m s-2
tidyrho1 = fullrho ; density, kg m-3
tidysrho1 = fullsrho ; 1 sigma uncertainty in density, kg m-3
endif ; if j eq 0 then begin

if j eq 1 then begin
tidyacc7 = fullacc
tidysacc7 = fullsacc
tidyrho7 = fullrho
tidysrho7 = fullsrho
endif ; if j eq 1 then begin

if j eq 2 then begin
tidyacc39 = fullacc
tidysacc39 = fullsacc
tidyrho39 = fullrho
tidysrho39 = fullsrho
endif ; if j eq 2 then begin

j=j+1
endwhile ; while j lt numaveragedaccs do begin

; I don't want to print a long string of zeroes
; for the pre-pass/post-pass portion of the dataset

aa = where(tidyacc1 ne 0.)
bb = where(tidyacc7 ne 0.)
cc = where(tidyacc39 ne 0.)
ff = [aa,bb,cc]
; lists all times when one or more datatype is good
; repeats values if more than one is good

gg = where(ff ge 0.)
if gg(0) eq -1 then begin
; if absolutely no datapoints are good
; then set everything to zero
t = [0.]
ronpass = t
zonpass = t
latonpass = t
lononpass = t
tidyacc1 = t
tidysacc1 = t
tidyrho1 = t
tidysrho1 = t
tidyacc7 = t
tidysacc7 = t
tidyrho7 = t
tidysrho7 = t
tidyacc39 = t
tidysacc39 = t
tidyrho39 = t
tidysrho39 = t
endif else begin ; if gg(0) eq -1 then begin
; if some datapoints are good
; then keep all datapoints between first good
; datapoint and last good datapoint
ff = ff(gg)
dd = min(ff)
ee = max(ff)
t = t(dd:ee)
ronpass = ronpass(dd:ee)
zonpass = zonpass(dd:ee)
latonpass = latonpass(dd:ee)
lononpass = lononpass(dd:ee)
tidyacc1 = tidyacc1(dd:ee)
tidysacc1 = tidysacc1(dd:ee)
tidyrho1 = tidyrho1(dd:ee)
tidysrho1 = tidysrho1(dd:ee)
tidyacc7 = tidyacc7(dd:ee)
tidysacc7 = tidysacc7(dd:ee)
tidyrho7 = tidyrho7(dd:ee)
tidysrho7 = tidysrho7(dd:ee)
tidyacc39 = tidyacc39(dd:ee)
tidysacc39 = tidysacc39(dd:ee)
tidyrho39 = tidyrho39(dd:ee)
tidysrho39 = tidysrho39(dd:ee)
endelse ; if gg(0) eq -1 then begin

;;;
; End of messy business of selecting the
; good datapoints
;;;

;;;;;;;;;
; DENSITY PROFILES HAVE BEEN DETERMINED
;;;;;;;;;

;;;
; Track periapsis densities
;;;

aa = where( abs(t) eq min(abs(t)) )
aa = aa(0)
ptime(ii) = t(aa)
pratptime(ii) = ronpass(aa)
prho1(ii) = tidyrho1(aa)
psrho1(ii) = tidysrho1(aa)
prho7(ii) = tidyrho7(aa)
psrho7(ii) = tidysrho7(aa)
prho39(ii) = tidyrho39(aa)
psrho39(ii) = tidysrho39(aa)



;;;
; Write out density profile
;;;

outfileprofile = 'p' + strorb(ii) + 'prof.txt'

openw, lun1, outdir+outfileprofile, /get_lun
formatprof = '(E12.5,2(",",E12.5),14(",",E12.5))'
printf, lun1, ';;; start of header ;;;'
printf, lun1, '40 lines of header'
printf, lun1, strcompress(string(n_elements(t)), /remove_all) + $
 ' lines of data'
printf, lun1, '17 columns of data separated by commas'
printf, lun1, 'Sources of data are JPL server MMONT1, file ' + $
 'M01-AB/NAV/Reconstruction_Data/Nav_reconstr_thru_P0338.xls'
printf, lun1, 'and pXXXacc.dat from Jim Murphy/PDS'
printf, lun1, format='(A60, F9.1)', 'orbit number (-): ', orb(ii)
printf, lun1, format='(A60, A30)', 'periapsis UTC date/time (-): ', utcdatetime(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis altitude above unknown reference surface (km): ', palt(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis altitude above 3397.2 km sphere (km): ', prefalt(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis radius (km): ', pradius(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis latitude (degN): ', plat(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis longitude (degE): ', plon(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis LST (hrs): ', plst(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis Ls (deg): ', pls(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis vrel (km s-1): ', pvrel(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis semi-major axis (km): ', psma(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis eccentricity (-): ', pecc(ii)
printf, lun1, format='(A60, E12.5)', 'periapsis inclination (deg): ', pinc(ii)
printf, lun1, format='(A60, E12.5)', 'spacecraft mass (kg): ', pmass(ii)
printf, lun1, format='(A60, E12.5)', 'period (hrs): ', period(ii)
printf, lun1, format='(A60, E12.5)', 'noise in raw acc measurements (m s-2): ', pnoise(ii)
printf, lun1, ' 1: time after periapsis (s)'
printf, lun1, ' 2: radial distance (km)'
printf, lun1, ' 3: altitude = radius - periapsis radius + periapsis altitude (km)'
printf, lun1, ' 4: raw acceleration, nominally 1-s average (m s-2)'
printf, lun1, ' 5: 1-sigma uncertainty in raw acceleration (m s-2)'
printf, lun1, ' 6: density derived from raw acceleration (kg m-3)'
printf, lun1, ' 7: 1-sigma uncertainty in density derived from raw acceleration (kg m-3)'
printf, lun1, ' 8: 7-point running mean of raw acceleration, nominally 7-s average (m s-2)'
printf, lun1, ' 9: 1-sigma uncertainty in 7-point mean acceleration (m s-2)'
printf, lun1, '10: density derived from 7-point mean acceleration (kg m-3)'
printf, lun1, '11: 1-sigma uncertainty in density derived from 7-point mean acceleration (kg m-3)'
printf, lun1, '12: 39-point running mean of raw acceleration, nominally 39-s average (m s-2)'
printf, lun1, '13: 1-sigma uncertainty in 39-point mean acceleration (m s-2)'
printf, lun1, '14: density derived from 39-point mean acceleration (kg m-3)'
printf, lun1, '15: 1-sigma uncertainty in density derived from 39-point mean acceleration (kg m-3)'
printf, lun1, '16: latitude (degN)'
printf, lun1, '17: longitude (degE)'

printf, lun1, ';;; end of header ;;;'
i=0
while i lt n_elements(t) do begin
printf, lun1, format=formatprof, $
 t(i), ronpass(i), zonpass(i), $
 tidyacc1(i), tidysacc1(i), tidyrho1(i), tidysrho1(i), $
 tidyacc7(i), tidysacc7(i), tidyrho7(i), tidysrho7(i), $
 tidyacc39(i), tidysacc39(i), tidyrho39(i), tidysrho39(i), $
 latonpass(i), lononpass(i)
i=i+1
endwhile
free_lun, lun1


;;;
; Start generating constant altitude data products for this orbit
; VERY dependent on your areoid
; Use 39 pt mean data type only
;;;

biglist = dblarr(n_elements(z0list)*2*nrep)
; will contain all results, all altitudes for this orbit
; useful for printing out results later

rho0list = dblarr(n_elements(z0list)*2)
; fitted density at each level for this orbit, kg m-3
srho0list = rho0list ; 1-sigma uncertainty in fitted density, kg m-3
dsh0list = rho0list ; fitted density scale height, km
sdsh0list = rho0list ; 1-sigma uncertainty in fitted density scale height, km
redchisqd0list = rho0list ; reduced chi-squared for fit
nfitlist = rho0list ; number of datapoints ued in fit

lat0list = rho0list ; latitude at appropriate altitude, degN
lon0list = rho0list ; longitude at appropriate altitude, degE


aa = where(t ge 0. and tidyrho39 gt 0.)
; select outbound part of density profile
if aa(0) eq -1 then begin
; if there aren't any good datapoints, set arrays to zero
rhoout = [0.] 
srhoout = rhoout
zout = rhoout
latout = rhoout
lonout = rhoout
endif else begin ; if aa(0) eq -1 then begin
rhoout = tidyrho39(aa) ; outbound density, kg m-3
srhoout = tidysrho39(aa) ; 1-sigma uncertainty in outbound density, kg m-3
zout = zonpass(aa) ; outbound altitude, km
latout = latonpass(aa) ; outbound latitude, degN
lonout = lononpass(aa) ; outbound longitude, degE
endelse ; if aa(0) eq -1 then begin

; Repeat for inbound
aa = where(t lt 0. and tidyrho39 gt 0.)
if aa(0) eq -1 then begin
rhoin = [0.]
srhoin = rhoin
zin = rhoin
latin = rhoin
lonin = rhoin
endif else begin ; if aa(0) eq -1 then begin
rhoin = tidyrho39(aa) ; inbound density, kg m-3
srhoin = tidysrho39(aa) ; 1-sigma uncertainty in inbound density, kg m-3
zin = zonpass(aa) ; inbound altitude, km
latin = latonpass(aa) ; inbound latitude, degN
lonin = lononpass(aa) ; inbound longitude, degN
endelse ; if aa(0) eq -1 then begin

;plot, tidyrho39, zonpass, /xlog, /ynozero, $
; xrange=[1e-10,1e-7], yrange=[100,150]

j=0
while j lt n_elements(z0list) do begin
z0 = z0list(j) ; relevant altitude level, km

; This bit provides a way to do inbound
; and outbound in the same loop 
; without needing to repeat everything twice
i=0
while i lt maxi do begin
if i eq 0 then begin
r = rhoin ; density, kg m-3
sr = srhoin ; 1-sigma uncertainty in density, kg m-3
z = zin ; altitude, km
la = latin ; latitude, degN
lo = lonin ; longitude, degE
endif else begin ; if i eq 0 then begin
r = rhoout
sr = srhoout
z = zout
la = latout
lo = lonout
endelse ; if i eq 0 then begin

aa = where(abs(z-z0) lt dz)
; select datapoints that lie between z0+dz and z0-dz
if aa(0) ne -1 then begin
if max(z(aa))-z0 ge altlimit and z0-min(z(aa)) ge altlimit then begin
; only proceed if we have data over a wide enough span

lat0 = interpol(la(aa), z(aa), z0)
; latitude at appropriate altitude, degN
lon0 = interpol(lo(aa), z(aa), z0)
; longitude at appropriate altitude, degE

y = alog(r(aa)) ; ln(rho)
sy = sr(aa) / abs(r(aa)) ; uncertainty in ln(rho)
x = z(aa)-z0 ; altitude relative to z0, km

result = poly_fit(x, y, 1, chisq=chisq, covar=covar, $
 measure_errors=sy, sigma=sigma, status=status, $
 yband=yband, yerror=yerror, yfit=yfit)

;oplot, exp(yfit), x+z0, color=255

rho0 = exp(result(0)) ; fitted density, kg m-3
srho0 = sigma(0) * rho0 ; 1-sigma uncertainty in fitted density, kg m-3
dsh0 = -1./result(1) ; fitted density scale height, km
sdsh0 = dsh0^2 * sigma(1) ; 1-sigma uncertainty in fitted density scale height, km
redchisqd0 = chisq / (n_elements(x)-2) ; reduced chi-squared from fit
nfit = n_elements(x) ; number of datapoints used in fit

if status eq 0 then begin
; only store results if IDL says that the fitting procedure worked OK
rho0list(maxi*j+i) = rho0
srho0list(maxi*j+i) = srho0
dsh0list(maxi*j+i) = dsh0
sdsh0list(maxi*j+i) = sdsh0
redchisqd0list(maxi*j+i) = redchisqd0
nfitlist(maxi*j+i) = nfit
lat0list(maxi*j+i) = lat0
lon0list(maxi*j+i) = lon0
biglist(nrep*(maxi*j+i):nrep*(maxi*j+i)+(nrep-1)) = [rho0, srho0, dsh0, sdsh0, redchisqd0, nfit, lat0, lon0]
endif ; if status eq 0 then begin

endif ; if max(z(aa))-z0 ge altlimit and z0-min(z(aa)) ge altlimit then begin
endif ; if aa(0) ne -1 then begin

i=i+1
endwhile ; while i lt maxi do begin
j=j+1
endwhile ; while j lt n_elements(z0list) do begin

; Organize results before printing them out
printline = [orb(ii), $
 palt(ii), prefalt(ii), pradius(ii), $
 plat(ii), plon(ii), plst(ii), pls(ii), $
 pvrel(ii), psma(ii), pecc(ii), pinc(ii), pmass(ii), $
 period(ii), pnoise(ii), biglist]
bigaltarray(ii,*) = printline(*)

; print results into waiting and open file
printf, lun2, format=formatalt, printline(0), utcdatetime(ii), printline(1:*)

endif ; if file_test(datadir+infile) eq 1 then begin
; Does raw acc datafile exist or not?

print, ii
ii = ii+1
endwhile ; while ii lt norb do begin

free_lun, lun2
; close the alt results file

;;;;;;;;;
; CONSTANT ALTITUDE PRODUCTS HAVE BEEN DETERMINED
;;;;;;;;;

;;;
; Print out constant altitude products
; in files like in130rho.dat
;;;

formatminialt = '(F9.1,",",A30,14(",",E12.5),1(",",E12.5,",",E12.5,",",E12.5,",",E12.5,",",E12.5,",",E12.5,",",E12.5,",",E12.5))'

j=0
while j lt n_elements(z0list) do begin

i=0
while i lt maxi do begin

openw, lun1, outdir+inout(i)+strz0list(j)+'res.txt', /get_lun
printf, lun1, ';;; start of header ;;;'
printf, lun1, '33 lines of header'
printf, lun1, strcompress(string(nrawdatafiles), /remove_all) + $
 ' lines of data'
printf, lun1, '24 columns of data separated by commas'
printf, lun1, 'Sources of data are JPL server MMONT1, file ' + $
 'M01-AB/NAV/Reconstruction_Data/Nav_reconstr_thru_P0338.xls'
printf, lun1, 'and pXXXacc.dat from Jim Murphy/PDS'
printf, lun1, inout(i) + 'bound leg of aerobraking passes'
printf, lun1, 'Altitude level = ' + strz0list(j) + ' km'
printf, lun1, '1: orbit number (-)'
printf, lun1, '2: periapsis UTC date/time (-)'
printf, lun1, '3: periapsis altitude above unknown reference surface (km)'
printf, lun1, '4: periapsis altitude above 3397.2 km sphere (km)'
printf, lun1, '5: periapsis radius (km)'
printf, lun1, '6: periapsis latitude (degN)'
printf, lun1, '7: periapsis longitude (degE)'
printf, lun1, '8: periapsis LST (hrs)'
printf, lun1, '9: periapsis Ls (deg)'
printf, lun1, '10: periapsis vrel (km s-1)'
printf, lun1, '11: periapsis semi-major axis (km)'
printf, lun1, '12: periapsis eccentricity (-)'
printf, lun1, '13: periapsis inclination (deg)'
printf, lun1, '14: spacecraft mass (kg)'
printf, lun1, '15: period (hrs)'
printf, lun1, '16: noise in raw acc measurements (m s-2)'
printf, lun1, '17: fitted density (kg m-3)'
printf, lun1, '18: 1-sigma uncertainty in fitted density (kg m-3)'
printf, lun1, '19: fitted density scale height (km)'
printf, lun1, '20: 1-sigma uncertainty in fitted density scale height (km)'
printf, lun1, '21: Reduced chi-squared value for fit (-)'
printf, lun1, '22: Number of data points used in fit (-)'
printf, lun1, '23: latitude at appropriate altitude (degN)'
printf, lun1, '24: longitude at appropriate altitude (degE)'
printf, lun1, ';;; end of header ;;;'


ii=0
while ii lt norb do begin
xx = bigaltarray(ii,0:14)
yy = bigaltarray(ii,15+nrep*(maxi*j+i):15+nrep*(maxi*j+i)+nrep-1)
printline = [xx(*), yy(*)]
if bigaltarray(ii,0) ne 0 then $
 printf, lun1, format=formatminialt, printline(0), utcdatetime(ii), printline(1:*)
; don't print out orbits that have no acc data
ii=ii+1
endwhile ; while ii lt norb do begin
free_lun, lun1

i=i+1
endwhile ; while i lt maxi do begin

j=j+1
endwhile ; while j lt n_elements(z0list) do begin


;;;
; Print out periapsis products
; in perires.dat
;;;

periresfile = outdir + 'perires.txt'

aa = where(pratptime ne 0.)
ngoodperi = n_elements(aa)

openw, lun1, periresfile, /get_lun
printf, lun1, ';;; start of header ;;;'
printf, lun1, '31 lines of header'
printf, lun1, strcompress(string(ngoodperi), /remove_all) + $
 ' lines of data'
printf, lun1, '24 columns of data separated by commas'
printf, lun1, 'Sources of data are JPL server MMONT1, file ' + $
 'M01-AB/NAV/Reconstruction_Data/Nav_reconstr_thru_P0338.xls'
printf, lun1, 'and pXXXacc.dat from Jim Murphy/PDS'
printf, lun1, '1: orbit number (-)'
printf, lun1, '2: periapsis UTC date/time (-)'
printf, lun1, '3: periapsis altitude above unknown reference surface (km)'
printf, lun1, '4: periapsis altitude above 3397.2 km sphere (km)'
printf, lun1, '5: periapsis radius (km)'
printf, lun1, '6: periapsis latitude (degN)'
printf, lun1, '7: periapsis longitude (degE)'
printf, lun1, '8: periapsis LST (hrs)'
printf, lun1, '9: periapsis Ls (deg)'
printf, lun1, '10: periapsis vrel (km s-1)'
printf, lun1, '11: periapsis semi-major axis (km)'
printf, lun1, '12: periapsis eccentricity (-)'
printf, lun1, '13: periapsis inclination (deg)'
printf, lun1, '14: spacecraft mass (kg)'
printf, lun1, '15: period (hrs)'
printf, lun1, '16: noise in raw acc measurements (m s-2)'
printf, lun1, '17: time after periapsis of closest datapoint (s)'
printf, lun1, '18: radius at this closet datapoint (km)'
printf, lun1, '19: density derived from raw acceleration (kg m-3)'
printf, lun1, '20: 1-sigma uncertainty in density derived from raw acceleration (kg m-3)'
printf, lun1, '21: density derived from 7-point mean acceleration (kg m-3)'
printf, lun1, '22: 1-sigma uncertainty in density derived from 7-point mean acceleration (kg m-3)'
printf, lun1, '23: density derived from 39-point mean acceleration (kg m-3)'
printf, lun1, '24: 1-sigma uncertainty in density derived from 39-point mean acceleration (kg m-3)'
printf, lun1, ';;; end of header ;;;'

formatperi = '(F9.1,",",A30,16(",",E12.5),6(",",E12.5))'

j=0
while j lt ngoodperi do begin
i = aa(j)
printf, lun1, format=formatperi, $
orb(i), utcdatetime(i), palt(i), prefalt(i), pradius(i), $
 plat(i), plon(i), plst(i), pls(i), $
 pvrel(i), psma(i), pecc(i), pinc(i), pmass(i), period(i), $
 pnoise(i), $
 ptime(i), pratptime(i), prho1(i), psrho1(i), $
 prho7(i), psrho7(i), prho39(i), psrho39(i)
j=j+1
endwhile
free_lun, lun1

stop
end