[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 
to Q(i+1) is 180 degrees apart from the rotation corresponding to Q(i), 
is much too far; the whole idea of smoothing doesn't make sense unless 
quaternions are very close to each other.  The quaternions must have 
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: 


   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 
   non-standard (in Fortran 77, that is), though I don't know of any 
   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 
   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

      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

      VADDG respectively find the sum of two 3-vectors and two general-
            dimension vectors

      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.

      VLCOM3 respectively find linear combinations of two or three 

   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 
   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
            qnorm = sqrt(qnorm)
            do i=1,4
               quats(i,j) = quats(i,j) / qnorm

   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
>>>>   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  
>>>> 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  This routine is called with 1-3 parameters:
>C	INFILE	= The name of the CK file to process.
>C	OUTFILE = The name of the output file.  If not passed, then
>C		  defaults to "cksmooth_out.bc".
>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  Example:
>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  Note that the output file must not exist when calling this program.
>      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