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