/*This file is TransCal.cpp, version 050324. 24 March 2005. Changed main and function get_qmatrix. Attitude transformations from quaternions are now interpolated. No interpolation over gaps greater than QLIMIT seconds. QLIMIT now defined as 16 sec (was 60 sec). Improved reporting of CHATT file errors. 5 Jan 2005, added calls to SPICE error routines so that SPICE error messages are written to report file as well as standard output. 10 Dec 2004, changed open_input_ff and open_output_ff. When run interactively, will try again if cannot open files. 28 June 2004, fixed bug in get_qmat. Uninitialized qtime caused crash if first data time < first CHATT time. 6 June 2004, added KSO Saturn coord sys. 10 Dec 2003. Erroneous calls to fswap, dswap, iswap replaced by FFSwapRecords. Includes Saturn Coordinate Systems too. Linked against CSPICE N0056 and FF_IGPP(12/01/03). March 2003 Added check on time of quaternion from CHATT file. Now if difference between data time and closest available quaternion is more than 60 sec, no coordinate transformation is calculated, and the data vector is not written to the output file. Previously, the last available quaternion was used for the s/c attitude, no matter how old it was with respect to the data time! End March/2003 comments. */ /* TransCal gets Cassini position and attitude and transforms a flatfile of calibrated magnetic field vectors from Cassini body-fixed to a user-selected coordinate system. Input times in flatfile are assumed SCLK(1958). The output flatfile contains times converted from SCLK to SCET, the transformed field vector, the field magnitude, and the XYZ position. The output times are seconds past 1/1/2000 12:00:00 TAI (= 11:59:28 UTC). */ #include #include #include #include #include #include #include #include extern "C" { #include "..\ff_user\time_igpp.h" #include "..\ff_user\ff_igpp.h" #include "..\spice_user\spiceusr.h" } /* In order to compile, the header file ff_igpp.h must be included. This file contains typedefs and function prototypes for the basic flatfile functions developed by UCLA/IGPP. The file ff_igpp.h "includes" the file time_igpp.h. This file contains function prototypes for the UCLA flatfile time functions. In order to link successfully, the files ff_igpp.obj and time_igpp.obj must be added to the TransCal project. These could be re-created by compiling ff_igpp.c and time_igpp.c, if necessary. All the above files are in the directory ff_user. */ /* TransCal uses many functions from the C version of SPICELIB. (The full CSPICE kit for the PC can be downloaded via anonymous ftp at naif.jpl.nasa.gov in directory /pub/naif/cspice. It is an executable of about 8MB that unpacks itself into 26MB of source, documentation, etc.) The top-level header file is spiceusr.h. This file includes other header files that contain the necessary typedefs and function prototypes that are needed to support the CSPICE application programming interface. In order for the linker to successfully resolve the references to the CSPICE functions, CSPICE.LIB must be added to the TransCal project, (Choose Project from the menu bar, then choose Add to Project, then choose Files, then in dialog box select Type of Files = Library. CSPICE.LIB is in spice_user.) */ #define FF_TIMEOFFSET 252460800 /*86400*(2+365*(1966-1958))*/ /* This offset will be added to input file times if epoch not in FFH */ #define NCOORDS 15 //Number of coordinate systems currently available #define EPOCH_IN Y1958 /*Epoch definition for input vector file and chatt*/ #define EPOCH_OUT J2000 /*Epoch definition for output flatfile*/ #define MAXKIND 10 /* Max Spice Kernel filenames of each kind */ #define QLIMIT 16 /* Limit (seconds) of quaternion interpolation */ // typedef struct { double time; float A1001; float A1002; float A1003; float A1004; float A1005; float A1006; float A1007; int STATUS; } ATTDATA; //record of CHATT (channelized attitude) flatfile //Function Prototypes defined later in this file FF_ID * open_in_transff (FILE *rpt); int get_transff ( FF_ID *idin, double *timeout, double b[3]); FF_ID * open_out_transff (FILE *rpt, FF_ID *idin, char * suffix1, char * suffix2); void write_transff (FF_ID *idout, double etime, double b[3], double pos[3]); void fftimeconv (double dtime, double *et, double *sclkdp); void get_qmatrix (FILE *rpt, double uclatime, double qmat[3][3], double *qtime, double *qnext, int *qerror); int get_ct ( double sclkdp, char * ref, double ct[3][3]); //void loadspicef (FILE *rpt); int loadkernels (int *numtext, int *numspk, int *numck, char textfiles[][72], char spkfiles[][72], char ckfiles[][72]); int get_attdata (FF_ID *idatt, ATTDATA *p); int leapsecs (double epochtai); FF_ID * open_input_ff (char * ffname); FF_ID * open_output_ff (char * ffname, int ncols, int recl); long get_nrows (FF_ID *id); int get_ncols (FF_ID *id); int get_recl (FF_ID *id); void exitspiceerr (FILE *rpt); void copy_info_abs (FF_ID *idin, FF_ID *idout, double time_first, double time_last); void trngse (double et, double tgse[3][3]); void trngsm (double et, double tgsm[3][3]); void trnrtn (double et, double trtn[3][3]); void trnjmxyz (double et, double tjmxyz[3][3]); void trnjmrtp (double et, double tjmrtp[3][3]); void trnsrtp (double et, double tsrtp[3][3]); void trnsunsat (double et, double tsunsat[3][3]); void trnkso (double et, double tkso[3][3]); //Global variables int swap_in = 1; //Input byte swapping is ON, but user can change int swap_out = 1; //Output byte swapping is ON, but user can change int epoch1958 = 0; /* Set to 1 if input time name includes "1958" */ double Jupiter_Dipole_Colat = 9.6; // Used by trnjmxyz and trnjmrtp double Jupiter_Dipole_Lon = 201.7; // CHANGE IF NECESSARY //MAIN void main () { //Variables FF_ID *idin; // flatfile pointer to input flatfile FF_ID *idout; //flatfile pointer to output flatfile char version[17] = "TransCal v050324"; char atime[32]; //text time string double time; //UCLA SCLK times double timefirst, timelast; //UCLA times converted to SCET double et, et1, et2; //Ephemeris Times char body[32]; //"399" for Earth, "299" for Venus, "10" for Sun int ibody; //Integer converted from 'body' string char att_opt[32] = " "; //"Q" for CHATT file attitudes, "C" for C-Kernel char line[80]; // Console input line double state[6]; //Position and velocity returned by SPKEZR double lt; //Light time returned by SPKEZR double pos[3]; //XYZ in J2000 double pout[3]; //XYZ or R/Lat/LON in selected Coord Sys double sclkdp; //Encoded SCLK time (ticks) int numtext, numspk, numck; /*Number of each kind of Spice kernel*/ char textfiles[MAXKIND][72]; /*Names of Spice text kernels*/ char spkfiles[MAXKIND][72]; /*Names of Spice SPK kernels*/ char ckfiles[MAXKIND][72]; /*Names of Spice CK kernels*/ int qerror = 0; double qtime; //Time of quaternion double qnext; //Time of next quaternion double qlimit = QLIMIT; //Limit to use quaternion time double qwarn = 0.; //Time of warning that quaternion time is past limit double qet; //qtime converted to ET double qdum; //dummy variable int icoord; // Number of Coordinate System chosen by user char coord[NCOORDS][32] = { "J2000", "ECLIPJ2000", "IAU_VENUS", "IAU_EARTH", "GSE", "GSM", "RTN", "CASSINI (no transformation)", "J3 (IAU_JUPITER)", "JMAG XYZ", "JMAG RTP", "KG (IAU_SATURN)", "SATURN RTP", "KSM (SUN-SATURN)", "KSO (SUN-SATURN_ORBIT_PLANE" }; char suffix[NCOORDS][7] = {"J2000", "ECLIP", "IAU_V", "IAU_E", "GSE", "GSM", "RTN", "S/C", "IAU_J", "JMXYZ", "JMRTP", "KG", "KRTP", "KSM", "KSO" }; /* Coord Sys Suffixes */ char suffixb[7]; /* Coord Sys Suffix for B Component Names*/ char suffixp[7]; /* Coord Sys Suffix for XYZ Position Names */ /* Position Coordinate System for RTN and S/C is J2000*/ int posindex[NCOORDS] = {0,1,2,3,4,5, 1,0, 8,9,9,11,11,13,14}; int ipos; //Position Coord Sys index double obliq = 23.43928982414818; //Obliquity for J2000 to ECLIPJ2000 double qmat[3][3]; //Trans from Cassini Body-fixed to J2000 double tmat[3][3]; //Trans from J2000 to selected Coord Sys double pmat[3][3]; /*Trans from J2000 to Position Coord Sys */ double ident[3][3] = { {1., 0., 0.,}, {0., 1., 0.}, {0., 0., 1.}}; //I double zmat[3][3] = { {0., 0., 0.}, {0., 0., 0.}, {0., 0., 0.}}; //ZERO double mout[3][3]; //Transformation matrix from Cassini to user-selected double b[3]; // Input B vector double bout[3]; //Output transformed B vector int foundc; //Found time in C Kernel; zero if not found int nread; //Result of reading input flatfile data int i; SpiceBoolean bfound; /* Body code found */ //Start cout << version << endl << endl; //Open Report File FILE *rpt = NULL; char report_def[256] = "..\\output\\TRANSFF.txt"; char report_name[256] = ""; cout << "Enter Report File Name: "; cin.getline (report_name, 255); if (strlen(report_name) == 0) strcpy (report_name, report_def); rpt = fopen (report_name, "w"); if (rpt == NULL) { cout << " Cannot open Report File " << report_name << endl; exit (1); } else { cout << "Opened Report File " << report_name << endl << endl; fprintf (rpt, "%s Report \n\n", version); } //Get Body cout << "Enter Reference Body ID " << endl; cout << " (399 for Earth, 699 for Saturn, 599 for Jupiter, 10 for Sun): "; cin.getline (body,32); fprintf (rpt, "Cassini Position is relative to Body %s \n", body); sscanf (body, "%d", &ibody); //Get Coordinate System cout << endl << "Available Coordinate Systems: " << endl; for (int j = 0; j= NCOORDS) { cout << "Invalid Coordinate System Number" << endl; cout << "Exit TransCal" << endl; fprintf (rpt, "Invalid Coordinate System Number %d\n", icoord); fprintf (rpt, "Exit TransCal"); exit (1); } ipos = posindex[icoord]; fprintf (rpt, "Position Coordinate System is %s \n\n", coord[ipos]); fprintf (rpt, "Vectors Will Be Transformed from Cassini body-fixed to %s \n\n", coord[icoord]); //Choose Attitude Source if (icoord != 7) { cout << endl << "Enter Q to use Quaternions from flatfile, C to use C-Kernel file: "; cin.getline (att_opt, 3); if (att_opt[0] == 'Q' || att_opt[0] == 'q') fprintf (rpt, "Attitudes will be calculated from Quaternion Flatfile \n\n"); if (att_opt[0] == 'C' || att_opt[0] == 'c') fprintf (rpt, "Attitudes will be calculated from C Kernel File(s) \n\n"); } cout << endl; //Set SPICE error option erract_c ("SET", 50, "RETURN"); errprt_c ("SET", 50, "NONE,SHORT,LONG,EXPLAIN"); //Load SPICE Kernel Files int nkernels = loadkernels(&numtext, &numspk, &numck, textfiles, spkfiles, ckfiles); if (failed_c()) exitspiceerr (rpt); cout << "Spice Kernels Loaded = " << nkernels << endl; fprintf (rpt, "%d Spice Text Kernel Files Loaded\n", numtext); for (i=0; i qlimit && fabs(qnext-time)>qlimit) || qnext-qtime > qlimit) { mequ_c (zmat, qmat); noatt++; if (qnext != qwarn) { fftimeconv (qtime, &qet, &qdum); timout_c (qet, "YYYY-DOYTHR:MN:SC.### (UTC)", 28,atime); fprintf (rpt, "Missing Attitude Info\n"); fprintf (rpt, " Previous CHATT Record = %s \n", atime); fftimeconv (qnext, &qet, &qdum); timout_c (qet, "YYYY-DOYTHR:MN:SC.### (UTC)", 28, atime); fprintf (rpt, " Current CHATT Record = %s \n\n", atime); qwarn = qnext; } } } if (qerror) break; if (att_opt[0] == 'C' || att_opt[0] == 'c') { foundc = get_ct (sclkdp, "J2000", qmat); if (failed_c()) exitspiceerr(rpt); if (!foundc) noatt++; } switch (icoord) { case 0 : /* J2000 */ mequ_c (ident, tmat); //J2000 - identity matrix mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 1 : /* ECLIPJ2000 */ rotate_c (obliq/deg, 1, tmat); //FROM J2000 to ECLIPJ2000 mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 2 : /* IAU_VENUS */ case 3 : /* IAU_EARTH */ tipbod_c ("J2000", ibody, et, tmat); //Venus and Earth mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 4 : /* GSE (Geocentric Solar Ecliptic) */ trngse (et, tmat); mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 5 : /* GSM (Geocentric Solar Magnetic) */ trngsm (et, tmat); mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 6 : /* RTN (Radial-Tangential-Normal, Sun Centered) */ trnrtn (et, tmat); rotate_c (obliq/deg, 1, pmat); /*Position in ECLIPJ2000 */ break; case 7 : /* Spacecraft fixed-body */ mequ_c (ident, pmat); /* Position in J2000 */ break; case 8 : /* IAU_JUPITER */ tipbod_c ("J2000", 599, et, tmat); mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 9 : /* JMAG XYZ (Jupiter Magnetic) */ trnjmxyz (et, tmat); mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 10 : /* JMAG RTP (Jupiter Magnetic R-Theta-Phi) */ trnjmrtp (et, tmat); trnjmxyz (et, pmat); /* Position in JMAG XYZ */ break; case 11 : /* IAU_SATURN */ tipbod_c ("J2000", 699, et, tmat); mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 12 : /* Saturn R-Theta-Phi */ trnsrtp (et, tmat); tipbod_c ("J2000", 699, et, pmat); /* Position in IAU_SATURN */ break; case 13 : /* Sun-Saturn */ trnsunsat (et, tmat); mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; case 14 : /* KSO, Sun, Saturn orbital plane */ trnkso (et, tmat); mequ_c (tmat, pmat); /* Position in same Coord Sys */ break; } if (failed_c()) exitspiceerr(rpt); mxv_c (pmat, pos, pout); //Transform Position from J2000 to CoordSys if (icoord != 7) { mxm_c (tmat, qmat, mout); //mout is FROM CasBody-fixed TO CoordSys mxv_c (mout, b, bout); // Transform mag vector to CoordSys } else mxv_c (ident, b, bout); // Or don't transform mag vector // Don't write if transformation matrix was zero because times not found if (bout[0] != 0. || bout[1] != 0. || bout[2] != 0.) { if (nwrit == 0) { et1 = et; //save first time written cout << endl << "Working."; } et2 = et; write_transff (idout, et, bout, pout); nwrit++; if (nwrit%20000 == 0) cout << "."; /* fprintf (rpt, "%s %11.1f %11.1f %11.1f", aatime, pout[0], pout[1],pout[2]); bmag = sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]); fprintf (rpt, " %10.2f %10.2f %10.2f %10.2f", bout[0], bout[1], bout[2], bmag); fprintf (rpt, "\n"); */ } } //end while //Finish timefirst = unitim_c (et1, "ET", "TAI"); timelast = unitim_c (et2, "ET", "TAI"); SetEpoch2JT(); int leap = leapsecs (timefirst); timefirst -= leap; timelast -= leap; copy_info_abs (idin, idout, timefirst, timelast); //Write TransCal Abstract char comment[72] = "#########"; ff_cput(idout->ffh, &comment[0]); sprintf (comment, "%-71s", version); ff_cput(idout->ffh, &comment[0]); sprintf (comment, "Input Flatfile was %-53s",idin->ffd->name); ff_cput(idout->ffh, &comment[0]); if (swap_out) sprintf (comment, "Output Byte Swapping was ON"); else sprintf (comment, "Output Byte Swapping was OFF"); ff_cput(idout->ffh, &comment[0]); sprintf (comment, "TransCal FFD times are TAI (International Atomic Time)"); ff_cput(idout->ffh, &comment[0]); sprintf (comment, "Reference Epoch is 2000 Jan 1 12:00:00 TAI (= 11:59:28 UTC)"); ff_cput(idout->ffh, &comment[0]); sprintf (comment, "TransCal header FIRST TIME and LAST TIME are UTC"); ff_cput(idout->ffh, &comment[0]); sprintf (comment, "Magnetic Field Coord Sys is %-43s", coord[icoord]); ff_cput(idout->ffh, &comment[0]); sprintf (comment, "Position Coord Sys is %-49s", coord[ipos]); ff_cput(idout->ffh, &comment[0]); bodc2n_c (ibody, 32, body, &bfound); sprintf (comment, "Position Reference Body is %-44s", body); ff_cput(idout->ffh, &comment[0]); sprintf (comment, "No Attitude information used"); if (att_opt[0] == 'Q' || att_opt[0] == 'q') sprintf (comment, "Spacecraft Attitude is from quaternions in CHATT file "); if (att_opt[0] == 'C' || att_opt[0] == 'c') sprintf (comment, "Spacecraft Attitude is from C Kernel "); ff_cput(idout->ffh, &comment[0]); ff_cput(idout->ffh, "Spice Kernel files:"); for (i=0; iffh, &textfiles[i][0]); for (i=0; iffh, &spkfiles[i][0]); for (i=0; iffh, &ckfiles[i][0]); ff_close(idin); ff_close(idout); fprintf (rpt, "\nNumber of Vectors Read: %d\n", nrecs); fprintf (rpt, "\nNumber of Vectors Written: %d\n", nwrit); if (noatt) fprintf (rpt, "Number of Vectors with no Attitude info: %d\n", noatt); fprintf (rpt, "End Report \n"); cout << endl; cout << "Number of Vectors Read: " << nrecs << endl; cout << "Number of Vectors Written: " << nwrit << endl; if (noatt) cout << "Missing Attitude information" << endl; if (qerror) cout << "Error reading CHATT file" << endl; cout << "End TRANSCAL" << endl; cout << endl << "Press Enter to Exit"; cin.getline(line,80); if (qerror) exit(1); }//End Main //__________________________________________________________________ //////////////////////////////////////////////////////////////////// FF_ID * open_in_transff (FILE *rpt) { FF_ID *idin = NULL; char flat_name[256] =""; //Name of input flatfile char line[80]; // console input line cout << endl << "Enter Input Flatfile Name (s/c coords): "; cin.getline (flat_name, 255); idin = open_input_ff (flat_name); if (idin == NULL) { fprintf (rpt, "Could not open %s\n", flat_name); exit(1); } /*end if*/ fprintf (rpt, "Opened Input Flatfile %s \n", flat_name); fprintf (rpt, " X, Y, Z components assumed in Columns 2, 3, 4 \n"); fprintf (rpt, " Input Flatfile Times assumed SCLK(1958) \n\n" ); int ncols = get_ncols (idin); if (ncols < 4) { fprintf (rpt, " Error on Input Flatfile: less than 4 columns\n"); cout << "Error on Input Flatfile: less than 4 columns"; exit (1); } cout << "Enter 1 to SWAP BYTES on INPUT Flatfile; 0 not to swap: "; cin.getline (line, 80); sscanf (line, "%d", &swap_in); if (swap_in) fprintf (rpt, "Input Byte Swapping is ON\n"); else fprintf (rpt, "Input Byte Swapping is OFF\n"); return (idin); } //End open_in_transff. //__________________________________________________________________ //////////////////////////////////////////////////////////////////// int get_transff ( FF_ID * idin, double *timeout, double b[3]) { //Returns result of reading flatfile data record. //Timeout = time just read; b = mag vector just read static int nrecs = 0; struct { double time; float stuff[255]; //Can't be any bigger???? } flat; char *p = (char *)♭ FILE *ffd = idin->ffd->fp; int recl = get_recl(idin); int xcol = 2; //Column number of X vector component int ycol = 3; //Y int zcol = 4; //Z int nread; nread = fread (&flat, recl, 1, ffd); if (swap_in) FFSwapRecords(p, recl); b[0] = flat.stuff[xcol-2]; b[1] = flat.stuff[ycol-2]; b[2] = flat.stuff[zcol-2]; nrecs ++; *timeout = flat.time; return (nread); } // End get_transff //____________________________________________________________________ ///////////////////////////////////////////////////////////////////// void get_qmatrix (FILE *rpt, double datatime, double qmat[3][3], double *qtime, double *qnext, int *qerror) { /* Rewritten 17 March 2005 to provide interpolation between rotation matrices from quaternions in consecutive CHATT records, Input ----- FILE *rpt is pointer to previously opened report file for error logging. double datatime = SCLK time of data field vector. Must be same time system as CHATT flatfile records. Output ------ double qmat[3][3] is rotation matrix FROM Cassini Body-Fixed TO J2000 at time = datatime. double *qtime is pointer to time of CHATT rec just before datatime. double *qnext is pointer to time of CHATT rec just after datatime. qtime < datatime <= qnext, except that qtime = first CHATT time if datatime is before time of first CHATT rec, and qnext = last CHATT time if datatime is after time of last CHATT rec. int *qerror is pointer to CHATT file status, nonzero if err or eof. Method ------ Interpolation method, following NAIF documentation Rotation.req. Let Q1 = rotation matrix from SC frame to J2000 at time t1. Let Q2 = rotation matrix from SC frame to J2000 at time t2. Want Q = rotation matrix from SC frame to J2000 at time t, t1 < t <= t2. Define R = product of Q2 times inverse of Q1. Note that for any vector v, R Q1 v = Q2 v. Calculate R: mxmt_c(Q2,Q1,R); Get rotation axis and angle of R: raxisa_c (R, axis, &angle); frac = fraction of angle rotated from t1 to t, frac = (t-t1)/(t2-t1). Now get rotation matrix for rotating frac*angle around axis: axisar_c(axis,frac*angle, RF); Then Q = product of RF times Q1: mxm_c(RF, Q1, Q); */ /* CHATT Flatfile */ static FF_ID *idatt = NULL; //Pointer to CHATT flatfile structure static ATTDATA att; // Storage for CHATT record static int nrecs = 0; //Number of records already read from CHATT data file char att_name[256]; //CHATT flatfile name static int nrows; //Number of rows in flatfile static int recl; //Number of bytes in flatfile record //Variables static double qmat1[3][3]; /* Q matrix from previous CHATT rec */ static double qmat2[3][3]; /* Q matrix from current CHATT rec */ static double qsave; /* time of previous CHATT rec */ double rmat[3][3]; /* Rotation matrix */ double fmat[3][3]; /* Fractional rotation matrix */ double axis[3]; /* Rotation axis */ double angle; /* Rotation angle */ double frac; /* fraction of rotation angle from t1 to t */ int result; /* Returned result of reading CHATT file */ double q[4]; /* Quaternion */ double qmag; /* Re-normalize quaternion */ int j; /* Loop index */ /*CHATT flatfile is not open..*/ if (idatt == NULL) { /*Get filename */ cout << endl << "Enter CHATT flatfile name: "; cin.getline (att_name, 255); /*Open CHATT flatfile */ idatt = open_input_ff(att_name); if (idatt == NULL) { fprintf (rpt, "Could not open %s\n", att_name); exit(1); } /* end if */ fprintf (rpt, "Opened CHATT Quaternion Flatfile %s\n\n", att_name); nrows = get_nrows(idatt); recl = get_recl(idatt); /*Read first record and calculate matrix from quaternion*/ result = get_attdata (idatt, &att); nrecs = 1; if (result <= 0) { *qerror = 1; fprintf (rpt, "Error on reading CHATT File, Record %d\n",nrecs); return; } /* fprintf (rpt, "CHATT TIME %20.3f\n", att.time); */ /* printf ("CHATT TIME %20.3f\n", att.time); */ q[0] = att.A1004; q[1] = att.A1001; q[2] = att.A1002; q[3] = att.A1003; qsave = att.time; qmag = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]); for (j=0; j<4; j++) q[j] /= qmag; q2m_c (q, qmat1); /*Read second record and calculate matrix */ result = get_attdata (idatt, &att); nrecs = 2; if (result <= 0) { *qerror = 1; fprintf (rpt, "Error on reading CHATT File, Record %d\n",nrecs); return; } /* fprintf (rpt, "CHATT TIME %20.3f\n", att.time); */ /* printf ("CHATT TIME %20.3f\n", att.time); */ q[0] = att.A1004; q[1] = att.A1001; q[2] = att.A1002; q[3] = att.A1003; qmag = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]); for (j=0; j<4; j++) q[j] /= qmag; q2m_c (q, qmat2); } /* end if */ /*CHATT flatfile is open */ while (datatime >= att.time && nrecs < nrows) { qsave = att.time; mequ_c (qmat2, qmat1); /* qmat1=qmat2 */ /*Get next CHATT record and calculate new Q matrix */ result = get_attdata (idatt, &att); nrecs ++; if (result <= 0) { *qerror = 1; fprintf (rpt, "Error on reading CHATT File, Record %d\n",nrecs); return; } if (nrecs==nrows) fprintf(rpt, "Last Record read from CHATT File\n"); q[0] = att.A1004; q[1] = att.A1001; q[2] = att.A1002; q[3] = att.A1003; qmag = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]); for (j=0; j<4; j++) q[j] /= qmag; q2m_c (q, qmat2); } // end while /* Interpolate */ mxmt_c (qmat2, qmat1, rmat); raxisa_c (rmat, axis, &angle); frac = 1.; if (att.time != qsave) frac = (datatime-qsave)/(att.time-qsave); if (frac < 0.) frac = 0.; if (frac > 1.) frac = 1.; axisar_c (axis, frac*angle, fmat); mxm_c (fmat, qmat1, qmat); *qnext = att.time; *qtime =qsave; return; } //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// void fftimeconv (double dtime, double *et, double *sclkdp) { // //Changed Nov 2003. Reconstructs SCLK and calls scencd to get sclkdp. //Must change FF_TIMEOFFSET to 0 if new flatfile times are Epoch 1958.0 //Returns Ephemeris Time corresponding to UCLA flatfile time "dtime", assumed SCLK; //Also returns "Encoded SCLK" to use with C Kernel files extern int epoch1958; /* 1 if input time col name contains "1958" */ int sc = -82; //Cassini double ethold; double sclkhold; double ftick; /* Fraction of tick */ int iclk1, iclk2; char clkstr[17]; iclk1 = (int) dtime; iclk2 = (int) (256*(dtime-iclk1)); ftick = 256*(dtime-iclk1) - iclk2; if (epoch1958 == 0) iclk1 += FF_TIMEOFFSET; sprintf (clkstr,"1/%10.10d:%3.3d", iclk1, iclk2); scencd_c (sc, clkstr, &sclkhold); sct2e_c (sc, sclkhold, ðold); *sclkdp = sclkhold; /* Add back fraction of tick */ ethold += ftick/256.; *et = ethold; return; } //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// int get_ct (double sclkdp, char * ref, double ct[3][3]) { //Gets transpose of C-matrix in reference frame "ref" at time given by // encoded SCLK "sclkdp" // Returns 1 if time found in C kernel, 0 if not found. char tolstr[] = "0:001"; // Minimum tolerance of 1 tick = 1/256 sec double tol; //Tolerance in SCLK ticks double clkout; //Returned time in SCLK ticks from ckgp double cmat[3][3]; //C matrix: from ref frame to s/c frame int sc = -82; //Cassini int inst = -82000; //Cassini body-fixed frame SpiceBoolean found; int i,j; //Start sctiks_c (sc, tolstr, &tol); ckgp_c (inst, sclkdp, tol, ref, cmat, &clkout, &found); if (found) xpose_c (cmat, ct); //Transpose cmat else { for (i=0; i<3; i++) { for (j=0; j<3; j++) ct[i][j] = 0.0; } } return ((int)found); } // End get_ct //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// int loadkernels (int *numtext, int *numspk, int *numck, char textfiles[][72], char spkfiles[][72], char ckfiles[][72]) { //Reads input file to get Spice filenames. // Loads Spice text and SPK and CK kernel files char inputname[255] = ""; //Name of Spice metafile with KERNELS_TO_LOAD cout << "Enter Name of Spice Metafile with KERNELS_TO_LOAD: " << endl; cin.getline (inputname,255); furnsh_c (inputname); // furnsh_c ("g:\\test03\\klist.txt"); #define TYPLEN 32 #define SRCLEN 128 char kertype[TYPLEN]; char srcname[SRCLEN]; SpiceInt tot = 0; SpiceInt which; SpiceInt handle; SpiceBoolean found; int sum =0; ktotal_c("TEXT", &tot); *numtext = tot; sum += tot; if (tot > MAXKIND) tot = MAXKIND; for (which=0; which < tot; which++) { kdata_c (which, "TEXT", 72, TYPLEN, SRCLEN, &textfiles[which][0], kertype, srcname, &handle, &found); } ktotal_c("SPK", &tot); *numspk = tot; sum += tot; if (tot > MAXKIND) tot = MAXKIND; for (which=0; which < tot; which++) { kdata_c (which, "SPK", 72, TYPLEN, SRCLEN, &spkfiles[which][0], kertype, srcname, &handle, &found); } ktotal_c("CK", &tot); *numck = tot; sum += tot; if (tot > MAXKIND) tot = MAXKIND; for (which=0; which < tot; which++) { kdata_c (which, "CK", 72, TYPLEN, SRCLEN, &ckfiles[which][0], kertype, srcname, &handle, &found); } return (sum); }//end loadkernels //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// FF_ID * open_input_ff (char * ffname) { /* 2 Dec 2004. If file open fails and standard input is console, user is asked whether to try again to enter input filename. */ //Opens input flatfile header and data files //Returns toplevel flatfile pointer. //Now returns instead of exiting if file open fails. -jw, june 2001. FF_ID *id; FF_O_PARAM param; param.status = FF_READ | FF_EXIST; param.ncols = 0; param.recl = 0; param.nbufs = 1; param.epoch = EPOCH_IN; id = ff_open (ffname, ¶m); int tryagain = _isatty(_fileno(stdin)); int ntries = 0; while (id == NULL && tryagain && ntries < 5) { cout << endl << "Cannot open " << ffname << endl; cout << "Enter 1 to try again, 0 to exit: "; cin >> tryagain; if (tryagain) { cout << "Enter Input Flatfile name: "; cin.get(); cin.getline (ffname, 255); id = ff_open (ffname, ¶m); ntries ++; } /* end if */ } /* end while */ return (id); } /* End open_input_ff() */ //------------------------------------------------------------------- FF_ID * open_output_ff (char * ffname, int ncols, int recl) { //Opens output flatfile header and data files //Returns toplevel flatfile pointer. //Now returns instead of exiting if file open fails. -jw, june 2001. FF_ID *id; FF_O_PARAM param; param.status = FF_CREATE | FF_WRITE | FF_EXIST; param.ncols = ncols; param.recl = recl; param.nbufs = 1; param.epoch = EPOCH_OUT; id = ff_open (ffname, ¶m); int tryagain = _isatty(_fileno(stdin)); int ntries = 0; while (id == NULL && tryagain && ntries < 7) { cout << endl << "Cannot open Output " << ffname << endl; cout << "Enter 1 to try again, 0 to exit: "; cin >> tryagain; if (tryagain) { cout << "Enter Output Flatfile name: "; cin.get(); cin.getline (ffname, 255); id = ff_open (ffname, ¶m); ntries ++; } /* end if */ } /* end while */ return (id); } /* End open_output_ff() */ //------------------------------------------------------------------- long get_nrows (FF_ID *id) { return (id ->ffd->nrows); } //------------------------------------------------------------------- int get_ncols (FF_ID *id) { return (id ->ffd->ncols); } //------------------------------------------------------------------- int get_recl (FF_ID *id) { return (id ->ffd->recl); } ///////////////////////////////////////////////////////////////////// int get_attdata (FF_ID *id, ATTDATA *p) { //Gets next data record from open CHATT flatfile //Swaps bytes if Global variable swap_in is non-zero. FILE *ffd = id->ffd->fp; //pointer to CHATT data flatfile ATTDATA natt; /*Storage for next chatt record*/ ATTDATA *pn = &natt; // char *pa = (char *)p; int recl = get_recl(id); int numread = fread (pn, recl, 1, ffd); if (numread > 0 && swap_in) FFSwapRecords ((char *)pn, recl); p->time = natt.time; p->A1001 = natt.A1001; p->A1002 = natt.A1002; p->A1003 = natt.A1003; p->A1004 = natt.A1004; return (numread); } //End get_attdata. //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// void trngse (double et, double tgse[3][3]) { // et is Ephemeris Time // Returns tgse, transformation from J2000 to GSE at time et double lt; //Storage for light time returned from spkezr double state_sun[6]; double sunhat[3]; //normalized double obliq = 23.439291; //Obliquity in J2000 double eclpole[3]; spkezr_c ("SUN", et, "J2000", "NONE", "EARTH", state_sun, <); vhat_c (state_sun, sunhat); eclpole[0] = 0.; eclpole[1] = -sin(obliq/dpr_c()); eclpole[2] = cos(obliq/dpr_c()); twovec_c (sunhat, 1, eclpole, 3, tgse); return; } // End trngse //_____________________________________________________________________ /////////////////////////////////////////////////////////////////////// void trngsm (double et, double tgsm[3][3]) { // et is Ephemeris Time // Returns tgsm, transformation from J2000 to GSM at time et // X GSM = XGSE = Earth-to-Sun // Z GSM lies in plane of X and Earth's Magnetic North Pole double state_sun[6]; // Earth to Sun in IAU_EARTH double sunhat[3]; //normalized Sun double lt; //light time returned from spkezr // double mag_colat = 11.2; //from EJS, 1980.0 // double mag_lon = -70.76; //from EJS, 1980.0 // Refer to http://nssdc.gsfc.nasa.gov/space/model/magnetos/tsygan.html double mag_colat = 10.547; //GEOPACK (subroutine RECALC), double mag_lon = -71.719; // IGRF1995 model for 1999,Day 230 // Refer to http://www.ngdc.noaa.gov/IAGA/wg8/wg8.htm for IGRF2000 model // double mag_colat = 10.476; //IGRF2000 for 1999,Day 230 // double mag_lon = -71.559; // double magpole[3]; // Mag Pole in IAU_EARTH double tipm[3][3]; //from J2000 to IAU_EARTH double mout[3][3]; // from IAU_EARTH to GSM sphrec_c (1.0, mag_colat/dpr_c(), mag_lon/dpr_c(), magpole); spkezr_c ("SUN", et, "IAU_EARTH", "NONE", "EARTH", state_sun, <); vhat_c (state_sun, sunhat); //Normalize big sun vector twovec_c (sunhat, 1, magpole, 3, mout); //from IAU_EARTH to GSM tipbod_c ("J2000", 399, et, tipm); //from J2000 to IAU_EARTH mxm_c (mout, tipm, tgsm); return; } //End tgsm //_____________________________________________________________________ /////////////////////////////////////////////////////////////////////// void trnrtn (double et, double trtn[3][3]) { // et is Ephemeris Time // Returns trtn, transformation from J2000 to RTN at time et double state_sun[6]; //Cassini to Sun in J2000 double sunpole[3]; double sunpole_ra = 286.13; //from pck00006 double sunpole_dec = 63.87; double lt; double r[3]; spkezr_c ("SUN", et, "J2000", "NONE", "-82", state_sun, <); vminus_c (state_sun, r); //r is Sun to Cassini latrec_c (1.0, sunpole_ra/dpr_c(), sunpole_dec/dpr_c(), sunpole); twovec_c (r, 1, sunpole, 3, trtn); return; } //End trtn //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// void trnjmxyz (double et, double tjmxyz[3][3]) { // et is Ephemeris Time // tjmxyz is returned transformation from J2000 to JMAGXYZ // JMAGXYZ: Z is Northward in direction of Jovian Dipole // X lies in System III X-Z plane // Y completes orthogonal right hand system // Get global dipole angles and convert to radians double dcolat = Jupiter_Dipole_Colat/dpr_c(); double dlon = Jupiter_Dipole_Lon/dpr_c(); double dipole[3]; double tipm[3][3]; //from J2000 to IAU_Jupiter double mout[3][3]; //from IAU_Jupiter to JMAGXYZ double yaxis[3] = {0., 1., 0.}; double one = 1.0; sphrec_c (one, dcolat, dlon, dipole); twovec_c (dipole, 3, yaxis, 2, mout); tipbod_c ("J2000", 599, et, tipm); mxm_c (mout, tipm, tjmxyz); return; } //End trnjmxyz //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// void trnjmrtp (double et, double tjmrtp[3][3]) { // et is Ephemeris Time // tjmrtp is returned transformation from J2000 to JMAGRTP // JMAGRTP: R is Jupiter-to-Cassini direction // T (Theta) lies in R - Jovian Dipole plane, positive Southward // P completes orthogonal right hand system // Get global dipole angles and convert to radians double dcolat = Jupiter_Dipole_Colat/dpr_c(); double dlon = Jupiter_Dipole_Lon/dpr_c(); double dipole[3]; double one = 1.0; double ndipole[3]; double tipm[3][3]; //from J2000 to IAU_Jupiter double mout[3][3]; //from IAU_Jupiter to JMAGRTP double state[6]; double lt; // light time returned by spkezr double rjup[3]; //Jupiter-to Cassini sphrec_c (one, dcolat, dlon, dipole); vminus_c (dipole, ndipole); spkezr_c ("599", et, "IAU_JUPITER", "NONE", "-82", state, <); vminus_c (state, rjup); twovec_c (rjup, 1, ndipole, 2, mout); tipbod_c ("J2000", 599, et, tipm); mxm_c (mout, tipm, tjmrtp); return; } //End trnjmrtp //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// void trnsrtp (double et, double tsrtp[3][3]) { // et is Ephemeris Time // tsrtp is returned transformation from J2000 to SAT_RTP // SAT_RTP: R is Saturn-to-Cassini direction // T (Theta) lies in R - Saturn_Axis plane, positive Southward // P (Phi) completes orthogonal right hand system double southpole[3] = {0., 0., -1.}; double tipsat[3][3]; //from J2000 to IAU_SATURN double mout[3][3]; //from IAU_SATURN to SAT_RTP double state[6]; double lt; // light time returned by spkezr double rsat[3]; //Saturn-to Cassini in IAU_SATURN spkezr_c ("699", et, "IAU_SATURN", "NONE", "-82", state, <); vminus_c (state, rsat); twovec_c (rsat, 1, southpole, 2, mout); tipbod_c ("J2000", 699, et, tipsat); mxm_c (mout, tipsat, tsrtp); return; } //End trnsrtp //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// void trnsunsat (double et, double tsunsat[3][3]) { // et is Ephemeris Time // tsunsat is returned transformation from J2000 to SUN-SATURN // SUN-SATURN: X is Saturn-to-Sun direction // Z lies in X - Saturn_Rotation_Axis plane, positive northward // Y completes orthogonal right hand system double axis[3] = {0., 0., 1.}; double axisj[3]; //Saturn Rotation Axis in J2000 double tipsat[3][3]; //from J2000 to IAU_SATURN double state[6]; double lt; // light time returned by spkezr spkezr_c ("SUN", et, "J2000", "NONE", "699", state, <); tipbod_c ("J2000", 699, et, tipsat); mtxv_c (tipsat, axis, axisj); twovec_c (state, 1, axisj, 3, tsunsat); return; } //End trnsunsat // ////////////////////////////////////////////////////////////////////// void trnkso (double et, double tkso[3][3]) { // et is Ephemeris Time // tkso is returned transformation from J2000 to KSO // KSO: X is Saturn-to-Sun direction // Z lies in X - O plane, // where O is north perpendicular to Sat Orbit Plane // Y completes orthogonal right hand system, // Y = -Velocity of Saturn, positive toward Dusk double x[3], y[3], z[3]; double state[6]; double lt; // light time returned by spkezr /*Get X, Saturn-to-Sun vector */ spkezr_c ("SUN", et, "J2000", "NONE", "699", state, <); y[0] = state[3]; y[1] = state[4]; y[2] = state[5]; vhat_c (state, x); vhat_c (y, y); // printf ("X %f %f %f \n", x[0],x[1],x[2]); // printf ("Y %f %f %f \n", state[3],state[4],state[5]); /* Y is Velocity of Sun wrt Saturn = -V of Saturn */ /* Y not exactly perpendicular to X, so use Z = X x Y */ ucrss_c (x, y, z); twovec_c (x, 1, z, 3, tkso); return; } //End trnkso //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// FF_ID * open_out_transff (FILE *rpt, FF_ID *idin, char * suffix1, char * suffix2) { //Opens an output flatfile for the transformed vectors and position extern int epoch1958; /* Will be set to 1 if input time column name contains "1958" */ FF_COL_DESC indesc; /* Input Column Descriptor */ FF_COL_DESC coldesc[] = { {1, "TIME_TAI", "SEC", "(source)", "T", 0}, {2, "BX_", "nT", "(source)", "R", 8}, {3, "BY_", "nT", "(source)", "R", 12}, {4, "BZ_", "nT", "(source)", "R", 16}, {5, "BTOTAL", "nT", "(source)", "R", 20}, {6, "X_", "km", "Spice SPK", "R", 24}, {7, "Y_", "km", "Spice SPK", "R", 28}, {8, "Z_", "km", "Spice SPK", "R", 32} }; FF_ID *idout = NULL; // output flatfile pointer int ncols = 8; int recl = 4 + 4*ncols; int j, ierror; /* Copy Source from input header to Column Descriptor*/ /* for Time and B Components */ for (j=0;j<4; j++) { indesc.ncol = j+1; ierror = ff_hget (idin->ffh, &indesc); if (ierror != 72) fprintf (rpt, "Error from FF_HGET on Input COL %d", j+1); strcpy (coldesc[j].source, indesc.source); /* Check if input time column name contains "1958" */ if (j==0 && strstr(indesc.name, "1958") != NULL) epoch1958 = 1; } /* Copy BZ Source to BTOTAL Source */ strcpy (coldesc[4].source, indesc.source); /* Append Suffix to B Component Names */ for (j=1; j<4; j++) strcat(coldesc[j].name, suffix1); /* Append Suffix to Position Component Names */ for (j=5; j<8; j++) strcat(coldesc[j].name, suffix2); //Get output flatfile name and open flatfile char out_name[256]; char line[80]; //console input line cout << endl << "Enter Output Flatfile Name: "; cin.getline (out_name, 255); idout = open_output_ff (out_name, ncols, recl); if (idout == NULL) { fprintf (rpt,"Could not open Output %s\n", out_name); exit(1); } // end if fprintf (rpt, "\nOpened Output Flatfile %s \n", out_name); cout << "Enter 1 to SWAP BYTES on OUTPUT Flatfile; 0 not to swap: "; cin.getline (line, 80); sscanf (line, "%d", &swap_out); if (swap_out) fprintf (rpt, "Output Byte Swapping is ON\n\n"); else fprintf (rpt, "Output Byte Swapping is OFF\n\n"); // Write column descriptors to flatfile header for (j=0; j<8; j++) { ff_hput (idout->ffh, &coldesc[j]); } return (idout); } //End open_transff //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// void write_transff (FF_ID *idout, double etime, double b[3], double pos[3]) { // Writes an output flatfile record : SCET, transformed vector, magnitude, // position FILE *ffd = idout->ffd->fp; struct { double time; float fval[7]; } out; // char atime[32]; //time string char *p = (char *)&out; int i; int num; int recl = 36; out.time = unitim_c (etime, "ET", "TAI"); // 1 JAN 2000 12:00:00 UTC = 32.0 TAI ! // 1 JAN 2000 11:59:28 UTC = 0.0 TAI for (i=0; i<3; i++) { out.fval[i] = (float) b[i]; out.fval[i+4] = (float) pos[i]; } out.fval[3] = (float) sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]); //Field magnitude if (swap_out) FFSwapRecords(p, recl); num = fwrite (&out, recl, 1, ffd); return; } //End write_transff //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// void copy_info_abs (FF_ID *idin, FF_ID *idout, double time_first, double time_last) { /*Copies info parameters from one flatfile header file to another, except that info.first_time and info.last_time are not copied but are set to time_first and time_last. Also copies the abstract, that is, all header records after the column descriptors that are not included in the info keywords. */ FF_H_INFO info; if (ff_hgetinfo(idin->ffh, &info) < 0) { cout << "Error from ff_hgetinfo.\nExit.\n"; exit(1); } info.first_time = time_first; info.last_time = time_last; if (ff_hputinfo (idout->ffh, &info) < 0) { cout << "error from ff_hputinfo.\nExit.\n"; exit(1); } if (ff_hcopyabstract (idin->ffh, idout->ffh) < 0) { cout << "Error from ff_copyabstract.\nExit.\n"; exit(1); } } // End copy_info_abs(). //____________________________________________________________________ ////////////////////////////////////////////////////////////////////// int leapsecs (double epochtai) { /* Returns number of leap seconds at epoch "epochtai". Input: epochtai is an epoch of type TAI (International Atomic Time); that is, seconds past 2000 Jan 1 12:00:00 TAI. (2000 Jan 1 12:00:00 TAI is 2000 Jan 1 11:59:28 UTC) For epochtai = 0., leapsecs returns the value 32. One may think of TAI as "ahead" of UTC by the value of leapsecs. Method: Leap second arrays are retrieved from the kernel pool variable "DELTET/DELTA_AT" using CSPICE function gdpool_C. A leap second kernel such as naif0007.tls must have been loaded. Reference: CSPICE function deltet.c or SPICE routine DELTET.FOR, */ #define MAXLEAP 200 SpiceDouble dleap[2*MAXLEAP]; SpiceInt nleap; SpiceBoolean found; int leap, i; gdpool_c ("DELTET/DELTA_AT", 0, 2*MAXLEAP, &nleap, dleap, &found); leap = (int) dleap[0] -1; for (i=0; i= dleap[i+1]) leap = (int) dleap[i]; } return (leap); } /*end leapsecs()*/ //__________________________________________________________________ //////////////////////////////////////////////////////////////////// void exitspiceerr (FILE *rpt) { /* Writes SPICE error messages to report file and exits */ #define LENSHORT 26 #define LENLONG 321 #define LENEXPLN 81 char shortmsg[LENSHORT]; char longmsg[LENLONG]; char explnmsg[LENEXPLN]; getmsg_c ("SHORT", LENSHORT, shortmsg); fprintf (rpt, "\nCSPICE Error: %s", shortmsg); getmsg_c ("EXPLAIN", LENEXPLN, explnmsg); fprintf (rpt, "\n%s", explnmsg); getmsg_c ("LONG", LENLONG, longmsg); fprintf (rpt, "\n%s", longmsg); fprintf (rpt, "\nEXIT Transcal"); _flushall(); exit (1); } /*End exitspicerr */ //////////////////////////////////////////////////////////////////////