[Spice_discussion] Changing resolution of CK files
Nat Bachman
Nathaniel.Bachman at jpl.nasa.gov
Tue Jul 25 19:46:26 PDT 2006
Hi Bill,
This is from Nat Bachman; I'm a member of the NAIF team.
I too have a few comments on the CK smoothing program you sent us. :)
Math issues
===========
Regarding quaternion smoothing: you can do it, but you must be
sure that you're combining consistent quaternions.
This is really analogous to averaging longitudes: you'd have
a problem if you averaged consecutive values of 359, 0 degrees and 1
degree, but
if you compute your average using 359 degrees, the closer choice to 359
of 0 and
360 degrees, and the closer choice to 360 degrees of 1 and 361 degrees,
you'll get a valid average ("smoothed") value. If you don't know a
priori that
the correct average value is closer to 360 degrees than 180 degrees, then
averaging is not valid.
Assuming you have a situation where smoothing is known to be valid,
for each quaternion Q(i+1) following quaternion Q(i) in the "box," you
should replace Q(i+1) with -Q(i+1) if -Q(i+1) is closer than Q(i+1)
(measured by the
4-dimensional Euclidean norm of the quaternion difference) to Q(i).
If the two choices are equidistant from Q(i), then the rotation
corresponding
to Q(i+1) is 180 degrees apart from the rotation corresponding to Q(i),
which
is much too far; the whole idea of smoothing doesn't make sense unless
consecutive
quaternions are very close to each other. The quaternions must have
identical
length (normally unit length) for this algorithm to work.
Averaging (close) Euler angles is valid, but this approach still presents
the same branch singularity issue.
Averaging rotation matrix elements is also valid (given the same
precondition that consecutive matrices are "close"), but you must make the
interpolated matrix back into a rotation matrix (having orthogonal
columns and
determinant == 1) before using the SPICE routine M2Q to convert the
interpolated matrix to a quaternion.
Concerning box bounds at segment boundaries: I'm a little worried
about this code: I may be mistaken, but it looks like averages are
being taken over boxes that truncated asymmetrically, and hence are not
centered at the time of interest. Unlike the case of polynomial fitting,
where this would usually be ok, it seems to me that for averaging, you
really will
end up with a displaced value if each box is not symmetric about the time
for which the average is to be computed. If the same were done with
a finite sequence of consecutive integers, using a nominal total box
width of 3,
you'd normally obtain an average of N at integer N, but at the left
boundary L the average would be L+0.5 (the box width here would be 2) and
at the right boundary R the average would be R-0.5.
Fortran issues
==============
The following comments are based on my assumption that you're using
Fortran 77. If you're using a more modern version of Fortran, some/all
of these comments might not be applicable. Still, FWIW:
1) I strongly recommend placing an IMPLICIT NONE statement right
under your PROGRAM statement. The Fortran bugs that arise from
misspelled variable names can be quite nasty and hard to find;
IMPLICIT NONE will catch these at compile time.
2) I recommend turning on all (reasonable) warning options when
you compile your code. Since Fortran 77 is so weak at detecting
problems at run time, one needs all the help one can get at
compile time.
3) There is in fact an undeclared variable "qnorm" whose type will
default to REAL. This will cause loss of precision (not necessarily
critical for your application) in your normalized quaternions.
4) Unlike C, floating point constants in Fortran are REAL unless
written with DOUBLE PRECISION syntax. For example, the double
precision constant 3.141592653589793D0 will get rounded to
single precision if you write it as follows:
3.141592653589793
Even if you store this constant in a double precision variable,
the stored value may be truncated to single precision; this behavior
is platform-dependent.
Also unlike C, you can obtain disastrous results if you pass a literal
single precision floating point constant into a routine expecting
a double precision input. No automatic type promotions occur.
I'm just noting this; I don't see such a problem in your code.
The initializations you have using integer 0 won't be a problem.
5) Double quote characters (") delimiting literal strings are an
extension to Fortran 77 that may be non-portable. Single quotes
are standard.
6) Use of lower case letters outside of literal strings is, strictly
speaking,
non-standard (in Fortran 77, that is), though I don't know of any
platforms
that still fail to support it.
SPICE issues
============
All of the following are for future reference. I see no problems with
your use of SPICE that would prevent the program from working properly.
1) The values of ND and NI for C-kernels are part of the published
interface:
you don't have to look them up from the file record. For CKs, ND is
2 and
NI is 6. From the CK Required Reading, here's the layout of a CK
segment descriptor:
-----------------------------------
DCD(1) | Initial SCLK |
-----------------------------------
DCD(2) | Final SCLK |
-----------------------------------
ICD(1) | Instrument |
------------------
ICD(2) | Reference |
------------------
ICD(3) | Data type |
------------------
ICD(4) | Rates Flag |
------------------
ICD(5) | Begin Address |
------------------
ICD(6) | End Address |
------------------
You can safely declare 'dc' with bound 2 and 'ic' with bound 6.
You can safely declare 'sum' with bound 5 ( this is ND + (NI+1)/2 )
and 'name' with length 40 characters (this is the number of bytes
occupied by 'sum').
IMO it's still a good idea, as you've done, to use symbolic names for
all of these dimensions. However, the values can be safely declared using
PARAMETER statements.
Our not providing these declarations in a public INCLUDE file is
a deficiency of the Fortran SPICE Toolkit.
2) To open the existing CK, all you need is the DAFOPR call.
3) To open the new CK, all you need is a call to CKOPN.
4) The code current simply copies the "segment identifier" 'name' from
the input segment to the output segment. Although this is far from
critical, your program might create new names indicating that these
data have been changed.
5) The call to DAFGH is unnecessary.
6) The call to DAFCLS for the input file is unnecessary. (For the
output file, the corresponding DAFCLS call, which you have,
is essential.)
7) The variables 'na' and 'values' are declared but never used.
8) There are many SPICE routines that could be used to simplify
the code (possibly at a very slight cost in terms of run-time
speed).
VDIST
VDISTG respectively find the Euclidean distance between two 3-vectors
and between two general-dimension vectors
VNORMG finds the norm of a general-dimension vector
VHATG normalizes ("unitizes") a general-dimension vector
MOVED transfers a sequence of elements from one double precision array
to another
VADD
VADDG respectively find the sum of two 3-vectors and two general-
dimension vectors
VSCL
VSCLG
VSCLIP respectively scale a 3-vector producing a distinct output,
scale a general-dimension vector producing a distinct output,
and scale a 3-vector in place.
VLCOM
VLCOM3 respectively find linear combinations of two or three
3-vectors.
There are many more routines that may allow you eliminate loop code.
SPICE routines having names starting with V are typically vector
manipulation routines. The SPICE "Most Useful/Most Used" routines
document
is a good place to find brief descriptions of these as well.
One such application would be transformation of this code
C Renormalize the quaternions.
qnorm = 0.
do i=1,4
qnorm = qnorm + quats(i,j)**2
enddo
qnorm = sqrt(qnorm)
do i=1,4
quats(i,j) = quats(i,j) / qnorm
enddo
to this:
C Renormalize the quaternions.
call vhatg ( quats(1,j), 4, qtemp )
call moved ( qtemp, 4, quats(1,j) )
Another deficiency of SPICE is that we don't have general-dimension,
in-place routines that allow
call vhat_gen_dimension_in_place ( 4, quats(1,j) )
You may see in various bits of SPICE documentation calls like
call vhatg( quats(1,j), 4, quats(1,j) )
but these are not allowed by the Fortran 77 standard, and in fact
they cause problems on some platforms, which we've observed in some
cases where optimization is used. As of the N0060 Fortran SPICE Toolkit,
we've started adding in-place vector manipulation routines such as VHATIP
to work around this Fortran limitation.
9) If you plan to port your code, you may want to use the SPICE
support library routine GETCML to fetch the command line arguments.
NAIF normally tells SPICE users to *not* use support library code,
but this routine is an exception: the interface and functionality
are stable, and it's supported on all platforms on which SPICE
runs. It gives you a more portable interface than you'd get
by making direct calls to platform-specific system routines.
The SPICE routine NEXTWD is convenient for picking off space-delimited
arguments from the string returned to your program by GETCML.
If you're at all curious about this, I can show you an example
of how the code would look.
Also, one question for you: given that you had to write this attitude
smoothing program, in which language would you prefer to have written it?
Best regards,
-Nat Bachman (JPL/NAIF)
Nathaniel.Bachman at jpl.nasa.gov
William Thompson wrote:
> Ed, Boris:
>
> Thanks for all your help. I've now got a working program. It's been
> so long since I did any serious programming in FORTRAN that I
> completely forgot about that 72 character limit. I also had some
> typos which were keeping me from calculating the averaged quaternions
> properly, which I've now corrected. Here's the finished program--I
> thought you'd be interested.
>
> I will use this to first smooth out the jitter from the CK files into
> an intermediate file, and then use CKSMRG to downsample to a more
> manageable size.
>
> Thanks again,
>
> Bill Thompson
>
>
> Ed Wright wrote:
>
>> To: William Thompson
>> From: Ed Wright
>>
>> On Jul 24, 2006, at 11:09 AM, William Thompson wrote:
>>
>>> Boris, Ed, Chuck:
>>>
>>> I've made some progress in reading and processing the CK files.
>>> I've been able to read the files, and print the contents out to a
>>> text file. However, I'm having trouble writing out the data. When
>>> I try to run my program (attached), I get the following error message.
>>>
>>>> argo:~/daf> cksmooth
>>>> Name of input file?
>>>> ahead_2006_303_01.ah.bc
>>>> 3.94467348E+11 3.94481664E+11 -234000 55924 1
>>>> STEREO Ahead spacecraft bus
>>>
>>>
>>>
>>>>
>>>> ======================================================================
>>>> ========== Toolkit version: N0060
>>>> SPICE(INVALIDNUMREC) --
>>>> 0 is an invalid number of pointing instances for type 3.
>>>> A traceback follows. The name of the highest level module is first.
>>>> CKW03
>>>> Oh, by the way: The SPICELIB error handling actions are
>>>> USER-TAILORABLE. You
>>>> can choose whether the Toolkit aborts or continues when errors
>>>> occur, which
>>>> error messages to output, and where to send the output. Please
>>>> read the ERROR
>>>> "Required Reading" file, or see the routines ERRACT, ERRDEV, and
>>>> ERRPRT.
>>>>
>>>> ======================================================================
>>>> ==========
>>>
>>>
>>>
>>> Can you suggest what the error is?
>>
>>
>>
>> As Boris mentioned, the CKW03 argument 'nints' should be a boolean.
>>
>> Another problem seems that your smoothing algorithm corrupts the
>> quaternion arrays, ||q|| = 1. A better technique for smoothing might
>> be to process the angular values represented by the quaternions
>> then convert the results back to quaternions.
>>
>> As always,
>> Ed Wright
>> ed.wright at jpl.nasa.gov
>> 1-818-354-0371
>>
>
>
>------------------------------------------------------------------------
>
>C This program applies a smoothing algorithm to a STEREO Type III
>C Attitude History (CK) file. It is used in conjunction with the
>C CKSMRG utility to produce a reduced-resolution attitude history file
>C representative of the average pointing.
>C
>C This routine is called with 1-3 parameters:
>C
>C INFILE = The name of the CK file to process.
>C
>C OUTFILE = The name of the output file. If not passed, then
>C defaults to "cksmooth_out.bc".
>C
>C NSMOOTH = A parameter related to the size of the boxcar average
>C smoothing function. The box will run between
>C I-NSMOOTH to I+NSMOOTH, so that the total size of the
>C box is 2*NSMOOTH+1. The default is NSMOOTH=150,
>C which for STEREO is equivalent to a box 5 minutes
>C wide.
>C
>C Example:
>C
>C rm -f ahead_2006_303_01_temp.ah.bc
>C cksmooth ahead_2006_303_01.ah.bc ahead_2006_303_01_temp.ah.bc
>C
>C Note that the output file must not exist when calling this program.
>C
>
> program cksmooth
>
> character*(128) infile, outfile
> character*(80) arch, type, ref, ssmooth
> character*(1000) name
> integer inhandle, outhandle, nd, ni, ic(250), nrec, na, nints
> integer i,j,k,navg
> logical found
> double precision dc(125), sum(125), values(10000)
> parameter (maxrec=1000000)
> double precision record(8),inarray(8,maxrec)
> double precision sclkdp(maxrec), quats(4,maxrec), avvs(3,maxrec)
>
>C Get the number of arguments. The first argument should be the name
>C of the input file. If supplied, the second argument will be the
>C name of the output file, and the third argument will be the window
>C size (+/-).
>
> nargs = iargc()
> call getarg (1, infile )
> outfile = "cksmooth_out.bc"
> if (nargs .ge. 2) call getarg (2, outfile)
> nsmooth = 150
> if (nargs .ge. 3) then
> call getarg (3, ssmooth)
> read (ssmooth,'(I5)') nsmooth
> endif
>
>C Get the architecture and type. Open the file, and get the values
>C of ND and NI.
>
> call getfat (infile, arch, type)
> call dafopr (infile, inhandle)
> call dafhsf (inhandle, nd, ni)
>
>C Open the new file to initialize it.
>
> call dafonw(outfile, type, nd, ni, "Smoothed file", 0, outhandle)
>
>C Begin forward search, and find the first array.
>
> call dafbfs (inhandle)
> call daffna (found)
>
>C Step through the arrays. Get the summary and unpack it. Get the
>C name.
>
> do while (found)
> call dafgs (sum)
> call dafus (sum, nd, ni, dc, ic)
> call dafgn (name)
>
>C Get the number of pointing instances. Step through the pointing
>C instances, and collect the array.
>
> call cknr03(inhandle, sum, nrec)
> do j=1,nrec
> call ckgr03(inhandle, sum, j, record)
> do i=1,8
> inarray(i,j) = record(i)
> enddo
> enddo
>
>C Separate out the data into spacecraft clock values, quaternions,
>C and angular velocities.
>
> do j=1,nrec
> sclkdp(j) = inarray(1,j)
>
> j1 = j-nsmooth
> if (j1 .lt. 1) j1 = 1
> j2 = j+nsmooth
> if (j2 .gt. nrec) j2 = nrec
> navg = j2-j1+1
>
>C Sum together the quaternions.
>
> do i=1,4
> quats(i,j) = 0
> do k=j1,j2
> quats(i,j) = quats(i,j) + inarray(i+1,k)
> enddo
> enddo
>
>C Renormalize the quaternions.
>
> qnorm = 0.
> do i=1,4
> qnorm = qnorm + quats(i,j)**2
> enddo
> qnorm = sqrt(qnorm)
> do i=1,4
> quats(i,j) = quats(i,j) / qnorm
> enddo
>
>C Average together the angular velocities.
>
> do i=1,3
> avvs(i,j) = 0
> do k=j1,j2
> avvs(i,j) = avvs(i,j) + inarray(i+5,k)
> enddo
> avvs(i,j) = avvs(i,j) / navg
> enddo
> enddo
>
>C Write out the smoothed pointing segment.
>
> nints = 1
> ref = "J2000"
> call ckw03(outhandle, dc(1), dc(2), ic(1), ref, nints, name,
> $ nrec, sclkdp, quats, avvs, 1, sclkdp(1))
>
>C Step to the next array
>
> call dafgh(inhandle)
> call daffna (found)
> end do
>
>C Close the input and output files.
>
> call dafcls(inhandle)
> call dafcls(outhandle)
>
> end
>
>
More information about the Spice_discussion
mailing list