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