JSOC_Documentation: proj/mag/synop/apps/dailysynframe

00001 
00007 // -----------------------------------------------------------
00008 #include <stdio.h>
00009 #include <stdlib.h>
00010 #include <time.h>
00011 #include "jsoc_main.h"
00012 #include <mkl_blas.h>
00013 #include <mkl_service.h>
00014 #include <mkl_lapack.h>
00015 #include <mkl_vml_functions.h>
00016 #include <omp.h>
00017 #include "fresize.h"
00018 #include "fstats.h"
00019 
00020 char *module_name = "dailysynframe_nrt";
00021 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status = %d \n", msg, status); return(status);}
00022 #define PARAMETER_ERROR(PNAME)
00023 #define     PI     4.0 * atan(1.0)
00024 #define QUAL_CHECK      (0xfffefb00)
00025 
00026 void frebinbox(float *image_in, float *image_out, int nx, int ny, int nbinx, int nbiny);
00027 
00028 ModuleArgs_t module_args[] =
00029   {
00030     {ARG_STRING, "in", "NOTSPECIFIED", "in"},
00031     {ARG_STRING, "out", "NOTSPECIFIED", "out"},
00032     {ARG_STRING, "synoptic", "NOTSPECIFIED", "synoptic"},
00033     {ARG_STRING, "smallsyn", "NOTSPECIFIED", ""},
00034     {ARG_STRING, "drmethod", "NOCORRDIFFROT", "drmethod"},
00035     {ARG_INT, "magresoln", "NOTSPECIFIED", "magresoln"},
00036     {ARG_INT, "synresoln", "NOTSPECIFIED", "synresoln"},
00037     {ARG_INT, "xx1", "-1", "xx1"},
00038     {ARG_INT, "yy1", "-1", "yy1"},
00039     {ARG_END}
00040   };
00041 
00042   float zgrid(int jph, int ith, int cmp, int rup, int dbl, 
00043        float phd[], float thd[], float lad[], float cth[], 
00044        float sth[], float csc[], float scs[]);
00045 
00046 int DoIt(void)
00047 {
00048   DRMS_RecordSet_t *inRS, *inRSfinal;
00049   DRMS_Record_t *inRec, *inRecfinal, *outRec = NULL, *smallRD; 
00050   DRMS_Segment_t *inSeg, *inSegfinal, *outSeg, *smalloutSeg;
00051   DRMS_Array_t *inArray, *inArrayfinal;
00052   TIME t_rec, t_rec0;
00053   TIME halfw = 7200.0; //half window = 2.0 hours
00054   char *t_window = "240m"; //hard-coded T-window
00055   char *inQueryfinal, *trec_str = NULL, *smallRecQuery;
00056   float crlt, crln;
00057 
00058   float clog0;
00059   float aa, bb, cc, dd, ee;
00060   int status = DRMS_SUCCESS, nrecs, irec, quality;
00061   int i, j, crn;
00062   int xdim_syn, ydim_syn, xmg, ymg;
00063   char sdatestime[100], sdatetmp[100], timetmp[100], timeprint[100];
00064   memset(sdatetmp, 0, sizeof(sdatetmp));
00065   memset(timetmp, 0, sizeof(timetmp));
00066   memset(timeprint, 0, sizeof(timeprint));
00067   memset(sdatestime, 0, sizeof(sdatestime));
00068 
00069   char *sdate, *stime, *inQuery, *outQuery, *synQuery;
00070   char *drmethod;
00071   int magresoln = params_get_int(&cmdparams, "magresoln");
00072   int synresoln = params_get_int(&cmdparams, "synresoln");
00073   int xx1 = params_get_int(&cmdparams, "xx1");
00074   int yy1 = params_get_int(&cmdparams, "yy1");
00075   int nbin = 5;
00076 
00077   inQuery = (char *)params_get_str(&cmdparams, "in");
00078   outQuery = (char *)params_get_str(&cmdparams, "out");
00079   synQuery = (char *)params_get_str(&cmdparams, "synoptic");
00080   smallRecQuery = (char *)cmdparams_get_str(&cmdparams, "smallsyn", &status);
00081   drmethod = (char *)params_get_str(&cmdparams, "drmethod");
00082   char historyofthemodule[2048]; // put history info into the data
00083   sprintf(historyofthemodule,"o2helio.c bug corrected -- July 2013");
00084 
00085      aa = 13.1988; bb = 0.0; cc = 0.0;
00086   if (strcmp(drmethod, "Meunier") == 0)  //  Meunier
00087     {
00088      aa = 13.562; bb = -2.04; cc = -1.4875;
00089     }
00090 
00091   if (strcmp(drmethod, "Phil") == 0)     // Phil, ApJ 241:811-819, 1980
00092    {
00093     aa = 2.917; bb = -0.40; cc = -0.40;  //in mu rad/s;sidereal
00094     dd = 0.202006;                       // 1 deg/day=0.202006 mu rad/s
00095     aa = aa/dd; bb = bb/dd; cc=cc/dd;    // in deg / day
00096     ee = 0.930505;                       // synodic/sidereal=0.930505
00097     aa = aa*ee; bb = bb*ee; cc = cc*ee;  // in synodic
00098    }                                     // end "phil"
00099 
00100   if (strcmp(drmethod, "Pevtsov") == 0)
00101     {
00102      aa = 13.51; bb = -1.72;  cc = -2.31;
00103     }
00104 
00105   if (strcmp(drmethod, "Snodgrass") == 0)
00106    {
00107     aa = 2.897; bb = -0.339; cc = -0.485;//in mu rad/s,sidereal
00108     dd = 0.202006;                       // 1 deg/day=0.202006 mu rad/s
00109     aa = aa/dd;  bb=bb/dd; cc=cc/dd;     // in deg / day
00110     ee = 0.930505;                       // synodic/sidereal=0.930505
00111     aa = aa*ee; bb = bb*ee; cc = cc*ee;  //in synodic
00112    }
00113 
00114   outRec = drms_create_record(drms_env, outQuery, DRMS_PERMANENT, &status);
00115   if (status) DIE("Output recordset not created");
00116   smallRD = drms_create_record(drms_env, smallRecQuery, DRMS_PERMANENT, &status);
00117   if (status) DIE("Output record not created");
00118 
00119   inRS = drms_open_records(drms_env, inQuery, &status);
00120   if (status || inRS->n == 0)
00121      DIE("No input data found -- no remapped files");
00122   inRec = inRS->records[0];
00123   t_rec = drms_getkey_time(inRec, "T_REC", &status);
00124   t_rec0 = t_rec - halfw;
00125 
00126   inQueryfinal = (char *)malloc(100 * sizeof(char));
00127   trec_str = (char *)malloc(30 * sizeof(char));
00128   sprint_time(trec_str, t_rec0, "TAI", 0);
00129   sprintf(inQueryfinal, "%s[%s/%s]", inRec->seriesinfo->seriesname, trec_str, t_window);
00130   printf("%s\n", inQueryfinal);
00131   drms_close_records(inRS, DRMS_FREE_RECORD);
00132   free(trec_str);
00133 
00134   inRSfinal = drms_open_records(drms_env, inQueryfinal, &status);
00135   if (status || inRSfinal->n == 0) DIE("No input data found -- files contain no data");
00136   nrecs = inRSfinal->n;
00137 
00138 // find the middle data
00139 
00140   int count = 0;
00141   int *recp;
00142   int nref, rec_cen;
00143   recp = (int *)malloc(nrecs * sizeof(int));
00144   for (i = 0; i < nrecs; i++)
00145     {
00146         inRecfinal = inRSfinal->records[i];
00147         inSeg = drms_segment_lookupnum(inRecfinal, 0);
00148         inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00149         if (status)
00150            {
00151               printf(" No data file found, status=%d\n", status);
00152               drms_free_array(inArray);
00153               continue;
00154             }
00155         quality = drms_getkey_int(inRecfinal, "QUALITY", &status);
00156         if (quality & QUAL_CHECK)
00157            {
00158              printf("SKIP: error getting keyword %s: iRec = %d, Bad QUALITY = 0x%08x\n",
00159             "QUALITY", i, quality);
00160             continue;
00161            }
00162 
00163         recp[count] = i;
00164         count += 1;
00165         drms_free_array(inArray);
00166     }
00167 
00168   if (count==0) DIE("No input remapped data found");
00169   rec_cen = (int)(count/2);
00170   nref = recp[rec_cen];
00171   printf("middle data id=%d\n", nref);
00172   inRecfinal = inRSfinal->records[nref];
00173   inSeg = drms_segment_lookupnum(inRecfinal, 0);
00174   inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00175   int naxis = inArray->naxis;
00176   xmg = inArray->axis[0]; ymg = inArray->axis[1];
00177   t_rec = drms_getkey_time(inRecfinal, "T_REC", &status);
00178   t_rec0 = drms_getkey_time(inRecfinal, "T_OBS", &status);
00179   crn = drms_getkey_int(inRecfinal, "CAR_ROT", &status);
00180   crlt = drms_getkey_float(inRecfinal, "CRLT_OBS", &status);
00181   crln = drms_getkey_float(inRecfinal, "CRLN_OBS", &status);
00182   clog0 = drms_getkey_float(inRecfinal, "CRVAL1", &status);
00183 
00184 printf("crn=%d, clog0=%f\n", crn, clog0);
00185 
00186 //  drms_copykey(outRec, inRecfinal, "CADENCE");
00187   drms_copykey(outRec, inRecfinal, "DATASIGN");
00188   drms_copykey(outRec, inRecfinal, "DSUN_OBS");
00189   drms_copykey(outRec, inRecfinal, "OBS_VR");
00190   drms_copykey(outRec, inRecfinal, "OBS_VW");
00191   drms_copykey(outRec, inRecfinal, "OBS_VN");
00192   drms_copykey(outRec, inRecfinal, "QUALITY");
00193   drms_copykey(outRec, inRecfinal, "MAPMMAX");
00194   drms_copykey(outRec, inRecfinal, "MAPLGMAX");
00195   drms_copykey(outRec, inRecfinal, "MAPLGMIN");
00196   drms_copykey(outRec, inRecfinal, "SINBDIVS");
00197   drms_copykey(outRec, inRecfinal, "MAPBMAX");
00198   drms_copykey(outRec, inRecfinal, "MAPRMAX");
00199   drms_copykey(outRec, inRecfinal, "LGSHIFT");
00200   drms_copykey(outRec, inRecfinal, "INTERPO");
00201   drms_copykey(outRec, inRecfinal, "MCORLEV");
00202   drms_copykey(outRec, inRecfinal, "MOFFSET");
00203   drms_copykey(outRec, inRecfinal, "CARSTRCH");
00204   drms_copykey(outRec, inRecfinal, "DIFROT_A");
00205   drms_copykey(outRec, inRecfinal, "DIFROT_B");
00206   drms_copykey(outRec, inRecfinal, "DIFROT_C");
00207   drms_copykey(outRec, inRecfinal, "INSTRUME");
00208   drms_copykey(outRec, inRecfinal, "BLD_VERS");
00209   drms_copykey(outRec, inRecfinal, "CALVER64");
00210 
00211 //  smallRD;
00212 int itmp;
00213 float ftmp;
00214 double dtmp;
00215 
00216   drms_copykey(smallRD, inRecfinal, "DATASIGN");
00217   drms_copykey(smallRD, inRecfinal, "DSUN_OBS");
00218   drms_copykey(smallRD, inRecfinal, "OBS_VR");
00219   drms_copykey(smallRD, inRecfinal, "OBS_VW");
00220   drms_copykey(smallRD, inRecfinal, "OBS_VN");
00221   drms_copykey(smallRD, inRecfinal, "QUALITY");
00222   drms_copykey(smallRD, inRecfinal, "MAPMMAX");
00223   drms_copykey(smallRD, inRecfinal, "MAPLGMAX");
00224   drms_copykey(smallRD, inRecfinal, "MAPLGMIN");
00225   drms_copykey(smallRD, inRecfinal, "SINBDIVS");
00226   drms_copykey(smallRD, inRecfinal, "MAPBMAX");
00227   drms_copykey(smallRD, inRecfinal, "MAPRMAX");
00228   drms_copykey(smallRD, inRecfinal, "LGSHIFT");
00229   drms_copykey(smallRD, inRecfinal, "INTERPO");
00230   drms_copykey(smallRD, inRecfinal, "MCORLEV");
00231   drms_copykey(smallRD, inRecfinal, "MOFFSET");
00232   drms_copykey(smallRD, inRecfinal, "CARSTRCH");
00233   drms_copykey(smallRD, inRecfinal, "DIFROT_A");
00234   drms_copykey(smallRD, inRecfinal, "DIFROT_B");
00235   drms_copykey(smallRD, inRecfinal, "DIFROT_C");
00236   drms_copykey(smallRD, inRecfinal, "INSTRUME");
00237   drms_copykey(smallRD, inRecfinal, "BLD_VERS");
00238   drms_copykey(smallRD, inRecfinal, "CALVER64");
00239   drms_free_array(inArray);
00240 
00241 // average remapped mags
00242 
00243   float *aveData;
00244 //  xmg = 1801; ymg = 1440;
00245   xdim_syn = 3600; ydim_syn = 1440;
00246   int inDims[2] = {xmg, ymg};
00247   int dxsz = 2 * inDims[0];     // jph in IDL, zgrid.pro
00248   int ith = inDims[1];
00249   int ppd = xdim_syn/360;       // pixels per degree
00250   int xbeg = 30;
00251   if (xx1 == -1) xx1 = 60;              // in degrees
00252   if (yy1 == -1) yy1 = 0;                  // in pixels
00253   int hwd = xx1;       // in degree
00254   int ii, jj;
00255   int smallDims[2], xout, yout;
00256   TIME tobs_total = 0.0, tobs_ave;
00257   xx1 *= ppd;          // in pixels
00258   xbeg *= ppd;
00259   aveData = (float *)malloc(xmg * ymg * sizeof(float));
00260   xout = xdim_syn/nbin; yout = ydim_syn/(nbin-1);
00261   smallDims[0] = xout; smallDims[1] = yout;
00262 
00263   for (i = 0; i < count; i++)
00264     {
00265         inRecfinal = inRSfinal->records[recp[i]];
00266         inSeg = drms_segment_lookupnum(inRecfinal, 0);
00267         inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00268         float *inData = (float *)inArray->data;
00269         int crnn = drms_getkey_int(inRecfinal, "CAR_ROT", &status);
00270         float clogn = drms_getkey_float(inRecfinal, "CRVAL1", &status);
00271         int xshift = (int)(ppd * ((clogn - clog0) - 360.0 * (crnn - crn)));
00272         TIME Tobs = drms_getkey_time(inRecfinal, "T_OBS", &status);
00273         tobs_total += Tobs;
00274 
00275 // printf("count=, record id=, indata[]=%d, %d,%f\n", count, recp[i], inData[1081*50 + 900]);
00276         for (jj = 0; jj < ymg; jj++)
00277             {
00278             for (ii = xbeg; ii < xmg - xbeg; ii++)
00279                 {
00280                   aveData[jj * xmg + ii] += inData[jj * xmg + ii - xshift];
00281                 }
00282             } 
00283         drms_free_array(inArray);
00284     }
00285 
00286    tobs_ave = tobs_total/count;
00287    for (jj = 0; jj < ymg; jj++)
00288        for (ii = 0; ii < xmg; ii++)
00289            aveData[jj * xmg + ii] /= count;
00290 
00291   drms_close_records(inRSfinal, DRMS_FREE_RECORD);
00292   float thd[ith], csc[ith], phd[dxsz];
00293   float lad[ith], sth[ith], cth[ith], scs[dxsz];
00294   double dtor = PI/180.;
00295   float constrate = 13.1988; // Carrington rotation rate in degree/day
00296   float ratelow, sinlat;
00297   DRMS_RecordSet_t *synRS = NULL;
00298   DRMS_Record_t *synRec;
00299   DRMS_Segment_t *synSeg;
00300   DRMS_Array_t *synArray, *supsynArray, *outArray, *smalloutArray;
00301   float *synData, *supsynData, *outData, *smalloutData;
00302   int i1, j1;
00303   int synleftst = ppd * hwd, supleftst = ppd * clog0;
00304   zgrid(dxsz, ith, 0, 0, 0, phd, thd, lad, cth, sth, csc, scs);
00305 
00306   snprintf(timetmp, sizeof(timetmp), "%s[%d/2]", synQuery, crn-2);
00307   snprintf(timeprint, sizeof(timeprint), "%s", synQuery);
00308   synQuery = timetmp;
00309   printf("inputname= %s\n", synQuery);
00310 
00311   synRS = drms_open_records(drms_env, synQuery, &status);
00312   if (status || synRS->n == 0) DIE("No input data found -- no synoptic charts");
00313                      // start combining the synoptic charts
00314   int nds = synRS->n;
00315   int supxDim = (nds + 1) * xdim_syn;
00316   int supsynDim[2] = {supxDim, ydim_syn};
00317   int synDim[2] = {xdim_syn, ydim_syn};
00318   supsynArray = drms_array_create(DRMS_TYPE_FLOAT, 2, supsynDim, NULL, &status);
00319   supsynData = supsynArray->data;
00320   outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, synDim, NULL, &status);
00321   outData = outArray->data;
00322   smalloutArray = drms_array_create(DRMS_TYPE_FLOAT, 2, smallDims, NULL, &status);
00323   smalloutData = smalloutArray->data;
00324 
00325 printf("synopN=%d\n", nds);
00326   for (i = 0; i < nds; i++)
00327     {
00328        synRec = synRS->records[i];
00329        synSeg = drms_segment_lookupnum(synRec, 0);
00330        synArray = drms_segment_read(synSeg, DRMS_TYPE_FLOAT, &status);
00331        synData = synArray->data;
00332        int ii = ((nds + 1) - 1 - i) * xdim_syn;
00333        for (j1 = 0; j1 < ydim_syn; j1++)
00334           for (i1 = 0; i1 < xdim_syn; i1++)
00335             {
00336                supsynData[supxDim * j1 + ii + i1] = synData[xdim_syn * j1 + i1];
00337             }
00338      }
00339 
00340   drms_free_array(synArray);
00341   drms_close_records(synRS, DRMS_FREE_RECORD);
00342 
00343   snprintf(timetmp, sizeof(timetmp), "%s[%d]", timeprint, crn);
00344   synQuery = timetmp;
00345   printf("inputname= %s\n", synQuery);
00346 
00347   synRS = drms_open_records(drms_env, synQuery, &status);
00348   if (status || synRS->n == 0) DIE("No input data found -- no synoptic data files");
00349 
00350        synRec = synRS->records[0];
00351        synSeg = drms_segment_lookupnum(synRec, 0);
00352        synArray = drms_segment_read(synSeg, DRMS_TYPE_FLOAT, &status);
00353        synData = synArray->data;
00354        ii = 0;
00355        for (j1 = 0; j1 < ydim_syn; j1++)
00356           for (i1 = 0; i1 < xdim_syn; i1++)
00357             {
00358                supsynData[supxDim * j1 + ii + i1] = synData[xdim_syn * j1 + i1];
00359             }
00360 
00361                    // end of the combination: the super synoptic map--supsynData
00362 
00363                    // start to generate the right portion of synchronic map
00364                    // the differential rotation is taken into account.
00365    for (j = 0; j < ydim_syn/2; j++)
00366      {
00367         sinlat = sin(lad[j]*dtor);
00368         ratelow = aa + bb * sinlat * sinlat + cc * pow(sinlat, 4); // in degree/day
00369 
00370         for (i = synleftst; i < xdim_syn; i++)
00371         {
00372              float lon = (i - synleftst) * constrate/ratelow;
00373              float delta_x = lon - (int)lon;
00374              int lonpixel = (int)(supleftst + lon);
00375              float x1 = supsynData[j * supxDim + lonpixel];
00376              float x2 = supsynData[j * supxDim + lonpixel + 1];
00377              outData[j*xdim_syn + i] = (1.0 - delta_x) * x1 + delta_x * x2;
00378              jj = ydim_syn - 1 - j;
00379              x1 = supsynData[jj * supxDim + lonpixel];
00380              x2 = supsynData[jj * supxDim + lonpixel + 1];
00381              outData[jj*xdim_syn + i] = (1.0 - delta_x) * x1 + delta_x * x2;
00382         }
00383       }
00384 
00385 
00386     int magleft = ppd * (90 - hwd);
00387     for (j = 0; j < ydim_syn; j++)
00388         for (i = 0; i < 2*synleftst; i++)
00389             {
00390                 outData[j * xdim_syn + i] = aveData[j * inDims[0] + (int)magleft + i];
00391             }
00392 
00393   double statMin, statMax, statMedn, statMean, statSig, statSkew, statKurt;
00394   int statNgood;
00395   if (fstats(xdim_syn*ydim_syn, outData, &statMin, &statMax, &statMedn, &statMean, &statSig,
00396       &statSkew, &statKurt, &statNgood)) printf("\n Statistics computation failed\n");
00397 
00398   frebinbox(outData, smalloutData, xdim_syn, ydim_syn, nbin, nbin-1);
00399   double smallstatMin, smallstatMax, smallstatMedn, smallstatMean, smallstatSig, smallstatSkew, smallstatKurt;
00400   int smallstatNgood;
00401   if (fstats(xout*yout, smalloutData, &smallstatMin, &smallstatMax, &smallstatMedn, &smallstatMean, &smallstatSig,
00402       &smallstatSkew, &smallstatKurt, &smallstatNgood)) printf("\n Statistics computation failed\n");
00403 
00404 //    outRec = drms_create_record(drms_env, outQuery, DRMS_PERMANENT, &status);
00405 //    if (status) DIE("Output recordset not created");
00406 
00407     drms_setkey_time(outRec, "T_REC", t_rec);
00408     trec_str = (char *)malloc(30 * sizeof(char));
00409     sprint_time(trec_str, tobs_ave, "TAI", 0);
00410     drms_setkey_time(outRec, "T_OBS", tobs_ave);
00411     drms_setkey_int(outRec, "CAR_ROT", crn);
00412     drms_setkey_float(outRec, "CRLT_OBS", crlt);
00413     drms_setkey_float(outRec, "CRLN_OBS", crln);
00414     drms_setkey_float(outRec, "CADENCE", 24.0 * 60.0 * 60.0);
00415     drms_setkey_float(outRec, "CROTA2", 0.0);
00416     drms_setkey_string(outRec, "WCSNAME", "Carrington Heliographic");
00417     drms_setkey_string(outRec, "HISTORY", historyofthemodule);
00418 
00419     double loncen = xdim_syn/2 + 0.0;
00420     drms_setkey_double(outRec, "CRPIX1", loncen);
00421         // origin is at the left corner of the first pixel
00422     double latcen = ydim_syn/2 + 0.0;
00423     drms_setkey_double(outRec, "CRPIX2", latcen);
00424         // origin is at the left corner of the first pixel
00425     double lonstep = -360.0/xdim_syn;
00426     drms_setkey_double(outRec, "CDELT1", lonstep);
00427     double latstep = 2.0/ydim_syn;
00428     drms_setkey_double(outRec, "CDELT2", latstep);
00429     double lonfirst = 360.0 * (crn  - 1) - clog0 + hwd;
00430                          // longitude for the last pixel is counted at the center of the pixel.
00431     drms_setkey_double(outRec, "LON_FRST", lonfirst);
00432     double lonlast = 360.0 * crn - clog0 + hwd - 360.0/xdim_syn;
00433                          // longitude for the first pixel is counted at the center of the pixel.
00434     drms_setkey_double(outRec, "LON_LAST", lonlast);
00435     double loncenter = lonfirst + 180.0 - 0.5 * 360.0/xdim_syn;
00436     drms_setkey_double(outRec, "CRVAL1", loncenter);
00437     drms_setkey_double(outRec, "CARRTIME", loncenter);
00438     drms_setkey_double(outRec, "LON_STEP", lonstep);
00439 
00440 // synoptic map info
00441     float hwnwidth = drms_getkey_float(synRec, "HWNWIDTH", &status);
00442     drms_setkey_float(outRec, "HWNWIDTH", hwnwidth);
00443     float eqpoints = drms_getkey_float(synRec, "EQPOINTS", &status);
00444     drms_setkey_float(outRec, "EQPOINTS", eqpoints);
00445     float syndro_a = drms_getkey_float(synRec, "DIFROT_A", &status);
00446     drms_setkey_float(outRec, "OSYNDR_A", syndro_a);
00447     float syndro_b = drms_getkey_float(synRec, "DIFROT_B", &status);
00448     drms_setkey_float(outRec, "OSYNDR_B", syndro_b);
00449     float syndro_c = drms_getkey_float(synRec, "DIFROT_C", &status);
00450     drms_setkey_float(outRec, "OSYNDR_C", syndro_c);
00451 // image statistics
00452     i = xdim_syn*ydim_syn;
00453     drms_setkey_int(outRec, "TOTVALS", i);
00454     drms_setkey_int(outRec, "DATAVALS", statNgood);
00455     i = xdim_syn*ydim_syn-statNgood;
00456     drms_setkey_int(outRec, "MISSVALS", i);
00457     drms_setkey_double(outRec, "DATAMIN", statMin);
00458     drms_setkey_double(outRec, "DATAMAX", statMax);
00459     drms_setkey_double(outRec, "DATAMEDN", statMedn);
00460     drms_setkey_double(outRec, "DATAMEAN", statMean);
00461     drms_setkey_double(outRec, "DATARMS", statSig);
00462     drms_setkey_double(outRec, "DATASKEW", statSkew);
00463     drms_setkey_double(outRec, "DATAKURT", statKurt);
00464 // frame info.  
00465     drms_setkey_string(outRec, "FRTIMWDN", t_window); 
00466     drms_setkey_string(outRec, "SYNDRORA", drmethod);
00467     drms_setkey_int(outRec, "FRAVEPNT", count);
00468     drms_setkey_float(outRec, "FRWINDOW", 2.0 * hwd);
00469     drms_setkey_float(outRec, "SYNDRO_A", aa);
00470     drms_setkey_float(outRec, "SYNDRO_B", bb);
00471     drms_setkey_float(outRec, "SYNDRO_C", cc);
00472 
00473     drms_keyword_setdate(outRec);
00474 
00475        printf("       WRITING OUTPUT\n");
00476 
00477     outSeg = drms_segment_lookupnum(outRec, 0);
00478     outArray->parent_segment = outSeg;
00479     status = drms_segment_write(outSeg, outArray, 0);
00480     if (status) DIE("problem writing file");
00481 
00482 // writting the small size map 
00483     drms_setkey_time(smallRD, "T_REC", t_rec);
00484     drms_setkey_time(smallRD, "T_OBS", tobs_ave);
00485     drms_setkey_int(smallRD, "CAR_ROT", crn);
00486     drms_setkey_float(smallRD, "CRLT_OBS", crlt);
00487     drms_setkey_float(smallRD, "CRLN_OBS", crln);
00488     drms_setkey_float(smallRD, "CADENCE", 24.0 * 60.0 * 60.0);
00489     drms_setkey_float(smallRD, "CROTA2", 0.0);
00490     drms_setkey_string(smallRD, "WCSNAME", "Carrington Heliographic");
00491     drms_setkey_string(smallRD, "HISTORY", historyofthemodule);
00492 
00493     loncen = xout/2 + 0.0;
00494     drms_setkey_double(smallRD, "CRPIX1", loncen);
00495     latcen = yout/2 + 0.0;
00496     drms_setkey_double(smallRD, "CRPIX2", latcen);
00497     lonstep = -360.0/xout;
00498     drms_setkey_double(smallRD, "CDELT1", lonstep);
00499     latstep = 2.0/yout;
00500     drms_setkey_double(smallRD, "CDELT2", latstep);
00501     lonfirst = 360.0 * (crn  - 1) - clog0 + hwd;
00502     drms_setkey_double(smallRD, "LON_FRST", lonfirst);
00503     lonlast = 360.0 * crn - clog0 + hwd - 360.0/xout;
00504     drms_setkey_double(smallRD, "LON_LAST", lonlast);
00505     loncenter = lonfirst + 180.0 - 0.5 * 360.0/xout;
00506     drms_setkey_double(smallRD, "CRVAL1", loncenter);
00507     drms_setkey_double(smallRD, "CARRTIME", loncenter);
00508     drms_setkey_double(smallRD, "LON_STEP", lonstep);
00509 // synoptic map info
00510     hwnwidth = drms_getkey_float(synRec, "HWNWIDTH", &status);
00511     drms_setkey_float(smallRD, "HWNWIDTH", hwnwidth);
00512     eqpoints = drms_getkey_float(synRec, "EQPOINTS", &status);
00513     drms_setkey_float(smallRD, "EQPOINTS", eqpoints);
00514     syndro_a = drms_getkey_float(synRec, "DIFROT_A", &status);
00515     drms_setkey_float(smallRD, "OSYNDR_A", syndro_a);
00516     syndro_b = drms_getkey_float(synRec, "DIFROT_B", &status);
00517     drms_setkey_float(smallRD, "OSYNDR_B", syndro_b);
00518     syndro_c = drms_getkey_float(synRec, "DIFROT_C", &status);
00519     drms_setkey_float(smallRD, "OSYNDR_C", syndro_c);
00520 // image statistics
00521     i = xout*yout;
00522     drms_setkey_int(smallRD, "TOTVALS", i);
00523     drms_setkey_int(smallRD, "DATAVALS", smallstatNgood);
00524     i = xout*yout-smallstatNgood;
00525     drms_setkey_int(smallRD, "MISSVALS", i);
00526     drms_setkey_double(smallRD, "DATAMIN", smallstatMin);
00527     drms_setkey_double(smallRD, "DATAMAX", smallstatMax);
00528     drms_setkey_double(smallRD, "DATAMEDN", smallstatMedn);
00529     drms_setkey_double(smallRD, "DATAMEAN", smallstatMean);
00530     drms_setkey_double(smallRD, "DATARMS", smallstatSig);
00531     drms_setkey_double(smallRD, "DATASKEW", smallstatSkew);
00532     drms_setkey_double(smallRD, "DATAKURT", smallstatKurt);
00533 // frame info.  
00534     drms_setkey_string(smallRD, "FRTIMWDN", t_window);
00535     drms_setkey_string(smallRD, "SYNDRORA", drmethod);
00536     drms_setkey_int(smallRD, "FRAVEPNT", count);
00537     drms_setkey_float(smallRD, "FRWINDOW", 2.0 * hwd);
00538     drms_setkey_float(smallRD, "SYNDRO_A", aa);
00539     drms_setkey_float(smallRD, "SYNDRO_B", bb);
00540     drms_setkey_float(smallRD, "SYNDRO_C", cc);
00541 
00542     drms_keyword_setdate(smallRD);
00543 
00544        printf("       WRITING OUTPUT--small-size version\n");
00545 
00546     smalloutSeg = drms_segment_lookupnum(smallRD, 0);
00547     smalloutArray->parent_segment = smalloutSeg;
00548     status = drms_segment_write(smalloutSeg, smalloutArray, 0);
00549     if (status) DIE("problem writing file");
00550 
00551     free(trec_str); free(recp); free(aveData);
00552     drms_free_array(smalloutArray);
00553     drms_free_array(outArray);
00554     drms_free_array(supsynArray);
00555     drms_free_array(synArray);
00556 
00557     drms_close_record(smallRD, DRMS_INSERT_RECORD);
00558     drms_close_record(outRec, DRMS_INSERT_RECORD);
00559     drms_close_records(synRS, DRMS_FREE_RECORD);
00560 //    drms_close_records(inRS, DRMS_FREE_RECORD);
00561   return 0;
00562 }        // end DoIt
00563 
00564 /*                
00565  *            zgrid()
00566  *
00567  */
00568 float zgrid(int jph, int ith, int cmp, int rup, int dbl, 
00569    float phd[], float thd[], float lad[], float cth[], 
00570    float sth[], float csc[], float scs[])
00571 {
00572 
00573   double dtor = PI/180.;
00574   int i, j;
00575 
00576   float dcth = 2.0/ith;
00577   float dph = 360./jph;
00578   float cpc[jph], thr, dphr;
00579   float lthd[ith], lphd[jph], lscs[jph];
00580   float llad[ith], lsth[ith], lcth[ith];
00581 
00582 #if 0
00583   if (dbl == 1)  // checking condition "dbl"
00584      {
00585       double thd[ith], phd[jph], csc[ith], cpc[jph];
00586       double dph = 360./jph;
00587       double dcth = 2.0/ith;
00588       double lad[ith], sth[ith], cth[ith], thr;
00589       double lthd[ith], lphd[jph], lscs[jph];
00590       double llad[ith], lsth[ith], lcth[ith];
00591       double dphr;
00592       float scs[jph];
00593      }  // end if "dbl" 
00594 #endif
00595 
00596   for (i =0; i < ith; i++)
00597      {
00598       cth[i] = (i + 0.5) * dcth - 1.0;   // from south to north
00599       thr = acos(cth[i]);
00600       sth[i] = sin(thr);
00601       thd[i] = thr/dtor;
00602       lad[i] = 90-thd[i];
00603      }
00604 
00605   for (j = 0; j < jph; j++)
00606      phd[j] = (j + 0.5) * dph;          // from left to right
00607 
00608   if (cmp == 1)                   // checking condition "cmp"
00609       {
00610        for (i = 0; i < ith; i++)
00611          csc[i] = 1/sth[i];
00612        for (j = 0; j < jph; j++)
00613          {
00614           dphr = (phd[j] - cmp) * dtor;
00615           cpc[j] = cos(dphr);
00616           scs[j] = 1/cpc[j];
00617          }
00618       }   //end condition: cmp = 1
00619 
00620   if (rup == 1)   // condition "rup"
00621       {
00622        for (i = 0; i < ith; i++)
00623          {
00624           thd[i] = lthd[ith - i];
00625           lad[i] = llad[ith - 1];
00626           cth[i] = lcth[ith - i];
00627           sth[i] = lsth[ith - 1];
00628          }
00629        for (j = 0; j < jph; j++)
00630          {
00631           phd[j] = lphd[jph-j];
00632           scs[j] = lscs[jph-j];
00633          }
00634       } // end condition "rup"
00635   return 0;
00636 } // end function
00637 
00638 /* ################## Wrapper for Jesper's code ################## */
00639 
00640 void frebin(float *image_in, float *image_out, int nx, int ny, int nbin)
00641 {
00642   struct fresize_struct fresizes;
00643   int nxout, nyout;
00644   int nlead = nx;
00645 
00646   nxout = nx / nbin; nyout = ny / nbin;
00647   init_fresize_gaussian(&fresizes, (nbin / 2) * 2, (nbin / 2) * 2, nbin);
00648   fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, (nbin / 2) * 2, (nbin / 2) * 2, DRMS_MISSING_FLOAT);
00649   free_fresize(&fresizes);
00650 
00651 }
00652 
00653 void frebinbox(float *image_in, float *image_out, int nx, int ny, int nbinx, int nbiny)
00654 {
00655   int nxout, nyout;
00656   int ii, jj, i, j;
00657   nxout = nx / nbinx; nyout = ny / nbiny;
00658   for (j = 0; j < nyout; j++) {
00659     int yOff, jy;
00660     jy = j * nbiny;
00661     for (i = 0; i < nxout; i++) {
00662       int ix;
00663       ix = i * nbinx;
00664       float aveval = 0.0;
00665       int number = 0;
00666       for (jj = 0; jj < nbiny; jj++) {
00667         yOff = (jy + jj) * nx;
00668         for (ii = 0; ii < nbinx; ii++) {
00669           int iData = yOff + ix + ii;
00670           if (!(isnan(image_in[iData]))) {
00671             aveval += image_in[iData];
00672             number += 1;
00673           }
00674         }
00675       }
00676     image_out[j*nxout+i] = aveval/number;
00677     }
00678   }
00679 }
00680 
00681 // ************ END ********** END ******** END ********** END ***********
proj/mag/synop/apps/dailysynframe_nrt.c
