Changeset 41710

Timestamp:

Jul 10, 2021, 2:40:08 PM (5 years ago)

Author:

eugene

Message:

add fits output to mpcorb_predict (moment)

Location:

trunk/Ohana/src/tools

Files:

• Makefile (modified) (1 diff)
• src/mpcorb_predict.c (modified) (4 diffs)

trunk/Ohana/src/tools/Makefile

@@ -76 +76 @@
 $(DESTLIB)/libsla.a:
         echo "missing SLALIB libsla.a: install in $(DESTLIB)"
+        @exit 2
 
 $(DESTINC)/slalib.h:
         echo "missing SLALIB slalib.h: install in $(DESTINC)"
+        @exit 2
 
-mpcorb_predict: $(DESTBIN)/mpcorb_predict $(DESTLIB)/libsla.a $(DESTINC)/slalib.h
+mpcorb_predict: $(DESTLIB)/libsla.a $(DESTINC)/slalib.h $(DESTBIN)/mpcorb_predict 
+
+mpcorb: mpcorb_predict
 
 %.clean:

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

@@ -38 +38 @@
 Planets *mpcorb_read_fits (char *filename, int *nplanets);
 void     mpcorb_save_fits (char *filename, Planets *planets, int Nplanets);
+void     mpcorb_refs_fits (char *filename, Planets *planets, int Nplanets);
 
 void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc);
@@ -151 +152 @@
   char  *output   = argv[8];
 
-  FILE *fout = fopen (output, "w");
-  if (!fout) Shutdown ("ERROR: unable to open file for output %s\n", output);
+  // FILE *fout = fopen (output, "w");
+  // if (!fout) Shutdown ("ERROR: unable to open file for output %s\n", output);
 
   int Nplanets = 0;
   Planets *planets = mpcorb_read_fits (filename, &Nplanets);
   if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
+
+  int NplanetsSave = 0;
+  int NPLANETSSAVE = 1000;
+  ALLOCATE_PTR (planetsSave, Planets, NPLANETSSAVE);
 
   // transform all objects and identify those in the target region:
@@ -173 +178 @@
     // if it is in the region, use the more accurate calculation
     mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE);
-    fprintf (fout, "%8s %12.6f %12.6f\n", planets[i].ID, Robs, Dobs);
+
+    // fprintf (fout, "%8s %12.6f %12.6f\n", planets[i].ID, Robs, Dobs);
+
+    // save this prediction, with R,D in both @min & @max
+    planetsSave[NplanetsSave] = planets[i];
+    planetsSave[NplanetsSave].mjdMin = mjdObs;
+    planetsSave[NplanetsSave].RatMin = Robs;
+    planetsSave[NplanetsSave].DatMin = Dobs;
+    planetsSave[NplanetsSave].mjdMax = mjdObs;
+    planetsSave[NplanetsSave].RatMax = Robs;
+    planetsSave[NplanetsSave].DatMax = Dobs;
+    
+    NplanetsSave++;
+    CHECK_REALLOCATE (planetsSave, Planets, NPLANETSSAVE, NplanetsSave, 1000);
   }
-  fclose (fout);
   fprintf (stderr, "\n");
+
+  mpcorb_refs_fits (output, planetsSave, NplanetsSave);
 
   exit (0);
@@ -183 +202 @@
 // STATUS is value expected for success
 # define CHECK_STATUS(STATUS,MSG,...) if (!(STATUS)) { Shutdown (MSG, __VA_ARGS__); }
+
+// write minor planet positions (Rref,Dref,Mref) to a FITS table
+void mpcorb_refs_fits (char *filename, Planets *planets, int Nplanets) {
+
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  gfits_init_header (&header);
+  header.extend = TRUE;
+  gfits_create_header (&header);
+  gfits_create_matrix (&header, &matrix);
+
+  gfits_create_table_header (&theader, "BINTABLE", "DATA");
+
+  gfits_define_bintable_column (&theader, "8A", "ID",   "name", "none",        1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "Rref", "RA",   "degrees",     1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "Dref", "DEC",  "degrees",     1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "Mref", "Mag",  "magnitudes",  1.0, 0.0);
+
+  // generate the output array that carries the data
+  gfits_create_table (&theader, &ftable);
+
+  // create intermediate storage arrays
+  ALLOCATE_PTR (ID,                     char, Nplanets*8);
+  ALLOCATE_PTR (Rref,                 double, Nplanets);
+  ALLOCATE_PTR (Dref,                 double, Nplanets);
+  ALLOCATE_PTR (Mref,                 double, Nplanets);
+
+  // set intermediate storage arrays
+  for (int i = 0; i < Nplanets; i++) {
+    Rref[i]        = planets[i].RatMin;
+    Dref[i]        = planets[i].DatMin;
+    Mref[i]        = 16.0; // need to add this to prediction
+    memcpy(&ID[i*8], planets[i].ID, 8);
+  }
+
+  // add the columns to the output array
+  gfits_set_bintable_column (&theader, &ftable, "ID",   ID,   Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "Rref", Rref, Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "Dref", Dref, Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "Mref", Mref, Nplanets);
+
+  // free intermediate storage arrays
+  free (ID);
+  free (Rref);
+  free (Dref);
+  free (Mref);
+
+  // open the output file
+  FILE *f = fopen (filename, "w");
+  if (!f) Shutdown ("ERROR: cannot open mpcorb file for output %s\n", filename);
+
+  int status;
+  status = gfits_fwrite_header  (f, &header);
+  CHECK_STATUS (status, "ERROR: cannot write header for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_matrix  (f, &matrix);
+  CHECK_STATUS (status, "ERROR: cannot write matrix for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_Theader (f, &theader);
+  CHECK_STATUS (status, "ERROR: cannot write table header for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_table  (f, &ftable);
+  CHECK_STATUS (status, "ERROR: cannot write table data for mpcorbs %s\n", filename);
+
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  gfits_free_header (&theader);
+  gfits_free_table (&ftable);
+
+  int fd = fileno (f);
+
+  status = fflush (f);
+  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
+
+  status = fsync (fd);
+  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
+
+  status = fclose (f);
+  CHECK_STATUS (!status, "ERROR: problem closing mpcorbs file %s\n", filename);
+
+  return;
+}
 
 // write minor planet orbital elements (plus (R,D,T)[min,max]) to a FITS table