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