IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41714


Ignore:
Timestamp:
Jul 13, 2021, 8:52:42 AM (5 years ago)
Author:
eugene
Message:

add apparent magnitude calculation (based on phase curve for Mercury)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/tools/src/mpcorb_predict.c

    r41713 r41714  
    2626  double eccentricity;
    2727  double semimajor_a;
     28  double Hmag;
     29  int    Nobs;  // number of observations
     30  int    Nopp;  // number of oppositions
    2831  double mjdMin;
    2932  double mjdMax;
    … …  
    3942  double Robs;
    4043  double Dobs;
     44  double Mobs; // apparent mag
    4145  double dist;
    4246} PlanetDatum;
    4347
    44 int Shutdown (char *format, ...);
    45 int mpcorb_parseline (char *line, Planets *planet, int Nmax);
     48void mpcorb_trange (int argc, char **argv);
     49void mpcorb_moment (int argc, char **argv);
     50
     51PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc);
     52
    4653Planets *mpcorb_read_text (char *filename, int *nplanets);
    4754Planets *mpcorb_read_fits (char *filename, int *nplanets);
    … …  
    4956void     mpcorb_refs_fits (char *filename, PlanetDatum *planets, int Nplanets);
    5057
    51 PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc);
    52 
    53 void mpcorb_trange (int argc, char **argv);
    54 void mpcorb_moment (int argc, char **argv);
     58int Shutdown (char *format, ...);
     59int mpcorb_parseline (char *line, Planets *planet, int Nmax);
     60
     61void solar_coords (double mjdObs, double *Rsun, double *Dsun, double *distSun);
    5562
    5663# define PS1_LONGITUDE -156.2559041668965 /* degrees East */
    … …  
    5865# define PS1_ALTITUDE  3067.7 /* meters */
    5966
     67# define SKIP_BAD FALSE
     68
     69// save solar position for apparent mag calc
     70static double R_SUN = NAN;
     71static double D_SUN = NAN;
     72static double DIST_SUN = NAN;
     73
    6074int main (int argc, char **argv) {
     75
     76  // XXX test
     77  if (TESTING) {
     78    double mjdObs   = atof (argv[1]);
     79    char  *filename = argv[2];
     80 
     81    // save the coordinates of the sun:
     82    solar_coords (mjdObs, &R_SUN, &D_SUN, &DIST_SUN);
     83   
     84    // load the planets
     85    int Nplanets = 0;
     86    Planets *planets = mpcorb_read_text (filename, &Nplanets);
     87    for (int i = 0; i < 10; i++) {
     88      // if it is in the region, use the more accurate calculation
     89      mpcorb_predict (&planets[i], mjdObs, TRUE);
     90    }
     91    exit (0);
     92  }
    6193
    6294  if ((argc != 9) && (argc != 10)) goto usage;
    … …  
    100132  for (int i = 0; i < Nplanets; i++) {
    101133    if (i % 1000 == 0) fprintf (stderr, ".");
     134
     135    // XXX make these quality cuts user options
     136    if (SKIP_BAD && ((planets[i].Nobs < 20) || (planets[i].Nopp < 2))) continue;
    102137
    103138    // is the object in the region for the start of the period?
    … …  
    153188  char  *output   = argv[8];
    154189
     190  // save the coordinates of the sun:
     191  solar_coords (mjdObs, &R_SUN, &D_SUN, &DIST_SUN);
     192
    155193  int Nplanets = 0;
    156194  Planets *planets = mpcorb_read_fits (filename, &Nplanets);
    … …  
    188226
    189227  exit (0);
    190 }
    191 
    192 // STATUS is value expected for success
    193 # define CHECK_STATUS(STATUS,MSG,...) if (!(STATUS)) { Shutdown (MSG, __VA_ARGS__); }
    194 
    195 // write minor planet positions (Rref,Dref,Mref) to a FITS table
    196 void mpcorb_refs_fits (char *filename, PlanetDatum *planets, int Nplanets) {
    197 
    198   Header header;
    199   Header theader;
    200   Matrix matrix;
    201   FTable ftable;
    202 
    203   gfits_init_header (&header);
    204   header.extend = TRUE;
    205   gfits_create_header (&header);
    206   gfits_create_matrix (&header, &matrix);
    207 
    208   gfits_create_table_header (&theader, "BINTABLE", "DATA");
    209 
    210   gfits_define_bintable_column (&theader, "8A", "ID",   "name",    "none",        1.0, 0.0);
    211   gfits_define_bintable_column (&theader,  "D", "Rref", "RA",      "degrees",     1.0, 0.0);
    212   gfits_define_bintable_column (&theader,  "D", "Dref", "DEC",     "degrees",     1.0, 0.0);
    213   gfits_define_bintable_column (&theader,  "D", "dist", "DIST",    "AU",          1.0, 0.0);
    214   gfits_define_bintable_column (&theader,  "D", "Mref", "Mag",     "magnitudes",  1.0, 0.0);
    215 
    216   // generate the output array that carries the data
    217   gfits_create_table (&theader, &ftable);
    218 
    219   // create intermediate storage arrays
    220   ALLOCATE_PTR (ID,                     char, Nplanets*8);
    221   ALLOCATE_PTR (Rref,                 double, Nplanets);
    222   ALLOCATE_PTR (Dref,                 double, Nplanets);
    223   ALLOCATE_PTR (dist,                 double, Nplanets);
    224   ALLOCATE_PTR (Mref,                 double, Nplanets);
    225 
    226   // set intermediate storage arrays
    227   for (int i = 0; i < Nplanets; i++) {
    228     Rref[i]        = planets[i].Robs;
    229     Dref[i]        = planets[i].Dobs;
    230     dist[i]        = planets[i].dist;
    231     Mref[i]        = 16.0; // need to add this to prediction
    232     memcpy(&ID[i*8], planets[i].ID, 8);
    233   }
    234 
    235   // add the columns to the output array
    236   gfits_set_bintable_column (&theader, &ftable, "ID",   ID,   Nplanets);
    237   gfits_set_bintable_column (&theader, &ftable, "Rref", Rref, Nplanets);
    238   gfits_set_bintable_column (&theader, &ftable, "Dref", Dref, Nplanets);
    239   gfits_set_bintable_column (&theader, &ftable, "dist", dist, Nplanets);
    240   gfits_set_bintable_column (&theader, &ftable, "Mref", Mref, Nplanets);
    241 
    242   // free intermediate storage arrays
    243   free (ID);
    244   free (Rref);
    245   free (Dref);
    246   free (dist);
    247   free (Mref);
    248 
    249   // open the output file
    250   FILE *f = fopen (filename, "w");
    251   if (!f) Shutdown ("ERROR: cannot open mpcorb file for output %s\n", filename);
    252 
    253   int status;
    254   status = gfits_fwrite_header  (f, &header);
    255   CHECK_STATUS (status, "ERROR: cannot write header for mpcorbs %s\n", filename);
    256 
    257   status = gfits_fwrite_matrix  (f, &matrix);
    258   CHECK_STATUS (status, "ERROR: cannot write matrix for mpcorbs %s\n", filename);
    259 
    260   status = gfits_fwrite_Theader (f, &theader);
    261   CHECK_STATUS (status, "ERROR: cannot write table header for mpcorbs %s\n", filename);
    262 
    263   status = gfits_fwrite_table  (f, &ftable);
    264   CHECK_STATUS (status, "ERROR: cannot write table data for mpcorbs %s\n", filename);
    265 
    266   gfits_free_header (&header);
    267   gfits_free_matrix (&matrix);
    268   gfits_free_header (&theader);
    269   gfits_free_table (&ftable);
    270 
    271   int fd = fileno (f);
    272 
    273   status = fflush (f);
    274   CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
    275 
    276   status = fsync (fd);
    277   CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
    278 
    279   status = fclose (f);
    280   CHECK_STATUS (!status, "ERROR: problem closing mpcorbs file %s\n", filename);
    281 
    282   return;
    283 }
    284 
    285 // write minor planet orbital elements (plus (R,D,T)[min,max]) to a FITS table
    286 void mpcorb_save_fits (char *filename, Planets *planets, int Nplanets) {
    287 
    288   Header header;
    289   Header theader;
    290   Matrix matrix;
    291   FTable ftable;
    292 
    293   gfits_init_header (&header);
    294   header.extend = TRUE;
    295   gfits_create_header (&header);
    296   gfits_create_matrix (&header, &matrix);
    297 
    298   gfits_create_table_header (&theader, "BINTABLE", "MPCORB");
    299 
    300   gfits_define_bintable_column (&theader, "8A", "ID",           "name",         "none",     1.0, 0.0);
    301   gfits_define_bintable_column (&theader,  "D", "EPOCH",        "epoch",        "MJD",      1.0, 0.0);
    302   gfits_define_bintable_column (&theader,  "D", "MEAN_ANOMALY", "mean_anomaly", "degrees",  1.0, 0.0);
    303   gfits_define_bintable_column (&theader,  "D", "PERIHELION",   "perihelion",   "degrees",  1.0, 0.0);
    304   gfits_define_bintable_column (&theader,  "D", "ASCEND_NODE",  "ascend_node",  "degrees",  1.0, 0.0);
    305   gfits_define_bintable_column (&theader,  "D", "INCLINATION",  "inclination",  "degrees",  1.0, 0.0);
    306   gfits_define_bintable_column (&theader,  "D", "ECCENTRICITY", "eccentricity", "unitless", 1.0, 0.0);
    307   gfits_define_bintable_column (&theader,  "D", "SEMIMAJOR_A",  "semimajor_a",  "AU",       1.0, 0.0);
    308   gfits_define_bintable_column (&theader,  "D", "MJD_MIN",      "mjdMin",       "MJD",      1.0, 0.0);
    309   gfits_define_bintable_column (&theader,  "D", "MJD_MAX",      "mjdMax",       "MJD",      1.0, 0.0);
    310   gfits_define_bintable_column (&theader,  "D", "RAT_MIN",      "RatMin",       "degrees",  1.0, 0.0);
    311   gfits_define_bintable_column (&theader,  "D", "RAT_MAX",      "RatMax",       "degrees",  1.0, 0.0);
    312   gfits_define_bintable_column (&theader,  "D", "DAT_MIN",      "DatMin",       "degrees",  1.0, 0.0);
    313   gfits_define_bintable_column (&theader,  "D", "DAT_MAX",      "DatMax",       "degrees",  1.0, 0.0);
    314 
    315   // generate the output array that carries the data
    316   gfits_create_table (&theader, &ftable);
    317 
    318   // create intermediate storage arrays
    319   ALLOCATE_PTR (ID,                     char,   Nplanets*8);
    320   ALLOCATE_PTR (epoch,                  double, Nplanets);
    321   ALLOCATE_PTR (mean_anomaly,           double, Nplanets);
    322   ALLOCATE_PTR (perihelion,             double, Nplanets);
    323   ALLOCATE_PTR (ascend_node,            double, Nplanets);
    324   ALLOCATE_PTR (inclination,            double, Nplanets);
    325   ALLOCATE_PTR (eccentricity,           double, Nplanets);
    326   ALLOCATE_PTR (semimajor_a,            double, Nplanets);
    327   ALLOCATE_PTR (mjdMin,                 double, Nplanets);
    328   ALLOCATE_PTR (mjdMax,                 double, Nplanets);
    329   ALLOCATE_PTR (RatMin,                 double, Nplanets);
    330   ALLOCATE_PTR (RatMax,                 double, Nplanets);
    331   ALLOCATE_PTR (DatMin,                 double, Nplanets);
    332   ALLOCATE_PTR (DatMax,                 double, Nplanets);
    333 
    334   // set intermediate storage arrays
    335   for (int i = 0; i < Nplanets; i++) {
    336     epoch[i]           = planets[i].epoch;
    337     mean_anomaly[i]    = planets[i].mean_anomaly;
    338     perihelion[i]      = planets[i].perihelion;
    339     ascend_node[i]     = planets[i].ascend_node;
    340     inclination[i]     = planets[i].inclination;
    341     eccentricity[i]    = planets[i].eccentricity;
    342     semimajor_a[i]     = planets[i].semimajor_a;
    343     mjdMin[i]          = planets[i].mjdMin;
    344     mjdMax[i]          = planets[i].mjdMax;
    345     RatMin[i]          = planets[i].RatMin;
    346     RatMax[i]          = planets[i].RatMax;
    347     DatMin[i]          = planets[i].DatMin;
    348     DatMax[i]          = planets[i].DatMax;
    349 
    350     memcpy(&ID[i*8], planets[i].ID, 8);
    351   }
    352 
    353   // add the columns to the output array
    354   gfits_set_bintable_column (&theader, &ftable, "ID",           ID,           Nplanets);
    355   gfits_set_bintable_column (&theader, &ftable, "EPOCH",        epoch,        Nplanets);
    356   gfits_set_bintable_column (&theader, &ftable, "MEAN_ANOMALY", mean_anomaly, Nplanets);
    357   gfits_set_bintable_column (&theader, &ftable, "PERIHELION",   perihelion,   Nplanets);
    358   gfits_set_bintable_column (&theader, &ftable, "ASCEND_NODE",  ascend_node,  Nplanets);
    359   gfits_set_bintable_column (&theader, &ftable, "INCLINATION",  inclination,  Nplanets);
    360   gfits_set_bintable_column (&theader, &ftable, "ECCENTRICITY", eccentricity, Nplanets);
    361   gfits_set_bintable_column (&theader, &ftable, "SEMIMAJOR_A",  semimajor_a,  Nplanets);
    362   gfits_set_bintable_column (&theader, &ftable, "MJD_MIN",      mjdMin,       Nplanets);
    363   gfits_set_bintable_column (&theader, &ftable, "MJD_MAX",      mjdMax,       Nplanets);
    364   gfits_set_bintable_column (&theader, &ftable, "RAT_MIN",      RatMin,       Nplanets);
    365   gfits_set_bintable_column (&theader, &ftable, "RAT_MAX",      RatMax,       Nplanets);
    366   gfits_set_bintable_column (&theader, &ftable, "DAT_MIN",      DatMin,       Nplanets);
    367   gfits_set_bintable_column (&theader, &ftable, "DAT_MAX",      DatMax,       Nplanets);
    368 
    369   // free intermediate storage arrays
    370   free (ID);
    371   free (epoch);
    372   free (mean_anomaly);
    373   free (perihelion);
    374   free (ascend_node);
    375   free (inclination);
    376   free (eccentricity);
    377   free (semimajor_a);
    378   free (mjdMin);
    379   free (mjdMax);
    380   free (RatMin);
    381   free (RatMax);
    382   free (DatMin);
    383   free (DatMax);
    384 
    385   // open the output file
    386   FILE *f = fopen (filename, "w");
    387   if (!f) Shutdown ("ERROR: cannot open mpcorb file for output %s\n", filename);
    388 
    389   int status;
    390   status = gfits_fwrite_header  (f, &header);
    391   CHECK_STATUS (status, "ERROR: cannot write header for mpcorbs %s\n", filename);
    392 
    393   status = gfits_fwrite_matrix  (f, &matrix);
    394   CHECK_STATUS (status, "ERROR: cannot write matrix for mpcorbs %s\n", filename);
    395 
    396   status = gfits_fwrite_Theader (f, &theader);
    397   CHECK_STATUS (status, "ERROR: cannot write table header for mpcorbs %s\n", filename);
    398 
    399   status = gfits_fwrite_table  (f, &ftable);
    400   CHECK_STATUS (status, "ERROR: cannot write table data for mpcorbs %s\n", filename);
    401 
    402   gfits_free_header (&header);
    403   gfits_free_matrix (&matrix);
    404   gfits_free_header (&theader);
    405   gfits_free_table (&ftable);
    406 
    407   int fd = fileno (f);
    408 
    409   status = fflush (f);
    410   CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
    411 
    412   status = fsync (fd);
    413   CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
    414 
    415   status = fclose (f);
    416   CHECK_STATUS (!status, "ERROR: problem closing mpcorbs file %s\n", filename);
    417 
    418   return;
    419 }
    420 
    421 # define GET_COLUMN(OUT,NAME,TYPE) \
    422   TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
    423   myAssert (!strcmp(type, #TYPE), "wrong column type");
    424 
    425 // load planet data from a fits table
    426 Planets *mpcorb_read_fits (char *filename, int *nplanets) {
    427 
    428   int Ncol;
    429   off_t Nrow;
    430   Header header;
    431   Header theader;
    432   Matrix matrix;
    433   FTable ftable;
    434 
    435   FILE *f = fopen (filename, "r");
    436   if (!f) Shutdown ("ERROR: cannot open mpcorb fits file %s\n", filename);
    437 
    438   // load in PHU segment (ignore)
    439   if (!gfits_fread_header (f, &header)) Shutdown ("can't read mpcorb fits file header %s\n", filename);
    440   if (!gfits_fread_matrix (f, &matrix, &header)) Shutdown ("can't read mpcorb fits file matrix %s\n", filename);
    441   ftable.header = &theader;
    442 
    443   // load data for this header
    444   if (!gfits_load_header (f, &theader)) Shutdown ("can't read mpcorb fits file table header %s\n", filename);
    445   if (!gfits_fread_ftable_data (f, &ftable, FALSE)) Shutdown ("can't read mpcorb fits file table data %s\n", filename);
    446 
    447   // we read the entire block of data, then extract the columns, then set the Planets structure values.
    448   // I free the FITS table data after extracting the colums to avoid having 3 copies in memory.
    449 
    450   char type[16]; // used by the functions called by GET_COLUMN
    451 
    452   GET_COLUMN (ID,           "ID",            char);
    453   GET_COLUMN (epoch,        "EPOCH",         double);
    454   GET_COLUMN (mean_anomaly, "MEAN_ANOMALY",  double);
    455   GET_COLUMN (perihelion,   "PERIHELION",    double);
    456   GET_COLUMN (ascend_node,  "ASCEND_NODE",   double);
    457   GET_COLUMN (inclination,  "INCLINATION",   double);
    458   GET_COLUMN (eccentricity, "ECCENTRICITY",  double);
    459   GET_COLUMN (semimajor_a,  "SEMIMAJOR_A",   double);
    460   GET_COLUMN (mjdMin,       "MJD_MIN",       double);
    461   GET_COLUMN (mjdMax,       "MJD_MAX",       double);
    462   GET_COLUMN (RatMin,       "RAT_MIN",       double);
    463   GET_COLUMN (RatMax,       "RAT_MAX",       double);
    464   GET_COLUMN (DatMin,       "DAT_MIN",       double);
    465   GET_COLUMN (DatMax,       "DAT_MAX",       double);
    466 
    467   // free the memory associated with the FITS files
    468   gfits_free_header (&header);
    469   gfits_free_matrix (&matrix);
    470   gfits_free_header (&theader);
    471   gfits_free_table (&ftable);
    472 
    473   ALLOCATE_PTR (planets, Planets, Nrow);
    474   for (int i = 0; i < Nrow; i++) {
    475     planets[i].epoch        = epoch[i];
    476     planets[i].mean_anomaly = mean_anomaly[i];
    477     planets[i].perihelion   = perihelion[i];
    478     planets[i].ascend_node  = ascend_node[i];
    479     planets[i].inclination  = inclination[i];
    480     planets[i].eccentricity = eccentricity[i];
    481     planets[i].semimajor_a  = semimajor_a[i];
    482     planets[i].mjdMin       = mjdMin[i];
    483     planets[i].mjdMax       = mjdMax[i];
    484     planets[i].RatMin       = RatMin[i];
    485     planets[i].RatMax       = RatMax[i];
    486     planets[i].DatMin       = DatMin[i];
    487     planets[i].DatMax       = DatMax[i];
    488 
    489     memcpy(planets[i].ID, &ID[i*8], 8);
    490     planets[i].ID[7] = 0;
    491   }
    492   fprintf (stderr, "loaded data for %lld planets\n", (long long) Nrow);
    493 
    494   free (ID          );
    495   free (epoch       );
    496   free (mean_anomaly);
    497   free (perihelion  );
    498   free (ascend_node );
    499   free (inclination );
    500   free (eccentricity);
    501   free (semimajor_a );
    502   free (mjdMin      );
    503   free (mjdMax      );
    504   free (RatMin      );
    505   free (RatMax      );
    506   free (DatMin      );
    507   free (DatMax      );
    508 
    509   *nplanets = Nrow;
    510   return planets;
    511228}
    512229
    … …  
    603320  myPlanet.Robs = Rmean * DEG_RAD;
    604321  myPlanet.Dobs = Dmean * DEG_RAD;
     322  myPlanet.Mobs = NAN;
     323
     324  // calculate the apparent magnitude from the H magnitude:
     325  if (FullCalc) {
     326    // find the angle between the sun and the object
     327    double dRmean = Rmean - R_SUN;
     328    double Chi = sin(D_SUN)*sin(Dmean) + cos(D_SUN)*cos(Dmean)*cos(dRmean);
     329    double Psi1 = cos(Dmean)*sin(dRmean);
     330    double Psi2 = cos(D_SUN)*sin(Dmean) - sin(D_SUN)*cos(Dmean)*cos(dRmean);
     331    double Psi = sqrt(Psi1*Psi1 + Psi2*Psi2);
     332    double phi = atan2 (Psi, Chi);
     333
     334    // place phi in the range 0 360
     335    if ((phi) <    0.0) (phi) += M_PI;
     336    if ((phi) > 2*M_PI) (phi) -= M_PI;
     337
     338    double DsunTgt = sqrt (SQ(DIST_SUN) + SQ(dist) - 2.0*DIST_SUN*dist*cos(phi));
     339
     340    double phaseArg = (SQ(dist) + SQ(DsunTgt) - SQ(DIST_SUN)) / (2.0*DsunTgt*dist);
     341    double phase = DEG_RAD*acos(phaseArg);
     342    double phcen = phase / 100.0;
     343
     344    // mercury phase curve from Hilton 2005
     345    myPlanet.Mobs = planet->Hmag + 5*log10(dist) + 5*log10(DsunTgt) + 4.98*phcen - 4.88*phcen*phcen + 3.02*(phcen*phcen*phcen);
     346
     347    if (TESTING) fprintf (stderr, "R,D: %12.6f %12.6f : SunDist: %6.4f : phi %6.4f : phase %6.4f : Mapp %6.4f\n",
     348             myPlanet.Robs, myPlanet.Dobs, DsunTgt, phi*DEG_RAD, phase, myPlanet.Mobs);
     349  }
     350
    605351  return myPlanet;
    606352}
    607353
    608 # define MPC_DOUBLE(FIELD,START,END) {FIELD = strtod(&line[START-1], NULL);}
    609 
    610 // part the data in a single line of the MPCORB.DAT file
     354void solar_coords (double mjdObs, double *Rsun, double *Dsun, double *distSun) {
     355
     356# define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
     357# define J1900 15019.5       /* Modified Julian at 1900 */
     358
     359  /* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
     360
     361  double n = mjdObs - J2000;          // day number relative to standard epoch
     362  double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration)
     363  double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly
     364
     365  double lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees
     366  double beta   = 0.0;                                   // approx latitude
     367  double dSun   = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU
     368
     369  // year since 1900
     370  double year = (mjdObs - J1900) / 365.25;
     371
     372  double T = year / 100.0;
     373
     374  double phi = RAD_DEG*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
     375  double Xo = 0.0;
     376  double xo = 0.0;
     377
     378  // pre-calculated constants:
     379  double sin_phi_cos_Xo = sin(phi)*cos(Xo);
     380  double sin_phi_sin_Xo = sin(phi)*sin(Xo);
     381  double cos_phi        = cos(phi);
     382  double cos_phi_cos_Xo = cos(phi)*cos(Xo);
     383  double cos_phi_sin_Xo = cos(phi)*sin(Xo);
     384  double sin_phi        = sin(phi);
     385  double cos_Xo         = cos(Xo);
     386  double sin_Xo         = sin(Xo);
     387
     388  double X = lambda*RAD_DEG;
     389  double Y = beta  *RAD_DEG;
     390
     391  double cosY = cos(Y);
     392  double sinY = sin(Y);
     393 
     394  double cosX = cos(X);
     395  double sinX = sin(X);
     396
     397  // recast with constants extracted:
     398  double sin_y = cosY*sinX*sin_phi_cos_Xo - cosY*cosX*sin_phi_sin_Xo + sinY*cos_phi;
     399  double cos_y = sqrt (1 - sin_y*sin_y);
     400
     401  double sin_x = cosY*sinX*cos_phi_cos_Xo - cosY*cosX*cos_phi_sin_Xo - sinY*sin_phi;
     402  double cos_x = cosY*cosX*cos_Xo + cosY*sinX*sin_Xo;
     403     
     404  // atan2 returns -pi : +pi
     405  double RsunRad = atan2 (sin_x, cos_x) + xo;
     406  if ((RsunRad) <    0.0) (RsunRad) += M_PI;
     407  if ((RsunRad) > 2*M_PI) (RsunRad) -= M_PI;
     408
     409  // should be in range -pi/2 : +pi/2
     410  double DsunRad = atan2 (sin_y, cos_y);
     411
     412  double epochYear = 2000.0 + (mjdObs - J2000) / 365.25; // J2000 = 51544.5
     413
     414  slaPreces ("FK5", epochYear, 2000.000000, &RsunRad, &DsunRad);
     415
     416  // *Rsun = RsunRad * DEG_RAD;
     417  // *Dsun = DsunRad * DEG_RAD;
     418
     419  // keep these in radians for calculations in mpcorb_predict
     420  *Rsun = RsunRad;
     421  *Dsun = DsunRad;
     422  *distSun = dSun;
     423}
     424
     425// STATUS is value expected for success
     426# define CHECK_STATUS(STATUS,MSG,...) if (!(STATUS)) { Shutdown (MSG, __VA_ARGS__); }
     427
     428// write minor planet positions (Rref,Dref,Mref) to a FITS table
     429void mpcorb_refs_fits (char *filename, PlanetDatum *planets, int Nplanets) {
     430
     431  Header header;
     432  Header theader;
     433  Matrix matrix;
     434  FTable ftable;
     435
     436  gfits_init_header (&header);
     437  header.extend = TRUE;
     438  gfits_create_header (&header);
     439  gfits_create_matrix (&header, &matrix);
     440
     441  gfits_create_table_header (&theader, "BINTABLE", "DATA");
     442
     443  gfits_define_bintable_column (&theader, "8A", "ID",   "name",    "none",        1.0, 0.0);
     444  gfits_define_bintable_column (&theader,  "D", "Rref", "RA",      "degrees",     1.0, 0.0);
     445  gfits_define_bintable_column (&theader,  "D", "Dref", "DEC",     "degrees",     1.0, 0.0);
     446  gfits_define_bintable_column (&theader,  "D", "Mref", "Mag",     "magnitudes",  1.0, 0.0);
     447  gfits_define_bintable_column (&theader,  "D", "dist", "DIST",    "AU",          1.0, 0.0);
     448
     449  // generate the output array that carries the data
     450  gfits_create_table (&theader, &ftable);
     451
     452  // create intermediate storage arrays
     453  ALLOCATE_PTR (ID,                     char, Nplanets*8);
     454  ALLOCATE_PTR (Rref,                 double, Nplanets);
     455  ALLOCATE_PTR (Dref,                 double, Nplanets);
     456  ALLOCATE_PTR (Mref,                 double, Nplanets);
     457  ALLOCATE_PTR (dist,                 double, Nplanets);
     458
     459  // set intermediate storage arrays
     460  for (int i = 0; i < Nplanets; i++) {
     461    Rref[i]        = planets[i].Robs;
     462    Dref[i]        = planets[i].Dobs;
     463    dist[i]        = planets[i].dist;
     464    Mref[i]        = planets[i].Mobs;
     465    memcpy(&ID[i*8], planets[i].ID, 8);
     466  }
     467
     468  // add the columns to the output array
     469  gfits_set_bintable_column (&theader, &ftable, "ID",   ID,   Nplanets);
     470  gfits_set_bintable_column (&theader, &ftable, "Rref", Rref, Nplanets);
     471  gfits_set_bintable_column (&theader, &ftable, "Dref", Dref, Nplanets);
     472  gfits_set_bintable_column (&theader, &ftable, "Mref", Mref, Nplanets);
     473  gfits_set_bintable_column (&theader, &ftable, "dist", dist, Nplanets);
     474
     475  // free intermediate storage arrays
     476  free (ID);
     477  free (Rref);
     478  free (Dref);
     479  free (Mref);
     480  free (dist);
     481
     482  // open the output file
     483  FILE *f = fopen (filename, "w");
     484  if (!f) Shutdown ("ERROR: cannot open mpcorb file for output %s\n", filename);
     485
     486  int status;
     487  status = gfits_fwrite_header  (f, &header);
     488  CHECK_STATUS (status, "ERROR: cannot write header for mpcorbs %s\n", filename);
     489
     490  status = gfits_fwrite_matrix  (f, &matrix);
     491  CHECK_STATUS (status, "ERROR: cannot write matrix for mpcorbs %s\n", filename);
     492
     493  status = gfits_fwrite_Theader (f, &theader);
     494  CHECK_STATUS (status, "ERROR: cannot write table header for mpcorbs %s\n", filename);
     495
     496  status = gfits_fwrite_table  (f, &ftable);
     497  CHECK_STATUS (status, "ERROR: cannot write table data for mpcorbs %s\n", filename);
     498
     499  gfits_free_header (&header);
     500  gfits_free_matrix (&matrix);
     501  gfits_free_header (&theader);
     502  gfits_free_table (&ftable);
     503
     504  int fd = fileno (f);
     505
     506  status = fflush (f);
     507  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
     508
     509  status = fsync (fd);
     510  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
     511
     512  status = fclose (f);
     513  CHECK_STATUS (!status, "ERROR: problem closing mpcorbs file %s\n", filename);
     514
     515  return;
     516}
     517
     518// write minor planet orbital elements (plus (R,D,T)[min,max]) to a FITS table
     519void mpcorb_save_fits (char *filename, Planets *planets, int Nplanets) {
     520
     521  Header header;
     522  Header theader;
     523  Matrix matrix;
     524  FTable ftable;
     525
     526  gfits_init_header (&header);
     527  header.extend = TRUE;
     528  gfits_create_header (&header);
     529  gfits_create_matrix (&header, &matrix);
     530
     531  gfits_create_table_header (&theader, "BINTABLE", "MPCORB");
     532
     533  gfits_define_bintable_column (&theader, "8A", "ID",           "name",         "none",     1.0, 0.0);
     534  gfits_define_bintable_column (&theader,  "D", "EPOCH",        "epoch",        "MJD",      1.0, 0.0);
     535  gfits_define_bintable_column (&theader,  "D", "MEAN_ANOMALY", "mean_anomaly", "degrees",  1.0, 0.0);
     536  gfits_define_bintable_column (&theader,  "D", "PERIHELION",   "perihelion",   "degrees",  1.0, 0.0);
     537  gfits_define_bintable_column (&theader,  "D", "ASCEND_NODE",  "ascend_node",  "degrees",  1.0, 0.0);
     538  gfits_define_bintable_column (&theader,  "D", "INCLINATION",  "inclination",  "degrees",  1.0, 0.0);
     539  gfits_define_bintable_column (&theader,  "D", "ECCENTRICITY", "eccentricity", "unitless", 1.0, 0.0);
     540  gfits_define_bintable_column (&theader,  "D", "SEMIMAJOR_A",  "semimajor_a",  "AU",       1.0, 0.0);
     541  gfits_define_bintable_column (&theader,  "D", "MJD_MIN",      "mjdMin",       "MJD",      1.0, 0.0);
     542  gfits_define_bintable_column (&theader,  "D", "MJD_MAX",      "mjdMax",       "MJD",      1.0, 0.0);
     543  gfits_define_bintable_column (&theader,  "D", "RAT_MIN",      "RatMin",       "degrees",  1.0, 0.0);
     544  gfits_define_bintable_column (&theader,  "D", "RAT_MAX",      "RatMax",       "degrees",  1.0, 0.0);
     545  gfits_define_bintable_column (&theader,  "D", "DAT_MIN",      "DatMin",       "degrees",  1.0, 0.0);
     546  gfits_define_bintable_column (&theader,  "D", "DAT_MAX",      "DatMax",       "degrees",  1.0, 0.0);
     547  gfits_define_bintable_column (&theader,  "D", "HMAG",         "Hmag",         "mag",      1.0, 0.0);
     548
     549  // generate the output array that carries the data
     550  gfits_create_table (&theader, &ftable);
     551
     552  // create intermediate storage arrays
     553  ALLOCATE_PTR (ID,                     char,   Nplanets*8);
     554  ALLOCATE_PTR (epoch,                  double, Nplanets);
     555  ALLOCATE_PTR (mean_anomaly,           double, Nplanets);
     556  ALLOCATE_PTR (perihelion,             double, Nplanets);
     557  ALLOCATE_PTR (ascend_node,            double, Nplanets);
     558  ALLOCATE_PTR (inclination,            double, Nplanets);
     559  ALLOCATE_PTR (eccentricity,           double, Nplanets);
     560  ALLOCATE_PTR (semimajor_a,            double, Nplanets);
     561  ALLOCATE_PTR (mjdMin,                 double, Nplanets);
     562  ALLOCATE_PTR (mjdMax,                 double, Nplanets);
     563  ALLOCATE_PTR (RatMin,                 double, Nplanets);
     564  ALLOCATE_PTR (RatMax,                 double, Nplanets);
     565  ALLOCATE_PTR (DatMin,                 double, Nplanets);
     566  ALLOCATE_PTR (DatMax,                 double, Nplanets);
     567  ALLOCATE_PTR (Hmag,                   double, Nplanets);
     568
     569  // set intermediate storage arrays
     570  for (int i = 0; i < Nplanets; i++) {
     571    epoch[i]           = planets[i].epoch;
     572    mean_anomaly[i]    = planets[i].mean_anomaly;
     573    perihelion[i]      = planets[i].perihelion;
     574    ascend_node[i]     = planets[i].ascend_node;
     575    inclination[i]     = planets[i].inclination;
     576    eccentricity[i]    = planets[i].eccentricity;
     577    semimajor_a[i]     = planets[i].semimajor_a;
     578    mjdMin[i]          = planets[i].mjdMin;
     579    mjdMax[i]          = planets[i].mjdMax;
     580    RatMin[i]          = planets[i].RatMin;
     581    RatMax[i]          = planets[i].RatMax;
     582    DatMin[i]          = planets[i].DatMin;
     583    DatMax[i]          = planets[i].DatMax;
     584    Hmag[i]            = planets[i].Hmag;
     585
     586    memcpy(&ID[i*8], planets[i].ID, 8);
     587  }
     588
     589  // add the columns to the output array
     590  gfits_set_bintable_column (&theader, &ftable, "ID",           ID,           Nplanets);
     591  gfits_set_bintable_column (&theader, &ftable, "EPOCH",        epoch,        Nplanets);
     592  gfits_set_bintable_column (&theader, &ftable, "MEAN_ANOMALY", mean_anomaly, Nplanets);
     593  gfits_set_bintable_column (&theader, &ftable, "PERIHELION",   perihelion,   Nplanets);
     594  gfits_set_bintable_column (&theader, &ftable, "ASCEND_NODE",  ascend_node,  Nplanets);
     595  gfits_set_bintable_column (&theader, &ftable, "INCLINATION",  inclination,  Nplanets);
     596  gfits_set_bintable_column (&theader, &ftable, "ECCENTRICITY", eccentricity, Nplanets);
     597  gfits_set_bintable_column (&theader, &ftable, "SEMIMAJOR_A",  semimajor_a,  Nplanets);
     598  gfits_set_bintable_column (&theader, &ftable, "MJD_MIN",      mjdMin,       Nplanets);
     599  gfits_set_bintable_column (&theader, &ftable, "MJD_MAX",      mjdMax,       Nplanets);
     600  gfits_set_bintable_column (&theader, &ftable, "RAT_MIN",      RatMin,       Nplanets);
     601  gfits_set_bintable_column (&theader, &ftable, "RAT_MAX",      RatMax,       Nplanets);
     602  gfits_set_bintable_column (&theader, &ftable, "DAT_MIN",      DatMin,       Nplanets);
     603  gfits_set_bintable_column (&theader, &ftable, "DAT_MAX",      DatMax,       Nplanets);
     604  gfits_set_bintable_column (&theader, &ftable, "HMAG",         Hmag,         Nplanets);
     605
     606  // free intermediate storage arrays
     607  free (ID);
     608  free (epoch);
     609  free (mean_anomaly);
     610  free (perihelion);
     611  free (ascend_node);
     612  free (inclination);
     613  free (eccentricity);
     614  free (semimajor_a);
     615  free (mjdMin);
     616  free (mjdMax);
     617  free (RatMin);
     618  free (RatMax);
     619  free (DatMin);
     620  free (DatMax);
     621  free (Hmag);
     622
     623  // open the output file
     624  FILE *f = fopen (filename, "w");
     625  if (!f) Shutdown ("ERROR: cannot open mpcorb file for output %s\n", filename);
     626
     627  int status;
     628  status = gfits_fwrite_header  (f, &header);
     629  CHECK_STATUS (status, "ERROR: cannot write header for mpcorbs %s\n", filename);
     630
     631  status = gfits_fwrite_matrix  (f, &matrix);
     632  CHECK_STATUS (status, "ERROR: cannot write matrix for mpcorbs %s\n", filename);
     633
     634  status = gfits_fwrite_Theader (f, &theader);
     635  CHECK_STATUS (status, "ERROR: cannot write table header for mpcorbs %s\n", filename);
     636
     637  status = gfits_fwrite_table  (f, &ftable);
     638  CHECK_STATUS (status, "ERROR: cannot write table data for mpcorbs %s\n", filename);
     639
     640  gfits_free_header (&header);
     641  gfits_free_matrix (&matrix);
     642  gfits_free_header (&theader);
     643  gfits_free_table (&ftable);
     644
     645  int fd = fileno (f);
     646
     647  status = fflush (f);
     648  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
     649
     650  status = fsync (fd);
     651  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
     652
     653  status = fclose (f);
     654  CHECK_STATUS (!status, "ERROR: problem closing mpcorbs file %s\n", filename);
     655
     656  return;
     657}
     658
     659# define GET_COLUMN(OUT,NAME,TYPE) \
     660  TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
     661  myAssert (!strcmp(type, #TYPE), "wrong column type");
     662
     663// load planet data from a fits table
     664Planets *mpcorb_read_fits (char *filename, int *nplanets) {
     665
     666  int Ncol;
     667  off_t Nrow;
     668  Header header;
     669  Header theader;
     670  Matrix matrix;
     671  FTable ftable;
     672
     673  FILE *f = fopen (filename, "r");
     674  if (!f) Shutdown ("ERROR: cannot open mpcorb fits file %s\n", filename);
     675
     676  // load in PHU segment (ignore)
     677  if (!gfits_fread_header (f, &header)) Shutdown ("can't read mpcorb fits file header %s\n", filename);
     678  if (!gfits_fread_matrix (f, &matrix, &header)) Shutdown ("can't read mpcorb fits file matrix %s\n", filename);
     679  ftable.header = &theader;
     680
     681  // load data for this header
     682  if (!gfits_load_header (f, &theader)) Shutdown ("can't read mpcorb fits file table header %s\n", filename);
     683  if (!gfits_fread_ftable_data (f, &ftable, FALSE)) Shutdown ("can't read mpcorb fits file table data %s\n", filename);
     684
     685  // we read the entire block of data, then extract the columns, then set the Planets structure values.
     686  // I free the FITS table data after extracting the colums to avoid having 3 copies in memory.
     687
     688  char type[16]; // used by the functions called by GET_COLUMN
     689
     690  GET_COLUMN (ID,           "ID",            char);
     691  GET_COLUMN (epoch,        "EPOCH",         double);
     692  GET_COLUMN (mean_anomaly, "MEAN_ANOMALY",  double);
     693  GET_COLUMN (perihelion,   "PERIHELION",    double);
     694  GET_COLUMN (ascend_node,  "ASCEND_NODE",   double);
     695  GET_COLUMN (inclination,  "INCLINATION",   double);
     696  GET_COLUMN (eccentricity, "ECCENTRICITY",  double);
     697  GET_COLUMN (semimajor_a,  "SEMIMAJOR_A",   double);
     698  GET_COLUMN (mjdMin,       "MJD_MIN",       double);
     699  GET_COLUMN (mjdMax,       "MJD_MAX",       double);
     700  GET_COLUMN (RatMin,       "RAT_MIN",       double);
     701  GET_COLUMN (RatMax,       "RAT_MAX",       double);
     702  GET_COLUMN (DatMin,       "DAT_MIN",       double);
     703  GET_COLUMN (DatMax,       "DAT_MAX",       double);
     704  GET_COLUMN (Hmag,         "HMAG",          double);
     705
     706  // free the memory associated with the FITS files
     707  gfits_free_header (&header);
     708  gfits_free_matrix (&matrix);
     709  gfits_free_header (&theader);
     710  gfits_free_table (&ftable);
     711
     712  ALLOCATE_PTR (planets, Planets, Nrow);
     713  for (int i = 0; i < Nrow; i++) {
     714    planets[i].epoch        = epoch[i];
     715    planets[i].mean_anomaly = mean_anomaly[i];
     716    planets[i].perihelion   = perihelion[i];
     717    planets[i].ascend_node  = ascend_node[i];
     718    planets[i].inclination  = inclination[i];
     719    planets[i].eccentricity = eccentricity[i];
     720    planets[i].semimajor_a  = semimajor_a[i];
     721    planets[i].mjdMin       = mjdMin[i];
     722    planets[i].mjdMax       = mjdMax[i];
     723    planets[i].RatMin       = RatMin[i];
     724    planets[i].RatMax       = RatMax[i];
     725    planets[i].DatMin       = DatMin[i];
     726    planets[i].DatMax       = DatMax[i];
     727    planets[i].Hmag         = Hmag[i];
     728
     729    memcpy(planets[i].ID, &ID[i*8], 8);
     730    planets[i].ID[7] = 0;
     731  }
     732  fprintf (stderr, "loaded data for %lld planets\n", (long long) Nrow);
     733
     734  free (ID          );
     735  free (epoch       );
     736  free (mean_anomaly);
     737  free (perihelion  );
     738  free (ascend_node );
     739  free (inclination );
     740  free (eccentricity);
     741  free (semimajor_a );
     742  free (mjdMin      );
     743  free (mjdMax      );
     744  free (RatMin      );
     745  free (RatMax      );
     746  free (DatMin      );
     747  free (DatMax      );
     748  free (Hmag        );
     749
     750  *nplanets = Nrow;
     751  return planets;
     752}
     753
     754# define MPC_GETDBL(FIELD,START,END) {FIELD = strtod(&line[START-1], &endptr); if (endptr == &line[START-1]) FIELD = NAN;}
     755# define MPC_GETINT(FIELD,START,END) {FIELD = strtol(&line[START-1], NULL, 10); }
     756
     757// parse the data in a single line of the MPCORB.DAT file
    611758int mpcorb_parseline (char *line, Planets *planet, int Nmax) {
    612759
    … …  
    626773  strncpy (planet->ID, line, 7); planet->ID[7] = 0;
    627774
    628   MPC_DOUBLE (planet->mean_anomaly, 27, 35);
    629   MPC_DOUBLE (planet->perihelion,   38, 46);
    630   MPC_DOUBLE (planet->ascend_node,  49, 57);
    631   MPC_DOUBLE (planet->inclination,  60, 68);
    632   MPC_DOUBLE (planet->eccentricity, 71, 79);
    633   MPC_DOUBLE (planet->semimajor_a,  93, 103);
     775  char *endptr;
     776  MPC_GETDBL (planet->Hmag,           9,  13);
     777  MPC_GETDBL (planet->mean_anomaly,  27,  35);
     778  MPC_GETDBL (planet->perihelion,    38,  46);
     779  MPC_GETDBL (planet->ascend_node,   49,  57);
     780  MPC_GETDBL (planet->inclination,   60,  68);
     781  MPC_GETDBL (planet->eccentricity,  71,  79);
     782  MPC_GETDBL (planet->semimajor_a,   93, 103);
     783  MPC_GETINT (planet->Nobs,         118, 122);
     784  MPC_GETINT (planet->Nopp,         124, 126);
    634785
    635786  // the epoch is stored in a compressed format:
Note: See TracChangeset for help on using the changeset viewer.