JSOC_Documentation: proj/mag/synop/apps/brblosdailysynframe.c Source File

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