JSOC_Documentation: proj/maps_avgs/apps/derot

00001 /* derot_mean.c */
00002 
00003 /* **************************************************************
00004  *
00005  * Module Name: derot_mean
00006  * Purpose:  Produces remapped, weight averages of input data sets
00007  * Programmer: John Beck    Beck@Sun.Stanford.EDU
00008  *
00009  * Notes: This is based on the DSDS module dr_mean (written by 
00010  * J. Beck & R. Bogart)
00011  *
00012  * Usage: derot_mean InData= OutData= [flags]
00013  * 
00014  * **************************************************************
00015  */
00016 
00100 #include <stdio.h>
00101 #include <stdlib.h>
00102 #include <strings.h>
00103 #include <string.h>
00104 #include <ctype.h>
00105 #include <time.h>
00106 #include <math.h>
00107 #include <limits.h>
00108 #include <complex.h>
00109 #include <sys/time.h>
00110 #include <jsoc_main.h>
00111 //#include <drms_types.h>
00112 #include "./derot_mean.h"
00113 #include "./cartography.c"
00114 #include "./ccinterp.c"
00115 //  #include "./heliographic_coords.c"
00116 #include <astro.h>
00117 #include <assert.h>
00118 
00119 char *module_name = "derot_mean";
00120 char *version_id  = "0.0.1";
00121 
00122 ModuleArgs_t module_args[] =
00123 {
00124   {ARG_STRING,  "in",      "",       "Input data series."},
00125   {ARG_STRING,  "out",     "",       "Output data series"},
00126   {ARG_TIME,    "start_t", "",       "First output record time "},     // 1
00127   {ARG_TIME,    "end_t",   "-1",     "Last output record time (or upper limit)"},     
00128   {ARG_INT,     "step_t",  "900",    "Steps between records in seconds"},
00129   {ARG_INT,     "cols",    "-1",     "Number of columns in output image"},  
00130   {ARG_INT,     "rows",    "-1",     "Number of rows"},  
00131   {ARG_FLOAT,   "a0",      " 2.851", "Solid Body Rotation rate (uRad/s)"},
00132   {ARG_FLOAT,   "a2",      "0.0",    "A2 rotation term (uRad/S)"},  // set to zero for test
00133   {ARG_FLOAT,   "a4",      "0.0",    "A4 rotation term (uRad/s)"},  // set to zero for test
00134 //  {ARG_FLOAT,   "a2",      "-0.343", "A2 rotation term (uRad/S)"},
00135 //  {ARG_FLOAT,   "a4",      "-0.474", "A4 rotation term (uRad/s)"},
00136   {ARG_FLOAT,   "merid_v", "0.0",    "Meridional Flow in m/s"},         // 9
00137   {ARG_STRING,  "wt_func", "hath",   "Weighting Filter Function Name"}, // 10
00138   {ARG_INT,     "wt_len",  "1860.0", "Width of Weighting Window in Seconds"},
00139   {ARG_FLOAT,   "wt_parm", "900.0",  "Weighting Filter Parameter"},
00140   {ARG_STRING,  "proj",    "ortho",  "Projection Name (MDI style)"},
00141   {ARG_FLOAT,   "crot",    "0.0",    "WCS Coord Rotation Angle (-pangle)"},   
00142   {ARG_FLOAT,   "crpix1",  "None",    "Reference Pixel (Columnwise)"},     
00143   {ARG_FLOAT,   "crpix2",  "None",    "Reference Pixel (Rowwise)"},   
00144   {ARG_FLOAT,   "crval1",  "",    "Coord at Ref Pixel (deg/arcsec)"},
00145   {ARG_FLOAT,   "crval2",  "",    "Coord at Ref Pixel (deg/arcsec)"},
00146   {ARG_FLOAT,   "cdelt1",  "None",    "Pixel size at Ref Pix (in appropriate units)"},
00147   {ARG_FLOAT,   "cdelt2",  "None",    "Pixel size at Ref Pix (in appropriate units)"},
00148   {ARG_STRING,  "ctype1",  "CRLN-SIN", "Coordinate Name for col axis"},
00149   {ARG_STRING,  "ctype2",  "CRLT-SIN", "Coordinate Name for row axis"}, // 22
00150   {ARG_FLOAT,   "ang_rsun", "960.",  "Solar Radius in arcsecs (if AERIAL)"},
00151   {ARG_FLOAT,   "dist",    "1.0",    "Obs-Solar Distance in AU (if AERIAL)"}, 
00152 //  {ARG_STRING,  "qualfile", "", "", ""},     Not used
00153   {ARG_FLAG,    "v",        "",      "Run in Verbose Mode"},
00154   {ARG_END}
00155 };
00156 
00157 ParmList parms;
00158 
00159 /* These functions involve reading, checking cmdline params & subsequent calcs */
00160 int    getcheck_cmd_parms(void);
00161 int    get_cparms(void);
00162 int    check_cparms(void);
00163 int    cparm_getcheck_float(char kword[], double *out);
00164 int    cparm_getcheck_int(char kword[], int *out);
00165 int    cparm_getcheck_time(char kword[], TIME *out);
00166 int    cparm_getcheck_str(char kword[], char out[]);
00167 int    wcs_map_proj(void);
00168 int    calc_weight_function(void);
00169 
00170 // These functions involve setting up output arrays & db
00171 int    setup_outdb(OutImgs **outimg, DRMS_RecordSet_t **outDB);
00172 int    check_outdb(DRMS_Record_t *rec);
00173 int    init_outimgs(OutImgs **outimg);
00174 int    calc_out_geom(OutImgs outimg);
00175 
00176 // These functions involve opening & testing input db for conformity
00177 int    open_check_indb(DRMS_RecordSet_t **inDB);
00178 int    get_query_str(char *qry);
00179 int    check_indb(DRMS_Record_t *rec);
00180 
00181 // These involve the actual derotation summation
00182 int    addin(DRMS_RecordSet_t *inDB, OutImgs *out);
00183 int    readin(DRMS_Record_t *irec, TIME *timg, double *Bo, double *Lo, 
00184           double *Rsun, double *phi,
00185               double *xc, double *yc, char ctype1[], char ctype2[]);
00186 int    rec_getcheck_double(DRMS_Record_t *rec, char kword[], double *out);
00187 int    rec_getcheck_int(DRMS_Record_t *rec, char kword[], int *out);
00188 int    rec_getcheck_str(DRMS_Record_t *rec, char kword[], char out[]);
00189 int    rec_getcheck_time(DRMS_Record_t *rec, char kword[], TIME *out);
00190 
00191 // These involve  writing out records
00192 int    write_out_dbs(OutImgs *outimgs, DRMS_RecordSet_t *outDB, DRMS_RecordSet_t *inDB);
00193 int    set_bscale_bzero(OutImgs *outimg);
00194 
00195 // Some Utility & cleanup functions
00196 double wcs_parm_unit(double parm, char *unit);
00197 int    free_outimgs(OutImgs *outimg);
00198 int    If_Err_Print(void);
00199 void   V_print(char *str);
00200 
00201 /* 
00202 **********************************************************************
00203 *
00204 *    THIS IS THE DOIT FUNCTION WHERE THE TOP CONTROL TAKES PLACE
00205 *
00206 *    Note, this is where all error writing occurs.  Errors are passed
00207 *    using the parms structure
00208 *
00209 **********************************************************************
00210 */
00211 
00212 int DoIt(void) 
00213 {
00214   int status, done;
00215   OutImgs *out_imgs;
00216   DRMS_RecordSet_t *inDB  = NULL;
00217   DRMS_RecordSet_t *outDB = NULL;
00218   
00219   V_print("\n===== Checking Commandline Arguments =====\n");
00220   if(getcheck_cmd_parms() !=  kMyMod_Success) // read command parms
00221     return (status = If_Err_Print());
00222 
00223   V_print("\n===== Checking Output Data Series =====\n");
00224   if(setup_outdb(&out_imgs, &outDB) != kMyMod_Success) // open & check output DB
00225     return (status = If_Err_Print());
00226 
00227   V_print("\n===== Checking Input Data Series =====\n");
00228   if(open_check_indb(&inDB) !=  kMyMod_Success ) // open & check  input DB
00229     return (status = If_Err_Print());
00230 
00231   V_print("\n===== Summing Derotated Mean =====\n");
00232   if(addin(inDB, out_imgs) !=  kMyMod_Success) // Add input images to output
00233     return (status = If_Err_Print());
00234 
00235   V_print("\n===== Writing Records to Output =====\n\n\n");
00236   if(write_out_dbs(out_imgs, outDB, inDB) !=  kMyMod_Success)  // Write out all output images
00237     return (status = If_Err_Print());
00238 
00239   V_print("Done.  derot_mean exiting normally\n");
00240   return kMyMod_Success;
00241 }
00242 
00243 
00244 
00245 /* 
00246  ******************************************************************* 
00247  * 
00248  *  Functions involving reading & calculating cmd-line arguments
00249  * 
00250  ******************************************************************** 
00251 */
00252 //GETCHECK_CMD_PARMS
00253 int  getcheck_cmd_parms(void)
00254 {
00255   int status;
00256 
00257   strcpy(parms.SigMsg," ");
00258   strcpy(parms.Msg," ");
00259   if((status = get_cparms()))
00260     return status;
00261 
00262   if((status = wcs_map_proj()))
00263     return status;
00264 
00265   if((status = check_cparms()))
00266     return status;
00267 
00268   if((status = calc_weight_function()))
00269     return status;
00270   
00271   V_print("OK\n");
00272   return kMyMod_Success;    // Exit getcheck_cmd_params normally
00273 }
00274 
00275 // GET_CPARMS
00276 int  get_cparms(void)
00277 {
00278   int status;
00279   /* Have to allocate memory for the strings */
00280   
00281   /* 
00282      Note: this only returns an error if the reading is an error
00283      it does not check that the values are valid.  That is done
00284      later as it requires more complicated logic 
00285   */
00286   if((status = cparm_getcheck_str("in",(parms.iser))))
00287     return status;
00288   if((status = cparm_getcheck_str("out",(parms.oser))))
00289     return status;
00290 
00291   if((status = cparm_getcheck_time("start_t",&(parms.Tstart))))
00292     return status;
00293   if((status = cparm_getcheck_time("end_t",&(parms.Tend))))
00294     return status;
00295   if((status = cparm_getcheck_int("step_t",&(parms.Tstep))))
00296     return status;
00297   if((status = cparm_getcheck_int("rows",&(parms.rows))))
00298     return status;
00299   if((status = cparm_getcheck_int("cols",&(parms.cols))))
00300     return status;
00301   if((status = cparm_getcheck_float("a0",&(parms.A0))))
00302     return status;
00303   if((status = cparm_getcheck_float("a2",&(parms.A2))))
00304     return status;
00305   if((status = cparm_getcheck_float("a4",&(parms.A4))))
00306     return status;
00307   if((status = cparm_getcheck_float("merid_v",&(parms.Meri_V))))
00308     return status;
00309 
00310   if((status = cparm_getcheck_str("wt_func",(parms.WtFunc))))
00311     return status;
00312   if((status = cparm_getcheck_int("wt_len",&(parms.WtLen))))
00313     return status;
00314   if((status = cparm_getcheck_float("wt_parm",&(parms.WtParm))))
00315     return status;
00316 
00317   if((status = cparm_getcheck_str("proj",(parms.projname))))
00318     return status;
00319   if((status = cparm_getcheck_str("ctype1",(parms.CTYPE1))))
00320     return status;
00321   if((status = cparm_getcheck_str("ctype2",(parms.CTYPE2))))
00322     return status;
00323   if((status = cparm_getcheck_float("crot",&(parms.CROTA2))))
00324     return status;
00325   if((status = cparm_getcheck_float("crpix1",&(parms.CRPIX1))))
00326     return status;
00327   if((status = cparm_getcheck_float("crpix2",&(parms.CRPIX2))))
00328     return status;
00329   if((status = cparm_getcheck_float("crval1",&(parms.CRVAL1))))
00330     return status;
00331   if((status = cparm_getcheck_float("crval2",&(parms.CRVAL2))))
00332     return status;
00333   if((status = cparm_getcheck_float("cdelt1",&(parms.CDELT1))))
00334     return status;
00335   if((status = cparm_getcheck_float("cdelt2",&(parms.CDELT2))))
00336     return status;
00337   if((status = cparm_getcheck_float("ang_rsun",&(parms.ang_rad))))
00338     return status;
00339   if((status = cparm_getcheck_float("dist",&(parms.dist))))
00340     return status;
00341 
00342   parms.verbose = cmdparams_isflagset(&cmdparams,"v");
00343 
00344   return kMyMod_Success;     // exit get_cparms normally
00345 }
00346 
00347 // CHECK_CPARMS
00348 int check_cparms(void)
00349 {
00350   char temp[100];
00351 
00352   if(parms.Tend < 0)                   //  if Tend is not set, set it to Tstart
00353     parms.Tend = parms.Tstart;
00354   else if(parms.Tend < parms.Tstart)   // First check that Tend >= Tstart
00355   {
00356     fprintf(stderr,"Error: Tend must be greater  or equal to Tstart\n");
00357     return (parms.status = kMyMod_ValErr);
00358   }
00359   if (parms.Tstep == 0)
00360   {
00361     fprintf(stderr,"Error: Tstep must be greater or equal to zero\n");
00362     return (parms.status = kMyMod_ValErr);
00363   }
00364   else if (parms.Tstep < 0)
00365   {
00366     V_print("Warning: Tstep < 0.  Absolute value taken\n");
00367     parms.Tstep = abs(parms.Tstep);
00368   }
00369 
00370   // Then verify that required geometry values are valid
00371   // (Note: CRPIXs are checked in check_outdb() where they can be set
00372   // to 0.5*cols or rows
00373 
00374   if( isnan(parms.CROTA2))
00375     parms.CROTA2 = 0.0;
00376 
00377   if( isnan(parms.CRVAL1) || isnan(parms.CRVAL2) )
00378   {
00379     fprintf(stderr,"Error: Missing or Invalid CRVALs\n");
00380     return (parms.status = kMyMod_ValErr);
00381   }
00382   if( isnan(parms.CDELT1) || isnan(parms.CDELT2) )
00383   {
00384     fprintf(stderr,"Error: Missing or Invalid CDELTs\n");
00385     return (parms.status = kMyMod_ValErr);
00386   }
00387   if( isnan(parms.A0) || isnan(parms.A2) || 
00388       isnan(parms.A4) || isnan(parms.Meri_V) )
00389   {
00390     fprintf(stderr,"Error: Invalid Rotation Coefs or V_meridional\n");
00391     return (parms.status = kMyMod_ValErr);
00392   }
00393 
00394   if(  isnan(parms.WtParm) || parms.WtLen < 1)
00395   {
00396     fprintf(stderr,"Error: Invalid Weighting Function Parameters\n");
00397     return (parms.status = kMyMod_ValErr);
00398   }
00399  
00400   // Convert CRVALs & CDELTs to appropriate (Projection specific) units
00401   if(parms.projcode == AERIAL)
00402   {
00403     if(parms.CRVAL1 != 0 || parms.CRVAL2 != 0)
00404     {
00405       sprintf(temp,"Warning: AERIAL map requires CRVALn to be Zero - setting to Zero\n");
00406       V_print(temp);
00407       parms.CRVAL1 = parms.CRVAL2 = 0.0;
00408     }
00409     if(parms.ang_rad < 0 && parms.dist < 0)
00410     {
00411       sprintf(parms.Msg,"%sWarning: AERIAL map requires ang_rsun or dist0\n",parms.Msg);
00412       return (parms.status = kMyMod_ValErr);
00413     }
00414     if(parms.ang_rad <= 0)
00415       parms.ang_rad = asin( 1.0/(parms.dist*214.93950));
00416     if(parms.dist <= 0)
00417       parms.dist = 1.0 / (sin(parms.ang_rad)*214.93950);    
00418 
00419     parms.CDELT1 *= (M_PI / 180.0 / 3600.0);   // Convert from arcsec to rad
00420     parms.CDELT2 *= (M_PI / 180.0 / 3600.0);   // Convert from arcsec to rad
00421 
00422     if(parms.rsun <= 0  && (parms.CDELT1 == parms.CDELT2) )
00423       parms.rsun = parms.ang_rad / parms.CDELT1;
00424   }
00425   else if (parms.projcode == LAMBERT || parms.projcode == ORTHOGRAPHIC ||  
00426        parms.projcode == STEREOGRAPHIC || parms.projcode == POSTEL || 
00427        parms.projcode == GNOMONIC || parms.projcode == RECTANGULAR)
00428   {
00429     parms.CRVAL1 *= (M_PI / 180.0);  // Convert from deg to rad
00430     parms.CRVAL2 *= (M_PI / 180.0);  // Convert from deg to rad
00431     parms.CDELT1 *= (M_PI / 180.0);  // Convert from deg to rad
00432     parms.CDELT2 *= (M_PI / 180.0);  // Convert from deg to rad
00433   }
00434   else if (parms.projcode == CYLEQA || parms.projcode == SINEQA ||  
00435        parms.projcode == MERCATOR || parms.projcode == CASSINI)
00436   {
00437     parms.CRVAL1 *= (M_PI / 180.0);  // Convert from deg to rad
00438     parms.CDELT1 *= (M_PI / 180.0);  // Convert from deg to rad
00439   }
00440 
00441   //  If projcode is CYLEQA, SINEQA, MERCATOR, CASSINI do not convert!!
00442   
00443   // Do some angle converstions & fill in some values
00444   parms.CROTA2  *= (M_PI / 180.0);  // Convert from deg to rad
00445   parms.pangle -1.0 * parms.CROTA2; // Change sign 
00446   parms.A0  =   parms.A0 / 1.0e6;   // convert from uRad/s to Rad/s
00447   parms.A2  =   parms.A2 / 1.0e6;   // convert from uRad/s to Rad/s
00448   parms.A4  =   parms.A4 / 1.0e6;   // convert from uRad/s to Rad/s
00449   parms.Meri_V  =   parms.Meri_V / 6.955e8;   // convert from m/s to Rad/s
00450   parms.size = parms.rows * parms.cols; // Calc number of pixels
00451   parms.Nmaps = ceil(((1 + parms.Tend - parms.Tstart)/parms.Tstep));
00452   if(parms.Nmaps  < 1) 
00453     parms.Nmaps = 1;
00454 
00455   return kMyMod_Success;     // exit check_cparms normally
00456 }
00457 
00458 // CPARM_GETCHECK_FLOAT
00459 int cparm_getcheck_float(char kword[], double *out)
00460 {
00461   int stat;
00462   double temp;
00463   temp = cmdparams_get_float(&cmdparams,kword,&stat);
00464   if(stat)
00465   {
00466     sprintf(parms.Msg,"%sError %d: could not read cmd param %s\n",parms.Msg,stat,kword);
00467     return (parms.status = stat);
00468   }
00469   (*out) = temp;
00470   return kMyMod_Success;     // exit cparm_getcheck_float normally
00471 
00472 }
00473 
00474 // CPARM_GETCHECK_INT
00475 int cparm_getcheck_int(char kword[], int *out)
00476 {
00477   int temp, stat;
00478   
00479   temp = cmdparams_get_int(&cmdparams,kword,&stat);
00480   if(stat)
00481   {
00482     sprintf(parms.Msg,"%sError %d: could not read cmd param %s\n",parms.Msg,stat,kword);
00483     return (parms.status = stat);
00484   }
00485   (*out) = temp;
00486   return kMyMod_Success;     // exit cparm_getcheck_int normally
00487   
00488 }
00489 
00490 // CPARM_GETCHECK_TIME
00491 int cparm_getcheck_time(char kword[], TIME *out)
00492 {
00493   int stat;
00494   TIME temp;
00495   temp = cmdparams_get_time(&cmdparams,kword,&stat);
00496   if(stat)
00497   {
00498     sprintf(parms.Msg,"%sError %d: could not read cmd param %s\n",parms.Msg,stat,kword);
00499     return (parms.status = stat);
00500   }
00501   (*out) = temp;
00502   return kMyMod_Success;     // exit cparm_getcheck_time normally
00503   
00504 }
00505 
00506 // CPARM_GETCHECK_STR
00507 int cparm_getcheck_str(char kword[], char out[])
00508 {
00509   int stat;
00510   const char *temp;
00511   temp = cmdparams_get_str(&cmdparams,kword,&stat);
00512   if(stat)
00513   {
00514     sprintf(parms.Msg,"%sError %d: could not read cmd param %s\n",parms.Msg,stat,kword);
00515     return (parms.status = stat);
00516   }
00517   strcpy(out,temp);
00518   return  kMyMod_Success;     // exit cparm_getcheck_str normally
00519   
00520 }
00521 
00522 // WCS_MAP_PROJ
00523 int wcs_map_proj(void)
00524 {
00525   int i, done;
00526   char comp1[12],comp2[12];
00527   static struct wcs_projnames 
00528   {
00529     int code;
00530     char CTYPE1[9];
00531     char CTYPE2[9];
00532     char projname[20];
00533   } wcs[]= {
00534     {RECTANGULAR,  "CRLN-CAR",  "CRLT-CAR", "Plate-Caree"},
00535     {RECTANGULAR,  "HGLN-CAR",  "HGLT-CAR", "Plate-Caree"},
00536     {CASSINI,      "CRLN-CAS",  "CRLT-CAS", "Cassini"},
00537     {CASSINI,      "HGLN-CAS",  "HGLT-CAS", "Cassini"},
00538     {MERCATOR,     "CRLN-MER",  "CRLT-MER", "Mercator"},
00539     {MERCATOR,     "HGLN-MER",  "HGLT-MER", "Mercator"},
00540     {CYLEQA,       "CRLN-CEA",  "CRLT-CEA", "CylEqA"},
00541     {CYLEQA,       "HGLN-CEA",  "HGLT-CEA", "CylEqA"},
00542     {SINEQA,       "CRLN-SFL",  "CRLT-SFL", "Sanson-Flamsteed"},
00543     {SINEQA,       "HGLN-SFL",  "HGLT-SFL", "Sanson-Flamsteed"},
00544     {GNOMONIC,     "CRLN-TAN",  "CRLT-TAN", "Gnomonic"},
00545     {GNOMONIC,     "HGLN-TAN",  "HGLT-TAN", "Gnomonic"},
00546     {POSTEL,       "CRLN-ARC",  "CRLT-ARC", "Postel"},
00547     {POSTEL,       "HGLN-ARC",  "HGLT-ARC", "Postel"},
00548     {STEREOGRAPHIC,"CRLN-STG",  "CRLT-STG", "Stereographic"},
00549     {STEREOGRAPHIC,"HGLN-STG",  "HGLT-STG", "Stereographic"},
00550     {ORTHOGRAPHIC, "CRLN-SIN",  "CRLT-SIN", "Orthographic"},
00551     {ORTHOGRAPHIC, "HGLN-SIN",  "HGLT-SIN", "Orthographic"},
00552     {LAMBERT,      "CRLN-ZEA",  "CRLT-ZEA", "Lambert"},
00553     {LAMBERT,      "HGLN-ZEA",  "HGLT-ZEA", "Lambert"},
00554     {AERIAL,       "HPLN-TAN",  "HPLT-TAN", "Aerial"},
00555     {-1, "END ", ""}
00556   };
00557   
00558   done = i = 0;
00559   while(strncasecmp(wcs[i].CTYPE1,"END ",4) && !done)  // First check by CTYPE
00560   {
00561     if(!strncasecmp(parms.CTYPE1,wcs[i].CTYPE1,8) &&
00562        !strncasecmp(parms.CTYPE2,wcs[i].CTYPE2,8))
00563     {
00564       parms.projcode = wcs[i].code;
00565       strcpy(parms.projname,wcs[i].projname);
00566 //      printf("PROJ_CODE: %s %s - %s\n",parms.CTYPE1,parms.CTYPE2,parms.projname);
00567       done = 1;
00568     }
00569     i++;
00570   }
00571   i = 0;                              // If not done, check by projname
00572   while(strncasecmp(wcs[i].CTYPE1,"END ",4) && !done) 
00573   {
00574     if(!strncasecmp(parms.projname,wcs[i].projname,3))
00575     {
00576       parms.projcode = wcs[i].code;
00577       strcpy(parms.CTYPE1,wcs[i].CTYPE1);
00578       strcpy(parms.CTYPE2,wcs[i].CTYPE2);
00579       printf("PROJ_NAME: %s - %s %s\n",parms.projname,parms.CTYPE1,parms.CTYPE2);
00580       done = 1;
00581     }
00582     i++;
00583   }
00584   if(done)
00585     return kMyMod_Success;     // exit wcs_map_proj normally
00586 
00587   sprintf(parms.Msg,"%sError: Undetermined Map Projection- %s, %s, %s,\n",
00588       parms.Msg, parms.CTYPE1, parms.CTYPE2, parms.projname);
00589   return (parms.status = kMyMod_ValErr);
00590 }
00591 
00592 
00593 // CALC_WEIGHT_FUNCTION
00594 int calc_weight_function(void)
00595 {
00596   int j, mid;
00597   double arg, Dmid, dt, Sum = 0.0;
00598   int npts = parms.WtLen;  
00599   char *filter = parms.WtFunc;
00600   double fwhm = parms.WtParm;
00601   double *Wt;
00602   
00603   if(( parms.Wt = (double *) malloc(sizeof(double)*npts))==NULL)
00604   {
00605     sprintf(parms.Msg,"%sError: Allocating Memory for weighting function",parms.Msg);
00606     return  (parms.status = kMyMod_MallocErr);  
00607   }
00608   Wt = parms.Wt;
00609   Dmid = 0.5*(npts - 1);
00610   mid = Dmid;
00611   
00612   if(!strncasecmp(filter,"box",3))
00613   {
00614     for(j = 0; j < npts; j++) 
00615       Wt[j] = 1.0;
00616   }
00617   else  if(!strncasecmp(filter,"tri",3))
00618   {
00619     dt = 1.0 / fwhm;
00620     for(j = 0; j <= mid; j++) 
00621     {
00622       arg = 1.0 - (Dmid - j) * dt;
00623       Wt[j] = Wt[npts-j-1] = (arg >= 0.0) ? arg: 0.0;
00624     }
00625   }
00626   else  if(!strncasecmp(filter,"sinc2",5))
00627   {
00628     dt = 2 * 0.442946 / fwhm;
00629     for(j = 0; j <= mid; j++) 
00630     {
00631       arg = (Dmid - j) * dt;
00632       if(arg != 0.0)
00633       {
00634     Wt[j] = Wt[npts - j-1] = sin(arg)/arg;
00635     Wt[npts-j-1] *= Wt[j];
00636     Wt[j] *= Wt[j];
00637       }
00638       else
00639     Wt[j] = Wt[npts - j-1] = 1.0;       
00640     }
00641   }
00642   else  if(!strncasecmp(filter,"sinc",4))
00643   {
00644     dt = 2 * 0.603355 / fwhm;
00645     for(j = 0; j <= mid; j++) 
00646     {
00647       arg = (mid - j) * dt;
00648       if(arg != 0.0)
00649     Wt[j] = Wt[npts - j-1] = sin(arg)/arg;
00650       else
00651     Wt[j] = Wt[npts - j-1] = 1.0;       
00652     }
00653   }
00654   else if(!strncasecmp(filter,"hath",4))
00655   {
00656     dt = exp(-2.0);
00657     for(j = 0; j < npts; j++) 
00658     {
00659       arg = (j - mid)/(fwhm);
00660       arg *= -0.5 * arg;
00661       Wt[j] = exp(arg) - dt*(3.0 + arg);
00662     }
00663   }
00664   else if(!strncasecmp(filter,"gauss",5))
00665   {
00666     dt = 2 * sqrt (log (2.0)) / fwhm;
00667     for(j = 0; j <= mid; j++) 
00668     {
00669       arg = (mid - j) * dt;
00670       arg *= -arg;
00671       Wt[j] = Wt[npts-j-1] = exp (arg);
00672     }
00673   }
00674   else
00675   {
00676     sprintf(parms.Msg,"%sUnknown Weight Function: %s\n",parms.Msg,filter);
00677     return (parms.status = kMyMod_WrongType);
00678   }
00679 /*
00680   Sum = 0.0;
00681   for(j = 0; j < npts; j++)
00682     Sum += Wt[j];
00683 
00684   if (Sum == 0.0)
00685   {
00686     sprintf(parms.Msg,"%sError Calculating Weight Function\n",parms.Msg);
00687     return (parms.status = kMyMod_InitErr);
00688   }
00689   for(j = 0; j < npts; j++)
00690     Wt[j] /= Sum;
00691 */
00692   return kMyMod_Success;     // exit calc_weight_func normally
00693 }
00694 
00695 
00696 /* 
00697 ******************************************************************** 
00698 * 
00699 *  Functions involving setting up output arrays and checking outDB
00700 * 
00701 ******************************************************************** 
00702 */
00703 // SETUP_OUTDB
00704 int  setup_outdb(OutImgs **outimg, DRMS_RecordSet_t **outdb)
00705 {
00706   int status;
00707 
00708   strcpy(parms.SigMsg," ");
00709   strcpy(parms.Msg," ");
00710   (*outdb) = drms_create_records(drms_env, parms.Nmaps,
00711                  parms.oser, DRMS_PERMANENT,&status);
00712 
00713   if(!(*outdb) || (*outdb)->n != parms.Nmaps)
00714   {
00715     sprintf(parms.Msg,"%sError creating records in series %s\n",
00716         parms.Msg,parms.oser);
00717     return (parms.status = kMyMod_Missing);
00718   }
00719   
00720   if((status = check_outdb((*outdb)->records[0])))
00721     return status;
00722 
00723   if((status = init_outimgs(outimg)))
00724     return status;
00725 //   drms_close_records((*outdb), DRMS_FREE_RECORD);
00726   V_print("OK\n");
00727   return  kMyMod_Success;     // exit setup_outdb normally
00728 }
00729 
00730 // CHECK_OUTDB
00731 
00732 int    check_outdb(DRMS_Record_t *rec)
00733 {
00734   int i, done, segcnt;
00735   int status = 0;
00736   char buff[200];
00737   TIME ttime;
00738   DRMS_Segment_t *record_segment;
00739   DRMS_Keyword_t *keywd;
00740   static struct drmean_keys 
00741   {
00742     int type;
00743     char name[20];
00744   } keys[]= {
00745     {DRMS_TYPE_TIME,    "TSTART"},
00746     {DRMS_TYPE_TIME,    "TEND"},
00747     {DRMS_TYPE_TIME,    "T_REC_STEP"},
00748     {DRMS_TYPE_TIME,    "T_REC"},
00749     {DRMS_TYPE_DOUBLE,  "A0"},
00750     {DRMS_TYPE_DOUBLE,  "A2"},         
00751     {DRMS_TYPE_DOUBLE,  "A4"},         
00752     {DRMS_TYPE_DOUBLE,  "MERI_V"},     
00753     {DRMS_TYPE_STRING,  "WTFUNC"}, 
00754     {DRMS_TYPE_DOUBLE,  "WTPARM"},     
00755     {DRMS_TYPE_INT,     "WTLEN"},      
00756     {DRMS_TYPE_STRING,  "PROJNAME"},  
00757     {DRMS_TYPE_DOUBLE,  "RSUN"},       
00758     {DRMS_TYPE_DOUBLE,  "DIST"},       
00759     {DRMS_TYPE_DOUBLE,  "ANG_RAD"},    
00760     {DRMS_TYPE_DOUBLE,  "PANGLE"},     
00761     {DRMS_TYPE_DOUBLE,  "CROTA2"},      
00762     {DRMS_TYPE_DOUBLE,  "CRLN"},       
00763     {DRMS_TYPE_DOUBLE,  "CRLT"},       
00764     {DRMS_TYPE_STRING,  "CTYPE1"},     
00765     {DRMS_TYPE_STRING,  "CTYPE2"},     
00766     {DRMS_TYPE_DOUBLE,  "CRPIX1"},     
00767     {DRMS_TYPE_DOUBLE,  "CRPIX2"},     
00768     {DRMS_TYPE_DOUBLE,  "CRVAL1"},     
00769     {DRMS_TYPE_DOUBLE,  "CRVAL2"},     
00770     {DRMS_TYPE_DOUBLE,  "CDELT1"},     
00771     {DRMS_TYPE_DOUBLE,  "CDELT2"},   
00772     {DRMS_TYPE_RAW,     "END "}
00773   };
00774   
00775   // First check to ensure that there are two segments (one for data
00776   // and the other for weights and that the dimensions are consistent
00777   
00778   if((segcnt = drms_record_numsegments (rec)) != 2)
00779     {
00780       sprintf(parms.Msg,"%sError: Bad Output DB - need 2 segments\n",parms.Msg);
00781       return (parms.status = kMyMod_Missing);
00782     }
00783   for(i = 0; i < segcnt; i++) 
00784   {
00785     record_segment = drms_segment_lookupnum (rec, i);
00786     if(record_segment->info->naxis != 2) 
00787     {
00788       sprintf(parms.Msg,"%sError: Bad Output DB Seg %d - not 2D\n",parms.Msg,i);
00789       return (parms.status = kMyMod_Missing);
00790     }
00791     if(record_segment->axis[0] != parms.cols ||
00792     record_segment->axis[1] != parms.rows ) 
00793     {
00794       parms.cols = record_segment->axis[0];
00795       parms.rows = record_segment->axis[1];
00796       parms.size = parms.rows * parms.cols;
00797       sprintf(buff,"Specified dimensions do not match segment definition. Set to %d,%d\n",
00798           parms.cols, parms.rows);
00799       V_print(buff);
00800     }
00801     if( isnan(parms.CRPIX1))
00802       parms.CRPIX1 = 0.5*parms.cols;
00803     if( isnan(parms.CRPIX2))
00804       parms.CRPIX2 = 0.5*parms.rows;
00805 
00806   }
00807   
00808   done = i = 0;
00809   while(strcasecmp(keys[i].name,"END ") && !done)
00810   {
00811     if(!(keywd = drms_keyword_lookup (rec, keys[i].name, 0))) 
00812     {
00813       sprintf(parms.Msg,"%sError: Bad Output DB - missing keyword %s\n",
00814           parms.Msg, keys[i].name);
00815       return (parms.status = kMyMod_Missing);
00816     }
00817     if(keywd->info->type != keys[i].type)
00818     {
00819       sprintf(parms.Msg,"%sError: Bad Output DB - keyword %s is wrong type\n",
00820           parms.Msg,keys[i].name);
00821       return (parms.status = kMyMod_WrongType);
00822     }
00823     i++;
00824   }
00825   
00826   // Check that T_REC_step matches step_t variable
00827   ttime = drms_getkey_time(rec, "T_REC_step", &status);
00828   if(status || ttime > parms.Tstep)
00829   {
00830     sprintf(parms.Msg,"%sError: STEP_T is less than T_REC_step\n",parms.Msg);
00831     return (parms.status = kMyMod_WrongType);
00832   }
00833   
00834   return  kMyMod_Success;   // exit  check_outdb normally
00835 }
00836   
00837 // INIT_OUTIMGS
00838 int    init_outimgs(OutImgs **outimg)
00839 {
00840   int i, status;
00841   char buff[200], tbuff[200];
00842   int size = parms.size;
00843   double Bz;
00844   OutImgs *out;
00845   int naxis = 2;
00846   int naxes[2] = {parms.cols, parms.rows}; 
00847 
00848   
00849   out = (OutImgs *)malloc(sizeof(OutImgs)*parms.Nmaps);
00850   status = 0;
00851   sprintf(buff,"%4dx%4d %s,%s CDELT=[%8.5f,%8.5f] CRPIX=[%8.2f,%8.2f]\n", 
00852       parms.cols,parms.rows, 
00853       parms.CTYPE1, parms.CTYPE2,
00854       parms.CDELT1, parms.CDELT2,
00855       parms.CRPIX1, parms.CRPIX2
00856     );
00857   V_print(buff);
00858   V_print("-----------------------\n");
00859 
00860   for(i =0; i < parms.Nmaps && !status; i++)
00861   {
00862     out[i].DatArray = drms_array_create(DRMS_TYPE_FLOAT, naxis, naxes, NULL, &status);
00863     out[i].DatArray->israw = 0;
00864     if(status)
00865     {
00866       sprintf(parms.Msg,"%sError: Creating output arrays",parms.Msg);
00867       return (parms.status = status);
00868     }
00869     out[i].dat = (float *)out[i].DatArray->data;
00870 
00871     out[i].WgtArray = drms_array_create(DRMS_TYPE_FLOAT, naxis, naxes, NULL, &status);
00872 //    out[i].WgtArray->israw = 0;
00873     if(status)
00874     {
00875       sprintf(parms.Msg,"%sError: Creating output arrays",parms.Msg);
00876       return (parms.status = status);
00877     }
00878     out[i].wgt = (float *)out[i].WgtArray->data;
00879     
00880     if((out[i].lon = (float *)malloc(sizeof(float)*size))==NULL)
00881       status = kMyMod_MallocErr;
00882     if((out[i].lat = (float *)malloc(sizeof(float)*size))==NULL)
00883       status = kMyMod_MallocErr;
00884     if((out[i].osun = (short *)malloc(sizeof(short)*size))==NULL)
00885       status = kMyMod_MallocErr;
00886 
00887     if(status)
00888     {
00889       sprintf(parms.Msg,"%sError: Allocating Memory for output arrays",
00890           parms.Msg);
00891       return (parms.status = status);
00892     }
00893     out[i].timctr = parms.Tstart + i * parms.Tstep;
00894     out[i].LN = parms.CRVAL1 - (i * parms.Tstep * parms.A0);
00895     out[i].LT = parms.CRVAL2 + i * parms.Tstep * parms.Meri_V;
00896     HeliographicLocation(out[i].timctr, &(out[i].crot), &(out[i].crln), &Bz);
00897     out[i].irecdt = parms.WtLen;
00898     out[i].irecno = -1;
00899     sprint_time(tbuff,out[i].timctr,"TAI",0);
00900     sprintf(buff,"[%d] %27s[%.19s]; LN,LT=[%8.4f,%8.4f]\n",i, 
00901         parms.oser, tbuff, out[i].LN,out[i].LT);
00902     V_print(buff);
00903     
00904     if((status = calc_out_geom(out[i])))
00905     {
00906       sprintf(parms.Msg,"%sError: Calculating Geometry for output image",parms.Msg);
00907       return (parms.status = status);
00908     }
00909   }
00910   (*outimg) = out; 
00911   return kMyMod_Success;     // exit init_outimgs normally
00912 }
00913 
00914 
00915 // CALC_OUT_GEOM
00916 int calc_out_geom(OutImgs outimg)
00917 {
00918   double x, y, xrot, yrot, lon, lat;
00919   int row, col, i, status;
00920   static int first = 1;
00921   static double sin_phi, cos_phi, x0, y0;
00922 
00923   if(first)
00924   {
00925     first = 0;
00926     cos_phi = cos (parms.pangle);
00927     sin_phi = sin (parms.pangle);
00928     y0 = parms.CRPIX2;
00929     x0 = parms.CRPIX1;
00930   }
00931   
00932   status = i = 0;
00933   for(row=0; row < parms.rows; row++)
00934   {
00935     y = (row-y0)*parms.CDELT2;
00936     for(col = 0; col < parms.cols; col++ ) 
00937     {
00938       i = col + row * parms.cols;
00939       x = (col-x0)*parms.CDELT1;
00940       xrot = (x * cos_phi - y * sin_phi);
00941       yrot = (y * cos_phi + x * sin_phi);
00942       if(parms.projcode == AERIAL)
00943       {     // AERIAL: call img2sphere instead of plane2sphere.  
00944     double RHO,SLT,CLT,SIG,MU,CHI;  // Dummy Vars for passing 
00945     outimg.osun[i] = img2sphere(xrot,yrot,parms.ang_rad,outimg.LT,
00946                     outimg.LN,0.0, &RHO,&lat,&lon,&SLT,&CLT,
00947                     &SIG,&MU,&CHI);
00948       } 
00949       else 
00950       {
00951     outimg.osun[i] = plane2sphere(xrot,yrot,outimg.LT,outimg.LN,&lat,&lon,
00952                       parms.projcode);
00953       }
00954       if(isnan(lat) || isnan(lon))
00955     outimg.osun[i] = kMyMod_ValErr;
00956       outimg.dat[i] = 0.0;
00957       outimg.wgt[i] = 0.0;
00958       outimg.lat[i] = lat;
00959       outimg.lon[i++] = lon;
00960     } 
00961   }
00962   return kMyMod_Success;     // exit geom_outimgs normally
00963 }
00964 
00965 
00966 
00967 /* 
00968  ******************************************************************* 
00969  * 
00970  *  Functions involving opening and  checking input DB
00971  * 
00972  ******************************************************************** 
00973 */
00974 // OPEN_CHECK_INDB
00975 int  open_check_indb(DRMS_RecordSet_t **indb)
00976 {
00977   int status;
00978   char temp[2000];
00979   char qry_str[DRMS_MAXQUERYLEN];
00980   DRMS_Record_t *iRec;
00981   
00982   strcpy(parms.SigMsg," ");
00983   strcpy(parms.Msg," ");
00984   /*  Need to check that query returned actual records */
00985   if(get_query_str(qry_str) != kMyMod_Success)
00986     return (parms.status = kMyMod_Missing);
00987 
00988   sprintf(temp,"%s\n",qry_str);
00989   V_print(temp);
00990 
00991   (*indb) = drms_open_records(drms_env, qry_str, &status);
00992 
00993   if(status)
00994   {
00995     sprintf(parms.Msg,"%sError %d: in opening input data records\n",
00996         parms.Msg, status);
00997     return (parms.status = status);
00998   }
00999   if((*indb)->n == 0)
01000   {
01001     sprintf(parms.Msg,"%sNo input data records -Check date/time\n",parms.Msg);
01002     return (parms.status = kMyMod_Missing);
01003   }
01004   iRec = (*indb)->records[0];
01005   if(iRec == NULL) 
01006     printf("iRec is Null\n");
01007   //  Now, run through and check the existence of required keywords
01008   if((status = check_indb(iRec)) != kMyMod_Success)
01009      return status;
01010      
01011   V_print("OK\n");
01012   return kMyMod_Success;     // exit open_check_indb normally
01013 }
01014 
01015 // GET_QUERY_STR
01016 int    get_query_str(char *qry)
01017 {
01018   const char *test;
01019   char *p;
01020   char tbuf1[100],tbuf2[100];
01021 
01022   test = cmdparams_get_str(&cmdparams, "in", NULL);
01023   if (strchr(test,'[') != NULL)
01024     strcpy(qry,test);
01025   else
01026   {
01027     sprint_time(tbuf1, parms.Tstart - (parms.WtLen)/2 - 60, "TAI", 0);
01028     sprint_time(tbuf2, parms.Tend + (parms.WtLen)/2 + 60, "TAI", 0);
01029     sprintf(qry,"%s[%s - %s]",test,tbuf1,tbuf2);
01030   }
01031   return kMyMod_Success;     // exit open_check_indb normally
01032 }
01033 
01034 // CHECK_INDB    checks the structure of input data series
01035 int    check_indb(DRMS_Record_t *rec)
01036 {
01037   int i, done, segcnt;
01038   int status = 0;
01039   DRMS_Keyword_t *keywd;
01040   static struct drmean_keys
01041   {
01042     int type;
01043     char name[20];
01044   } keys[]= {
01045 //    {DRMS_TYPE_DOUBLE,  "R_SUN"},       
01046     {DRMS_TYPE_DOUBLE,  "DSUN_OBS"},       
01047 //    {DRMS_TYPE_DOUBLE,  "ANG_RAD"},    
01048 //    {DRMS_TYPE_DOUBLE,  "PANGLE"},     
01049     {DRMS_TYPE_DOUBLE,  "CROTA2"},      
01050     {DRMS_TYPE_DOUBLE,  "CRLN_OBS"},       
01051     {DRMS_TYPE_DOUBLE,  "CRLT_OBS"},       
01052     {DRMS_TYPE_STRING,  "CTYPE1"},     
01053     {DRMS_TYPE_STRING,  "CTYPE2"},     
01054     {DRMS_TYPE_DOUBLE,  "CRPIX1"},     
01055     {DRMS_TYPE_DOUBLE,  "CRPIX2"},     
01056     {DRMS_TYPE_DOUBLE,  "CRVAL1"},     
01057     {DRMS_TYPE_DOUBLE,  "CRVAL2"},     
01058     {DRMS_TYPE_DOUBLE,  "CDELT1"},     
01059     {DRMS_TYPE_DOUBLE,  "CDELT2"},   
01060     {DRMS_TYPE_RAW,     "END "}
01061   };
01062   
01063   done = i = 0;
01064   while(strcasecmp(keys[i].name,"END "))
01065   {
01066     if(!(keywd = drms_keyword_lookup (rec, keys[i].name, 0))) 
01067     {
01068       sprintf(parms.Msg,"%sError: Bad Input DB - missing keyword: %s\n",
01069           parms.Msg, keys[i].name);
01070       return (parms.status = kMyMod_Missing);
01071     }
01072 //    if(keywd->info->type != keys[i].type)
01073 //   {
01074 //      sprintf(parms.Msg,"%sError: Bad Input DB - keyword %s is wrong type\n",
01075 //        parms.Msg, keys[i].name);
01076 //      return (parms.status = kMyMod_WrongType);
01077 //    }
01078     i++;
01079   }
01080 
01081   if(!(keywd = drms_keyword_lookup (rec, "R_SUN", 0))) 
01082   {
01083     if(!(keywd = drms_keyword_lookup (rec, "RSUN_OBS", 0))) 
01084     {
01085       sprintf(parms.Msg,"%sError: Input DB - missing R_SUN && RSUN_OBS one is required\n",
01086           parms.Msg);
01087       return (parms.status = kMyMod_Missing);
01088     }
01089   }
01090   
01091   return kMyMod_Success;     // exit check_indb normally
01092 }
01093 
01094 
01095 
01096 /* 
01097  ******************************************************************* 
01098  * 
01099  *  Derotation & summation  Functions 
01100  * 
01101  ******************************************************************** 
01102 */
01103 // ADDIN
01104 
01105 int    addin(DRMS_RecordSet_t *inDB, OutImgs *outimgs)
01106 {
01107   int i, nInRecs, status;
01108   int dt, cols, rows, m, n;
01109   int WtLen2 = parms.WtLen / 2;
01110   char repstr[DRMS_MAXQUERYLEN], temp[1000], ctype1[20], ctype2[20];
01111   float *in_data;
01112   double Bo, DT, Lo, SB2, lat, lon, phi, Rsun, xc, xx, yc, yy, TEST; 
01113   TIME timg;
01114   DRMS_Array_t *data_array = NULL;
01115   DRMS_Segment_t *dseg;
01116   DRMS_Keyword_t *keywd;
01117   DRMS_Record_t    *inRec = NULL;
01118   
01119   strcpy(parms.SigMsg," ");
01120   strcpy(parms.Msg," ");
01121 
01122   nInRecs = inDB->n;
01123   for(i =0; i < nInRecs; i++)
01124   {
01125     inRec = inDB->records[i];
01126     drms_sprint_rec_query(repstr,inRec);
01127     if((readin(inRec, &timg, &Bo, &Lo, &Rsun, &phi, 
01128            &xc, &yc, ctype1, ctype2) == kMyMod_Success))
01129     {
01130       sprintf(temp,"%5.1f %5.1f %5.1f %5.1f %5.1f ",Lo,Bo,xc,yc,Rsun);
01131       strcat(repstr,temp);
01132       dseg = drms_segment_lookupnum (inRec, 0);
01133       if(data_array != NULL)
01134       {
01135     drms_free_array(data_array);
01136     data_array = NULL;
01137       }
01138       data_array = drms_segment_read (dseg, DRMS_TYPE_FLOAT, &status);
01139       in_data = (float *)data_array->data;
01140       cols = data_array->axis[0];
01141       rows = data_array->axis[1];
01142       for(m = 0; m < parms.Nmaps; m++) 
01143       {
01144     dt = DT = timg - outimgs[m].timctr;
01145     if(abs(DT) <= WtLen2)
01146     {
01147       sprintf(temp,"+%d,",m);
01148       strcat(repstr,temp);
01149       if(abs(DT) < abs(outimgs[m].irecdt))
01150       {
01151         outimgs[m].irecdt = DT;
01152         outimgs[m].irecno = i;
01153       }
01154       for(n = 0; n < parms.size; n++) 
01155       {
01156         if(!outimgs[m].osun[n])         
01157         {
01158           lat = outimgs[m].lat[n] + parms.Meri_V * DT;
01159           lon = outimgs[m].lon[n] + parms.A0 * DT;
01160           if(parms.A2 != 0 || parms.A4 != 0)
01161           {
01162         SB2 = sin(lat) * sin(lat);
01163         lon -= (SB2 * parms.A2 + parms.A4 * SB2 * SB2) * DT;
01164           }
01165           sphere2img (lat, lon, Bo, Lo, &xx, &yy, xc, yc, 
01166               Rsun, phi, 0.0, 1.0, 0.0, 0.0); 
01167           if(!isnan(TEST = ccint2(in_data,cols,rows,xx,yy)))
01168           {
01169         outimgs[m].dat[n] += TEST * parms.Wt[dt + WtLen2];
01170         outimgs[m].wgt[n] += parms.Wt[dt + WtLen2];
01171           }
01172         }
01173       }
01174     }
01175       }
01176     }
01177     strcat(repstr,"\n");
01178     V_print(repstr);
01179   }
01180 
01181   for(m = 0; m < parms.Nmaps; m++) 
01182   {
01183     for(n = 0; n < parms.size; n++) 
01184     {
01185       if(outimgs[m].wgt[n] > 0.0)
01186     outimgs[m].dat[n] = outimgs[m].dat[n] / outimgs[m].wgt[n];
01187       else
01188     outimgs[m].dat[n] = 0.0/0.0;
01189     }
01190   }    
01191   return kMyMod_Success;     // exit addin normally
01192 }
01193 
01194 
01195 // READIN reads input image parameters
01196 
01197 int    readin(DRMS_Record_t *irec, TIME *timg, double *Bo, double *Lo, 
01198           double *Rsun, double *phi,
01199               double *xc, double *yc, char ctype1[], char ctype2[])
01200 {
01201   int status;
01202   double dsun, cdelt, rang;
01203   double R2D = 180.0 / M_PI;
01204 //  char temp1[100], tbuff[100];
01205   
01206   if(irec == NULL)
01207     return kMyMod_Missing;
01208 
01209   if((status = rec_getcheck_time(irec, "T_REC",timg)))
01210     if((status = rec_getcheck_time(irec, "T_OBS",timg)))
01211       return status;
01212   
01213   if(status = rec_getcheck_str(irec,"CTYPE1",ctype1))
01214     return status;
01215   if(status = rec_getcheck_str(irec,"CTYPE2",ctype2))
01216     return status;
01217   if(status = rec_getcheck_double(irec,"CRVAL1",Lo))
01218     return status;
01219   if(status = rec_getcheck_double(irec,"CRVAL2",Bo))
01220     return status;
01221   if(status = rec_getcheck_double(irec,"CRPIX1",xc))
01222     return status;
01223   if(status = rec_getcheck_double(irec,"CRPIX2",yc))
01224     return status;
01225   if(status = rec_getcheck_double(irec,"CROTA2",phi) ||
01226      drms_ismissing_double(*phi))
01227     (*phi) = 0.0;
01228   (*phi) = -1*(*phi) / R2D;    // Convert from degrees to Radians
01229 
01230   if(
01231      drms_ismissing_double(*Lo) || drms_ismissing_double(*Bo) ||
01232      drms_ismissing_double(*xc) || drms_ismissing_double(*yc) ||
01233      drms_ismissing_string(ctype1) || drms_ismissing_string(ctype2) ||
01234      drms_ismissing_time(*timg) )
01235     return kMyMod_Missing;
01236 
01237 /* To get the Solar Radius - in pixels instead read DSUN_OBS and
01238     convert.
01239   if(status = rec_getcheck_double(irec,"RSUN_OBS",Rsun))
01240     if(status = rec_getcheck_double(irec,"R_SUN",Rsun))
01241 */
01242   if(status = rec_getcheck_double(irec,"DSUN_OBS",&dsun))
01243     return status;
01244   rang = 3600.*57.295778 * asin(6.96e8/dsun);
01245   if(status = rec_getcheck_double(irec,"CDELT1",&cdelt))
01246     return status;
01247   (*Rsun) = rang / cdelt;
01248   
01249   return kMyMod_Success;     // exit readin normally
01250 }
01251 
01252 
01253 // REC_GETCHECK_DOUBLE
01254 int rec_getcheck_double(DRMS_Record_t *rec, char kword[], double *out)
01255 {
01256   int stat;
01257   char ctmp[1000];
01258   double temp;
01259 
01260   temp = drms_getkey_double(rec, kword,&stat);
01261   if(stat)
01262   {  
01263     drms_sprint_rec_query(ctmp,rec);
01264     sprintf(ctmp,"%s:  MISSING KEYWORD %s ",ctmp, kword);
01265     V_print(ctmp);
01266     return stat;
01267   }
01268   (*out) = temp;
01269   return kMyMod_Success;     // exit rec_getcheck_double normally
01270   
01271 }
01272 
01273 // REC_GETCHECK_INT
01274 int rec_getcheck_int(DRMS_Record_t *rec, char kword[], int *out)
01275 {
01276   int stat;
01277   char ctmp[1000];
01278   int temp;
01279 
01280   temp = drms_getkey_int(rec, kword,&stat);
01281   if(stat)
01282   {
01283     drms_sprint_rec_query(ctmp,rec);
01284     sprintf(ctmp,"%s:  MISSING KEYWORD %s ",ctmp, kword);
01285     V_print(ctmp);
01286     return stat;
01287   }
01288   (*out) = temp;
01289   return kMyMod_Success;     // exit rec_getcheck_int normally
01290   
01291 }
01292 // REC_GETCHECK_STR
01293 int rec_getcheck_str(DRMS_Record_t *rec, char kword[], char out[])
01294 {
01295   int stat;
01296   char ctmp[1000];
01297   const char *temp;
01298   temp = drms_getkey_string(rec,kword,&stat);
01299   if(stat)
01300   {
01301     drms_sprint_rec_query(ctmp,rec);
01302     sprintf(ctmp,"%s:  MISSING KEYWORD %s ",ctmp, kword);
01303     V_print(ctmp);
01304     return stat;
01305   }
01306   strcpy(out,temp);
01307   return  kMyMod_Success;     // exit rec_getcheck_str normally
01308   
01309 }
01310 
01311 // REC_GETCHECK_TIME
01312 int rec_getcheck_time(DRMS_Record_t *rec, char kword[], TIME *out)
01313 {
01314   int stat;
01315   char ctmp[1000];
01316   TIME temp;
01317 
01318   temp = drms_getkey_time(rec, kword,&stat);
01319   if(stat)
01320   {
01321     drms_sprint_rec_query(ctmp,rec);
01322     sprintf(ctmp,"%s:  MISSING KEYWORD %s ",ctmp, kword);
01323     V_print(ctmp);
01324     return stat;
01325   }
01326   (*out) = temp;
01327   return kMyMod_Success;     // exit rec_getcheck_time normally
01328   
01329 }
01330 
01331 /* 
01332  ******************************************************************* 
01333  * 
01334  *  Checking & Writing Out Img Functions
01335  * 
01336  ******************************************************************** 
01337 */
01338 // WRITE_OUT_DBS
01339 int    write_out_dbs(OutImgs *outimgs, DRMS_RecordSet_t *outDB, DRMS_RecordSet_t *inDB)
01340 {
01341   int i, nrecs, status;
01342   char cunit[20];
01343   char buff[200];
01344   DRMS_Record_t *orec, *irec;
01345   DRMS_Segment_t *oseg;
01346   double SCL, R2D, R2AS;
01347 //  DRMS_RecordSet_t *outdb;
01348 
01349   strcpy(parms.SigMsg," ");
01350   strcpy(parms.Msg," ");
01351 
01352 //  outdb = drms_create_records(drms_env, parms.Nmaps,
01353 //               parms.oser, DRMS_PERMANENT,&status);
01354   R2D = 180.0 / M_PI;
01355   R2AS = 3600.0 *180.0 / M_PI;
01356 
01357 // Have to re-do SCL factor to be consistent with input,
01358 // comp with check_cparms... the key CDELT2 is off!!
01359 
01360   nrecs = parms.Nmaps;
01361   if(parms.projcode == AERIAL)
01362   {
01363     strcpy(cunit,"arcsec");
01364     SCL = R2D * 3600.00;
01365   }
01366   else
01367   {
01368      strcpy(cunit,"deg");
01369      SCL = R2D;
01370   }
01371 
01372   for(i = 0; i < nrecs; i++)
01373   {
01374     orec = outDB->records[i];
01375 //    orec = drms_create_record(drms_env,parms.oser,DRMS_PERMANENT, &status);
01376     irec = inDB->records[outimgs[i].irecno];
01377     drms_copykeys(orec,irec,0, kDRMS_KeyClass_Explicit ); 
01378     drms_setkey_time(orec,"T_REC",outimgs[i].timctr);
01379     drms_setkey_double(orec,"A0",parms.A0*1.0e6);
01380     drms_setkey_double(orec,"A2",parms.A2*1.0e6);
01381     drms_setkey_double(orec,"A4",parms.A4*1.0e6);
01382     drms_setkey_double(orec,"MERI_V",parms.Meri_V);
01383     drms_setkey_string(orec,"WTFUNC",parms.WtFunc);
01384     drms_setkey_double(orec,"WTPARM",parms.WtParm);
01385     drms_setkey_int(orec,"WTLEN",parms.WtLen);
01386     drms_setkey_string(orec,"PROJNAME",parms.projname);
01387     drms_setkey_double(orec,"RSUN",parms.rsun);
01388     drms_setkey_double(orec,"DIST",parms.dist);
01389     drms_setkey_double(orec,"ANG_RAD",parms.ang_rad*R2AS);
01390     drms_setkey_double(orec,"PANGLE_RAD",parms.pangle*R2D);
01391     drms_setkey_double(orec,"CROTA2",parms.CROTA2*R2D);
01392     drms_setkey_double(orec,"CRLN",outimgs[i].LN*SCL);
01393     drms_setkey_double(orec,"CRLT",outimgs[i].LT*SCL);
01394     drms_setkey_string(orec,"CTYPE1",parms.CTYPE1);
01395     drms_setkey_string(orec,"CTYPE2",parms.CTYPE2);
01396     drms_setkey_double(orec,"CRPIX1",parms.CRPIX1);
01397     drms_setkey_double(orec,"CRPIX2",parms.CRPIX2);
01398     drms_setkey_double(orec,"CRVAL1",outimgs[i].LN*SCL);
01399     drms_setkey_double(orec,"CRVAL2",outimgs[i].LT*SCL);
01400     drms_setkey_double(orec,"CDELT1",parms.CDELT1*SCL);
01401     drms_setkey_double(orec,"CDELT2",parms.CDELT2*SCL);
01402     drms_setkey_string(orec,"CUNIT1",cunit);
01403     drms_setkey_string(orec,"CUNIT2",cunit);
01404     drms_setkey_double(orec,"CADENCE",parms.Tstep);
01405     drms_setkey_double(orec,"CRLN_OBS",outimgs[i].crln);
01406     drms_setkey_int(orec,"CAR_ROT",outimgs[i].crot);
01407 
01408     if((parms.status = set_bscale_bzero(&outimgs[i])) != kMyMod_Success)
01409     {
01410       sprintf(parms.Msg,"%sError: setting BSCALE & BZERO\n",parms.Msg);
01411       return (parms.status = kMyMod_InitErr);
01412     } 
01413     sprintf(buff,"[%02d] Data:   Max=%9.4f, Min=%9.4f, Bscale=%9.4f, Bzero=%9.4f\n",
01414         i, outimgs[i].Dmax, outimgs[i].Dmin, 
01415         outimgs[i].DatArray->bscale, outimgs[i].DatArray->bzero);
01416     V_print(buff);
01417     sprintf(buff,"[%02d] Weight: Max=%9f, Min=%9f, Bscale=%9f, Bzero=%9f\n",
01418         i, outimgs[i].Wmax, outimgs[i].Wmin, 
01419         outimgs[i].WgtArray->bscale, outimgs[i].WgtArray->bzero);
01420     V_print(buff);
01421 
01422     oseg = drms_segment_lookupnum (orec, 0);
01423     oseg->bzero = outimgs[i].DatArray->bzero;
01424     oseg->bscale = outimgs[i].DatArray->bscale;
01425     outimgs[i].DatArray->parent_segment = oseg;
01426     outimgs[i].DatArray->israw = 0;
01427     drms_segment_write(oseg,outimgs[i].DatArray,0);
01428     drms_free_array(outimgs[i].DatArray);
01429 
01430     oseg = drms_segment_lookupnum (orec, 1);
01431     oseg->bzero = outimgs[i].WgtArray->bzero;
01432     oseg->bscale = outimgs[i].WgtArray->bscale;
01433     outimgs[i].WgtArray->parent_segment = oseg;
01434     outimgs[i].WgtArray->israw = 0;
01435     drms_segment_write(oseg,outimgs[i].WgtArray,0);
01436     drms_free_array(outimgs[i].WgtArray);
01437 
01438   }
01439 
01440   drms_close_records(outDB, DRMS_INSERT_RECORD);
01441   status = free_outimgs(outimgs);
01442 
01443   V_print("OK\n");
01444   return kMyMod_Success;     // exit write_out_dbs normally
01445 }
01446 
01447 // SET_BSCALE_BZERO
01448 int    set_bscale_bzero(OutImgs *outimg)
01449 {
01450   int i;
01451   double TAPERNG = 3.0e4;
01452   char buff[200];
01453 
01454   outimg->DatArray->bscale = 1.0;
01455   outimg->DatArray->bzero = 0.0;
01456   outimg->WgtArray->bscale = 0.001;
01457   outimg->WgtArray->bzero = 0.0;
01458   outimg->Dmax = SHRT_MIN;
01459   outimg->Dmin = SHRT_MAX;
01460   outimg->Wmax = SHRT_MIN;
01461   outimg->Wmin = SHRT_MAX;
01462 
01463   for (i = 0; i < parms.size; i++)
01464   {
01465     if(!isnan(outimg->dat[i]))
01466     {
01467       if(outimg->Dmax < outimg->dat[i])
01468     outimg->Dmax = outimg->dat[i];
01469       if(outimg->Dmin > outimg->dat[i])
01470     outimg->Dmin = outimg->dat[i];
01471     }
01472     if(!isnan(outimg->wgt[i]))
01473     {
01474       if(outimg->Wmax < outimg->wgt[i])
01475     outimg->Wmax = outimg->wgt[i];
01476       if(outimg->Wmin > outimg->wgt[i])
01477     outimg->Wmin = outimg->wgt[i];
01478     }
01479   }
01480 
01481   if (outimg->Wmin < 0.0)
01482   { 
01483     outimg->WgtArray->bscale = (outimg->Wmax - outimg->Wmin)/TAPERNG;
01484     outimg->WgtArray->bzero = outimg->Wmin;
01485   }
01486   else
01487   {
01488     outimg->WgtArray->bscale = (outimg->Wmax)/TAPERNG;
01489     outimg->WgtArray->bzero = 0.0;
01490   }
01491   
01492   if(outimg->Dmax >= TAPERNG || outimg->Dmin <= -1*TAPERNG || 
01493      abs(outimg->Dmax - outimg->Dmin) < 0.01*TAPERNG)
01494   {
01495     outimg->DatArray->bscale = (outimg->Dmax - outimg->Dmin)/TAPERNG/2.0;
01496     outimg->DatArray->bzero = 0.5*(outimg->Dmax + outimg->Dmin);
01497     if(abs(outimg->DatArray->bzero) < 0.08*(outimg->DatArray->bscale)*TAPERNG)
01498       outimg->DatArray->bzero = 0.0;
01499   }
01500 
01501   return kMyMod_Success;   // Exit set_bscale_bzero normally - setvals
01502 }
01503 
01504 
01505 
01506 /* 
01507  ******************************************************************* 
01508  * 
01509  *  Misc Utility Functions 
01510  * 
01511  ******************************************************************** 
01512 */
01513 
01514 
01515 // WCS_PARM_UNIT
01516 double wcs_parm_unit(double X, char *unit)
01517 {
01518   static double d2r =  M_PI / 180.0;
01519   static double m2r =  M_PI / 180.0 / 60.0;
01520   static double a2r =  M_PI / 180.0 / 3600.0;
01521   
01522   if(!strcmp(unit,"arcsec")) 
01523     return X * a2r;
01524   else if(!strcmp(unit,"arcmin"))
01525     return X * m2r;
01526   else if(!strcmp(unit,"deg")) 
01527     return X * d2r;
01528   else if(!strcmp(unit,"rad")) 
01529     return X;
01530   
01531   return 0.0/0.0;
01532 }
01533 
01534 // FREE_OUTIMGS
01535 int    free_outimgs(OutImgs *outimg)
01536 {
01537   int i;
01538 
01539   for(i =0; i < parms.Nmaps; i++)
01540   {
01541     free(outimg[i].lon);
01542     free(outimg[i].lat);
01543     free(outimg[i].osun);
01544   }
01545   free(outimg);
01546   return kMyMod_Success;     // exit free_outimgs normally
01547 }
01548 
01549 
01550 // IF_ERR_PRINT
01551 int If_Err_Print(void)
01552 {
01553   if(parms.status)
01554   {
01555     fprintf(stderr, "Error %d: %s\n",parms.status,parms.Msg);
01556     fflush(stderr);
01557   }
01558   return parms.status;
01559 }
01560 
01561 //   V_PRINT
01562 void   V_print(char *str)
01563 {
01564   if (!parms.verbose && !(parms.verbose = cmdparams_isflagset(&cmdparams,"v")))
01565     return;
01566   printf("%s",str);
01567   fflush(stdout);
01568   
01569 }
01570 
01571 
01572 
01573 /*
01574  *  Revision History
01575  *  v 0.0  09.12.07     John Beck      created this file
01576  *
01577  */
01578 
01579
proj/maps_avgs/apps/derot_mean.c
