Megapov+CSPICE documentation

I have added a feature to my version of Megapov which allows access to a limited subset of JPL Spice directly from the Povray scene description language. Here is the sad part: It was written for Project Blinn, but Project Blinn was (may have been) lost in the Great Hard Drive Crash of 2012. So, I am going to use it on the MSL movie I am (still) working on.

#furnsh

Yes, the spelling is different -- blame the six-character limit in fortran.

This furnishes a kernel, see furnsh_c for details.  In particular, the POV part searches the POV path for the kernel, then changes the current working directory to that of the kernel and furnishes it from there. This way any relative paths in the kernel are relative to where the kernel is. The argument to #furnsh is a POV string expression, so it may be calculated in POV script. This is exactly like the argument for an #include directive.

This fully supports any kind of kernel, including particularly metakernels, for which the path stuff above is intended.

The POV part calls kinfo_c() to check if the given kernel is already loaded, and if it is, furnsh_c() is not called. This is a limited form of persistence to handle the case of the animation loop calling #furnsh thousands of times in the same process.

The main POV support for this is added to tokenize.cpp in function Furnish();

#furnsh "generic.tm" //load generic kernels like lsk, de4xx, etc.
#furnsh "msl.tm" //load kernels directly related to MSL mission

gdpool()

The gdpool() function takes a string and a number, and calls gdpool_c() to look up the kernel constant with the given name and return the given element. An optional third argument is used to determine whether the constant was found.

The kernel constant pool is a map of strings to arrays of numbers (or strings, but that's a job for gcpool_c() which isn't supported yet). Some of these constants are native to Spice, while others come in incidentally with kernels, and some kernels are loaded specifically for their constants. Examples of constants are planet sizes and masses.

Because I haven't figured out how to create an array with a function, gdpool() takes the name of the function and the index in the constant array to get.

#furnsh "pck00010.tpc" //planetary constants kernel
#declare found=0;
#declare MarsR=gdpool("BODY499_RADII",0,found);

#if(!found)
  #declare MarsR=..some backup value...
#end

str2et()

This takes a string in any format that Spice understands and calculates the ephemeris time coordinate (Teph, roughly TDB) which is the number of seconds since the J2000 epoch. Teph ticks at the same average rate as clocks on the surface of the Earth, but doesn't include any of the relativity effects of being on the surface of the Earth. TDT is ephemeris time as measured from the surface of the earth, and TAI is an atomic clock model of this with a fixed offset. TDB is different from TDT by never more than a few milliseconds.

No kernels of any kind are needed to calculate this.

#declare ET=str2et("2012-07-24 12:34:56UTC");

sce2c()

Given a spacecraft number and an ephemeris time, get the count of sclk ticks on that time. You need a sclk kernel for the appropriate spacecraft to do this.

#furnsh "msl.tm" //metakernel pulls in the right sclk kernel

#declare ET=str2et("2012-07-24 12:34:56UTC");

#declare sclk=sce2c(-76,ET);

spkezr()

This is the heart of the matter -- the reason we use Spice. This function calculates the position and velocity of something, relative to something else, in a certain reference frame. As long as all the appropriate kernels are loaded, Spice automatically handles all the vagaries of coordinate transformation, including transforming between rotating frames, including properly handling velocity in that case.

Arguments are:

1. Target - Target name or number, but must be a string
2. ET - Teph in seconds since J2000 epoch, at observer location if doing light-time stuff
3. Ref - Reference frame to report the vector in
4. Abcor - aberration correction to apply
5. Obs - Observer name or number, but must be a string
6. Vel - Optional. If a pre-declared vector variable, vector is changed to velocity part of state vector

This returns a vector representing the position vector of the target relative to the observer (observer is at origin, axes are parallel to Ref axes). No correction for left-right handeness is done -- the x component of the result is the x component of the position as reported by Spice, likewise y and z. Same note for velocity.


#furnsh "generic.tm" //load generic kernels like lsk, de4xx, etc.

#furnsh "msl.tm" //metakernel pulls in the right sclk kernel

#declare ET=str2et("2012-07-24 12:34:56UTC");

#declare vel=<0,0,0>;

#declare pos=spkezr("MSL",ET,"J2000","LT+S","399",vel);

pxform()

Here is the other reason we use Spice -- orientation. This function calculates the quaternion needed to transform a position vector in one frame to a position vector in a different frame. This ignores the origins of the frames, even in cases where the frames have origins, so the assumed origin is the same for both to and from frames.

Arguments are:

1. From - name or number of the 'from' frame, but must be a string
2. To - name or number of the 'to' frame, but again must be a string
3. ET - Teph in seconds since J2000 epoch

pxform_c() returns a matrix, but there is no native type in POV to handle a matrix, and I still haven't figured out how to do arrays, so the POV code uses Spice m2q_c() to change the matrix into a quaternion, and returns that quaternion as a 4-component vector -- the vector imaginary part is stored in the return .x, .y, and .z components, and this works like you think it does. The scalar real part is stored in .t

#furnsh "generic.tm" //load generic kernels like lsk, de4xx, etc.

#furnsh "msl.tm" //metakernel pulls in the right sclk kernel

#declare ET=str2et("2012-08-06 05:10:00UTC");

#declare q=pxform("MSL_ROVER","J2000",ET);

ckgp()

This is almost like pxform() except that it works in sclk time instead of ephemeris time. Having said that, I'm not sure why ckgp() is needed, and why I took the time to include it.

Arguments are:

1. Instrument number
2. SclkDP - spacecraft clock tick count in DP format, the same form returned by sce2c().
3. tol - Tolerance in SCLK ticks
4. Ref - name or number of reference frame, but must be a string

ckgp_c() returns a matrix, but there is no native type in POV to handle a matrix, and I still haven't figured out how to do arrays, so the POV code uses Spice m2q_c() to change the matrix into a quaternion, and returns that quaternion as a 4-component vector -- the vector imaginary part is stored in the return .x, .y, and .z components, and this works like you think it does. The scalar real part is stored in .t


#furnsh "generic.tm" //load generic kernels like lsk, de4xx, etc.

#furnsh "msl.tm" //metakernel pulls in the right sclk kernel

#declare ET=str2et("2012-08-06 05:10:00UTC");
#declare SCLK=sce2c(-76,ET);

#declare q=pxform(-76000,SCLK,10,"J2000");

Quaternions

Ha! Quaternions come with 16 conventions built-in!

Really. Does the scalar part come first or last? Is the scalar the negative or positive cosine? Do you conjugate first or last? (Good-night everybody!) Is the direct transformation reference-to-body or vice-versa? Each permutation is possible, and I am sure that each permutation has been used.

Quaternions as used in Megapov+Spice follow these conventions:

1. All components are named. The vector part goes in components labeled .x, .y, and .z, while the scalar part goes into .t (Now there are three competing standards). When typing a quaternion constant in vector form, the first three components are the vector part .x, .y, and .z and the fourth component is the scalar part .t, so <1,2,3,4> represents a quaternion with a vector part <1,2,3> and a scalar part 4.
2. Consider an orientation expressed in axis-angle form, where the frames are related by a right-handed rotation around a particular unit-vector v by a particular angle a. In Megapov+Spice, that orientation is expressed as a quaternion <v.x*sin(a/2),v.x*sin(a/2),v.x*sin(a/2),cos(a/2)>.
3. To transform a vector, first make that vector a quaternion by appending a zero scalar to it, so to transform vector v, the corresponding quaternion is <v.x,v.y,v.z,0>. Now the quaternion transformation formula is v_to=q*v_from*q', where * represents quaternion multiplication and q' represents the quaternion conjugate of q, where the vector part is replaced by its negative and the scalar part stays the same
4. Quaternions created by pxform are used to transform a vector in the "from" frame into a vector in the "to" frame, and it is up to the user to decide which is which, and put them in the right order in a call to pxform. The formula in part 3 does not care which is which.

Quaternions are not naturally supported in Megapov, so we write some macros to add the support: We write macros when we can, and C++ code when we have to.

First, quaternion multiplication. The stone bridge equations i²=j²=k²=ijk=-1, along with the anti-comutativity of multiplying basis vectors, is one true constant across all quaternion implementations, the one thing that makes a quaternion a quaternion. We use the formula mentioned on Wikipedia to do the operation with vector cross and dot products, rather than the error-prone 16 multiplications and 12 adds where Jeppesen's law of Signs is sure to bite you.

#macro QuatMult(Q1,Q2)
  #local v1=<Q1.x,Q1.y,Q1.z>;
  #local r1=Q1.t;
  #local v2=<Q2.x,Q2.y,Q2.z>;
  #local r2=Q2.t;
  #local r=r1*r2-vdot(v1,v2);
  #local vv=r1*v2+r2*v1+vcross(v1,v2);
  (<vv.x,vv.y,vv.z,r>)
#end

Now we build quaternion rotation on top of it:

#macro QuatVectRotate(Q,V)
  #local Q1=QuatMult(Q,<V.x,V.y,V.z,0>);
  #local Q2=QuatMult(Q1,<-Q.x,-Q.y,-Q.z,Q.t>);
  (<Q2.x,Q2.y,Q2.z>)
#end

Matrices

Finally, we want to build a matrix<> as understood by POV. We will include the quaternion rotation, then a translation, so that we only need one matrix, because I have had bad experiences with using more than one matrix, or matrix and other kinds of transformation, for the same object. First, a word about homogeneous coordinates:

Normally, in 3D space, you use a 3x3 matrix to represent any arbitrary affine transformation which maps the origin to the origin. You can scale, shear, rotate, or do any combination. But, you can't translate. The origin maps to the origin.

But, there is a clever way around this. Use a 4D vector instead, and make the 4th coordinate 1. This opens up a 4x4 transformation matrix. By doing the right thing with the extra cells, we can translate as well as the other affine transformations. Further, the matrices even compose properly, so we can combine translates with arbitrary other transformations into a single matrix. Here's the math:

[V'x] [c11 c12 c13 Tx][Vx]

[V'y]=[c21 c22 c23 Ty][Vy]

[V'z] [c31 c32 c33 Tz][Vz]

[ 1 ] [ 0   0   0  1 ][1 ]

POV supports this, but it uses the transpose of that matrix, so you have to specify this:

matrix <c11,c21,c31,

        c12,c22,c32,
        c13,c23,c33,
        Tx ,Ty ,Tz>

We use the Spice formula, transposed, with the extra translation term, as follows:

#macro QuatTrans(Q,T)
  #local q0=Q.t;
  #local q1=Q.x;
  #local q2=Q.y;
  #local q3=Q.z;
  matrix <1-2*(q2*q2+q3*q3),2*(q1*q2+q0*q3)  ,2*(q1*q3-q0*q2),
          2*(q1*q2-q0*q3)  ,1-2*(q1*q1+q3*q3),2*(q2*q3+q0*q1),
          2*(q1*q3+q0*q2)  ,2*(q2*q3-q0*q1)  ,1-2*(q1*q1+q2*q2),
          T.x              ,T.y              ,T.z>
#end
 
So with this, my animation of the MSL entry sequence is done except for the art. I have a parachute model and entry capsule model with separable heatshield, but I also need the six entry balance masses, the descent stage, and the rover itself. The last, I have converted from the old scad model I had laying around.