JSOC_Documentation: proj/rings/apps/imginfo.c Source File

00001         /*  The following code is adapted from SOI functions/keywords.c  */
00002 #define PLATFORM_UNKNOWN    (0)
00003 #define PLATFORM_UNRECOGNIZED   (1)
00004 #define PLATFORM_SOHO       (2)
00005 #define PLATFORM_GONGPLUS   (3)
00006 #define PLATFORM_MWO60      (4)
00007 #define PLATFORM_BBSO       (5)
00008 #define PLATFORM_TRACE      (6)
00009 #define PLATFORM_SPOLE_JSO  (10)
00010 #define PLATFORM_GONG       (30)
00011 #define PLATFORM_OBSPM      (40)
00012 
00013 #define INSTRUMENT_UNKNOWN  (0)
00014 #define INSTRUMENT_UNRECOGNIZED (1)
00015 #define INSTRUMENT_SOHO_MDI (10)
00016 #define INSTRUMENT_SOHO_EIT (11)
00017 #define INSTRUMENT_GONG_TD  (20)
00018 #define INSTRUMENT_GONG_CT  (21)
00019 #define INSTRUMENT_GONG_TC  (22)
00020 #define INSTRUMENT_GONG_BB  (23)
00021 #define INSTRUMENT_GONG_ML  (24)
00022 #define INSTRUMENT_GONG_LE  (25)
00023 #define INSTRUMENT_GONG_UD  (26)
00024 #define INSTRUMENT_GONG_MERGE   (29)
00025 #define INSTRUMENT_MWO60_MOF    (30)
00026 #define INSTRUMENT_BBSO_SINGER  (40)
00027 #define INSTRUMENT_TRACE    (50)
00028 #define INSTRUMENT_MOTH     (60)
00029 #define INSTRUMENT_OBSPM_SPHG   (70)
00030 
00031 #define NO_DATA_DICT    (0x0001)
00032 #define NO_SEMIDIAMETER (0x0002)
00033 #define NO_XSCALE   (0x0004)
00034 #define NO_YSCALE   (0x0008)
00035 #define NO_XCENTERLOC   (0x0010)
00036 #define NO_YCENTERLOC   (0x0020)
00037 #define NO_HELIO_LATC   (0x0040)
00038 #define NO_HELIO_LONC   (0x0080)
00039 #define NO_HELIO_PA (0x0100)
00040 #define NO_XUNITS   (0x0200)
00041 #define NO_YUNITS   (0x0400)
00042 #define NO_OBSERVER_LAT (0x0002)
00043 #define NO_OBSERVER_LON (0x0004)
00044 
00045 #define KEYSCOPE_VARIABLE   (0x80000000)
00046 #define LOCALHS_IMGINFO_VERSION ("1.0")
00047 
00048 /*
00049  *  process image info (attitude, plate scale, distortion)
00050  *  stub function based on SOI solar_image_params()
00051  */
00052 typedef struct paramdef {
00053   double scale;
00054   double offset;
00055   double defval;
00056   unsigned int statusbit;
00057   char name[32];
00058 } ParamDef;
00059 
00060 static double lookup (DRMS_Record_t *rec, ParamDef key, int *status) {
00061   double value = key.defval;
00062   int lookupstat = 0;
00063 
00064   value = drms_getkey_double (rec, key.name, &lookupstat);
00065   value = value * key.scale + key.offset;
00066   if (lookupstat) *status |= key.statusbit;
00067   if (isnan (value)) *status |= key.statusbit;
00068   return value;
00069 }
00070 
00071 static char *lookup_str (DRMS_Record_t *rec, ParamDef key, int *status) {
00072   DRMS_Keyword_t *keywd;
00073   int lstat;
00074   char *value;
00075 
00076   value = drms_getkey_string (rec, key.name, &lstat);
00077   if (lstat) *status |= key.statusbit;
00078                          /*  cadence should be constant  */
00079   if ((keywd = drms_keyword_lookup (rec, key.name, 1))) {
00080     if (keywd->info->recscope != 1) *status |= KEYSCOPE_VARIABLE;
00081   } else *status |= key.statusbit;
00082   return value;
00083 }
00084 
00085 static int solar_image_info (DRMS_Record_t *img, double *xscl, double *yscl,
00086     double *ctrx, double *ctry, double *apsd, const char *rsun_key,
00087     const char *apsd_key, double *pang, double *ellipse_e, double *ellipse_pa,
00088     int *x_invrt, int *y_invrt, int *need_ephem, int AIPS_convention) {
00089 /*
00090  *  Provides the following values from the DRMS record:
00091  *    xscl  scale in the image column direction (arc-sec/pixel)
00092  *    yscl  scale in the image row direction (arc-sec/pixel)
00093  *    ctrx  (virtual) fractional pixel column of the center of the
00094  *  solar disc
00095  *    ctry  (virtual) fractional pixel row of the center of the
00096  *  solar disc
00097  *    apsd  apparent semi-diameter (semimajor-axis) of the solar disc, in
00098  *  pixel units
00099  *    pang  position angle of solar north relative to image vertical (y-axis,
00100  *  [0,0] -> [0,1]), measured westward (clockwise), in radians
00101  *    eecc  eccentricity of best-fit ellipse describing limb
00102  *    eang  position angle of best-fit ellipse describing limb, relative
00103  *  to direction of solar north, measured westward (clockwise), in radians
00104  *    xinv  0 if image is direct, 1 if flipped by columns
00105  *    yinv  0 if image is direct, 1 if flipped by rows
00106  *  If AIPS_convention is true (!=0), it is assumed that the input keywords
00107  *    representing position and ellipse angles are measured westward
00108  *    (clockwise) relative to their nominal axes; otherwise they are measured
00109  *    eastward (counter-clockwise) relative to their nominal axes.
00110  *  NO!
00111  *  The following data types are supported: SOHO-MDI, GONG+, Mt. Wilson MOF,
00112  *    SOHO-EIT, TRACE, BBSO Ha
00113  *  NO!
00114  *  If the data are not recognizably of one of these types, the function
00115  *    returns NO_DATA_DICT; if one or more required keywords are missing
00116  *    the function returns a status mask indicating which values could not
00117  *    be filled.  No flag is set if the image ellipse parameters are not
00118  *    available (unless they are required for other parameters), but the
00119  *    eccentricity and pericenter angle are quietly set to 0.
00120  */
00121   enum param {
00122     XSCL, YSCL, XUNI, YUNI, LATC, LONC, CTRX, CTRY, PANG, APSD, RSUN,
00123     ESMA, ESMI, EECC, EANG, XASP, YASP, PCT
00124   };
00125   enum known_plat {
00126     UNKNOWN
00127   };
00128   static ParamDef param[PCT];
00129   static double raddeg =  M_PI / 180.0;
00130   static double degrad = 180.0 / M_PI;
00131   double ella, ellb;
00132   int n, status = 0;
00133   static int scale_avail, xinv_type, yinv_type;
00134   static int hdrtype = UNKNOWN, lasthdr = UNKNOWN - 1;
00135   char *strval;
00136 /*
00137  *  Set up the appropriate dictionary for interpretation of keywords
00138  */
00139   if (lasthdr != hdrtype) {
00140     if (lasthdr >= UNKNOWN)
00141       fprintf (stderr,
00142       "Warning from solar_image_info(): record keywords may change\n");
00143     for (n = 0; n < PCT; n++) {
00144       sprintf (param[n].name, "No Keyword");
00145       param[n].scale = 1.0;
00146       param[n].offset = 0.0;
00147       param[n].defval = NAN;
00148       param[n].statusbit = 0;
00149     }
00150     param[RSUN].statusbit = NO_SEMIDIAMETER;
00151     param[APSD].statusbit = NO_SEMIDIAMETER;
00152     param[XSCL].statusbit = NO_XSCALE;
00153     param[YSCL].statusbit = NO_YSCALE;
00154     param[XUNI].statusbit = NO_XUNITS;
00155     param[YUNI].statusbit = NO_YUNITS;
00156     param[CTRX].statusbit = NO_XCENTERLOC;
00157     param[CTRY].statusbit = NO_YCENTERLOC;
00158     param[CTRX].statusbit = NO_XCENTERLOC;
00159     param[CTRY].statusbit = NO_YCENTERLOC;
00160     param[LATC].statusbit = NO_HELIO_LATC;
00161     param[LONC].statusbit = NO_HELIO_LONC;
00162     param[PANG].statusbit = NO_HELIO_PA;
00163     param[EECC].defval = 0.0;
00164     param[EANG].defval = 0.0;
00165     param[XASP].defval = 1.0;
00166     param[YASP].defval = 1.0;
00167 
00168     switch (hdrtype) {
00169       default:                        /*  Assume WCS HPLN/T  */
00170                      /*  WITH CERTAIN MDI SPECIFIC ONES  */
00171     scale_avail = 1;
00172     xinv_type = yinv_type = 0;
00173     sprintf (param[XUNI].name, "CUNIT1");
00174     sprintf (param[YUNI].name, "CUNIT2");
00175     sprintf (param[XSCL].name, "CDELT1");
00176     sprintf (param[YSCL].name, "CDELT2");
00177     sprintf (param[CTRX].name, "CRPIX1");
00178     param[CTRX].offset = -1.0;
00179     sprintf (param[CTRY].name, "CRPIX2");
00180     param[CTRY].offset = -1.0;
00181     *need_ephem = 0;
00182     strval = lookup_str (img, param[XUNI], &status);
00183     if (!(status & NO_XUNITS)) {
00184       if (!strcmp (strval, "arcsec")) param[XSCL].scale = 1.0;
00185       else if (!strcmp (strval, "arcmin")) param[XSCL].scale = 1.0 / 60.0;
00186       else if (!strcmp (strval, "deg")) param[XSCL].scale = 1.0 / 3600.0;
00187       else if (!strcmp (strval, "mas")) param[XSCL].scale = 1000.0;
00188       else if (!strcmp (strval, "rad")) param[XSCL].scale = degrad * 3600.0;
00189 /*
00190       need_units = status & KEYSCOPE_VARIABLE;
00191 */
00192     }
00193     if (strval) free (strval);
00194     strval = lookup_str (img, param[YUNI], &status);
00195     if (!(status & NO_YUNITS)) {
00196       if (!strcmp (strval, "arcsec")) param[YSCL].scale = 1.0;
00197       else if (!strcmp (strval, "arcmin")) param[YSCL].scale = 1.0 / 60.0;
00198       else if (!strcmp (strval, "deg")) param[YSCL].scale = 1.0 / 3600.0;
00199       else if (!strcmp (strval, "mas")) param[YSCL].scale = 1000.0;
00200       else if (!strcmp (strval, "rad")) param[YSCL].scale = degrad * 3600.0;
00201 /*
00202       need_units = status & KEYSCOPE_VARIABLE;
00203 */
00204     }
00205     if (strval) free (strval);
00206    /*  the following are appropriate for MDI, but not strictly based on WCS  */
00207 /*
00208     sprintf (param[RSUN].name, "R_SUN");
00209     sprintf (param[APSD].name, "OBS_ASD");
00210 */
00211     strncpy (param[RSUN].name, rsun_key, 31);
00212     strncpy (param[APSD].name, apsd_key, 31);
00213     sprintf (param[PANG].name, "CROTA2");
00214     param[PANG].scale = -raddeg;
00215     if (AIPS_convention) param[PANG].scale *= -1;
00216     sprintf (param[ESMA].name, "S_MAJOR");
00217     sprintf (param[ESMI].name, "S_MINOR");
00218     sprintf (param[EANG].name, "S_ANGLE");
00219     param[EANG].scale = -raddeg;
00220     if (AIPS_convention) param[EANG].scale *= -1;
00221     }
00222   }
00223   lasthdr = hdrtype;
00224                     /*  Plate info: image scale, distortion  */
00225   *apsd = lookup (img, param[RSUN], &status);
00226   if (scale_avail) {
00227     *xscl = lookup (img, param[XSCL], &status);
00228     *yscl = lookup (img, param[YSCL], &status);
00229     if (status & NO_SEMIDIAMETER) {
00230       status &= ~NO_SEMIDIAMETER;
00231       *apsd = lookup (img, param[APSD], &status);
00232       if (status & NO_SEMIDIAMETER) {
00233         *need_ephem = 1;
00234       } else {
00235     if (!(status & (NO_XSCALE | NO_YSCALE))) {
00236       *apsd /= (*xscl <= *yscl) ? *xscl : *yscl;
00237       status &= ~NO_SEMIDIAMETER;
00238     }
00239       }
00240     }
00241   }
00242   ella = lookup (img, param[ESMA], &status);
00243   ellb = lookup (img, param[ESMI], &status);
00244   *ellipse_e = sqrt ((ella - ellb) * (ella + ellb)) / ella;
00245   *ellipse_pa = lookup (img, param[EANG], &status);
00246                     /*  Pointing (attitude: image location)  */
00247   *ctrx = lookup (img, param[CTRX], &status);
00248   *ctry = lookup (img, param[CTRY], &status);
00249                       /*  Position angle of solar north  */
00250   *pang = lookup (img, param[PANG], &status);
00251                                                       /*  Image orientation  */
00252   *x_invrt = xinv_type;
00253   *y_invrt = yinv_type;
00254 
00255   return status;
00256 }
00257 
00258 /*
00259  *  Revision History (all mods by Rick Bogart unless otherwise noted)
00260  *
00261  *  09.10.06        version that went into first release under proj/rings
00262  *  09.12.03        fixed two icc11 compiler warnings; added keyword name
00263  *      passing to solar_image_info for selected cases
00264  *  10.04.01        fixed up processing for semidiameter when apsd_key
00265  *      is missing
00266  *  10.08.19        make sure status is initialized to 0
00267  *  10.09.03        fixed initialization of default values for ellipse
00268  *      elements
00269  *  10.10.05        fixed define of keyscope_variable (removed ;)
00270  *  12.04.17        changed return value of pang to conform to Thompson
00271  *      (solar WCS) convention rather than the AIPS convention for
00272  *      CROTA2; added argument AIPS_convention to solar_image_info to
00273  *      be set if the CROTA2 keyword conforms to the AIPS convention
00274  *  14.10.04        added AIPS-convention correction for ellipse angle
00275  *  15.06.19        corrected erroneous overwrite of apsd_key
00276  *          added definition of code version number
00277  */
00278