pro readodyaccresultsformro ; Paul Withers, 2006.01.07 ; Center for Space Physics ; Read in the results of makeodyaccresultsformro.pro ; What do you want to read? ; 0: ancinfo.txt, ancilliary information ; 1: allalt.txt, the complete constant altitude dataset ; 2: in100.txt, results for all orbits at one altitude level ; 3: pXXX.txt, a specific profile ; 4: perires.txt, periapsis results ; if 2, inbound or outbound and 100, 110, 120, 130, or 140 km? ; if 3, what orbit? ; Enter your answers below ; Names and details of arrays containing data ; are given just after data are ingested from ; file within the appropriate filetype=x loop filetype = 1 ; set filetype = 0, 1, 2, 3, 4 based on above decision print, 'filetype = ', filetype filedir = 'outdir/' if filetype eq 0 then infile = filedir + 'ancinfo.txt' if filetype eq 1 then infile = filedir + 'allalt.txt' if filetype eq 2 then begin inout = 'in' ; set inout to 'in' or 'out' zlevel = '110' ; set zlevel to '100', '110', '120', '130', or '140' infile = filedir + inout + zlevel + 'res.txt' endif if filetype eq 3 then begin orbnum = '020' ; set orbnum to '076' or appropriate 3-character string infile = filedir + 'p' + orbnum + 'prof.txt' endif if filetype eq 4 then infile = 'outdir/perires.txt' print, 'infile = ', infile if file_test(infile) eq 0 then begin print, 'Your file does not exist!' stop endif ;;; ; Determine size of header ; and data table ;;; openr, lun1, infile, /get_lun junk='' readf, lun1, junk readf, lun1, junk junk1 = str_sep(junk, ' ') nheader = double(junk1(0)) ; number of lines for header readf, lun1, junk junk1 = str_sep(junk, ' ') nlines = double(junk1(0)) ; number of lines for data readf, lun1, junk junk1 = str_sep(junk, ' ') ncols = double(junk1(0)) ; number of columns for data free_lun, lun1 if filetype eq 0 then begin ancdataarr = dblarr(nlines, ncols-1) ancpdatetime = strarr(nlines) ; array containing almost all data ; plus string array for date/time ; Read in data from file openr, lun1, infile, /get_lun junk='' i=0 while i lt nheader do begin readf, lun1, junk i=i+1 endwhile i=0 while i lt nlines do begin readf, lun1, junk junk1 = str_sep(junk, ',') ancdataarr(i,0) = junk1(0) ancpdatetime(i) = strtrim(junk1(1),2) ancdataarr(i,1:*) = junk1(2:*) i=i+1 endwhile free_lun, lun1 ; ancpdatetime = periapsis UTC date/time ; ancdataarr(*,0) = orbit number (-) ; ancdataarr(*,1) = periapsis altitude above unknown reference surface (km) ; ancdataarr(*,2) = periapsis altitude above 3397.2 km sphere (km) ; ancdataarr(*,3) = periapsis radius (km) ; ancdataarr(*,4) = periapsis latitude (degN) ; ancdataarr(*,5) = periapsis longitude (degE) ; ancdataarr(*,6) = periapsis LST (hrs) ; ancdataarr(*,7) = periapsis Ls (deg) ; ancdataarr(*,8) = periapsis vrel (km s-1) ; ancdataarr(*,9) = periapsis semi-major axis (km) ; ancdataarr(*,10) = periapsis eccentricity (-) ; ancdataarr(*,11) = periapsis inclination (deg) ; ancdataarr(*,12) = spacecraft mass (kg) ; ancdataarr(*,13) = period (hrs) endif ; if filetype eq 0 if filetype eq 1 then begin allaltarr = dblarr(nlines, ncols-1) allaltpdatetime = strarr(nlines) ; array containing almost all data ; plus string array for date/time ; Read in data from file openr, lun1, infile, /get_lun junk='' i=0 while i lt nheader do begin readf, lun1, junk i=i+1 endwhile i=0 while i lt nlines do begin readf, lun1, junk junk1 = str_sep(junk, ',') allaltarr(i,0) = junk1(0) allaltpdatetime(i) = strtrim(junk1(1),2) allaltarr(i,1:*) = junk1(2:*) i=i+1 endwhile free_lun, lun1 ; allaltpdatetime = periapsis UTC date/time ; allaltarr(*,0:13) = ancdataarr(*,0:13), see filetype=0 ; allaltarr(*,14) = noise in raw acc measurement, m s-2 ; allaltarr(*,15:22) = 100 km inbound ; allaltarr(*,23:30) = 100 km outbound ; allaltarr(*,31:38) = 110 km inbound ; allaltarr(*,39:46) = 110 km outbound ; allaltarr(*,47:54) = 120 km inbound ; allaltarr(*,55:62) = 120 km outbound ; allaltarr(*,63:70) = 130 km inbound ; allaltarr(*,71:78) = 130 km outbound ; allaltarr(*,79:86) = 140 km inbound ; allaltarr(*,87:94) = 140 km outbound ; Each set of eight columns is: ; fitted density (kg m-3) ; 1-sigma uncertainty in fitted density (kg m-3) ; fitted density scale height (km) ; 1-sigma uncertainty in fitted density scale height (km) ; Reduced chi-squared value for fit (-) ; Number of data points used in fit (-) ; latitude at appropriate altitude (degN) ; longitude at appropriate altitude (degE) ;;; ; Read in odya_0001 data as well ;;; odya0001file = '~/odya_0001/data/110kmin.tab' ; filename for relevant ODYA_0001 datafile iallalt = 31 + 8*0. ; index in allaltarr for density ; at appropriate in/out and alt level ; for odya0001file ; Read in data from file i=0 junk='' openr, lun1, odya0001file, /get_lun while not eof(lun1) do begin readf, lun1, junk i=i+1 endwhile ; while not eof(lun1) do begin free_lun, lun1 nodya0001lines = i nodya0001cols = 13 ; number of lines and columns ; of data in odya0001file odya0001block = $ dblarr(nodya0001lines, nodya0001cols-1) ; big array for all data from ; odya0001file - except date/time i=0 junk = dblarr(nodya0001cols) openr, lun1, odya0001file, /get_lun while i lt nodya0001lines do begin readf, lun1, junk odya0001block(i,0) = junk(0) odya0001block(i,1:*) = junk(2:*) ; throw away date/time string in junk(1) i=i+1 endwhile free_lun, lun1 ; Quantities in odya0001block(*,i) for each i ; 0 1 2 3 4 5 ; orb, alt, lat, lon, LST, SZA, ; ---, km, degN, degE, hr, deg, ; 6 7 8 9 10 11 ; rho, srho, T, H, sH, rms ; kg km-3, kg km-3, K, km, km, kg km-3 ; Compare ODYA_0001 densities to my densities orb1 = allaltarr(*,0) orb2 = odya0001block(*,0) rho1 = allaltarr(*,iallalt)*1e9 ; kg km-3 rho2 = odya0001block(*,6) ; kg km-3 ; Move onto common reference for orbit number neworb = dindgen(338)+1 newrho1 = dblarr(n_elements(neworb)) newrho2 = newrho1 i=0 while i lt n_elements(neworb) do begin aa = where(orb1 eq neworb(i)) if aa(0) ne -1 then newrho1(i) = rho1(aa(0)) aa = where(orb2 eq neworb(i)) if aa(0) ne -1 then newrho2(i) = rho2(aa(0)) i=i+1 endwhile ; Keep only orbits that have results fo ; both sources of results aa = where(newrho1 ne 0. and newrho2 ne 0.) orb = neworb(aa) r1 = newrho1(aa) r2 = newrho2(aa) ; Is comparison any good? plot, orb, (r1-r2)/r1, psym=1 endif ; if filetype eq 1 then begin if filetype eq 2 then begin constaltarr = dblarr(nlines, ncols-1) constaltpdatetime = strarr(nlines) ; array containing almost all data ; plus string array for date/time ; Read in data from file openr, lun1, infile, /get_lun junk='' i=0 while i lt nheader do begin readf, lun1, junk i=i+1 endwhile i=0 while i lt nlines do begin readf, lun1, junk junk1 = str_sep(junk, ',') constaltarr(i,0) = junk1(0) constaltpdatetime(i) = strtrim(junk1(1),2) constaltarr(i,1:*) = junk1(2:*) i=i+1 endwhile free_lun, lun1 ; constaltpdatetime = periapsis UTC date/time ; constaltarr(*,0:14) = allaltarr(*,0:14), see filetype=1 ; constaltarr(*,15) = fitted density (kg m-3) ; constaltarr(*,16) = 1-sigma uncertainty in fitted density (kg m-3) ; constaltarr(*,17) = fitted density scale height (km) ; constaltarr(*,18) = 1-sigma uncertainty in fitted density scale height (km) ; constaltarr(*,19) = Reduced chi-squared value for fit (-) ; constaltarr(*,20) = Number of data points used in fit (-) ; constaltarr(*,21) = latitude at appropriate altitude (degN) ; constaltarr(*,22) = longitude at appropriate altitude (degE) endif ; if filetype eq 2 then begin if filetype eq 3 then begin profarr = dblarr(nlines, ncols) nusefulheader=15 profheader = dblarr(nusefulheader) profpdatetime = '' ; Store information from header ; don't just discard it ; useful size of header, ; ie nusefulheader, is set manually openr, lun1, infile, /get_lun junk='' i=0 while i lt nheader do begin readf, lun1, junk ; Ugly way to read in the useful bits of ; the header, including string for date/time if i ge 6 and i le 21 then begin if i eq 6 then $ profheader(0) = double(strmid(junk, 60, 80)) if i eq 7 then $ profpdatetime = strmid(junk, 60, 80) if i ge 8 then $ profheader(i-7) = double(strmid(junk, 60, 80)) endif ; Examine infile to find out what ; each element in profheader is profpdatetime=strtrim(profpdatetime,2) i=i+1 endwhile ; Now move beyond header to rest of data i=0 while i lt nlines do begin readf, lun1, junk junk1 = str_sep(junk, ',') profarr(i,*) = junk1(*) i=i+1 endwhile free_lun, lun1 ; profarr(*,0) = time after periapsis (s) ; profarr(*,1) = radial distance (km) ; profarr(*,2) = altitude = radius - periapsis radius + periapsis altitude (km) ; profarr(*,3) = raw acceleration, nominally 1-s average (m s-2) ; profarr(*,4) = 1-sigma uncertainty in raw acceleration (m s-2) ; profarr(*,5) = density derived from raw acceleration (kg m-3) ; profarr(*,6) = 1-sigma uncertainty in density derived from raw acceleration (kg m-3) ; profarr(*,7) = 7-point running mean of raw acceleration, nominally 7-s average (m s-2) ; profarr(*,8) = 1-sigma uncertainty in 7-point mean acceleration (m s-2) ; profarr(*,9) = density derived from 7-point mean acceleration (kg m-3) ; profarr(*,10) = 1-sigma uncertainty in density derived from 7-point mean acceleration (kg m-3) ; profarr(*,11) = 39-point running mean of raw acceleration, nominally 39-s average (m s-2) ; profarr(*,12) = 1-sigma uncertainty in 39-point mean acceleration (m s-2) ; profarr(*,13) = density derived from 39-point mean acceleration (kg m-3) ; profarr(*,14) = 1-sigma uncertainty in density derived from 39-point mean acceleration (kg m-3) ; profarr(*,15) = latitude (degN) ; profarr(*,16) = longitude (degE) endif ; if filetype eq 3 then begin if filetype eq 4 then begin peridataarr = dblarr(nlines, ncols-1) peripdatetime = strarr(nlines) ; array containing almost all data ; plus string array for date/time ; Read in data from file openr, lun1, infile, /get_lun junk='' i=0 while i lt nheader do begin readf, lun1, junk i=i+1 endwhile i=0 while i lt nlines do begin readf, lun1, junk junk1 = str_sep(junk, ',') peridataarr(i,0) = junk1(0) peripdatetime(i) = strtrim(junk1(1),2) peridataarr(i,1:*) = junk1(2:*) i=i+1 endwhile free_lun, lun1 ; peridataarr(*,0) = orbit number (-) ; peridataarr(*,1) = periapsis altitude above unknown reference surface (km) ; peridataarr(*,2) = periapsis altitude above 3397.2 km sphere (km) ; peridataarr(*,3) = periapsis radius (km) ; peridataarr(*,4) = periapsis latitude (degN) ; peridataarr(*,5) = periapsis longitude (degE) ; peridataarr(*,6) = periapsis LST (hrs) ; peridataarr(*,7) = periapsis Ls (deg) ; peridataarr(*,8) = periapsis vrel (km s-1) ; peridataarr(*,9) = periapsis semi-major axis (km) ; peridataarr(*,10) = periapsis eccentricity (-) ; peridataarr(*,11) = periapsis inclination (deg) ; peridataarr(*,12) = spacecraft mass (kg) ; peridataarr(*,13) = period (hrs) ; peridataarr(*,14) = noise in raw acc measurements (m s-2) ; peridataarr(*,15) = time after periapsis of closest datapoint (s) ; peridataarr(*,16) = radius at this closet datapoint (km) ; peridataarr(*,17) = density derived from raw acceleration (kg m-3) ; peridataarr(*,18) = 1-sigma uncertainty in density derived from raw acceleration (kg m-3) ; peridataarr(*,19) = density derived from 7-point mean acceleration (kg m-3) ; peridataarr(*,20) = 1-sigma uncertainty in density derived from 7-point mean acceleration (kg m-3) ; peridataarr(*,21) = density derived from 39-point mean acceleration (kg m-3) ; peridataarr(*,22) = 1-sigma uncertainty in density derived from 39-point mean acceleration (kg m-3) ;;; ; Read ODY AAG QLR archives ;;; odyaagqlrdir = 'odyaagrework/data/' nqlr = 0 openr, lun1, odyaagqlrdir + 'num.dat', /get_lun readf, lun1, nqlr free_lun, lun1 odyaagqlrfiles = ['orb', 'pht', 'plat', 'plon', $ 'plst', 'period', 'prho', 'sigma_prho'] + '.dat' ; list of the ODY AAG QLR files to be ingested nodyaagqlrfiles = n_elements(odyaagqlrfiles) odyaagqlrarray = dblarr(nqlr, nodyaagqlrfiles) ; entries odyaagqlrarray(*,i) correspond ; to data type in odyaagqlrfiles(i) ; 0: orbit number (-) ; 1: periapsis altitude (km) ; 2: periapsis latitude (degN) ; 3: periapsis longitude (degE) ; 4: periapsis LST (hrs) ; 5: period (hrs) ; 6: periapsis density (kg km-3) ; 7: 1-sigma uncertainty in periapsis density (kg km-3) junkarray = dblarr(nqlr) i=0 while i lt nodyaagqlrfiles do begin openr, lun1, odyaagqlrdir + odyaagqlrfiles(i) readf, lun1, junkarray odyaagqlrarray(*,i) = junkarray free_lun, lun1 i=i+1 endwhile aa = where(odyaagqlrarray(*,3) lt 0.) if aa(0) ne -1 then $ odyaagqlrarray(aa,3) = odyaagqlrarray(aa,3)+360. ; make longitude range 0-360, not -180 to +360 ; Compare periapsis densities orb1 = peridataarr(*,0) orb2 = odyaagqlrarray(*,0) rho1 = peridataarr(*,20)*1e9 ; kg km-3 rho2 = odyaagqlrarray(*,6) ; kg km-3 neworb = dindgen(338)+1 newrho1 = dblarr(338) newrho2 = newrho1 i=0 while i lt n_elements(neworb) do begin aa = where(orb1 eq neworb(i)) if aa(0) ne -1 then newrho1(i) = rho1(aa(0)) aa = where(orb2 eq neworb(i)) if aa(0) ne -1 then newrho2(i) = rho2(aa(0)) i=i+1 endwhile aa = where(newrho1 ne 0. and newrho2 ne 0.) orb = neworb(aa) r1 = newrho1(aa) r2 = newrho2(aa) plot, orb, (r1-r2)/r1, psym=1 xx = r1/r2-1. junk = moment(xx) print, junk(0), sqrt(junk(1)) endif ; if filetype eq 4 then begin stop end