which shows the annual modulations quite nicely.
These are my notes on how to do some basic analysis of the JPL ATDF format, which is used for doppler and range tracking of deep space spacecraft.
The files contain records with information about doppler and range values, and also about the configuration of the uplink and downlink antennas during the tracking passes. The goal is usually to extract this information and reduce it for the purposes of orbit determination.
I should say first that I am not an expert at this analysis, but rather an interested person persuing a hobby. I will give some code in IDL that can get you a long way.
Outside of JPL itself there is little information on how this data is accessed or reduced. I have found a few URLs which provide extremely valuable information about the formatting and interpretation of the data in the files.
The format itself is documented in TRK-2-25, a JPL DSN document.
The NSSDC has a hardcopy of this document, but it is also available
online at:
TRK-2-25 Format Specification See, for example, http://pds-geosciences.wustl.edu/geodata/mpf-m-rss-1-5-radiotrack-v1/mprs_0001/document/ http://pdssbn.astro.umd.edu/volume/hal_1001/document/trk_2_25.asc
A fantastic link is the following article by Morabito & Asmar, which
appeared in JPL's own refereed journal.
Morabito & Asmar "Radio-Science Performance Analysis Software" http://tda.jpl.nasa.gov/progress_report/42-120/120B.pdf
It describes the basic processes involved in reducing Doppler data, including how to convert from doppler "counts" to a sky frequency. There is a good review of how precession, nutation and relativistic effects are handled. This can be important for high precision work because the observing platform on the earth is moving with respect to an inertial frame.
As documented in TRK-2-25, the ATDF physical record is 2016 32-bit words in length (=8064 bytes), and each physical record consists of 28 72-bit-word (=288 byte) logical records.
The NSSDC staff loaded staged two tapes on a VMS machine and I downloaded to Unix. Unix does not have a concept of record blocking, and there appears to be a small artifact introduced into the files. The artifact is an additional byte is inserted at the *physical* record boundaries.
In other words, every 8065th byte is a record marker and is not part of the ATDF stream. The easiest way to deal with this is to strip out the marker.
Here is IDL code to read the file and remove the record marker
bytes:
openr, 50, 'pioneer10.fl1' nbytes = (fstat(50)).size ;; Number of bytes in file nphrec = nbytes / 8065L ;; Number of physical records (+1 marker byte) bb = bytarr(8065, nphrec) ;; Read the records from file readu, 50, bb close, 50 bb = bb(0:8063,*) ;; Strip out the marker byte
Next, the physical records can be reformatted into logical records, knowing that each logical record is 288 bytes long.
nrec = n_elements(bb)/288 bb = reform(bb, 288, nrec, /overwrite)
The OVERWRITE keyword saves memory. In the case of the file being larger than core memory, obviously one would want to split the file into chunks.
With the final REFORM operation above, the byte array BB is now divided into logical records of 288 bytes each.
The first two logical records contain meta-information about the dataset (see TRK-2-25). The first record stores information about when the tape file was created, and an indicator that it is ATDF data.
The second record stores the start and stop time of the data, as well as the spacecraft identification number and transponder frequency.
The TRK-2-25 document divides the 288-byte logical record up into bit-packed fields. In some cases the fields appear to be aligned on 36-bit boundaries, which would appear to be indicative of the Univac heritage which generated the data.
Thus, any software which reads the data will need to be adept at shifting and masking bits. I have included a small IDL routine called BITSTRACT which does this (see DOWNLOADS).
There appear to be no floating point numbers stored in an ATDF file. Instead, numbers are stored as decimal integers. Where fractional values are needed, the document describes the multiplier used. [ For example column 30, Sampler Time, is the sampler time in seconds multiplied by 100. Thus, this column can keep up to 0.01 seconds of precision. ]
Where large dynamic range is needed in a value, the field is divided into a "high" and "low" part (ie, most significant and least significant digits), and stored as two separate fields. [ See, for example, column 31-32, Doppler Count. ] The documentation establishes predetermined multipliers for both parts.
While some fields are longer than 32 bits, there are no significant digits beyond the 31st. I conclude this based on the use of "sign" bits which extend the value to 36 bits. The number of sign bits is never less than 5. Thus, all of the fields are appropriate to be stored in 32-bit integers (or smaller in many cases).
There are some values listed as "two's complement." This is the normal notation for negative numbers in today's computers anyway. The only care that needs to be taken is when a value has fewer than 8 bits, in which case the sign bits must be extended to fill the full byte. BITSTRACT does this.
The first thing to do is to remove the first two records to a different array, since they are header and not tracking records.
bbhead1 = bb(*,0) ;; File header bbhead2 = bb(*,1) ;; Transponder header bb = bb(*,2:*) ;; Remove first two records from main array
Next, the data can be read using the IDL procedure named BITSTRACT. This routine accepts the raw byte array, a description of each field's data, and the desired data.
A call looks something like this:
bitstract, bb, tkform, 4, day
where
Upon return, the variable DAY will contain the values for field number 4 in each record, which is hopefully the day of the year.
The TKFORM variable describes each field, according to which bit
positions it is stored in, the data type, and whether the field is
signed or not. These fields are already defined in the
"trk225fmt.pro" file, which I am supplying. Each field is described
by a structure. For example, field number 4 is:
{TRK_BIT_FIELD, 4, 85, 100, 0, 16, 'UINT', 'DAY_NUM'}, $
Where each of the fields in this structure is defined by:
hformat = { TRK_BIT_FIELD, $
ITEM: 0L, $ ;; Item number
START: 0L, $ ;; Bit position, starting at 1
STOP: 0L, $ ;; Input length in bits
SIGNLENGTH: 0L,$ ;; Length of sign field in bits
OUTLENGTH: 0L, $ ;; Convert to output length in bits (MSB lost)
TYPE: '', $ ;; Data type to convert to
NAME: '' $ ;; Name of field
}
In this case, field 4, day of year, extends from bit 85 to bit 100, is unsigned, and can be stored in a 16-bit unsigned integer. The NAME field ('DAY_NUM') string gives a descriptive name which can be used mnemonically, or in automated processing of the files.
Using BITSTRACT one can extract all of the columns of the ATDF file into IDL variables. In fact, that is what the BITVERT and BITREAD procedures will do (see below).
I have a set of routines to convert a binary ATDF file to ASCII text. The resulting text representation is about a factor of 3 bigger than the original binary file. However, compressing using gzip or zip can reduce the size by a factor of 25!
The main routine is BITVERT, which coordinates the reading, converting, and writing of the data. Two forms of ASCII output are created, one of which is a direct verbatim copy of the integer values, and a second form is converted to more useful units. See the download section.
Please note that the "processed" ASCII files do not contain a complete set of data. Many of the columns have been dropped because they seemed irrelevant for the Pioneer 10 data at hand. You will need to customize BITVCONV if you want additional fields.
At this point you probably want to get some real data to work on. Here are some starting points.
The date is stored in fields 3-7. ASSUMPTION: the time is in UTC.
Also, the year field only stores the last two digits, so it is Y2K
ambiguous. Here is code to read the date information and turn it into
a Julian date:
bitstract, bb, tkform, 3, year bitstract, bb, tkform, 4, day bitstract, bb, tkform, 5, hour bitstract, bb, tkform, 6, minute bitstract, bb, tkform, 7, second jdy = julday(1d,1d,1900+year) - 0.5d ;; Convert to midnight jd = jdy + (day-1) + (hour + (minute + second/60d)/60d)/24d
This formulation ignores leap seconds, which can be important for high precision work.
The doppler "count" is a measure of the integrated doppler
frequency. It measures the number of beat cycles of the downlinked
signal, compared to a reference oscillator. The number of counts can
be extracted with this:
bitstract, bb, tkform, 31, doppler_cnt_hp bitstract, bb, tkform, 32, doppler_cnt_lp dcnt = doppler_cnt_lp/1000d + doppler_cnt_hp*10000d
It is converted to a doppler frequency by computing a 1st difference and dividing by the integration time, which is the 1st difference of the clock times.
drate = (dcnt(1:*)-dcnt)/(jd(1:*)-jd)/86400d
This doppler rate is referred to the reference frequency in a
non-intuitive way, which depends on the uplink and downlink signal
frequency bands. To compute the "sky frequency", which is the actual
downlink signal frequency as it enters the antenna, one needs both the
reference frequency and the "bias" value:
bitstract, bb, tkform, 20, doppler_bias bitstract, bb, tkform, 40, doppler_ref doppler_ref = doppler_ref/10d ;; Divisor per trk-2-25
Please see Morabito & Asmar for the formula for converting S-band
data, and other bands as well. Here is how it is computed for
S-band:
sdoppler = (doppler_bias GE 0)*2L - 1 fsky = 96d*240d/221d*doppler_ref - sdoppler*(drate - 1d6)
The first expression computes the doppler bias sign value, which is negative if the doppler bias is negative. The second formula is Morabito & Asmar's equation for the sky frequency at the downlink receiver.
For other bands, I believe that the following monograph would be useful, but I have just recently found it (and the web page access appears to be sporadic):
T. D. Moyer, *Formulation for Observed and Computed Values of Deep Space Network (DSN) Data Types*, JPL Publication 00-07, Jet Propulsion Laboratory, Pasadena, California, 2000. http://descanso.jpl.nasa.gov/Monograph/mono.cfm
The sky frequency is the starting point for the reduction of doppler data.
There appear to be a lot of noisy points, which contain garbage or clearly wrong values. For both of the files, I saw about 10-12% bad data in each. However, this is not unexpected, especially for a spacecraft (Pioneer 10) which has essentially left the solar system.
To remove the noise, I made a set of successive filters. The main
filter is this:
mask = (abs(doppler_bias) LE 1 AND jd GE 0 AND $
doppler_ref GE 2197d4 AND doppler_ref LE 2200d4 AND $
dcnt(1:*) GT dcnt AND jd(1:*)-jd LT .1 AND $
abs(fsky-fmed0) LT .6d6)
which means:
where fmed0 = median(fsky).
This set of filtering criteria ensures that we are only looking at contiguous data with no gaps. The value of 0.6d6 is chosen so that the frequency modulation due to the annual motion of the earth is permitted, but large-scale noisy points are removed.
The good points can be collected using the IDL WHERE() function,
wh = where(mask)
And a second pass of filtering can be done. Here I did a sliding MEDIAN filter, which removes a lot of the garbage.
fmed1 = fsky fmed1(wh) = median(fsky(wh),500)
I chose the value of 500 somewhat arbitrarily, but hopefully so that
it included a few days worth of data. This filter is short enough
timescale, however, to remove the earth's annual modulation, which
leaves only the daily modulations, which are only +/- 0.5 km/s (~ +/-
10 kHz). Thus, we can use this to filter further:
mask = mask AND (abs(fsky-fmed1) LT 10000d) wh = where(mask)
At this point I did an even tighter filter by fitting a sinusoid to the data, but this requires the MPFIT & MPFITEXPR functions from my web page.
p0 = [-40715d,227616d,1.648D-02,2.13D+00] expr = 'p(0) + p(1)*cos(p(2)*x+p(3))' pp = mpfitexpr(expr, jd(wh), fsky(wh)-fmed0, 1000d, p0)
These coefficients will obviously be different for every observation and spacecraft, so this is just a toy example.
The resulting offset-subtracted doppler frequencies look like this:
which shows the annual modulations quite nicely.
Further analysis requires more information, much of it beyond the scope of what I can describe in this note. The Morabito & Asmar paper has a very good summary of the types of tasks involved.
Users will need to obtain the reference frequencies of the uplink and downlink transmitters. These can be gotten from the records with data type "6", which is the "ramp" data type (see field 12).
The uplink transmitter can be configured to send a ramped doppler signal. In principle this is stored in field 112, "Programmed Frequency Ramp Rate." However, I have found that only some of the Pioneer 10 data has a ramp value. Furthermore, the documentation describes field 112 as a ramp *rate* and yet the units are Hz*1000000, which do not appar to be the correct units. I believe the correct units are Hz/s.
For the uplink frequency, it appears that fields 113-114, programmed start frequency, should be used. This is the frequency of a programmable digital oscillator (DCO), which is then fed to the transmitter/exciter and multiplied by 96.
The documentation I have says that the DCO frequency should be near 24 MHz, however, the values printed in fields 113-114 appear to be about twice that value, or 48 MHz. To make the output be S-band, ~2.2 GHz, the multiplier should be 48 and not 96.
I have some reduced data from Anderson et al's reduction of the Pioneer data, and they use the DCO frequency, and not the value in field 116, which is the transmitter/exciter frequency. Part of the reason to avoid field 116 is that it is low precision.
It appears you can distinguish between uplink and downlink transmitters by the values in fields 11 and 64 (and perhaps 25 & 26). Both are stored in the same file.
None of the data in the Pioneer 10 files seems to have angle or "DRVID" data.
Virtually all of the Pioneer 10 data are low rate doppler, so almost none of the auxiliary fields, fields 42-59, are filled.
Field 30 is sampler time, but I do not understand how this correlates with the integration time.
There are almost no range data points in the Pioneer 10 data, so I cannot comment on those fields.
Fields 76 & 77, slipped cycles and doppler noise, in principle describe the quality of the signal, but I am not sure how to interpret these values.
Raw data
(is unfortunately no longer available on-line, after my old computer system was shut down)