JSOC_Documentation: proj/mag/synop/apps/brblosdailysynframe

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