#include <stdio.h>
#include <string.h>
#include "SpiceUsr.h"

#define  LSK     "naif0006.tls"
#define  DESPK   "de405s.bsp"
#define  SCSPK   "sirtf.bsp"

SpiceDouble  ra;
SpiceDouble  dec;
SpiceDouble  axdef[3];
SpiceInt     indexa;
SpiceDouble  plndef[3];
SpiceInt     indexp;
SpiceDouble  mat[3][3];
SpiceDouble  q[4];

SpiceDouble  lt;
SpiceDouble  epoch;
SpiceDouble  state[6];

SpiceInt     handle;
SpiceInt     obs;
SpiceInt     targ;

SpiceChar    utc_epoch[32];
SpiceChar    frame[32];
SpiceChar    corr[5];

void main()
{

   /*
   Load SPICE kernel files required for computation:
      LSK         Leapseconds Kernel file;
      DESPK       Planetary Ephemeris SPK file;
      SCSPK       Spacecraft SPK file.     
   */
   ldpool_c ( LSK );
   spklef_c ( DESPK, &handle );
   spklef_c ( SCSPK, &handle );

   /*
   Convert ra & dec to rectangular coordines. For simplicity assume 
   that ra & dec are given in J2000 inertial frame so result vector 
   doesn't need to be rotated to this frame. (Angles are converted to
   radians by multiplying by CSPICE function dpr_c returning number
   of radians per degree.)
   */
   ra = 120.0 * rpd_c();
   dec = 45.0 * rpd_c();
   radrec_c ( 1.0, ra, dec, axdef );
   
   /*
   Compute position of the Sun as seen fom the s/c in J2000 
   inertial frame. spkez_c input values are:
   
      targ        NAIF ID code of a target body (10 for Sun);
      
      obs         NAIF ID code of observer (-999 for SIRTF; Note 
                  that this is NOT offical ID for the project but
                  we have to use it because it is used in the SIRTF
                  SPK file);
                  
      utc_epoch   UTC epoch for which state will be computed (this 
                  time is converted to Ephemeris Time, ET, using 
                  str2et_c routine);
                  
      frame       reference frame relative to which state will be 
                  computed ("J2000" for J2000 inertial frame);

      corr        correction to be applied to the computed state 
                  ("lt" for the light time correction);

   */
   targ = 10;
   obs  = -999;
   strcpy ( frame, "J2000" );
   strcpy ( corr, "lt" );

   strcpy ( utc_epoch, "2003 MAR 11 16:07:25" );
   str2et_c ( utc_epoch, &epoch );

   spkez_c  ( targ, epoch, frame, corr, obs, state, &lt );
   
   /*
   Construct transformation matrix by calling twovec. Assume that 
   ra & dec direction defines Z axis and direction to the Sun
   defines ZX plane.
   */
   indexa = 3; 
   indexp = 1;
   twovec_c ( axdef, indexa, state, indexp, mat );
             
   /*
   Convert matrix to quaternion. Output quaternion is in SPICE form:
   
      "q is a 4-dimensional vector. If r rotates vectors by an angle
      of R radians about a unit vector a, where R is in [0, pi], then
      if h = R/2, 
      
            q = ( cos(h), sin(h)a ,  sin(h)a ,  sin(h)a )."
                                 1          2          3 
   */
   m2q_c( mat, q );
   printf ( "quaternion %f %f %f %f\n", q[0], q[1], q[2], q[3] );
   
}