JSOC_Documentation: proj/mag/synop/apps/synop-imginfo.c Source File

00001 #include "drms.h"
00002 #include "cartography.h"
00003 
00004 // There were a bunch of defines and a struct here. However, these defines and the struct
00005 // were already defined in cartography.h.
00006 static double lookup (DRMS_Record_t *rec, ParamDef key, int *status) {
00007   double value = key.defval;
00008   int lookupstat = 0;
00009 
00010   value = drms_getkey_double (rec, key.name, &lookupstat);
00011   value = value * key.scale + key.offset;
00012   if (lookupstat) *status |= key.statusbit;
00013   if (isnan (value)) *status |= key.statusbit;
00014   return value;
00015 }
00016 
00017 static char *lookup_str (DRMS_Record_t *rec, ParamDef key, int *status) {
00018   DRMS_Keyword_t *keywd;
00019   char *value;;
00020 
00021   value = drms_getkey_string (rec, key.name, status);
00022   if (*status) *status |= key.statusbit;
00023                          /*  cadence should be constant  */
00024   if ((keywd = drms_keyword_lookup (rec, key.name, 1))) {
00025     if (keywd->info->recscope != 1) *status |= KEYSCOPE_VARIABLE;
00026   } else *status |= key.statusbit;
00027   return value;
00028 }
00029 
00030 int synop_solar_image_info (DRMS_Record_t *img, double *xscl, double *yscl,
00031     double *ctrx, double *ctry, double *apsd, const char *rsun_key,
00032     const char *apsd_key, double *pang, double *ellipse_e, double *ellipse_pa,
00033     int *x_invrt, int *y_invrt, int *need_ephem, int AIPS_convention) {
00034 /*
00035  *  Provides the following values from the DRMS record:
00036  *    xscl  scale in the image column direction (arc-sec/pixel)
00037  *    yscl  scale in the image row direction (arc-sec/pixel)
00038  *    ctrx  (virtual) fractional pixel column of the center of the
00039  *  solar disc
00040  *    ctry  (virtual) fractional pixel row of the center of the
00041  *  solar disc
00042  *    NO latc  heliographic latitude of disc center, in radians
00043  *    NO lonc  heliographic longitude of disc center, in radians
00044  *    apsd  apparent semi-diameter (semimajor-axis) of the solar disc, in
00045  *  pixel units
00046  *    pang  position angle of solar north relative to image vertical
00047  *  ([0,0] -> [0,1]), measured westward (clockwise), in radians
00048  *    eecc  eccentricity of best-fit ellipse describing limb
00049  *    eang  position angle of best-fit ellipse describing limb, relative
00050  *  to ?, in radians
00051  *    xinv  0 if image is direct, 1 if flipped by columns
00052  *    yinv  0 if image is direct, 1 if flipped by rows
00053  *  NO!
00054  *  The following data types are supported: SOHO-MDI, GONG+, Mt. Wilson MOF,
00055  *    SOHO-EIT, TRACE, BBSO Ha
00056  *  NO!
00057  *  If the data are not recognizably of one of these types, the function
00058  *    returns NO_DATA_DICT; if one or more required keywords are missing
00059  *    the function returns a status mask indicating which values could not
00060  *    be filled.  No flag is set if the image ellipse parameters are not
00061  *    available (unless they are required for other parameters), but the
00062  *    eccentricity and pericenter angle are quietly set to 0.
00063  */
00064   enum param {
00065     XSCL, YSCL, XUNI, YUNI, LATC, LONC, CTRX, CTRY, PANG, APSD, RSUN,
00066     ESMA, ESMI, EECC, EANG, XASP, YASP, PCT
00067   };
00068   enum known_plat {
00069     UNKNOWN
00070   };
00071   static ParamDef param[PCT];
00072   static double raddeg =  M_PI / 180.0;
00073   static double degrad = 180.0 / M_PI;
00074   double ella, ellb;
00075   int n, status = 0;
00076   static int scale_avail, xinv_type, yinv_type;
00077   static int hdrtype = UNKNOWN, lasthdr = UNKNOWN - 1;
00078   char *strval;
00079 /*
00080  *  Set up the appropriate dictionary for interpretation of keywords
00081  */
00082   if (lasthdr != hdrtype) {
00083     if (lasthdr >= UNKNOWN)
00084       fprintf (stderr,
00085       "Warning from solar_image_info(): record keywords may change\n");
00086     for (n = 0; n < PCT; n++) {
00087       sprintf (param[n].name, "No Keyword");
00088       param[n].scale = 1.0;
00089       param[n].offset = 0.0;
00090       param[n].defval = NAN;
00091       param[n].statusbit = 0;
00092     }
00093     param[RSUN].statusbit = NO_SEMIDIAMETER;
00094     param[APSD].statusbit = NO_SEMIDIAMETER;
00095     param[XSCL].statusbit = NO_XSCALE;
00096     param[YSCL].statusbit = NO_YSCALE;
00097     param[XUNI].statusbit = NO_XUNITS;
00098     param[YUNI].statusbit = NO_YUNITS;
00099     param[CTRX].statusbit = NO_XCENTERLOC;
00100     param[CTRY].statusbit = NO_YCENTERLOC;
00101     param[CTRX].statusbit = NO_XCENTERLOC;
00102     param[CTRY].statusbit = NO_YCENTERLOC;
00103     param[LATC].statusbit = NO_HELIO_LATC;
00104     param[LONC].statusbit = NO_HELIO_LONC;
00105     param[PANG].statusbit = NO_HELIO_PA;
00106     param[EECC].defval = 0.0;
00107     param[EANG].defval = 0.0;
00108     param[XASP].defval = 1.0;
00109     param[YASP].defval = 1.0;
00110 
00111     switch (hdrtype) {
00112       default:                        /*  Assume WCS HPLN/T  */
00113                      /*  WITH CERTAIN MDI SPECIFIC ONES  */
00114     scale_avail = 1;
00115     xinv_type = yinv_type = 0;
00116     sprintf (param[XUNI].name, "CUNIT1");
00117     sprintf (param[YUNI].name, "CUNIT2");
00118     sprintf (param[XSCL].name, "CDELT1");
00119     sprintf (param[YSCL].name, "CDELT2");
00120     sprintf (param[CTRX].name, "CRPIX1");
00121     param[CTRX].offset = -1.0;
00122     sprintf (param[CTRY].name, "CRPIX2");
00123     param[CTRY].offset = -1.0;
00124     *need_ephem = 0;
00125     strval = lookup_str (img, param[XUNI], &status);
00126     if (!(status & NO_XUNITS)) {
00127       if (!strcmp (strval, "arcsec")) param[XSCL].scale = 1.0;
00128       else if (!strcmp (strval, "arcmin")) param[XSCL].scale = 1.0 / 60.0;
00129       else if (!strcmp (strval, "deg")) param[XSCL].scale = 1.0 / 3600.0;
00130       else if (!strcmp (strval, "mas")) param[XSCL].scale = 1000.0;
00131       else if (!strcmp (strval, "rad")) param[XSCL].scale = degrad * 3600.0;
00132 /*
00133       need_units = status & KEYSCOPE_VARIABLE;
00134 */
00135     }
00136     if (strval) free (strval);
00137     strval = lookup_str (img, param[YUNI], &status);
00138     if (!(status & NO_YUNITS)) {
00139       if (!strcmp (strval, "arcsec")) param[YSCL].scale = 1.0;
00140       else if (!strcmp (strval, "arcmin")) param[YSCL].scale = 1.0 / 60.0;
00141       else if (!strcmp (strval, "deg")) param[YSCL].scale = 1.0 / 3600.0;
00142       else if (!strcmp (strval, "mas")) param[YSCL].scale = 1000.0;
00143       else if (!strcmp (strval, "rad")) param[YSCL].scale = degrad * 3600.0;
00144 /*
00145       need_units = status & KEYSCOPE_VARIABLE;
00146 */
00147     }
00148     if (strval) free (strval);
00149    /*  the following are appropriate for MDI, but not strictly based on WCS  */
00150 /*
00151     sprintf (param[RSUN].name, "R_SUN");
00152     sprintf (param[APSD].name, "OBS_ASD");
00153 */
00154     strncpy (param[RSUN].name, rsun_key, 31);
00155     strncpy (param[APSD].name, apsd_key, 31);
00156     sprintf (param[APSD].name, "OBS_ASD");
00157     sprintf (param[PANG].name, "CROTA2");
00158     param[PANG].scale = -raddeg;
00159     if (AIPS_convention) param[PANG].scale *= -1;
00160     sprintf (param[ESMA].name, "S_MAJOR");
00161     sprintf (param[ESMI].name, "S_MINOR");
00162     sprintf (param[EANG].name, "S_ANGLE");
00163     param[EANG].scale = raddeg;
00164     }
00165   }
00166   lasthdr = hdrtype;
00167                     /*  Plate info: image scale, distortion  */
00168   *apsd = lookup (img, param[RSUN], &status);
00169   if (scale_avail) {
00170     *xscl = lookup (img, param[XSCL], &status);
00171     *yscl = lookup (img, param[YSCL], &status);
00172     if (status & NO_SEMIDIAMETER) {
00173       status &= ~NO_SEMIDIAMETER;
00174       *apsd = lookup (img, param[APSD], &status);
00175       if (status & NO_SEMIDIAMETER) {
00176         *need_ephem = 1;
00177       } else {
00178     if (!(status & (NO_XSCALE | NO_YSCALE))) {
00179       *apsd /= (*xscl <= *yscl) ? *xscl : *yscl;
00180       status &= ~NO_SEMIDIAMETER;
00181     }
00182       }
00183     }
00184   }
00185   ella = lookup (img, param[ESMA], &status);
00186   ellb = lookup (img, param[ESMI], &status);
00187   *ellipse_e = sqrt ((ella - ellb) * (ella + ellb)) / ella;
00188   *ellipse_pa = lookup (img, param[EANG], &status);
00189                     /*  Pointing (attitude: image location)  */
00190   *ctrx = lookup (img, param[CTRX], &status);
00191   *ctry = lookup (img, param[CTRY], &status);
00192                       /*  Position angle of solar north  */
00193   *pang = lookup (img, param[PANG], &status);
00194                                                       /*  Image orientation  */
00195   *x_invrt = xinv_type;
00196   *y_invrt = yinv_type;
00197 
00198   return status;
00199 }
00200 
00201 /*
00202  *  Revision History (all mods by Rick Bogart unless otherwise noted)
00203  *  
00204  *  09.10.06        version that went into first release under proj/rings
00205  *  09.12.03        fixed two icc11 compiler warnings; added keyword name
00206  *      passing to solar_image_info for selected cases
00207  *  10.04.01        fixed up processing for semidiameter when apsd_key
00208  *      is missing
00209  *  10.08.19        make sure status is initialized to 0
00210  *  10.09.03        fixed initialization of default values for ellipse
00211  *      elements
00212  *  10.10.05        fixed define of keyscope_variable (removed ;)
00213  *  12.04.17        changed return value of pang to conform to Thompson
00214  *      (solar WCS) convention rather than the AIPS convention for
00215  *      CROTA2; added argument AIPS_convention to solar_image_info to
00216  *      be set if the CROTA2 keyword conforms to the AIPS convention
00217  * Feb-18-2016      This was modified by Yang for use with two modules (I do not know the nature of the changes). Copied and enamed to synop-imginfo.c. For use only by these two modules - ART
00218  */
00219