00001
00002
00003 #include <sys/types.h>
00004 #include <sys/stat.h>
00005 #include <unistd.h>
00006 #include <stdlib.h>
00007 #include <string.h>
00008 #include <time.h>
00009 #include "jsoc_main.h"
00010 #include "drms.h"
00011 #include "fitsexport.h"
00012 #include "drms_names.h"
00013
00014 #define NOT_SPECIFIED "***Not Specified***"
00015 #define DIE(msg) {fprintf(stderr,"$$$$ %s: %s\n", module_name, msg); return 1;}
00016
00017 int limbcompute(float *x, int nx, int ny, float xcguess, float ycguess, float rguess, float rrange,
00018 int limbmode, int useprevflag, float fwhm);
00019
00020 float rsun_offset(int wavelength);
00021
00022 ModuleArgs_t module_args[] =
00023 {
00024 {ARG_STRING, "dsinp", NOT_SPECIFIED, "Series name w/optional record spec"},
00025 {ARG_STRING, "dsout", NOT_SPECIFIED, "Output series"},
00026 {ARG_STRING, "segment", "0", "segment number"},
00027 {ARG_STRING, "sum", "1", "number of images to sum to reduce noise"},
00028 {ARG_FLAG, "h", "0", "Print usage message and quit"},
00029 {ARG_FLAG, "v", "0", "verbose flag"},
00030 {ARG_END}
00031 };
00032
00033 extern float sdisk_xc, sdisk_yc, sdisk_r;
00034 char *module_name = "limbfit_aia";
00035 int verbose = 0;
00036
00037 int nice_intro(int help)
00038 {
00039 int usage = cmdparams_get_int(&cmdparams, "h", NULL) != 0;
00040 verbose = cmdparams_get_int(&cmdparams, "v", NULL) != 0;
00041 if (usage || help) {
00042 fprintf(stderr, "limbfit_aia {-h} {-v} dsinp=series_record_spec dsout=out_series\n"
00043 " -h: print this message\n"
00044 " -v: verbose\n"
00045 "dsinp=<recordset query> as <series>{[record specifier]} - required\n"
00046 "segment=<segment number> default is 0\n"
00047 "sum=<number of images to sum> to reduce noise\n"
00048 );
00049 return(1);
00050 }
00051 return(0);
00052 }
00053
00054 int DoIt ()
00055 {
00056 char *dsinp, *outpathbase, *date_obs, outfilename[512];
00057 int limbmode=0, useprevflag=0;
00058 int fsn, good_rec, irec, isum=0, nfits=0, nrecs, segment, status=0, sum, wl;
00059 int nx, ny;
00060 int yr, mo, da, hr, mn, sc, ss;
00061 float rguess, xcguess, ycguess, rrange=100.0, fwhm=3.0, *sumarr=NULL;
00062 DRMS_Record_t *inprec;
00063 DRMS_RecordSet_t *inprs;
00064 DRMS_Segment_t *inpseg;
00065 DRMS_Array_t *inparr=NULL;
00066 struct stat sbuf;
00067 TIME t_obs;
00068 int i, ut;
00069
00070 if (nice_intro(0)) return(0);
00071 dsinp = strdup(cmdparams_get_str(&cmdparams, "dsinp", NULL));
00072 segment = cmdparams_get_int(&cmdparams, "segment", NULL);
00073 if (strcmp(dsinp, NOT_SPECIFIED)==0) DIE("in argument is required");
00074 sum = cmdparams_get_int(&cmdparams, "sum", NULL);
00075 inprs = drms_open_records(drms_env, dsinp, &status);
00076 if (status) DIE("cant open recordset query");
00077 drms_stage_records(inprs, 1, 0);
00078 nrecs = inprs->n;
00079 fprintf(stderr, "%d records, sum: %d\n", nrecs, sum);
00080 for (irec=0; irec<nrecs; irec++) {
00081 char *filename, *inpsegname, outpath[512];
00082 inprec = inprs->records[irec];
00083 inpseg = drms_segment_lookupnum(inprec, segment);
00084 inpsegname = inpseg->info->name;
00085 filename = inpseg->filename;
00086 if(strlen(filename)) {
00087 if (verbose) printf("rec %d, protocol %d, '%s', '%s'\n",
00088 irec, inpseg->info->protocol, inpsegname, filename);
00089 wl = drms_getkey_int(inprec, "WAVELNTH", &status);
00090 fsn = drms_getkey_int(inprec, "FSN", &status);
00091 rguess = drms_getkey_float(inprec, "R_SUN", &status);
00092 rguess += rsun_offset(wl);
00093 xcguess = drms_getkey_float(inprec, "CRPIX1", &status);
00094 ycguess = drms_getkey_float(inprec, "CRPIX2", &status);
00095 t_obs = drms_getkey_time(inprec, "T_OBS", &status);
00096 date_obs = drms_getkey_string(inprec, "T_OBS", &status);
00097 if ((t_obs<0.0) || (fsn == 0x1c001c00)) continue;
00098 limbmode = 0;
00099 if (wl == 304) limbmode = 3;
00100 else if (wl > 2000) limbmode = 2;
00101 else if (wl > 1000) limbmode = 1;
00102 inparr = drms_segment_read(inpseg, DRMS_TYPE_FLOAT, &status);
00103 if (status) DIE("drms_segment_read failed!");
00104 nx = inparr->axis[0]; ny = inparr->axis[1];
00105 if (isum) { float *values = inparr->data;
00106 for (i=0; i<nx*ny; i++) sumarr[i] += values[i];
00107 } else {
00108 if (sum == 1) sumarr = inparr->data;
00109 else sumarr = (float*) calloc(nx*ny, sizeof(float));
00110 memcpy(sumarr, inparr->data, nx*ny*sizeof(float));
00111 }
00112 isum++;
00113 if (isum == sum) {
00114 limbcompute(sumarr, nx, ny, xcguess, ycguess, rguess, rrange,
00115 limbmode, useprevflag, fwhm);
00116 if (sum>1) { free(sumarr); sumarr = NULL; }
00117 nfits++; isum = 0;
00118 ut = t_obs + 220924763;
00119 printf("%.2f\t%.2f\t%.2f\t%d\t%s\t%.2f\n",
00120 sdisk_xc, sdisk_yc, sdisk_r, ut, date_obs, rguess);
00121 }
00122 if (inparr) drms_free_array(inparr);
00123 }
00124 }
00125
00126 return status;
00127 }
