JSOC_Documentation: proj/util/apps/hg

00001 
00117 #include "jsoc.h"
00118 #include "jsoc_main.h"
00119 #include "astro.h"
00120 #include "fstats.h"
00121 #include "atoinc.h"
00122 #include <math.h>
00123 
00124 // void HeliographicLocation(TIME t, int *crot, double *L, double *B);
00125 // TIME HeliographicTime(int crot, double L);
00126 
00127 char *module_name = "hg_patch";
00128 
00129 #define DIE(msg) {fprintf(stderr,"%s  Status=%d\n",msg, status); return(status?status:1);}
00130 #define DIE2(msg,val) {fprintf(stderr,"%s %s,  Status=%d\n",msg, val,  status); return(status?status:1);}
00131 #define TEST_PARAM(param) {if (status) DIE2("Required keyword missing: ", param);}
00132 
00133 #define     Deg2Rad    (M_PI/180.0)
00134 #define     Rad2arcsec (3600.0/Deg2Rad)
00135 #define     Rad2Deg    (180.0/M_PI)
00136 
00137 ModuleArgs_t module_args[] =
00138   {
00139     {ARG_STRING, "in", "NOTSPECIFIED", "input series. e.g 'mdi.fd_M_96m_lev18'"},
00140     {ARG_STRING, "out", "NOTSPECIFIED", "output seies. e.g. 'su_bala.extractreg'"},
00141     {ARG_STRING, "log", "NOTSPECIFIED", "output log file of records made"},
00142     {ARG_STRING, "requestid", "NOTSPECIFIED", "RequestID if hg_patch call originated in an export request."},
00143     {ARG_INT, "car_rot", "-1", "Carrington Rotation when the region crosses CM"},
00144     {ARG_FLOAT, "width", "0", "width of box in degrees of longitude"},
00145     {ARG_FLOAT, "height", "0", "height of box in degrees of latitude when it corsses CM"},
00146     {ARG_STRING, "boxunits", "NOTSPECIFIED", "units of patch, 'pixels', 'arcsecs', or 'degrees'"},
00147     {ARG_TIME, "t_start", "JD_0", "Start time, defaults to time at 90E"},
00148     {ARG_TIME, "t_stop", "JD_0", "End time, defauolts to 90W"},
00149     {ARG_TIME, "cadence", "NOTSPECIFIED", "Cadence of product, defaults to input cadence."},
00150     {ARG_TIME, "t_ref", "JD_0", "Time for which x and y apply, implies ref image."},
00151     {ARG_STRING, "locunits", "NOTSPECIFIED", "Location units in 'pixels', 'arcsec', or 'degrees'"},
00152     {ARG_STRING, "where", "NOTSPECIFIED", "Additional 'where' clause if needed"},
00153     {ARG_FLOAT, "x", "0", "Location of extract box center"},
00154     {ARG_FLOAT, "y", "0", "Location of extract box center"},
00155     {ARG_FLAG, "t", "0", "Disable Carrington rate tracking"},
00156     {ARG_END}
00157   };
00158 
00159 int img2sphere (double x, double y, double ang_r, double latc, double lonc,
00160         double pa, double *rho, double *lat, double *lon, double *sinlat,
00161         double *coslat, double *sig, double *mu, double *chi);
00162 int sphere2img (double lat, double lon, double latc, double lonc,
00163         double *x, double *y, double xcenter, double ycenter,
00164         double rsun, double peff, double ecc, double chi,
00165         int xinvrt, int yinvrt);
00166 
00167 // experimental make recordset list at cadence
00168 char *get_input_recset(DRMS_Env_t *env, char *in, TIME cadence);
00169 
00170 // Macros for WCS transformations.  assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
00171 // and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
00172 // PIX_X and PIX_Y are CCD pixel addresses, WX and WY are arc-sec W and N on the Sun from disk center.
00173 // These are in units where the first pixel is 1 not 0.
00174 #define PIX_X(wx,wy) ((((wx-crvalx)*cosa + (wy-crvaly)*sina)/cdelt)+crpix1)
00175 #define PIX_Y(wx,wy) ((((wy-crvaly)*cosa - (wx-crvalx)*sina)/cdelt)+crpix2)
00176 #define WX(pix_x,pix_y) (((pix_x-crpix1)*cosa - (pix_y-crpix2)*sina)*cdelt+crvalx)
00177 #define WY(pix_x,pix_y) (((pix_y-crpix2)*cosa + (pix_x-crpix1)*sina)*cdelt+crvaly)
00178 
00179 #define BOXBAD 0
00180 #define BOXARCSEC 1
00181 #define BOXPIXELS 2
00182 #define BOXDEGREE 3
00183 
00184 #define LOCBAD 0
00185 #define LOCARCSEC 1
00186 #define LOCPIXELS 2
00187 #define LOCSTONY 3
00188 #define LOCCARR 4
00189 
00190 #define TIMES_RECSET 0   // recordset spec contains time range
00191 #define TIMES_GIVEN 1    // explicit t_start and t_stop provided
00192 #define TIMES_IMPLICIT 2 // times to be deduced for disk transit of specified box
00193 int DoIt(void)
00194 {
00195   CmdParams_t *params = &cmdparams;
00196   DRMS_RecordSet_t *inRS, *outRS;
00197   DRMS_Record_t *inRec, *outRec, *inTemplate, *outTemplate;
00198   DRMS_Segment_t *inSeg, *outSeg;
00199   DRMS_Array_t *inArray, *outArray;
00200   int i, ii, status = DRMS_SUCCESS, nrecs; 
00201   int  irec;
00202   int  OK_recs = 0;
00203   double center_x, center_y, crpix1, crpix2, x0, y0; 
00204   double rsun_ref, dsun_obs, rsun, rsun_rad, rsunpix;
00205   double crln_obs, crlt_obs;
00206   double crln_obs_rad, crlt_obs_rad;
00207   double crln_rad, crlt_rad;
00208   double pa_rad;
00209   double cdelt;
00210   double urx, ury, llx, lly;
00211   int inAxis[2];
00212   int this_car_rot;
00213   char *ctype1, *ctype2;
00214   TIME t_rec;
00215   double box_x, box_y;
00216   double crln, crlt;
00217   char outseries[DRMS_MAXSERIESNAMELEN];
00218   char inseries[DRMS_MAXSERIESNAMELEN];
00219   char inQuery[DRMS_MAXQUERYLEN];
00220   char in[DRMS_MAXQUERYLEN];
00221 
00222   const char *ingiven = params_get_str(params, "in");
00223   char *inparam;
00224   char *lbracket;
00225   char *moreQuery = NULL;
00226   const char *outparam = params_get_str(params, "out");
00227   const char *logfile = params_get_str(params, "log");
00228   const char *requestid = params_get_str(params, "requestid");
00229   const char *where = params_get_str(params, "where");
00230   TIME t_start = params_get_time(params, "t_start");
00231   TIME t_stop = params_get_time(params, "t_stop");
00232   const char *cadence_str = params_get_str(params, "cadence");
00233   double cadence;
00234   TIME t_ref = params_get_time(params, "t_ref");
00235   double width = params_get_double(params, "width");
00236   double height = params_get_double(params, "height");
00237   int car_rot = params_get_int(params, "car_rot");
00238   double x = params_get_double(params, "x");
00239   double y = params_get_double(params, "y");
00240   const char *boxunits = params_get_str(params, "boxunits");
00241   const char *locunits = params_get_str(params, "locunits");
00242   int NoTrack = params_isflagset(params, "t");
00243   TIME tNotSpecified = sscan_time("JD_0");
00244   int do_reftime = 0;  // Target specified by reftime vs Carr rotation mode.
00245   double crvalx = 0.0;
00246   double crvaly = 0.0;
00247   double crota, sina, cosa;
00248   double width_arcsec, height_arcsec;
00249   double pa, deltlong;
00250   double center_x_first, center_y_first;  // used for NoTrack option
00251   int firstimage = 1;
00252   int boxtype, loctype;
00253   char *timekeyname;
00254   TIME t_step, t_epoch;
00255   int npkeys;
00256   int times_source = TIMES_RECSET;
00257   char timebuf[20];
00258   char *in_filename = NULL;
00259 
00260   FILE *log = NULL;
00261 
00262   if (strncasecmp(locunits, "arcsec", 6) == 0) loctype = LOCARCSEC;
00263   else if (strncasecmp(locunits, "pixels", 3) == 0) loctype = LOCPIXELS;
00264   else if (strncasecmp(locunits, "stony", 5) == 0) loctype = LOCSTONY;
00265   else if (strncasecmp(locunits, "carrlong", 4) == 0) loctype = LOCCARR;
00266   else loctype = LOCBAD;
00267 
00268   if (strncasecmp(boxunits, "arcsec", 6) == 0) boxtype = BOXARCSEC;
00269   else if (strncasecmp(boxunits, "pixels", 3) == 0) boxtype = BOXPIXELS;
00270   else if (strncasecmp(boxunits, "degrees", 3) == 0) boxtype = BOXDEGREE;
00271   else loctype = BOXBAD;
00272 
00273   if (loctype == LOCCARR)
00274     {
00275     if (car_rot < 0) DIE("Carrington rotation number must be provided for locunits=carrlong");
00276     if (NoTrack) DIE("Can not use locunits=carrlong for no-tracking mode");
00277     do_reftime = 0;
00278     }
00279   else 
00280     {
00281     if (t_ref < 0) DIE("t_ref must be provided for locunits!=carrlong");
00282     do_reftime = 1;
00283     }
00284 
00285   if (loctype==LOCBAD)
00286     DIE2("Location units not detected, locunits", locunits)
00287   if (boxtype==BOXBAD)
00288     DIE2("patch dimensions not understood,", boxunits);
00289 
00290   if (strcasecmp(logfile, "NOTSPECIFIED") != 0)
00291     {
00292     log = fopen(logfile, "w");
00293     if (!log) DIE2("Can not create log file.",logfile);
00294     }
00295 
00296   if (strcasecmp(where, "NOTSPECIFIED") == 0)
00297     {
00298     where = "";
00299     }
00300   else
00301     {
00302     char wherework[4096];
00303     sprintf(wherework,"[? %s ?]", where);
00304     where = strdup(wherework);
00305 fprintf(stderr,"where is %s\n", where);
00306     }
00307 
00308   inparam = strdup(ingiven);
00309   if (strcasecmp(inparam, "NOTSPECIFIED") == 0) DIE("Input series must be specified.");
00310   lbracket = index(inparam, '[');
00311   // first, get input series names.
00312   if (lbracket)
00313     {
00314     int n = lbracket - inparam;
00315     strncpy(inseries, inparam, n);
00316     inseries[n] = '\0';
00317     if (strncmp(lbracket, "[$]", 3) == 0) // Special case, discard the explicit last-record spec to allow exports via jsoc_fetch.
00318       {
00319       moreQuery = lbracket + 3;
00320       *lbracket = '\0';
00321       lbracket = NULL;
00322       if (t_start == tNotSpecified || t_stop == tNotSpecified)
00323           times_source = TIMES_IMPLICIT;
00324       else
00325           times_source = TIMES_GIVEN;
00326       }
00327     else
00328       moreQuery = lbracket;
00329     }
00330   else
00331     {
00332     strcpy(inseries, inparam);
00333     if (t_start == tNotSpecified || t_stop == tNotSpecified)
00334         times_source = TIMES_IMPLICIT;
00335     else
00336         times_source = TIMES_GIVEN;
00337     }
00338     
00339 fprintf(stderr,"general params processed.\n");
00340 // XXXXXXXXXXX Get box location in Carrington Coords for all location types XXXXXXXXXXXXXX
00341   if (do_reftime) // Arc-sec, pixel, or Stonyhurst specification, get ref image information
00342     { // the image for ref_time is supposed to be present in the series.
00343     char t_ref_text[100];
00344     sprint_at(t_ref_text, t_ref-7200);
00345     sprintf(in, "%s[%s/4h][? QUALITY >=0 ?]", inseries, t_ref_text);
00346     inRS = drms_open_records(drms_env, in, &status); if (status || inRS->n == 0) DIE2("No input data found within 2-hours of t_ref",in);
00347     int irec, nrecs = inRS->n;
00348     TIME tdiff = 10000;
00349     for (irec=0; irec<nrecs; irec++) // find record close to t_ref
00350       {
00351       TIME newdiff;
00352       inRec = inRS->records[irec];
00353       if (drms_getkey_int(inRec,"QUALITY",NULL) < 0)
00354         continue;
00355       if ((newdiff=(fabs(drms_getkey_time(inRec,"T_OBS",NULL) - t_ref))) > tdiff)
00356         break;
00357       tdiff = newdiff;
00358       }
00359     inRec = inRS->records[irec];
00360     this_car_rot = drms_getkey_int(inRec, "CAR_ROT", &status); TEST_PARAM("CAR_ROT");
00361     crlt_obs = drms_getkey_double(inRec, "CRLT_OBS", &status); TEST_PARAM("CRLT_OBS");
00362     crln_obs = drms_getkey_double(inRec, "CRLN_OBS", &status); TEST_PARAM("CRLN_OBS");
00363     crpix1 = drms_getkey_double(inRec, "CRPIX1", &status); TEST_PARAM("CRPIX1");
00364     x0 = crpix1 - 1;
00365     crpix2 = drms_getkey_double(inRec, "CRPIX2", &status); TEST_PARAM("CRPIX2");
00366     y0 = crpix2 - 1;
00367     crvalx = drms_getkey_double(inRec, "CRVAL1", &status); TEST_PARAM("CRVAL1");
00368     crvaly = drms_getkey_double(inRec, "CRVAL2", &status); TEST_PARAM("CRVAL2");
00369     crota = drms_getkey_double(inRec, "CROTA2", &status); TEST_PARAM("CROTA2");
00370     pa = -crota;
00371     crlt_obs_rad = Deg2Rad * crlt_obs;
00372     crln_obs_rad = Deg2Rad * crln_obs;
00373     pa_rad = Deg2Rad * pa;
00374     sina = sin(crota*Deg2Rad);
00375     cosa = cos(crota*Deg2Rad);
00376     ctype1 = strdup(drms_getkey_string(inRec, "CTYPE1", &status)); TEST_PARAM("CTYPE1");
00377       if (strcmp(ctype1, "HPLN-TAN") != 0) DIE2("CTYPE1 not HPLN-TAN as required, is: ", ctype1);
00378     ctype2 = strdup(drms_getkey_string(inRec, "CTYPE2", &status)); TEST_PARAM("CTYPE2");
00379       if (strcmp(ctype2, "HPLT-TAN") != 0) DIE2("CTYPE2 not HPLT-TAN as required, is: ", ctype2);
00380     rsun_ref = drms_getkey_double(inRec, "RSUN_REF", &status);
00381       if (status) rsun_ref = 6.96e8;
00382     dsun_obs = drms_getkey_double(inRec, "DSUN_OBS", &status); TEST_PARAM("DSUN_OBS");
00383     rsun_rad = asin(rsun_ref/dsun_obs); 
00384     rsun = rsun_rad*Rad2arcsec; 
00385     cdelt = drms_getkey_double(inRec, "CDELT1", &status); TEST_PARAM("CDELT1");
00386     rsunpix = rsun/cdelt;
00387     if (loctype == LOCPIXELS || loctype==LOCARCSEC)
00388       { /* box ref location is in arcsec - find Carrington equivalent */
00389       if (loctype==LOCARCSEC)
00390         {
00391         center_x = PIX_X(x,y) - 1 - x0;
00392         center_y = PIX_Y(x,y) - 1 - y0;
00393         }
00394       else
00395         {
00396         center_x = x - 1 - x0;
00397         center_y = y - 1 - y0;
00398         }
00399       inSeg = drms_segment_lookupnum(inRec, 0);
00400       inAxis[0] = inSeg->axis[0];
00401       inAxis[1] = inSeg->axis[1];
00402       if ( img2sphere(center_x/rsunpix, center_y/rsunpix, rsun_rad, crlt_obs_rad, crln_obs_rad, pa_rad,
00403            NULL, &crlt_rad, &crln_rad, NULL, NULL, NULL, NULL, NULL) < 0)
00404           DIE("Starting location is off the solar disk.");
00405       crln = crln_rad * Rad2Deg;
00406       crlt = crlt_rad * Rad2Deg;
00407       }
00408     else /* loctype==LOCSTONY */
00409       {
00410       crlt = y;
00411       crln = crln_obs + x;
00412       }
00413     deltlong = crln - crln_obs;
00414     car_rot = this_car_rot;
00415     if (deltlong >= 360.0)
00416       { car_rot--; crln -= 360.0; }
00417     else if (deltlong < 0)
00418       { car_rot++; crln += 360.0; }
00419     drms_close_records(inRS, DRMS_FREE_RECORD);
00420 fprintf(stderr,"reftime params processed.\n");
00421     }
00422   else // Carrington specification
00423     {
00424     crln = x;
00425     crlt = y;
00426     if (crln < 0) DIE("Box longitude must be specified.");
00427     if (crlt < -990) DIE("box latitude must be specified.");
00428     }
00429   crln_rad = crln * Deg2Rad;
00430   crlt_rad = crlt * Deg2Rad;
00431 
00432   // Now we have crln, crlt, crln_rad, crlt_rad, car_rot for the target box. 
00433 
00434   // Get implied time limits from car_rot and crln
00435   t_ref = HeliographicTime(car_rot, crln);
00436   if (times_source == TIMES_IMPLICIT)
00437       {
00438       if (t_start == tNotSpecified)
00439         t_start = HeliographicTime(car_rot, crln + 90);
00440       if (t_stop == tNotSpecified)
00441         t_stop = HeliographicTime(car_rot, crln - 90);
00442       }
00443 fprintf(stderr,"box location params processed.\n");
00444 // XXXXXXXXXXXXXXXXX End of get target box location
00445 
00446 // XXXXXXXXXXXXXXXX get input seriesname and output seriesname
00447 // first, get input and output series names.
00448 
00449   inTemplate = drms_template_record(drms_env, inseries, &status);
00450   if (status || !inTemplate) DIE2("Input series can not be found: ", inseries);
00451 
00452   if (strcasecmp(outparam, "NOTSPECIFIED") == 0)
00453     {
00454     strncpy(outseries, inseries, DRMS_MAXSERIESNAMELEN);
00455     strncat(outseries, "_hgpatch", DRMS_MAXSERIESNAMELEN);
00456     }
00457   else
00458    strncpy(outseries, outparam, DRMS_MAXSERIESNAMELEN);
00459 
00460   // Now, make sure output series exists and get template record.
00461   outTemplate = drms_template_record(drms_env, outseries, &status);
00462   if (status || !outTemplate) DIE2("Output series can not be found: ", outseries);
00463 
00464   // Now find the prime time keyword name
00465   npkeys = inTemplate->seriesinfo->pidx_num;
00466   timekeyname = NULL;
00467   if (npkeys > 0)
00468     {
00469     int i;
00470     for (i=0; i<npkeys; i++)
00471         {
00472         DRMS_Keyword_t *pkey, *skey;
00473         pkey = inTemplate->seriesinfo->pidx_keywords[i];
00474         if (pkey->info->recscope > 1)
00475            pkey = drms_keyword_slotfromindex(pkey);
00476         if (pkey->info->type != DRMS_TYPE_TIME)
00477           continue;
00478     if(i > 0) DIE("Input series must have TIME keyword first, for now...");
00479         timekeyname = pkey->info->name;
00480         t_step = drms_keyword_getdouble(drms_keyword_stepfromslot(pkey), &status);
00481     if (status) DIE("problem getting t_step");
00482         t_epoch = drms_keyword_getdouble(drms_keyword_epochfromslot(pkey), &status);
00483     if (status) DIE("problem getting t_epoch");
00484         }
00485     }
00486   else
00487     DIE("Must have time prime key");
00488 
00489 fprintf(stderr,"prime key params processed.\n");
00490   // Get cadence for output data.  The output series should be slotted on a multiple of
00491   // the input series, multiple may be 1.
00492   
00493   if (strcasecmp(cadence_str, "NOTSPECIFIED") != 0)
00494     {
00495     int ratio;
00496     double initial_cadence;
00497     initial_cadence = atoinc((char *)cadence_str);
00498     ratio = round(initial_cadence/t_step);
00499     cadence = ratio * t_step;
00500     if (cadence != initial_cadence)
00501       fprintf(stderr,"Cadence rounded from %f to %f, since t_step == %f\n", initial_cadence, cadence, t_step);
00502     }
00503   else
00504     cadence = t_step;
00505 
00506   // Finally, get input recordset spec
00507   if (lbracket && times_source == TIMES_RECSET) // RecordSet query specified
00508     strncpy(in, inparam, DRMS_MAXQUERYLEN);
00509   else // need to generate query from limit times
00510     { 
00511     char t_start_text[100], t_stop_text[100], cadence_text[100];
00512     if (cadence > t_step)
00513       {
00514       sprintf(cadence_text,"@%fs", cadence);
00515       if (times_source == TIMES_IMPLICIT) // round start and stop to cadence slots for auto limits
00516         {
00517         int nsteps;
00518         nsteps = round((t_start - t_epoch)/cadence);
00519         t_start = t_epoch + nsteps * cadence;
00520         nsteps = round((t_stop - t_epoch)/cadence);
00521         t_stop = t_epoch + nsteps * cadence;
00522         }
00523       }
00524     else
00525       cadence_text[0] = '\0';
00526     // strncpy(inseries, inparam, DRMS_MAXSERIESNAMELEN);
00527     sprint_at(t_start_text, t_start);
00528     sprint_at(t_stop_text, t_stop);
00529     if (strncmp(inseries,"aia.lev1",8)==0 && t_step == 1.0 && cadence > 1.0) // special case for AIA slots
00530         {
00531         // experimental, AIA is not tseq slots so get vector of times and convert to list of records
00532         // rounded to nearest slots.  Must put this list in a temp file since may be big.
00533         sprintf(in, "%s[%s-%s]%s%s", inseries, t_start_text, t_stop_text, where , (moreQuery ? moreQuery : ""));
00534         in_filename = get_input_recset(drms_env, in, cadence);
00535         if (!in_filename) DIE("Cant make AIA cadence recordset list file");
00536         sprintf(in, "@%s", in_filename);
00537         }
00538     else // normal case
00539         sprintf(in, "%s[%s-%s%s]%s%s", inseries, t_start_text, t_stop_text, cadence_text, where, (moreQuery ? moreQuery : ""));
00540     }
00541 fprintf(stderr,"Query is %s\n",in);
00542 
00543 // XXXXXXXXXXXXX End of get input and output series information
00544 
00545 // XXXXXXXXXXXXX Now read data and extract boxes.
00546 
00547   inRS = drms_open_records(drms_env, in, &status);
00548   if (status || inRS->n == 0)
00549            DIE("No input data found");
00550   nrecs = inRS->n;
00551 fprintf(stderr,"Query finds %d records\n", nrecs);
00552 
00553   // extract patches from each record
00554   for (irec = 0; irec < nrecs; irec ++)
00555     {
00556     inRec = inRS->records[irec];
00557     if (status || !inRec) DIE("Record read failed.");
00558     TIME trec = drms_getkey_time(inRec, timekeyname, &status); TEST_PARAM(timekeyname);
00559     if (t_start != tNotSpecified && trec < t_start)
00560       {
00561 fprintf(stderr,"Rec %d, skip, trec < tstart\n",irec);
00562       continue;
00563       }
00564     if (t_stop != tNotSpecified && trec > t_stop)
00565       break;
00566     // skip records without images.
00567     if (drms_keyword_lookup(inRec, "QUALITY", 1))
00568       {
00569       int quality = drms_getkey_int(inRec, "QUALITY", &status);
00570       if (quality < 0)
00571          {
00572 fprintf(stderr,"Rec %d, skip, QUALITY=%0x\n",irec, quality);
00573          continue;
00574          }
00575       }
00576     
00577     // get coordinate information for this image
00578     this_car_rot = drms_getkey_int(inRec, "CAR_ROT", &status); TEST_PARAM("CAR_ROT");
00579     crlt_obs = drms_getkey_double(inRec, "CRLT_OBS", &status); TEST_PARAM("CRLT_OBS");
00580     crln_obs = drms_getkey_double(inRec, "CRLN_OBS", &status); TEST_PARAM("CRLN_OBS");
00581     crpix1 = drms_getkey_double(inRec, "CRPIX1", &status); TEST_PARAM("CRPIX1");
00582     x0 = crpix1 - 1;
00583     crpix2 = drms_getkey_double(inRec, "CRPIX2", &status); TEST_PARAM("CRPIX2");
00584     y0 = crpix2 - 1;
00585     pa = -drms_getkey_double(inRec, "CROTA2", &status); TEST_PARAM("CROTA2");
00586     // convert to radians
00587     crlt_obs_rad = Deg2Rad*crlt_obs;;
00588     crln_obs_rad = Deg2Rad*crln_obs;;
00589     pa_rad = Deg2Rad*pa;
00590 
00591 //  XXXXXXXXXXXXXXXXX get coordinate mapping info using keywords from first record
00592     if (firstimage)
00593       {
00594 fprintf(stderr,"start exam of first image params.\n");
00595       firstimage = 0;
00596       if (loctype==LOCCARR)
00597         {
00598         ctype1 = strdup(drms_getkey_string(inRec, "CTYPE1", &status)); TEST_PARAM("CTYPE1");
00599         if (strcmp(ctype1, "HPLN-TAN") != 0) DIE2("CTYPE1 not HPLN-TAN as required, is: ", ctype1);
00600         ctype2 = strdup(drms_getkey_string(inRec, "CTYPE2", &status)); TEST_PARAM("CTYPE2");
00601         if (strcmp(ctype2, "HPLT-TAN") != 0) DIE2("CTYPE2 not HPLT-TAN as required, is: ", ctype2);
00602         rsun_ref = drms_getkey_double(inRec, "RSUN_REF", &status);
00603         if (status) rsun_ref = 6.96e8;
00604         dsun_obs = drms_getkey_double(inRec, "DSUN_OBS", &status); TEST_PARAM("DSUN_OBS");
00605         rsun = asin(rsun_ref/dsun_obs)*Rad2arcsec; 
00606         cdelt = drms_getkey_double(inRec, "CDELT1", &status); TEST_PARAM("CDELT1");
00607         // in principle use deltlong to get time to use for correcting crlt_obs, ignore for now
00608         // get ur, ll coords of box at CM.
00609         rsunpix = rsun/cdelt;
00610         }
00611       if (boxtype == BOXDEGREE)
00612         {
00613         int nx, ny;
00614         // get corners of sample box at zero longitude wrt center.
00615         sphere2img(Deg2Rad*(crlt+height/2), Deg2Rad*(width/2), crlt_obs_rad, 0.0,
00616           &urx, &ury, 0.0, 0.0, rsunpix, pa_rad, 0, 0, 0, 0);
00617         sphere2img(Deg2Rad*(crlt-height/2), -Deg2Rad*(width/2), crlt_obs_rad, 0.0,
00618           &llx, &lly, 0.0, 0.0, rsunpix, pa_rad, 0, 0, 0, 0);
00619         nx = urx-llx;
00620         ny = ury-lly;
00621         // fix for 180 flipped coords.
00622         if (nx < 0) nx = -nx;
00623         if (ny < 0) ny = -ny;
00624         llx = -nx/2;  lly = -ny/2;  urx = nx/2; ury = ny/2;
00625         urx = round(urx); ury = round(ury); llx = round(llx); lly = round(lly);
00626         // make size odd, increase smaller distance from center
00627         nx = urx - llx + 1;  ny = ury - lly + 1;
00628         if (!(nx&1)) {if (urx > -llx) llx-=1; else urx+=1;}
00629         if (!(ny&1)) {if (ury > -lly) lly-=1; else ury+=1;}
00630         width_arcsec = nx * cdelt;
00631         height_arcsec = nx * cdelt;
00632         }
00633       else if (boxtype==BOXARCSEC)
00634         {
00635         int nx, ny;
00636         urx = width/(2.0*cdelt); llx = -urx; ury = height/(2.0*cdelt); lly = -ury;
00637         urx = round(urx); ury = round(ury); llx = round(llx); lly = round(lly);
00638         nx = urx - llx + 1;  ny = ury - lly + 1;
00639         width_arcsec = nx * cdelt;
00640         height_arcsec = nx * cdelt;
00641         }
00642       else /* boxtype == BOXPIXELS */
00643         {
00644         int nx, ny;
00645         int iurx, iury, illx, illy;
00646         iurx = round(width/2.0); illx = -iurx; iury = round(height/2.0); illy = -iury;
00647         if (((int)(round(width)) & 1) == 0 && ((iurx-illx+1) & 1) == 1) iurx--;
00648         if (((int)(round(height)) & 1) == 0 && ((iury-illy+1) & 1) == 1) iury--;
00649         urx = iurx; ury = iury; llx = illx; lly = illy;
00650         nx = urx - llx + 1;  ny = ury - lly + 1;
00651         width_arcsec = nx * cdelt;
00652         height_arcsec = nx * cdelt;
00653         }
00654       if (NoTrack)
00655           {
00656           sphere2img(crlt_rad, crln_rad, crlt_obs_rad, crln_obs_rad, &center_x, &center_y, x0, y0, rsunpix, pa_rad, 0, 0, 0, 0);
00657           center_x_first = center_x - x0;
00658           center_y_first = center_y - y0;
00659           }
00660 fprintf(stderr,"finish exam of first image params.\n");
00661       }
00662 
00663     if (NoTrack)
00664         {
00665         center_x = center_x_first + x0;
00666         center_y = center_y_first + y0;
00667         }
00668     else
00669         {
00670         sphere2img(crlt_rad, crln_rad, crlt_obs_rad, crln_obs_rad, &center_x, &center_y, x0, y0, rsunpix, pa_rad, 0, 0, 0, 0);
00671         }
00672     int x1 = round(center_x + llx);
00673     int y1 = round(center_y + lly);
00674     int x2 = round(center_x + urx);
00675     int y2 = round(center_y + ury);
00676     crpix1 = 1 + x0 - x1;
00677     crpix2 = 1 + y0 - y1;
00678 
00679     int start1[2] = {x1, y1};
00680     int end1[2] = {x2, y2};
00681 
00682     inSeg = drms_segment_lookupnum(inRec, 0);
00683     inAxis[0] = inSeg->axis[0];
00684     inAxis[1] = inSeg->axis[1];
00685 
00686     if (x1 >= inAxis[0] || y1 >= inAxis[1] || x2 < 0 || y2 < 0)
00687       {
00688 fprintf(stderr, "Rec %d, slice outside image, (x1,y1)=(%d,%d), (x2,y2)=(%d,%d)\n",irec,x1,y1,x2,y2);
00689       continue; // slice outside image
00690       }
00691 
00692     if (x1>=0 && y1>=0 && x2<inAxis[0] && y2<inAxis[1])
00693       {
00694       status = 0;
00695       outArray = drms_segment_readslice(inSeg, DRMS_TYPE_FLOAT, start1, end1, &status);
00696       if (status) DIE("Cant read input record");
00697       }
00698     else
00699       {
00700       int dims[2] = {x2-x1+1,y2-y1+1};
00701       int start2[2], end2[2];
00702       int x,y;
00703       outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, NULL, &status);
00704       drms_array2missing(outArray);
00705       start2[0] = x1 < 0 ? 0 : x1;
00706       start2[1] = y1 < 0 ? 0 : y1;
00707       end2[0] = x2 >= inAxis[0] ? inAxis[0]-1 : x2;
00708       end2[1] = y2 >= inAxis[1] ? inAxis[1]-1 : y2;
00709       status = 0;
00710       inArray = drms_segment_readslice(inSeg, DRMS_TYPE_FLOAT, start2, end2, &status);
00711       if (status) DIE("Cant read input record");
00712       int start3[2] = {start2[0]-start1[0],start2[1]-start1[1]}; // fetched slice offset in outarray
00713       int end3[2] = {end2[0]-end1[0],end2[1]-end1[1]};           // fetched slice offset in outarray
00714       int n3x = end2[0] - start2[0] + 1;
00715       int n3y = end2[1] - start2[1] + 1;
00716       for (y=0; y< n3y; y++)
00717         for (x=0; x < n3x; x++)
00718           *((float *)outArray->data + ((start3[1]+y)*dims[0] + x + start3[0])) =
00719             *((float *)inArray->data + (y*n3x + x));
00720       drms_free_array(inArray);
00721       }
00722     // Handle special case where |pa| == 180
00723     if (fabs(fabs(pa)-180.0) < 1.0)
00724       {
00725       // Tracking should be OK for center, but patch rotated.  Flip on both axes.
00726       float *data = (float *)outArray->data;
00727       int i,j;
00728       int nx = outArray->axis[0];
00729       int ny = outArray->axis[1];
00730       int midrow = (ny)/2;
00731       int midcol = (nx)/2;
00732       float val;
00733       for (j=0; j<midrow; j++)
00734         for (i=0; i<nx; i++)
00735           {
00736           val = data[j*nx + i];
00737           data[j*nx + i] = data[(ny - 1 - j)*nx + nx - 1 - i];
00738       data[(ny - 1 - j)*nx + nx - 1 - i] = val;
00739       }
00740       if ((nx & 1) != 0)
00741         {
00742         for (i=0; i<midcol; i++)
00743           {
00744           val = data[midrow*nx + i];
00745           data[midrow*nx + i] = data[midrow*nx + nx - 1 - i];
00746       data[midrow*nx + nx - 1 - i] = val;
00747           }
00748         }
00749       pa -= 180.0;
00750       crpix1 = 1 + x2 - x0;
00751       crpix2 = 1 + y2 - y0;
00752       cosa = 1.0; sina = 0.0;
00753       }
00754 
00755 /*
00756  *               writing the extracted region data file
00757 */
00758     outRS = drms_create_records(drms_env, 1, outseries, DRMS_PERMANENT, &status);
00759     if (status) {fprintf(stderr,"Output series is %s, ",outseries); DIE("Cant make outout record");}
00760     outRec = outRS->records[0];
00761     drms_copykeys(outRec, inRec, 0, kDRMS_KeyClass_Explicit);
00762     outSeg = drms_segment_lookupnum(outRec, 0);
00763     outArray->bzero = outSeg->bzero;
00764     outArray->bscale = outSeg->bscale;
00765     set_statistics(outSeg, outArray, 1);
00766     drms_setkey_float(outRec, "XCEN", WX((outArray->axis[0]+1)/2, (outArray->axis[1]+1)/2));
00767     drms_setkey_float(outRec, "YCEN", WY((outArray->axis[0]+1)/2, (outArray->axis[1]+1)/2));
00768     status = drms_segment_write(outSeg, outArray, 0);
00769     if (status) DIE("problem writing file");
00770     drms_free_array(outArray);
00771     drms_sprint_rec_query(inQuery, inRec);
00772     drms_setkey_string(outRec, "SOURCE", inQuery);
00773     drms_setkey_time(outRec, "T_REC", trec);
00774     drms_setkey_double(outRec, "CRPIX1", crpix1);
00775     drms_setkey_double(outRec, "CRPIX2", crpix2);
00776     drms_setkey_double(outRec, "CRVAL1", 0.0);
00777     drms_setkey_double(outRec, "CRVAL2", 0.0);
00778     drms_setkey_double(outRec, "CRDELT1", cdelt);
00779     drms_setkey_double(outRec, "CRDELT2", cdelt);
00780     drms_setkey_double(outRec, "CROTA2", 0.0);
00781     drms_setkey_string(outRec, "CONTENT", "Tracked Extracted Patches, made by hg_patch");
00782     drms_setkey_string(outRec, "COMMENTS", "Patches");
00783     drms_setkey_string(outRec, "RequestID", requestid);
00784     drms_setkey_string(outRec, "HGBOXUNITS", boxunits);
00785     drms_setkey_string(outRec, "HGLOCUNITS", locunits);
00786     drms_setkey_float(outRec, "HGWIDE", width);
00787     drms_setkey_float(outRec, "HGHIGH", height);
00788     drms_setkey_float(outRec, "HGASWIDE", width_arcsec);
00789     drms_setkey_float(outRec, "HGASHIGH", height_arcsec);
00790     drms_setkey_float(outRec, "HGCRLN", crln);
00791     drms_setkey_float(outRec, "HGCRLT", crlt);
00792     drms_setkey_int(outRec, "HGCARROT", car_rot);
00793     drms_setkey_time(outRec, "HGTSTART", t_start);
00794     drms_setkey_time(outRec, "HGTSTOP", t_stop);
00795     drms_setkey_time(outRec, "DATE", time(0) + UNIX_EPOCH);
00796     drms_setkey_string(outRec, "HGQUERY", in);
00797     // drms_copykey(outRec, inRec, "EXPTIME");
00798     // drms_copykey(outRec, inRec, "QUALITY");
00799     // drms_copykey(outRec, inRec, "DSUN_OBS");
00800     // drms_copykey(outRec, inRec, "CRLN_OBS");
00801     // drms_copykey(outRec, inRec, "CRLT_OBS");
00802     // drms_copykey(outRec, inRec, "OBS_VR");
00803     // drms_copykey(outRec, inRec, "OBS_VW");
00804     // drms_copykey(outRec, inRec, "OBS_VN");
00805     // drms_copykey(outRec, inRec, "CAR_ROT");
00806     // drms_copykey(outRec, inRec, "T_OBS");
00807     // drms_copykey(outRec, inRec, "DATE__OBS");
00808     // drms_copykey(outRec, inRec, "WAVELNTH");
00809 fprintf(stderr,"Rec %d done OK\n", irec);
00810 
00811     OK_recs += 1;
00812 
00813     if (log)
00814       {
00815       drms_fprint_rec_query(log, outRec);
00816       fprintf(log, "\n");
00817       }
00818     drms_close_records(outRS, DRMS_INSERT_RECORD);
00819     } 
00820 
00821     drms_close_records(inRS, DRMS_FREE_RECORD); 
00822     if (in_filename)
00823       {
00824       unlink(in_filename);
00825       }
00826 
00827   if (OK_recs)
00828     printf("  DONE, %d records made! \n \n",OK_recs);
00829   else
00830     printf("  DONE but no useful records created! \n \n");
00831 
00832   if (log)
00833     fclose(log);
00834   return DRMS_SUCCESS;
00835 }    // 
00836 
00837 /*
00838  *
00839  *                     FUNCTIONS
00840  *
00841  *          img2sphere()
00842  *
00843  *  Written: Rick Bogart 
00844  *
00845 */
00846 int img2sphere (double x, double y, double ang_r, double latc, double lonc,
00847     double pa, double *p_rho, double *p_lat, double *p_lon, double *p_sinlat,
00848     double *p_coslat, double *p_sig, double *p_mu, double *p_chi)
00849   {
00850     /*
00851  *  Taken from cartography.c, written by Rick Bogart
00852  *  Map projected coordinates (x, y) to (lon, lat) and (rho | sig, chi)
00853  *  
00854  *  Arguments:
00855  *    x }           Plate locations, in same units as the image radius, relative
00856  *    y }               to the image center
00857  *    ang_r         Apparent semi-diameter of sun (angular radius of sun at
00858  *                      the observer's tangent line)
00859  *    latc          Latitude of disc center, uncorrected for light travel time
00860  *    lonc          Longitude of disc center
00861  *    pa            Position angle of solar north on image, measured eastward
00862  *                      from north (sky coordinates)
00863  *  Return values:
00864  *    rho           Angle point:sun_center:observer
00865  *    lon           Heliographic longitude
00866  *    lat           Heliographic latitude
00867  *    sinlat        sine of heliographic latitude
00868  *    coslat        cosine of heliographic latitude
00869  *    sig           Angle point:observer:sun_center
00870  *    mu            cosine of angle between the point:observer line and the
00871  *                      local normal
00872  *    chi           Position angle on image measured westward from solar
00873  *                      north
00874  *
00875  *  All angles are in radians.
00876  *  Return value is 1 if point is outside solar radius (in which case the
00877  *    heliographic coordinates and mu are meaningless), 0 otherwise.
00878  *  It is assumed that the image is direct; the x or y coordinates require a
00879  *    sign change if the image is inverted.
00880  *
00881  * modified to ignore setting return values for NULL pointers.
00882  *
00883  */
00884 
00885 
00886   double w_rho, w_lat, w_lon, w_sinlat, w_coslat, w_sig, w_mu, w_chi;
00887   double *rho= &w_rho, *lat= &w_lat, *lon= &w_lon, *sinlat= &w_sinlat, *coslat= &w_coslat, *sig= &w_sig, *mu= &w_mu, *chi= &w_chi;
00888   static double ang_r0 = 0.0, sinang_r = 0.0, tanang_r = 0.0;
00889   static double latc0 = 0.0, coslatc = 1.0, sinlatc = 0.0;
00890   double cosr, sinr, sinlon, sinsig;
00891 
00892   if (p_rho) rho=p_rho;
00893   if (p_lat) lat=p_lat;
00894   if (p_lon) lon=p_lon;
00895   if (p_sinlat) sinlat=p_sinlat;
00896   if (p_coslat) coslat=p_coslat;
00897   if (p_sig) sig=p_sig;
00898   if (p_mu) mu=p_mu;
00899   if (p_chi) chi=p_chi;
00900 
00901   if (ang_r != ang_r0) {
00902       sinang_r = sin(ang_r);
00903       tanang_r = tan(ang_r);
00904       ang_r0 = ang_r;
00905     }
00906 
00907   if (latc != latc0) {
00908       sinlatc = sin(latc);
00909       coslatc = cos(latc);
00910       latc0 = latc;
00911     }
00912 
00913   // position angle
00914   *chi = atan2 (x, y) + pa;
00915   while (*chi > 2 * M_PI) *chi -= 2 * M_PI;
00916   while (*chi < 0.0) *chi += 2 * M_PI;
00917        /*  Camera curvature correction, no small angle approximations  */
00918   *sig = atan (hypot (x, y) * tanang_r);
00919   sinsig = sin (*sig);
00920   *rho = asin (sinsig / sinang_r) - *sig;
00921 if (*sig > ang_r) return (-1);
00922   *mu = cos (*rho + *sig);
00923   sinr = sin (*rho);
00924   cosr = cos (*rho);
00925   *sinlat = sinlatc * cosr + coslatc * sinr * cos (*chi);
00926   *coslat = sqrt (1.0 - *sinlat * *sinlat);
00927   *lat = asin (*sinlat);
00928 
00929   sinlon = sinr * sin (*chi) / *coslat;
00930   *lon = asin (sinlon);
00931   if (cosr < (*sinlat * sinlatc)) *lon = M_PI - *lon;
00932   *lon += lonc;
00933   while (*lon < 0.0) *lon += 2 * M_PI;
00934   while (*lon >= 2 * M_PI) *lon -= 2 * M_PI;
00935   return (0);
00936 }
00937 /*
00938  *
00939  *          sphere2img()
00940  *
00941  *  Written: Rick Bogart 
00942  *
00943 */
00944 
00945 int sphere2img (double lat, double lon, double latc, double lonc,
00946         double *x, double *y, double xcenter, double ycenter,
00947         double rsun, double peff, double ecc, double chi,
00948         int xinvrt, int yinvrt) {
00949 /*
00950  *  Taken from cartography.c, written by Rick Bogart
00951  *  Perform a mapping from heliographic coordinates latitude and longitude
00952  *    (in radians) to plate location on an image of the sun.  The plate
00953  *    location is in units of the image radius and is given relative to
00954  *    the image center.  The function returns 1 if the requested point is
00955  *    on the far side (>90 deg from disc center), 0 otherwise.
00956  *
00957  *  Arguments:
00958  *      lat         Latitude (in radians)
00959  *      lon         Longitude (in radians)
00960  *      latc        Heliographic latitude of the disc center (in radians)
00961  *      lonc        Heliographic longitude of the disc center (in radians)
00962  *      *x }        Plate locations, in same units as the image radius, NOT relative
00963  *      *y }               to the image center
00964  *      xcenter }   Plate locations of the image center, in units of the
00965  *      ycenter }     image radius, and measured from an arbitrary origin
00966  *                    (presumably the plate center or a corner)
00967  *      rsun        Apparent semi-diameter of the solar disc, in plate
00968  *                    coordinates
00969  *      peff        Position angle of the heliographic pole, measured
00970  *                    eastward from north, relative to the north direction
00971  *                    on the plate, in radians
00972  *      ecc         Eccentricity of the fit ellipse presumed due to image
00973  *                    distortion (no distortion in direction of major axis)
00974  *      chi         Position angle of the major axis of the fit ellipse, 
00975  *                    measure eastward from north,  relative to the north
00976  *                    direction on the plate, in radians (ignored if ecc = 0)
00977  *      xinvrt}     Flag parameters: if not equal to 0, the respective
00978  *      yinvrt}       coordinates on the image x and y are inverted
00979  *
00980  *  The heliographic coordinates are first mapped into the polar coordinates
00981  *    in an orthographic projection centered at the appropriate location and
00982  *    oriented with north in direction of increasing y and west in direction
00983  *    of increasing x.  The radial coordinate is corrected for foreshortening
00984  *    due to the finite distance to the Sun. If the eccentricity of the fit
00985  *    ellipse is non-zero the coordinate of the mapped point is proportionately
00986  *    reduced in the direction parallel to the minor axis.
00987  *
00988  *  Bugs:
00989  *    The finite distance correction uses a fixed apparent semi-diameter
00990  *    of 16'01'' appropriate to 1.0 AU.  In principle the plate radius could
00991  *    be used, but this would require the plate scale to be supplied and the
00992  *    correction would probably be erroneous and in any case negligible.
00993  *
00994  *    The ellipsoidal correction has not been tested very thoroughly.
00995  *
00996  *    The return value is based on a test which does not take foreshortening
00997  *    into account.
00998  */
00999   static double sin_asd = 0.004660, cos_asd = 0.99998914;
01000                                                    /*  appropriate to 1 AU  */
01001   static double last_latc = 0.0, cos_latc = 1.0, sin_latc = 0.0;
01002   double r, cos_cang, xr, yr;
01003   double sin_lat, cos_lat, cos_lat_lon, cospa, sinpa;
01004   double squash, cchi, schi, c2chi, s2chi, xp, yp;
01005   int hemisphere;
01006 
01007   if (latc != last_latc) {
01008       sin_latc = sin (latc);
01009       cos_latc = cos (latc);
01010       last_latc = latc;
01011     }
01012   sin_lat     = sin (lat);
01013   cos_lat     = cos (lat);
01014   cos_lat_lon = cos_lat * cos (lon - lonc);
01015 
01016   cos_cang   = sin_lat * sin_latc + cos_latc * cos_lat_lon;
01017   hemisphere = (cos_cang < 0.0) ? 1 : 0;
01018   r = rsun * cos_asd / (1.0 - cos_cang * sin_asd);
01019   xr = r * cos_lat * sin (lon - lonc);
01020   yr = r * (sin_lat * cos_latc - sin_latc * cos_lat_lon);
01021                         /*  Change sign for inverted images  */
01022   if (xinvrt) xr *= -1.0;
01023   if (yinvrt) yr *= -1.0;
01024                         /*  Correction for ellipsoidal squashing of image  */
01025   if (ecc > 0.0 && ecc < 1.0) {
01026     squash = sqrt (1.0 - ecc * ecc);
01027     cchi = cos (chi);
01028     schi = sin (chi);
01029     s2chi = schi * schi;
01030     c2chi = 1.0 - s2chi;
01031     xp = xr * (s2chi + squash * c2chi) - yr * (1.0 - squash) * schi * cchi;
01032     yp = yr * (c2chi + squash * s2chi) - xr * (1.0 - squash) * schi * cchi;
01033     xr = xp;
01034     yr = yp;
01035   }
01036 
01037   cospa = cos (peff);
01038   sinpa = sin (peff);
01039   *x = xr * cospa - yr * sinpa;
01040   *y = xr * sinpa + yr * cospa;
01041 
01042   *y += ycenter;
01043   *x += xcenter;
01044 
01045   return (hemisphere);
01046 }
01047 
01048 
01049 // Generate explicit recordset list of closest good record to desired grid
01050 // First get vector of times and quality
01051 // Then if vector is not OK, quit.
01052 // then: make temp file to hold recordset list
01053 //       start with first time to define desired grid, 
01054 //       make array of desired times.
01055 //       make empty array of recnums
01056 //       search vector for good images nearest desired times
01057 //       for each found time, write record query
01058 char *get_input_recset(DRMS_Env_t *drms_env, char *in, TIME cadence)
01059   {
01060   DRMS_Array_t *data;
01061   TIME t_start, t_stop, t_now, t_want, t_diff, this_t_diff;
01062   int status;
01063   int nrecs, irec;
01064   int nslots, islot;
01065   long long *recnums;
01066   TIME *t_this, half = cadence/2.0;
01067   double *drecnum, *dquality;
01068   int quality;
01069   long long recnum;
01070   char keylist[DRMS_MAXQUERYLEN];
01071   static char filename[100];
01072   char *tmpdir;
01073   FILE *tmpfile;
01074   char *lbracket;
01075   char seriesname[DRMS_MAXQUERYLEN];
01076 
01077   sprintf(keylist, "T_OBS,QUALITY,recnum");
01078   data = drms_record_getvector(drms_env, in, keylist, DRMS_TYPE_DOUBLE, 0, &status);
01079   if (!data || status)
01080     {
01081     fprintf(stderr, "getkey_vector failed status=%d\n", status);
01082     return(NULL);
01083     }
01084   nrecs = data->axis[1];
01085   irec = 0;
01086   t_this = (TIME *)data->data;
01087   dquality = (double *)data->data + 1*nrecs;
01088   drecnum = (double *)data->data + 2*nrecs;
01089   t_start = t_this[0];
01090   t_stop = t_this[nrecs-1];
01091   nslots = (t_stop - t_start + cadence/2)/cadence;
01092   recnums = (long long *)malloc(nslots*sizeof(long long));
01093   for (islot=0; islot<nslots; islot++)
01094     recnums[islot] = 0;
01095   islot = 0;
01096   t_want = t_start;
01097   t_diff = 1.0e9; 
01098   for (irec = 0; irec<nrecs; irec++)
01099       {
01100       t_now = t_this[irec];
01101       quality = (int)dquality[irec] & 0xFFFFFFFF;
01102       recnum = (long long)drecnum[irec];
01103       this_t_diff = fabs(t_now - t_want);
01104       if (quality < 0)
01105         continue;
01106       if (t_now <= (t_want-half))
01107         continue;
01108       while (t_now > (t_want+half))
01109         {
01110         islot++;
01111         if (islot >= nslots)
01112            break;
01113         t_want = t_start + cadence * islot;
01114         this_t_diff = fabs(t_now - t_want);
01115         t_diff = 1.0e8;
01116         }
01117       if (this_t_diff <= t_diff)
01118         recnums[islot] = recnum;
01119       t_diff = fabs(t_now - t_want);
01120       }
01121   if (islot+1 < nslots)
01122     nslots = islot+1;  // take what we got.
01123   strcpy(seriesname, in);
01124   lbracket = index(seriesname,'[');
01125   if (lbracket) *lbracket = '\0';
01126   tmpdir = getenv("TMPDIR");
01127   if (!tmpdir) tmpdir = "/tmp";
01128   sprintf(filename, "%s/hg_patchXXXXXX", tmpdir);
01129   mkstemp(filename);
01130   tmpfile = fopen(filename,"w");
01131   for (islot=0; islot<nslots; islot++)
01132     if (recnums[islot])
01133       fprintf(tmpfile, "%s[:#%lld]\n", seriesname, recnums[islot]);
01134   fclose(tmpfile);
01135   free(recnums);
01136   drms_free_array(data);
01137   return(filename);
01138   }
proj/util/apps/hg_patch.c
