JSOC_Documentation: proj/libs/projection/imginfo.c Source File

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