pro orbs4and6 ; Paul Withers, 2007.08.14 ; Center for Space Physics, Boston University ; Use SPICE to process ODY ACC data spawn, 'date', startdate forward_function getrho3, getpt ; file:///deimos00/spice/icy/doc/html/cspice/index.html ; file:///deimos00/spice/icy/doc/html/icy/index.html ;;; Begin - Ingest MOLA areoid ;;; dlat = 0.25 ; degrees per MOLA pixel dlon = dlat fulllat = 180. ; range in degrees of latitude fulllon = 360. ; range in degrees of longitude nmolalat = fulllat/dlat nmolalon = fulllon/dlon junk=intarr(nmolalon,nmolalat) ; ftp://pds-geosciences.wustl.edu/mgs-m-mola-5-megdr-l3-v1/ ; mgsl_300x/meg004/mega90n000cb.img openr, lun1, 'mega90n000cb.img', /get_lun ; Areoid radius at the center of the 0.25 by 0.25 degree area, ; after subtracting 3,396,000 meters. Units are m readu, lun1, junk free_lun, lun1 molaareoid = swap_endian(junk) molaareoid = double(molaareoid) rrefareoid = 3396000L ; areoid offset in metres molalat = dblarr(nmolalon,nmolalat) ; lon, lat molalon = molalat linelat = dindgen(nmolalat) * dlat + dlat/2. - fulllat/2. linelat= reverse(linelat) linelon = dindgen(nmolalon) * dlon + dlon/2. i=0 while i lt nmolalon do begin molalat(i,*) = linelat(*) molalon(i,*) = linelon(i) i=i+1 endwhile i=0 while i lt nmolalat do begin molalat(*,i) = linelat(i) molalon(*,i) = linelon(*) i=i+1 endwhile ;;; End - Ingest MOLA areoid ;;; minimolalat = dblarr(180,360) minimolalon = minimolalat minimolaz = minimolalat i=0 while i lt 180 do begin j=0 while j lt 360 do begin minimolaz(i,j) = mean(molaareoid(4.*j:4.*j+3,4.*i:4.*i+3) ) minimolalat(i,j) = mean(molalat(4.*j:4.*j+3,4.*i:4.*i+3) ) minimolalon(i,j) = mean(molalon(4.*j:4.*j+3,4.*i:4.*i+3) ) ;minimolalat(i,j) = i-89.5 ;minimolalon(i,j) = j+0.5 ;aa = where(molalat ge minimolalat(i,j)-0.5 and $ ; molalat lt minimolalat(i,j)+0.5 and $ ; molalon ge minimolalon(i,j)-0.5 and $ ; molalon lt minimolalon(i,j)+0.5 ) ;if aa(0) eq -1 then stop ;minimolaz(i,j) = mean(molaareoid(aa)) j=j+1 endwhile print,i i=i+1 endwhile ;;; Begin - Furnish necessary SPICE files ;;; cspice_furnsh, 'naif0008.tls' ; leap seconds cspice_furnsh, 'de414.bsp' ; SPK for planets cspice_furnsh, 'orb1_sclkscet_00129.tsc' ; SCLK/time conversions cspice_furnsh, 'pck00008.tpc' ; planetary rotational states cspice_furnsh, 'm01_ab.bsp' ; ODY state (position/velocity) cspice_furnsh, 'm01_v28.tf' ; definition of various ODY sct frames cspice_furnsh, 'm01_sc_ab0110.bc' ; ODY orientation and ang velocity for YYYY=2001, MM=10 cspice_furnsh, 'm01_sc_ab0111.bc' ; 2001, 11 cspice_furnsh, 'm01_sc_ab0112.bc' ; 2001, 12 cspice_furnsh, 'm01_sc_ab0201.bc' ; 2002, 01 cspice_furnsh, 'm01_sc_ab0202.bc' ; 2002, 02 ;;; End - Furnish necessary SPICE files ;;; ;;; Begin - Set some SPICE constants ;;; sctint = -53L sctstr = '-53' ; ODY = -53 ; deimos:/deimos00/spice/icy/exe> brief ../odypsspice/m01_ab.bsp ssbint = 0L ssbstr = '0' ; Solar system barycentre sunint = 10L sunstr = '10' ; The Sun marsint = 499L marsstr = '499' ; Mars frame='j2000' ; Use J2000 frame as default choice abcorr='none' ; what type of aberration correction? method='intercept' ; needed to find subsolar point rmars = double(3400.) ; approx radius of Mars, km ;;; End - Set some SPICE constants ;;; ;;; Begin - Identify available lo and hi rate files ;;; spawn, 'ls -1 ../rawdata/', filearrlo spawn, 'ls -1 ../exmatlab/*acc*', filearrhi ;filearrlo = filearrlo(1:*) ;; discard P004 because no useful density profiles obtained ; No, keep it for this special case orbstrlo = strmid(filearrlo,1,3) orbstrhi = strmid(filearrhi,13,3) orbnumlo = double(orbstrlo) orbnumhi = double(orbstrhi) nfileslo = n_elements(filearrlo) nfileshi = n_elements(filearrhi) junk = [orbnumlo, orbnumhi] orbnumboth = junk[uniq(junk, sort(junk))] norbboth = n_elements(orbnumboth) nfilesboth = norbboth filearrhiboth = strarr(norbboth) filearrloboth = filearrhiboth orbstrhiboth = filearrhiboth orbstrloboth = filearrhiboth hiexists = dblarr(norbboth) loexists = hiexists orbnumhiboth = hiexists orbnumloboth = hiexists i=0 while i lt n_elements(orbnumboth) do begin aa = where(orbnumhi eq orbnumboth(i)) if aa(0) ne -1 and orbnumboth(i) ne 137. then begin ; 4s of hirate P137 data exists, but is not useful hiexists(i) = 1. filearrhiboth(i) = filearrhi(aa(0)) orbstrhiboth(i) = orbstrhi(aa(0)) orbnumhiboth(i) = orbnumhi(aa(0)) endif aa = where(orbnumlo eq orbnumboth(i)) if aa(0) ne -1 then begin loexists(i) = 1. filearrloboth(i) = filearrlo(aa(0)) orbstrloboth(i) = orbstrlo(aa(0)) orbnumloboth(i) = orbnumlo(aa(0)) endif i=i+1 endwhile ;;; End - Identify available lo and hi rate files ;;; ;;; ;;; Begin - Loop for individual ACC file ;;; ;;; j=0 while j lt 2 do begin ; P004 and P006 only if hiexists(j) eq 1 then accrate = 'hi' if hiexists(j) ne 1 and loexists(j) eq 1 then accrate = 'lo' ;;; Begin - Ingest the ACC file ;;; if accrate eq 'lo' then begin orbstr = orbstrloboth(j) ;orbstr = strmid(filearr(j),1,3) ; Three character string for the orbit number ; Likely to be superceded by a search in the file listing at = strarr(1) ; String for measurement time, eg 02/008-09:05:20.087 ; at = UTC time ; labelled as "SCET 3" in files ; If assumed to be UTC, then amax for P076 ; occurs close to periapsis ; so can't be ET, which would introduce a ~60s offset ax = dblarr(1) ; Array for measured x-axis acceleration, m/s2 ay = ax ; Array for measured y-axis acceleration, m/s2 az = ax ; Array for measured z-axis acceleration, m/s2 junk='' accfile = '../rawdata/' + filearrloboth(j) ; Data file ;accfile = 'rawdata/' + filearr(j) ; Data file openr, lun1, accfile, /get_lun readf, lun1, junk readf, lun1, junk readf, lun1, junk readf, lun1, junk ; Read past the four lines of header i=0 while not eof(lun1) do begin readf, lun1, junk junkarr = str_sep(strtrim(strcompress(junk)), ' ') at = [at, junkarr(0)] ax = [ax, double(junkarr(1))] ay = [ay, double(junkarr(2))] az = [az, double(junkarr(3))] i=i+1 endwhile free_lun, lun1 ; Discard the place-holding first element in some arrays at = at(1:*) ax = ax(1:*) ay = ay(1:*) az = az(1:*) nacc = n_elements(at) ; number of data points ; now fudge time into SPICE-friendly format newat = strarr(nacc) ; Measurement times, string ; newat = UTC time, SPICE recognizes the format ; automatically (the central T) ; Since at is UTC as well, just need to reorder string a bit i=0 while i lt nacc do begin junkarr1 = str_sep(at(i), '-') junkarr2 = str_sep(junkarr1(0), '/') newat(i) = '20' + junkarr2(0) + '-' + junkarr2(1) + $ 'T' + junkarr1(1) i=i+1 endwhile endif ; if rate eq 'lo' then begin if accrate eq 'hi' then begin orbstr = orbstrhiboth(j) ;orbstr = strmid(filearr(j),10,3) ; Three character string for the orbit number ; Likely to be superceded by a search in the file listing ax = dblarr(1) ; Array for measured x-axis acceleration, m/s2 ay = ax ; Array for measured y-axis acceleration, m/s2 az = ax ; Array for measured z-axis acceleration, m/s2 junk='' accfile = filearrhiboth(j) ; Data file ;accfile = filearr(j) ; Data file openr, lun1, accfile, /get_lun i=0 while not eof(lun1) do begin readf, lun1, junk junkarr = str_sep(junk, ',') ax = [ax, double(junkarr(0))] ay = [ay, double(junkarr(1))] az = [az, double(junkarr(2))] i=i+1 endwhile free_lun, lun1 ax = ax(1:*) ay = ay(1:*) az = az(1:*) nacc = n_elements(ax) ; number of data points blah=0 if blah ne 0 then begin doyfile = '../exmatlab/p' + orbstr + 'doy.txt' newat = strarr(nacc) ; Measurement times, string ; Is this ET or UTC ??????? Not known ; There are more times than acc measurements ; also more times than J2000 times ; A mystery openr, lun1, doyfile, /get_lun i=0 ;while not eof(lun1) do begin while i lt nacc do begin readf, lun1, junk junkarr = str_sep(junk, ',') yyyy = junkarr(0) doy = junkarr(1) hh = junkarr(2) mm = junkarr(3) ss = junkarr(4) xx = str_sep(ss, '.') while strlen(doy) lt 3 do doy = '0'+doy while strlen(hh) lt 2 do hh = '0'+hh while strlen(mm) lt 2 do mm = '0'+mm while strlen(xx(0)) lt 2 do xx(0) = '0'+xx(0) while strlen(xx(1)) lt 3 do xx(1) = xx(1)+'0' ss = xx(0)+'.'+xx(1) newat(i) = yyyy+'-'+doy+'T'+hh+':' + $ mm+':'+ss ;while strlen(newat(i)) lt 21 do newat(i) = newat(i) + '0' if strlen(newat(i)) ne 21 then stop ;scet = '2002-010T19:43:21.61' i=i+1 endwhile free_lun, lun1 oldat = newat endif ; if blah ne 0 tfile = '../exmatlab/p' + orbstr + 't.txt' newj2000t = dblarr(nacc) ; Measurement times, string ; newj2000t = UTC seconds since J2000 ; Found by looking at max(ay) vs periapsis for P076 ; if this is ET, then max(ay) is offset by ~60s ; from periapsis ; if this is UTC, then max(ay) is close to periapsis ; consistent with Tolson figure and expectations openr, lun1, tfile, /get_lun readf, lun1, newj2000t free_lun, lun1 cspice_deltet, newj2000t, 'UTC', delta shiftedj2000t = newj2000t + delta ; shiftedj2000t = ET seconds since J2000 cspice_et2utc, shiftedj2000t, 'ISOD', 3, newat ; newat = UTC time, SPICE recognizes the format ; automatically (the central T) ; Only newat is used in rest of the code endif ;;; End - Ingest the ACC file ;;; ayold=ay ;;; Begin - Make storage arrays ;;; rarray = dblarr(nacc) ; Radial distance between sct and Mars CM, km latarray = rarray ; Planetocentric latitude of sct, degN lonarray = rarray ; Planetocentric longitude of sct, degE lstarray = rarray ; LST of sct, hrs szaarray = rarray ; solar zenith angle of sct, deg lsarray = rarray ; Ls at this time, deg vrelxarray = rarray ; Sct-atmosphere relative velocity vector, component along sct +x axis, km/s vrelyarray = rarray ; component along sct +y axis, km/s vrelzarray = rarray ; component along sct +z axis, km/s vrelmagarray = rarray ; magnitude of sct-atmosphere relative velocity vector, km/s etarray = rarray foundarray = rarray wxarray = rarray wyarray = rarray wzarray = rarray scetarray = newat ; scetarray = UTC time, SPICE recognizes the format ; automatically (the central T) ;;; End - Make storage arrays ;;; ;;; ;;; Begin - The big SPICE loop ;;; ;;; i=0 while i lt nacc do begin ;scet = '2002-010T19:43:21.61' ; shows a format that works scet = newat(i) ; scet = UTC time, SPICE recognizes the format ; automatically (the central T) ; Convert to ET cspice_str2et, scet, et ; et = ET seconds since J2000 ;print, FORMAT='(F20.8)', et ; Note how accuracy of time is much better than 1 second. ; Find position vectors of Mars, sct, and Sun in J2000 frame ; with origin at SSB, units are km cspice_spkgps, marsint, et, frame, ssbint, marsssb, ltimemarsssb cspice_spkgps, sctint, et, frame, ssbint, sctssb, ltimesctssb cspice_spkgps, sunint, et, frame, ssbint, sunssb, ltimesunssb ; marsssb, sctssb, sunssb are position vectors wrt ssb sctmars = sctssb - marsssb ; position vector of sct wrt Mars in J2000 frame, km cspice_tipbod, frame, marsint, et, tipm ; prepare to transform position ; from J2000 frame to Mars body-fixed frame cspice_mxv, tipm, sctmars, bfsct ; bfsct is sct position wrt Mars in Mars body-fixed cartesian frame, units km cspice_reclat, bfsct, sctradius, sctlongitude, sctlatitude sctlongitude = sctlongitude * 180./!pi sctlatitude = sctlatitude * 180./!pi ; Obtain radius (km), planetocentric latitude (degN), ; and planetocentric longitude (degE) of sct wrt Mars ls = cspice_lspcn(marsstr, et, abcorr) * 180./!pi ; Find Ls for Mars at this time, deg cspice_et2lst, et, marsint, !pi/180.*sctlongitude, 'planetocentric', $ hr, mn, sc, time, ampm lst = hr + (mn + sc/double(60.))/double(60.) ; Find LST of sct at this time, hrs marssun = marsssb - sunssb ; position vector of Mars wrt Sun in J2000 frame, km sza = cspice_vsep(double(-1.)*marssun, sctmars)*double(180.)/!pi ; Use Mars-Sun and sct-Mars position vectors to find SZA, deg ;print, sctradius, sctlongitude, sctlatitude ;print, ls, lst, sza cspice_spkgeo, sctint, et, frame, marsint, sctstatemars, ltimesctmars sctvelmars = sctstatemars(3:*) ; Velocity of sct wrt Mars in J2000 frame, km/s marsnorthpole_bf = [0.,0.,1.] * rmars ; Cartesian coordinates of north pole ; of Mars wrt Mars CM in Mars body-fixed coordinates, km cspice_pxform, 'iau_mars', frame, et, rotate cspice_mxv, rotate, marsnorthpole_bf, marsnorthpole_j2000 ; Cartesian coordinates of north pole ; of Mars wrt Mars CM in J2000, km marsnorthpole_j2000_norm = cspice_vnorm(marsnorthpole_j2000) marsnorthpole_j2000_unitvec = marsnorthpole_j2000 / marsnorthpole_j2000_norm ; Cartesian unit vector for position of north pole ; of Mars wrt Mars CM in J2000 omega = double(2.*!pi / (24.6229 *60.*60.)) ; print, '(rad s-1) omega = ', omega ; http://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html ; Mars sidereal rot rate in rad s-1 ; Add in a +/- 100 m/s zonal wind ;omega = omega - 100./3400e3 * cos(sctlatitude*!pi/180.) cspice_vcrss, omega * marsnorthpole_j2000_unitvec, sctmars, vwind ; Velocity of Mars atmosphere at sct wrt Mars CM due to solid body rotation ; in J2000 frame, km/s ; Add in a +/- 100 m/s meridional wind ;zonalvec = vwind ;radialvec = sctmars ;cspice_vcrss, zonalvec, radialvec, meridionalvec ;umeridionalvec = meridionalvec/norm(meridionalvec) ;vmeridionalwind = 0.1 * umeridionalvec ; 100 m/s = 0.1 km/s ;vwind = vwind - vmeridionalwind vrel = sctvelmars - vwind ; Velocity of sct wrt atmosphere in J2000 frame, km/s vrelmag = norm(vrel) ; Speed of sct wrt atmosphere, km/s ; Calculated here for all orbits ; Later calculated again for all orbits except P058, P059 odysctframeidstr = 'M01_SPACECRAFT' ; SPICE label for spacecraft reference frame odysctframeidint = -53000L if orbstr ne '058' and orbstr ne '059' then begin ; Don't do this for orbits P058, P059 because ; C-kernels don't exist, so no attitude (v wrt atm, alpha) ; Fudge a solution instead in else loop cspice_pxform, frame, odysctframeidstr, et, rotate ; Prepare to transform velocity vector of spacecraft wrt atmosphere ; from J2000 to sct frame ; from J2000 ; to ODY sct frame cspice_mxv, rotate, vrel, vrelinsctframe ; vrelinsctframe is velocity vector of spacecraft wrt atmosphere in sct frame cspice_unorm, vrelinsctframe, unitvrelinsctframe, vrelmag ; Unit vector and magnitude for vrelinsctframe cspice_sce2s, sctint, et, sclkch ; Convert ET into character representation of ; SCLK time (not necessarily integral number of ticks) cspice_scencd, sctint, sclkch, sclkdp ; Convert string/character SCLK time into encoded double precision SCLK time ;cspice_sce2t, sctint, et, sclkdp & print, format='(F30.8)', sclkdp ;cspice_sce2t, sctint, et+1, sclkdp & print, format='(F30.8)', sclkdp ; Determine that 1 second of ET time = 256 SCLK ticks tol = 128. ; tolerance = 128 ticks = 0.5s ; used to give angular rate data some slop, I think cspice_ckgpav, odysctframeidint, sclkdp, tol*10., frame, $ cmat, av, clkout, found cspice_pxform, frame, odysctframeidstr, et, rotate ; Prepare to transform angular velocity vector of spacecraft wrt atmosphere ; from J2000 ; to ODY sct frame cspice_mxv, rotate, av, avinsctframe ; avinsctframe is angular velocity vector of spacecraft in sct frame av = avinsctframe endif else begin ; This else loop is for P058, P059 only ; Make velocity-related values zero at present ; Nominal stable attitude will be used to correct them later ; vrelmag was calculated for all orbits before the ; "if orbstr ne '058' and orbstr ne '059' then begin" loop began ; so no need to change it vrelinsctframe = dblarr(3) ; this is changed later found=0 ; this is not changed later av = dblarr(3) ; this is not changed later and leads to angular velocity=0 ; hence angular velocity corrections are not possible endelse ; What is orientation of sct frame? ; http://www.igpp.ucla.edu/cache1/ODMA_0001/CATALOG/INST.CAT ; http://pds-geosciences.wustl.edu/missions/odyssey/catalog/inst_marie.txt ; aerobraking velocity is nominally in -y direction ; http://pdsproto.jpl.nasa.gov/catalog/host/Results.CFM?resultsselbox=ODY ; Shows +y axis originating in sct, then going out through solar arrays ; this document doesn't state aerobraking direction ; Takashima and Wilmoth AIAA2002-4809 Figure 1 supports the above ; Tolson JSR 2005 seems to have -y and -z labelled as +ve ; OK, proceed etarray(i) = et rarray(i) = sctradius ; km latarray(i) = sctlatitude ; degN lonarray(i) = sctlongitude ; degE lstarray(i) = lst ; hrs szaarray(i) = sza ; deg lsarray(i) = ls ; deg vrelxarray(i) = vrelinsctframe(0) ; km/s vrelyarray(i) = vrelinsctframe(1) ; km/s vrelzarray(i) = vrelinsctframe(2) ; km/s vrelmagarray(i) = vrelmag ; km/s foundarray(i) = found wxarray(i) = av(0) wyarray(i) = av(1) wzarray(i) = av(2) ; So wxarray, etc, are zero for P058, P059 i=i+1 endwhile ;;; ;;; End - The big SPICE loop ;;; ;;; ;;; Begin - Find areoid and tfromperi ;;; aa = where(lonarray ge fulllon) if aa(0) ne -1 then lonarray(aa) = lonarray(aa)-fulllon aa = where(lonarray lt 0.) if aa(0) ne -1 then lonarray(aa) = lonarray(aa)+fulllon ; lon = 0 to 359.99999 ; lat = -90 to 90 areoidnp = mean(molaareoid(*,0)) areoidsp = mean(molaareoid(*,n_elements(linelat)-1)) newareoid = dblarr(nmolalon+2,nmolalat+2) newlat = newareoid newlon = newareoid newareoid(1:nmolalon,1:nmolalat) = molaareoid newlon(1:nmolalon,1:nmolalat) = molalon newlat(1:nmolalon,1:nmolalat) = molalat newlinelat = [2.*linelat(0)-linelat(1), linelat, $ 2.*linelat(nmolalat-1)-linelat(nmolalat-2)] newlinelon = [2.*linelon(0)-linelon(1), linelon, $ 2.*linelon(nmolalon-1)-linelon(nmolalon-2)] newlat(0,*) = newlinelat(*) newlat(nmolalon+1,*) = newlinelat(*) newlat(*,0) = newlinelat(0) newlat(*,nmolalat+1) = newlinelat(nmolalat+1) newlon(0,*) = newlinelon(0) newlon(nmolalon+1,*) = newlinelon(nmolalon+1) newlon(*,0) = newlinelon(*) newlon(*,nmolalat+1) = newlinelon(*) newareoid(0,*) = newareoid(nmolalon,*) newareoid(nmolalon+1,*) = newareoid(1,*) newareoid(*,0) = 2.*newareoid(*,1) - newareoid(*,2) newareoid(*,nmolalat+1) = 2.*newareoid(*,nmolalat) - newareoid(*,nmolalat-1) interplon = (nmolalon+1) / $ (newlinelon(nmolalon+1)-newlinelon(0)) * (lonarray - newlinelon(0)) interplat = (nmolalat+1) / $ (newlinelat(nmolalat+1) - newlinelat(0)) * (latarray - newlinelat(0)) areoidarray = interpolate( newareoid, interplon, interplat) zarray = rarray - (areoidarray + rrefareoid)/1e3 aa = where(zarray eq min(zarray)) tfromperiarray = etarray - etarray(aa(0)) ;;; End - Find areoid and tfromperi ;;; ;;; Begin - Fix P058, P059 ;;; if orbstr eq '058' or orbstr eq '059' then begin ; Now that tfromperiarray is known, ; I can interpolate speeds ; This is old approach ;vinterpol = interpol(vfudge,tfudge,tfromperiarray) ;vrelmagarray(*) = vinterpol(*) ; New approach is to use nominal stable attitude phinullrad = 0. thetanullrad = (-4.)*!pi/180. phiconversionrad = abs(vrelmagarray)*0. + phinullrad thetaconversionrad = abs(vrelmagarray)*0. + thetanullrad ux = cos(thetaconversionrad) * sin(thetaconversionrad) uy = cos(thetaconversionrad) * cos(thetaconversionrad) uz = (-1.)*sin(thetaconversionrad) ; from Takashima vrelxarray = vrelmagarray * ux * (-1.) ; u values based on atm rel to sct vrelyarray = vrelmagarray * uy * (-1.) ; u values based on atm rel to sct vrelzarray = vrelmagarray * uz * (-1.) ; u values based on atm rel to sct endif ;;; End - Fix P058, P059 ;;; alphadeg = 180./!pi*acos( double(-1.)*vrelyarray / vrelmagarray ) ; Angle of attack, angle between velocity vector and -y axis ; -1 is needed because I want to know angle between velocity vector and -y axis print, accfile ;;; Begin - Find rho ;;; sx = vrelxarray/vrelmagarray * (-1.) sy = vrelyarray/vrelmagarray * (-1.) sz = vrelzarray/vrelmagarray * (-1.) ; Attitude works with unit vectors sx,sy,sz (ux, etc, in Takashima) ; for velocity of atm wrt sct, hence factor of -1. ; Then eqns and frame definitions from Takashima papers (two of them) ; Plus SPICE ODY frame definition ASCIIimages ; Verified on P183 phideg = abs(sz)*0. aa = where(abs(sz) ne 1.) if aa(0) ne -1 then begin phideg(aa) = atan( sx(aa), sy(aa) ) * 180./!pi endif thetadeg = (-1.) * asin(sz) * 180./!pi ;;; End - Find rho ;;; ;;; Begin - Write output ;;; ; First print out what I think will be PDS table file aa = where(tfromperiarray eq 0.) print, orbstr, ' ', $ newat(aa(0)), rarray(aa(0)), zarray(aa(0)), latarray(aa(0)), $ lonarray(aa(0)), lstarray(aa(0)), szaarray(aa(0)), lsarray(aa(0)) ;;; End - Write output ;;; j=j+1 endwhile ;;; ;;; End - Loop for individual ACC file ;;; ;;; ;;; Begin - Unload SPICE files ;;; cspice_clpool ; Clear kernel pool cspice_unload, 'naif0008.tls' ; leap seconds cspice_unload, 'de414.bsp' ; SPK for planets cspice_unload, 'orb1_sclkscet_00129.tsc' ; SCLK/time conversions cspice_unload, 'pck00008.tpc' ; planetary rotational states cspice_unload, 'm01_ab.bsp' ; ODY state (position/velocity) cspice_unload, 'm01_v28.tf' ; definition of various ODY sct frames cspice_unload, 'm01_sc_ab0110.bc' ; ODY orientation and ang velocity for YYYY=2001, MM=10 cspice_unload, 'm01_sc_ab0111.bc' ; 2001, 11 cspice_unload, 'm01_sc_ab0112.bc' ; 2001, 12 cspice_unload, 'm01_sc_ab0201.bc' ; 2002, 01 cspice_unload, 'm01_sc_ab0202.bc' ; 2002, 02 ;;; End - Unload SPICE files ;;; print, startdate spawn, 'date' ; when did this slow code end? stop end ;;; ;;; End - makeodyaccprofiles.pro ;;; ;;;