[Spice_discussion] Changing resolution of CK files
William Thompson
William.T.Thompson.1 at gsfc.nasa.gov
Wed Jul 26 07:56:01 PDT 2006
Nat, Boris:
Thank you for all your detailed comments. They're extremely helpful, and I will
certainly try to incorporate them.
I was aware of the dangers of averaging angles together in the presence of
singularities, but I didn't appreciate that the same danger existed for
quaternions. That didn't show up in any of the test cases that I ran. I will
need to modify my smoothing approach, probably using Euler angles.
I was also thinking about strategies to deal with the asymmetric smoothing boxes
at the beginning and ending of the file, and I think I know how to deal with that.
Nat, to answer your last question, my language of choice would have been IDL. I
grew up on Fortran, but it's been almost exclusively IDL for the last 20 years.
Writing this short program really brought home how much I've forgotten.
Cheers, and thanks again,
Bill Thompson
Nat Bachman wrote:
>
> 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
>>
>>
>
>
--
William Thompson
NASA Goddard Space Flight Center
Code 612.1
Greenbelt, MD 20771
USA
301-286-2040
William.T.Thompson.1 at gsfc.nasa.gov
More information about the Spice_discussion
mailing list