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

00001 
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include <time.h>
00005 #include <math.h>
00006 #include <sys/time.h>
00007 #include <sys/resource.h>
00008 #include <limits.h>
00009 #include <sys/types.h>
00010 #include <sys/stat.h>
00011 #include "jsoc_main.h"
00012 #include "astro.h" 
00013 #include "drms_dsdsapi.h"
00014 #include "fstats.h"
00015 #include "heliographic_coords.c"
00016 #include "/home/wso/src/libastro.d/solephem.c"
00017 
00018 /* cmd-line parameters */
00019 #define kRecSetIn       "in"
00020 #define kSeriesOut      "out"
00021 #define kSegIn          "segin"
00022 #define kSegOut         "segout"
00023 #define kNOTSPECIFIED   "not specified"
00024 
00025 #define EQPOINTS               "EQPOINTS"
00026 #define HALFWIN                "HALFWIN"
00027 #define LOSFLAG                "l"
00028 #define FORCEFLAG              "f" /* force - accept data without checking version */
00029 #define DLOGFLAG               "d" /* debug log - verbose stdout */
00030 #define NOISESIGMA             "NOISESIGMA"  /* std dev of noise */
00031 #define MAXNOISEADJ            "MAXNOISEADJ" /* maximum adjustment due to increases
00032                           * with latitude of the noise level
00033                           * (radial data only) */
00034 #define MINOUTLIERPTS          "MINOUTLIERPTS" /* minimum number of points that must
00035                           exist before outliers can be 
00036                               discarded when calulating summary
00037                               statistics */
00038 #define DPC_OBSR               "DPC_OBSR"
00039 #define FD_Magnetogram_Sum     "FD_Magnetogram_Sum"
00040 #define FD_Magnetogram         "FD_Magnetogram"
00041 #define HWNWIDTH               "HWNWIDTH"
00042 
00043 /* keywords in the output data */
00044 #define kCARSTRCH              "CARSTRCH"
00045 #define kDIFROT_A              "DIFROT_A"
00046 #define kDIFROT_B              "DIFROT_B"
00047 #define kDIFROT_C              "DIFROT_C"
00048 #define CR                     "CR"
00049 
00050 const float kOutScale = 0.1;
00051 const int kOutBPP = 16;
00052 const float kNOISE_EQ = 10.0;
00053 
00054 /*
00055 enum MagAggregate_t
00056 {
00057      kOneMinuteValue = 0,
00058      kFiveMinuteValue = 1
00059 };
00060 
00061 float gMagAgg[] =
00062 {
00063    1.0,
00064    2.75
00065 };
00066 */
00067 
00068 char *module_name = "hmisynop";
00069 
00070 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
00071 
00072 /* Describes a column of magnetogram data.  Fields other than slot refer to 
00073  * the magnetogram in which this data originated.  For example, dist is the 
00074  * distance (in pixels) between the CM - center of the magnetogram and the 
00075  * magnetogram's column of data.
00076  */
00077 
00078 typedef struct MagCol_struct
00079 {
00080   float dist;     /* Distance, in magnetogram pixels, from CM to column. */
00081   float *datacol; /* Ptr to column of data elonging to one magnetogram */
00082   float equivPts; /* 1-minute noise-equivalent points of this column */
00083   int ds;         /* ds number, for id purposes */
00084 //  int sn;         /* sn number, for id purposes */
00085   int col;        /* original magnetogram column */
00086   char valid;     /* 1 if in use */
00087 } MagCol_t;
00088 
00089 int CalcSynCols(int start,
00090         int end,
00091         int incr,
00092         MagCol_t **smc,
00093         float *synop,
00094         int *wt,
00095         char *ww,
00096         float *epts,
00097         float *losSynop,
00098         int *len,
00099         float center,
00100         float nEquivPtsReq, 
00101         int los,
00102         int radialFound,
00103         float sensAdj,
00104         float noiseS, 
00105         float nsig,
00106         float maxNoiseAdj, 
00107         int minOutPts,
00108         int dlog);
00109 
00110 int magStats(float **val, int npts, double sum, int outThreshold, double *avg, float *med);
00111 void SortMagCols(MagCol_t *mc, int nelem);
00112 void FreeMagColsData(int start, int end, int incr, int *n, MagCol_t **mc);
00113 void m_sort(MagCol_t *mc, MagCol_t *mc_cpy, int sz);
00114 
00115 ModuleArgs_t module_args[] =
00116 {
00117    {ARG_STRING,  kRecSetIn, "", "Input data records."},
00118    {ARG_STRING,  kSeriesOut, "", "Output data series."},
00119    {ARG_STRING,  kSegIn, kNOTSPECIFIED, ""},
00120    {ARG_STRING,  kSegOut, kNOTSPECIFIED, ""},
00121    {ARG_INT, CR, "", "", ""},
00122    {ARG_FLOAT, "nsig", "3.0", "", ""},
00123    {ARG_INT, "MAPMMAX", "1800", "", ""},
00124    {ARG_INT, "SINBDIVS", "720", "", ""},
00125    {ARG_FLOAT, "MAPLGMAX", "90.0", "", ""},
00126    {ARG_FLOAT, "MAPLGMIN", "-90.0", "", ""},
00127    {ARG_FLOAT, EQPOINTS, "20.0", "", ""},
00128    {ARG_FLOAT, HALFWIN, "15.0", "", ""},
00129    {ARG_FLOAT, "center", "0.0", "", ""},
00130    {ARG_INT, "checkqual", "0", "", ""},
00131    {ARG_STRING, "qualmask", "0x0", ""},
00132    {ARG_FLAG, LOSFLAG, "0", "-1", "1"},
00133    {ARG_FLAG, FORCEFLAG, "0", "-1", "1"},
00134    {ARG_FLAG, DLOGFLAG, "0", "-1", "1"},
00135    {ARG_FLOAT, NOISESIGMA, "3.0", "", ""},
00136    {ARG_FLOAT, MAXNOISEADJ, "3.0", "", ""},
00137    {ARG_INT, MINOUTLIERPTS, "4", "", ""},
00138    {ARG_END}
00139 };
00140 
00141 int DoIt(void)
00142 {
00143 
00144   int newstat = 0;
00145   int status = DRMS_SUCCESS;
00146   float *synop;
00147   int *wt;
00148   char *ww;
00149   float *epts = NULL;
00150   float *losSynop = NULL;
00151   float *ps, *pm;
00152   char *odir;
00153   struct stat stbuf;
00154   float nsig;
00155   int cr, mapmmax, sinbdivs;
00156   int checkqual;
00157   unsigned int quality, qualvld, qualmask;
00158   float lgmax, lgmin;
00159   float center;
00160   float halfWindow;
00161   int len[2]; /* 0 - x pixels, 1 - y pixels */
00162   int ngood, idx;
00163   int started, ended;
00164   int firstidx, lastidx;
00165   int mapmidcol, synmidcol;
00166   int cols2right, cols2left;
00167   int mrc, mlc, src, slc;
00168   int mapcols, col, row, syncol;
00169   double maprightct, mapleftct;
00170   double synleftct, synrightct;
00171   double synstart; /* CT, in degrees, of the right-most pixel of the synoptic chart to output */
00172   double synend; /* CT, in degrees, of the left-most pixel of the synoptic chart to output */
00173   double synstep; /* degrees per pixel */
00174   double subSynopStart, subSynopEnd, synRange;
00175   int SyncolStart, SyncolEnd, nsynop;
00176   char tstr[64];
00177   int mag_num = 0;
00178   TIME tstart, tstop, trot, tearth;
00179   TIME tdiff, magtime;
00180   double r,b,l,vr,vn,vw;
00181   int y,m,d;
00182   TIME tmod;
00183   double bearth;
00184   double carrtime;
00185   char dsname[50400], key[50400], logfn[50400], synfn[50400], wtfn[50400], losFn[50400];
00186   char eptsfn[50400];
00187   int ds, nds, noutRec;
00188   int fsn, lsn;
00189   int i,j;
00190   double statMin, statMax, statMedn, statMean, statSig, statSkew, statKurt;
00191   int statNgood;
00192 
00193   struct {
00194       int recno;
00195       double mapct; /* Carrington time of center pixel of remapped image (in degrees) MINUS center. */
00196       double mapCM; /* Carrington time of CM of remapped image (in degrees). */
00197       float mapdev; /* deviation between center pixel and either end of the map (in degrees) */
00198       int ds;     /* file id */
00199       float tmin;
00200       float tmax;
00201       TIME tobs;
00202           } imrec[50400];
00203 
00204   int nStackMags;
00205   float nEquivPtsReq;
00206   MagCol_t **sortedMagCol = NULL; 
00207 
00208   long long modVers = soi_vers_num;
00209   long long imgVers = 0;
00210   int radialFound = 0;
00211   int losFound = 0;
00212   int carrStretch = 1;
00213   float diffrotA = DRMS_MISSING_FLOAT;
00214   float diffrotB = DRMS_MISSING_FLOAT;
00215   float diffrotC = DRMS_MISSING_FLOAT;
00216 
00217   const char *synopFileName = NULL;
00218   const char *wtFileName = NULL;
00219   const char *eptsFileName = NULL;
00220 
00221   int los;
00222   int force;
00223   int dlog;
00224   float noiseS = cmdparams_get_float(&cmdparams, NOISESIGMA, &status);      /* noise sigma */
00225   float maxNoiseAdj = cmdparams_get_float(&cmdparams, MAXNOISEADJ, &status); /* maximum adj due to increases
00226                             * with latitude of the noise level
00227                             * (radial data only) */
00228   int minOutPts = cmdparams_get_int(&cmdparams, MINOUTLIERPTS, &status);     /* minimum number of points that 
00229                             * must exist before outliers can be
00230                             * discarded when calulating summary
00231                             * statistics */
00232   int sensAdjDone = -1;
00233   float sensAdj;
00234 
00235   int nxtSyncol = -1; /* the first synop column that hasn't been calculated */
00236   int calcsynret = 0;
00237 
00238   cr = cmdparams_get_int(&cmdparams, CR, &status);
00239   nsig = cmdparams_get_float(&cmdparams, "nsig", &status);
00240   mapmmax = cmdparams_get_int(&cmdparams, "MAPMMAX", &status);
00241   sinbdivs = cmdparams_get_int(&cmdparams, "SINBDIVS", &status);
00242   checkqual = cmdparams_get_int(&cmdparams, "checkqual", &status);
00243   qualmask = strtoul(cmdparams_get_str(&cmdparams, "qualmask", &status), (char **) 0, 0);
00244   center = cmdparams_get_float(&cmdparams, "center", &status);
00245   halfWindow = cmdparams_get_float(&cmdparams, HALFWIN, &status);
00246   lgmax = cmdparams_get_float(&cmdparams, "MAPLGMAX", &status);
00247   lgmin = cmdparams_get_float(&cmdparams, "MAPLGMIN", &status);
00248   firstidx = 0;
00249 
00250   if (lgmax <= lgmin) 
00251     {
00252       printf("MAPLGMAX must be greater than MAPLGMIN\n");
00253       DIE("MAKES_NO_SENSE");
00254     }
00255 
00256   nEquivPtsReq = cmdparams_get_float(&cmdparams, EQPOINTS, &status);
00257 
00258 /* synoptic char dimensions*/
00259   len[0] = rint(360.0 / (lgmax - lgmin) * mapmmax);
00260   len[1] = 2 * sinbdivs;
00261 
00262 printf("%d, %d\n", len[0], len[1]);
00263 
00264   synstep = 360.0 / len[0];
00265   synstart = (cr - 1) * 360.0 + 0.0 * synstep;
00266   synend = cr * 360.0 - 1.0 * synstep;
00267 
00268   nStackMags = rint(2 * halfWindow * 0.5); /* consecutive magnetograms are 
00269                       * shifted about 1 degree apart */
00270 /*
00271   tstart = meridian_crossing(360.0, cr);
00272   tstop = meridian_crossing(0.0, cr);
00273   trot = meridian_crossing(180.0, cr);
00274   tearth = local_earth_meridian_crossing(180.0, cr); 
00275   bearth = earth_B(tearth); 
00276 */
00277 
00278   tstart = HeliographicTime(cr, 360.0);
00279   tstop = HeliographicTime(cr, 0.0);
00280   trot = HeliographicTime(cr, 180.0);
00281   tearth = DRMS_MISSING_TIME;
00282   bearth = DRMS_MISSING_FLOAT;
00283 
00284   los = (cmdparams_isflagset(&cmdparams, LOSFLAG) != 0);
00285   force = (cmdparams_isflagset(&cmdparams, FORCEFLAG) != 0);
00286   dlog = (cmdparams_isflagset(&cmdparams, DLOGFLAG) != 0);
00287 
00288 printf("\n######### FIRST PASS #################\n");
00289 
00290 /* Determine how many of the input series are valid.  Store this in idx. */
00291 /* Also, determine the minimum and maximum CT covered by each image. */
00292 
00293   if (modVers < 0) modVers = -modVers;
00294     char *inRecQuery, *outRecQuery;
00295     inRecQuery = (char *)cmdparams_get_str(&cmdparams, kRecSetIn, &status);
00296     outRecQuery = (char *)cmdparams_get_str(&cmdparams, kSeriesOut, &status);
00297     DRMS_RecordSet_t *inRD, *outRD;
00298     DRMS_Record_t *outRec;
00299 
00300     inRD = drms_open_records(drms_env, inRecQuery, &status);
00301     if (status || inRD->n == 0)
00302         DIE("No input data found");
00303     nds = inRD->n;
00304     idx = 0;
00305 for (ds=0; ds < nds; ds++) 
00306   {
00307         double clong;
00308         double cmLong;
00309         double mapct = 0.0;
00310         double mapctnew;
00311         double mapCM = 0.0;
00312         double mapCMnew;
00313         double remapLgmax = 0.0;
00314         double remapLgmin = 0.0;
00315         int orot;
00316         char t_obs[64];
00317         char *imgVersStr = NULL;
00318         int rkey = 0;
00319         int csKey = 0;
00320         float csCoeffKey;
00321         int sensAdjKey;
00322         DRMS_Record_t *inRec;
00323         inRec = inRD->records[ds];
00324 
00325 /*
00326         if (checkqual) 
00327         {
00328         quality = drms_getkey_int(inRec, "QUALITY", &status);
00329         qualvld = drms_getkey_int(inRec, "QUAL_VLD", &status);
00330 
00331     if (drms_ismissing_int(quality) || drms_ismissing_int(qualvld)) 
00332       printf("  Can't find QUALITY in log file; quality check disabled\n"); 
00333       else if (quality & qualvld & qualmask) 
00334           {
00335           printf("  Bad QUALITY = 0x%08x; rejected\n", quality); 
00336           continue;
00337        }
00338         }
00339 
00340 */
00341 
00342 /*
00343       // Don't create synoptic charts from obsolete data.
00344       if (!force) // check data version (version 1.8) to avoid earlier version mags 
00345         {
00346      imgVers = 0;
00347      imgVersStr = drms_getkey_string(inRec, "BLDVER18", &status);
00348 
00349       if (imgVersStr)
00350         {
00351           char *endptr = NULL;
00352           imgVers = strtoll(imgVersStr, &endptr, 10);
00353           if (*endptr != '\0' || endptr == imgVersStr)
00354             {
00355                // Invalid build version 
00356                printf("  Invalid build version information.\n"); 
00357                printf("    Rejecting ds=%d.\n", ds); 
00358                continue;
00359             }
00360         }
00361         else
00362           {
00363             // Missing build number, bail on this image.
00364             printf("    Rejecting ds=%d \n", ds); 
00365             continue;
00366           }
00367 
00368         if (imgVers < 0) imgVers = -imgVers;
00369 
00370         if (imgVers < modVers)
00371           { 
00372             printf("  Attempt to use obsolete image (img lev18 vers %lld, mod lev18 vers %lld).\n", 
00373             imgVers, modVers);
00374             printf("    Rejecting ds=%d.\n", ds); 
00375             continue;
00376            }
00377           }
00378 
00379 */  
00380           /* Ensure that all images used are either radial values or 
00381             * line-of-sight values, not a mixture of the two. */
00382 
00383           rkey = drms_getkey_int(inRec, "FDRADIAL", &status);
00384           if (rkey > 0)
00385             {
00386            printf("  Found radial keyword %d, ds=%d.\n", rkey, ds); 
00387            if (losFound)
00388              {
00389                printf("  Attempt to use a mixture of radial and line-of-sight images.\n"); 
00390                printf("    Rejecting ds=%d.\n", ds); 
00391                continue;
00392              }
00393            else
00394              {
00395                radialFound = 1;
00396              }
00397              }
00398            else
00399              {
00400            if (radialFound)
00401              {
00402                printf("  Attempt to use a mixture of radial and line-of-sight images.\n"); 
00403                printf("    Rejecting ds=%d.\n", ds); 
00404                continue;
00405              }
00406            else
00407              {
00408                losFound = 1;
00409              }
00410              }
00411 
00412              csKey = drms_getkey_int(inRec, kCARSTRCH, &status);
00413              csKey = (drms_ismissing_int(csKey)) ? 0 : csKey;
00414 
00415              if (idx == 0)
00416                {
00417              carrStretch = csKey;
00418                }
00419              else
00420                {
00421              if (csKey != carrStretch)
00422                {
00423                  printf("  Attempt to mix carr-streched and non-carr-stretched images.\n");
00424                  printf("    Rejecting ds=%d.\n", ds); 
00425                  continue;
00426                }
00427                }
00428 
00429               if (carrStretch > 0)
00430                 {
00431               csCoeffKey = (float)drms_getkey_double(inRec, kDIFROT_A, &status);
00432               if (idx == 0)
00433                 {
00434                   diffrotA = (float)csCoeffKey;
00435                 }
00436               else
00437                 {
00438                   if (fabsf(csCoeffKey - diffrotA) > 0.001)
00439                     {
00440                   printf("  Attempt to use inconsistent carr stretch parameters.\n");
00441                   printf("    Rejecting ds=%d.\n", ds); 
00442                   continue;
00443                     }
00444                 }
00445 
00446                csCoeffKey = (float)drms_getkey_double(inRec, kDIFROT_B, &status);
00447                if (idx == 0)
00448                  {
00449                    diffrotB = (float)csCoeffKey;
00450                  }
00451                else
00452                  {
00453                    if (fabsf(csCoeffKey - diffrotB) > 0.001)
00454                      {
00455                    printf("  Attempt to use inconsistent carr stretch parameters.\n");
00456                    printf("    Rejecting ds=%d.\n", ds); 
00457                    continue;
00458                      }
00459                  }
00460 
00461                  csCoeffKey = (float)drms_getkey_double(inRec, kDIFROT_C, &status);
00462                  if (idx == 0)
00463                    {
00464                      diffrotC = (float)csCoeffKey;
00465                    }
00466                  else
00467                    {
00468                      if (fabsf(csCoeffKey - diffrotC) > 0.001)
00469                        {
00470                      printf("  Attempt to use inconsistent carr stretch parameters.\n");
00471                      printf("    Rejecting ds=%d.\n", ds); 
00472                      continue;
00473                         }
00474                    }
00475                    }
00476 
00477                       clong = drms_getkey_double(inRec, "CRVAL1", &status);
00478                       cmLong = drms_getkey_double(inRec, "CRLN_OBS", &status);
00479                       orot = drms_getkey_int(inRec, "CAR_ROT", &status);
00480                       if (drms_ismissing_float(clong) || drms_ismissing_int(orot)) 
00481                         {
00482                       printf("  Bad MAP_L0 or OBS_CR in header; skipped"); 
00483                       continue;
00484                         }
00485 
00486                       mapctnew = (double)(orot) * 360.0 - clong - center;
00487                       mapCMnew = (double)(orot) * 360.0 - cmLong - center;
00488                       if (mapCMnew < mapCM)
00489                         {
00490                       printf("  Source image (ds=%d) has a smaller CT than previous img; skipped.\n", ds);
00491                       continue;
00492                         }
00493                       else
00494                         {
00495                       mapct = mapctnew;
00496                       mapCM = mapCMnew;
00497                         }
00498 
00499                    remapLgmax = drms_getkey_double(inRec, "MAPLGMAX", &status);
00500                    remapLgmin = drms_getkey_double(inRec, "MAPLGMIN", &status);
00501                    imrec[idx].mapdev = (remapLgmax - remapLgmin) / 2.0;
00502                    imrec[idx].recno = drms_getkey_int(inRec, "I_DREC", &status);
00503                    imrec[idx].mapct = mapct;
00504                    imrec[idx].mapCM = mapCM;
00505                    imrec[idx].ds = ds;
00506                    strcpy(t_obs, drms_getkey_string(inRec, "T_OBS", &status));
00507                    imrec[idx].tobs = sscan_time(t_obs);
00508                    imrec[idx].tmin = MAX(mapCM - halfWindow, mapct + center - imrec[idx].mapdev);
00509                    imrec[idx].tmax = MIN(mapCM + halfWindow, mapct + center + imrec[idx].mapdev);
00510                    ++idx;
00511 
00512   } /* for each dataset */
00513   ngood = idx;
00514 
00515 /* Every image pixel has a Carrington time.  Calculate the CT range for each 
00516  * magnetogram.
00517  */
00518 //    ngood = idx;
00519 
00520 /*
00521     for (idx = 0; idx < ngood; ++idx) 
00522       {
00523          imrec[idx].tmin = MAX(imrec[idx].mapCM - halfWindow, imrec[idx].mapct + center - imrec[idx].mapdev);
00524          imrec[idx].tmax = MIN(imrec[idx].mapCM + halfWindow, imrec[idx].mapct + center + imrec[idx].mapdev);
00525     }
00526 */
00527 
00528 printf("\n######### SECOND PASS #################\n"); 
00529 
00530       if (!(wt = (int *) calloc(len[0],4)))
00531         DIE("MALLOC_FAILED");
00532       if (!(ww = (char *) malloc(len[0]*len[1]*2)))
00533         DIE("MALLOC_FAILED");
00534       if (!(epts = (float *) malloc(len[0] * len[1] * sizeof(float))))
00535         DIE("MALLOC_FAILED");
00536       if (!(synop = (float *) malloc(len[0]*len[1]*4)))
00537         DIE("MALLOC_FAILED");
00538       if (los && radialFound && !(losSynop = (float *)malloc(len[0] * len[1] * sizeof(float))))
00539           DIE("MALLOC_FAILED");
00540       sortedMagCol = (MagCol_t **)malloc(sizeof(MagCol_t *) * len[0]);
00541       if (!sortedMagCol)
00542           DIE("MALLOC_FAILED");
00543       memset(sortedMagCol, 0, sizeof(MagCol_t *) * len[0]);
00544 
00545 //    idx = 0;
00546     started = 0;
00547     ended = 0;
00548     nxtSyncol = len[0] - 1;
00549     synRange = 10.0 ; 
00550                     /* Longitude range in degree in synoptic chart;
00551                     divide the synoptic chart into pieces to avoid reading too much data */ 
00552     nsynop = rint(360.0/synRange);
00553 
00554 /* Step through successive magnetograms that overlap the synoptic chart's CT range.  
00555  * Stop as soon as the focal magnetogram ceases to overlap this range.
00556  */
00557 
00558 for (ds = 0; ds < nsynop; ds++)
00559   {
00560     subSynopStart = synstart + ds * synRange;
00561     subSynopEnd = synstart + (ds + 1) * synRange - synstep;
00562     SyncolStart = rint((synend - subSynopStart) / synstep);
00563     SyncolEnd = rint((synend - subSynopEnd) / synstep);
00564     printf("ds=%d\n", ds);
00565 
00566     for (idx = 0; idx < ngood; idx++) 
00567       {
00568         int recno;
00569         int col;
00570         float magtype = 0.0;
00571         char *mType = NULL;
00572         float equivPts = 0.0;
00573         char t_mag[64];
00574         TIME tmag;
00575         DRMS_Record_t *inRec;
00576         DRMS_Segment_t *inSeg = NULL;
00577         DRMS_Array_t *inArray;
00578         inRec = inRD->records[imrec[idx].ds];
00579  
00580         if (imrec[idx].tmax <= subSynopStart || imrec[idx].tmin >= subSynopEnd) continue;
00581 
00582         /* truncate the magnetogram's CT range at synstart or synend if necessary */
00583 //        maprightct = MAX(synstart, imrec[idx].tmin);
00584 //        mapleftct = MIN(synend, imrec[idx].tmax);
00585         maprightct = MAX(subSynopStart, imrec[idx].tmin);
00586         mapleftct = MIN(subSynopEnd, imrec[idx].tmax);
00587         mapcols = drms_getkey_int(inRec, "MAPMMAX", &status) + 1;
00588         mapmidcol = rint((mapcols - 1.0) / 2 + center / synstep);
00589         cols2right = rint((imrec[idx].mapct - maprightct) / synstep);
00590         cols2left = rint((imrec[idx].mapct - mapleftct) / synstep);
00591         synmidcol = rint((synend - imrec[idx].mapct) / synstep); 
00592                                                                /* NOT the middle column of the
00593                                 * synoptic chart, but essentially
00594                                 * the column of the synoptic chart
00595                                 * that corresponds to the middle
00596                                 * column of the remapped img. */
00597         mrc = mapmidcol + cols2right;
00598         mlc = mapmidcol + cols2left;
00599         /* src - right column of synoptic chart that matches with mrc */
00600         /* slc - left column of synoptic chart that matches with mlc */
00601         src = synmidcol + cols2right;
00602         slc = synmidcol + cols2left;
00603 
00604       /* For each column of the synoptic chart, push the corresponding input magnetogram's 
00605        * column of data to the stack at that column.  The stack is an array of 2D 
00606        * input images that will be averaged to make the synoptic chart.  The dimensions 
00607        * of these images match the final dimensions of the synoptic chart.  Upon 
00608        * completion of all input series, wt is the depth of the stack at each synoptic column. 
00609        */
00610  
00611       
00612 //        magtype = drms_getkey_float(inRec, "EXPTIME", &status);
00613 //        equivPts = gMagAgg[kOneMinuteValue];
00614 //        if (magtype == 300.0) equivPts = gMagAgg[kFiveMinuteValue];
00615         equivPts = 1.0;
00616         inSeg = drms_segment_lookupnum(inRec, 0);
00617         inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00618 
00619         if (status)
00620             {
00621             printf(" No data file found, status=%d\n", status);
00622             drms_free_array(inArray);
00623             continue;
00624             }
00625 
00626         for (col = mlc; col <= mrc; ++col) 
00627           {
00628         MagCol_t mMagCol;
00629 
00630             mMagCol.dist = ((col - mapmidcol) * synstep + imrec[idx].mapct) - imrec[idx].mapCM;
00631         mMagCol.datacol = (float *)malloc(sizeof(float) * len[1]);
00632         mMagCol.equivPts = equivPts;
00633         mMagCol.ds = imrec[idx].ds;
00634         mMagCol.col = col;
00635         mMagCol.valid = 1;
00636        
00637         pm = (float *)inArray->data + col;
00638         ps = mMagCol.datacol;
00639 
00640         for (row = 0; row < len[1]; ++row) 
00641           {
00642             *ps = *pm;
00643             pm += mapcols;
00644             ps++;
00645           }
00646 
00647         syncol = col + slc - mlc;
00648 
00649         if (!sortedMagCol[syncol])
00650           {
00651         sortedMagCol[syncol] = (MagCol_t *)malloc(sizeof(MagCol_t) * nStackMags);
00652         if (!sortedMagCol[syncol]) DIE("MALLOC_FAILED");
00653         memset(sortedMagCol[syncol], 0, sizeof(MagCol_t) * nStackMags);
00654           }
00655 
00656         if (wt[syncol] > 0 && (wt[syncol] % nStackMags == 0))
00657           {
00658             int frames = 1 + (wt[syncol] / nStackMags);
00659             MagCol_t *mc = (MagCol_t *)malloc(sizeof(MagCol_t) * frames * nStackMags);
00660             memcpy(mc, 
00661                sortedMagCol[syncol], 
00662                sizeof(MagCol_t) * (frames - 1) * nStackMags);
00663             free(sortedMagCol[syncol]);
00664             sortedMagCol[syncol] = mc;
00665           }
00666 
00667         sortedMagCol[syncol][wt[syncol]] = mMagCol;    
00668         wt[syncol]++;
00669 
00670        } /* for each series */
00671     drms_free_array(inArray);
00672     }
00673 
00674     calcsynret = CalcSynCols(SyncolStart,
00675                                 SyncolEnd,
00676                                 -1,
00677                                 sortedMagCol,
00678                                 synop,
00679                                 wt,
00680                                 ww,
00681                                 epts,
00682                                 losSynop,
00683                                 len,
00684                                 center,
00685                                 nEquivPtsReq,
00686                                 los,
00687                                 radialFound,
00688                                 sensAdj,
00689                                 noiseS,
00690                                 nsig,
00691                                 maxNoiseAdj,
00692                                 minOutPts,
00693                                 dlog);
00694 
00695     printf("ds=%d\n", ds);
00696     FreeMagColsData(SyncolStart, SyncolEnd, -1, wt, sortedMagCol);
00697 /*
00698     if (sortedMagCol)
00699       {
00700         for (i = SyncolStart; i >= SyncolEnd; i--)
00701           {
00702             if (sortedMagCol[i])
00703               {
00704                 free(sortedMagCol[i]);
00705               }
00706           }
00707        }
00708 
00709 */
00710     if (calcsynret) DIE(" Cannot be finished up!");
00711  }
00712 
00713     if (!ended) lastidx = ngood-1;
00714 
00715     tdiff = 1e99;
00716     for (idx=firstidx; idx<=lastidx; ++idx) 
00717         {
00718           if (fabs(imrec[idx].tobs - trot) < tdiff) 
00719             {
00720           tdiff = fabs(imrec[idx].tobs - trot);
00721           magtime = imrec[idx].tobs;
00722             }
00723         }
00724 
00725 // compute statistics
00726 
00727   if (fstats(len[0]*len[1], synop, &statMin, &statMax, &statMedn, &statMean, &statSig, 
00728       &statSkew, &statKurt, &statNgood)) printf("\n Statistics computation failed\n");
00729 
00730 /* start to archive the data */
00731 
00732     outRD = drms_create_records(drms_env, 1, outRecQuery, DRMS_PERMANENT, &status);
00733     if (status)
00734         DIE("Output record not created");
00735         DRMS_Segment_t *outSeg;
00736         DRMS_Array_t *outArray;
00737         int outDims[2] = {len[0],len[1]};
00738         outRec = outRD->records[0];
00739 
00740     // synop (first segment)
00741     outSeg = drms_segment_lookupnum(outRec, 0);
00742     outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, synop, &status);
00743     status = drms_segment_write(outSeg, outArray, 0);
00744         if (status)
00745           DIE("problem writing file");
00746         drms_free_array(outArray);
00747 
00748 
00749     // epts (second segment)
00750     outSeg = drms_segment_lookupnum(outRec, 1);
00751         outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, epts, &status);
00752     status = drms_segment_write(outSeg, outArray, 0);
00753         if (status)
00754           DIE("problem writing file");
00755         drms_free_array(outArray);
00756 
00757 /*
00758     // ww (third segment)
00759     outSeg = drms_segment_lookupnum(outRec, 2);
00760     outArray = drms_array_create(DRMS_TYPE_CHAR, 2, outDims, ww, &status);
00761     status = drms_segment_write(outSeg, outArray, 0);
00762         if (status)
00763           DIE("problem writing file");
00764         drms_free_array(outArray);
00765 */
00766 
00767 /* write the keywords */
00768 
00769 //for reserved keywords
00770 
00771         sprint_ut(tstr, CURRENT_SYSTEM_TIME);
00772         sscanf(tstr, "%d.%d.%d", &y, &m, &d);
00773         sprintf(tstr, "%04d-%02d-%02d", y, m, d);
00774         drms_setkey_string(outRec, "DATE", tstr);
00775         drms_setkey_string(outRec, "BUNIT", "Gauss");
00776 //        drms_setkey_int(outRec, "BLD_VERS", soi_version);
00777 
00778 //for image coordinate mapping keywords
00779 
00780         drms_setkey_string(outRec, "CTYPE1", "CRLN-CEA");
00781         drms_setkey_string(outRec, "CTYPE2", "CRLT-CEA");
00782         i=len[0]/2;
00783         drms_setkey_int(outRec, "CRPIX1", i);
00784         l=len[1]/2+0.5;
00785         drms_setkey_float(outRec, "CRPIX2", l);
00786         carrtime = cr*360.0 - 180.0;
00787         drms_setkey_float(outRec, "CRVAL1", carrtime);
00788         i=0.0;
00789         drms_setkey_float(outRec, "CRVAL2", i);
00790         l=-360.0/len[0];
00791         drms_setkey_float(outRec, "CDELT1", l);
00792         drms_setkey_float(outRec, "CDELT2", 1.0/sinbdivs);
00793         drms_setkey_string(outRec, "CUNIT1", "degree");
00794         drms_setkey_string(outRec, "CUNIT2", "Sine Latitude");
00795 
00796 //HMI observables keywords
00797         sprint_at(tstr, trot);
00798         drms_setkey_string(outRec, "T_REC", tstr);
00799 // image statistics
00800         i = len[0]*len[1];
00801         drms_setkey_int(outRec, "TOTVALS", i);
00802         drms_setkey_int(outRec, "DATAVALS", statNgood);
00803         i = len[0]*len[1]-statNgood;
00804         drms_setkey_int(outRec, "MISSVALS", i);
00805         drms_setkey_double(outRec, "DATAMIN", statMin);
00806         drms_setkey_double(outRec, "DATAMAX", statMax);
00807         drms_setkey_double(outRec, "DATAMEDN", statMedn);
00808         drms_setkey_double(outRec, "DATAMEAN", statMean);
00809         drms_setkey_double(outRec, "DATARMS", statSig);
00810         drms_setkey_double(outRec, "DATASKEW", statSkew);
00811         drms_setkey_double(outRec, "DATAKURT", statKurt);
00812 
00813 //Synoptic map keywords
00814         drms_setkey_string(outRec, "T_OBS", tstr);
00815         drms_setkey_string(outRec, "T_ROT", tstr);
00816         sprint_at(tstr, tstart);
00817         drms_setkey_string(outRec, "T_START", tstr);
00818         sprint_at(tstr, tstop);
00819         drms_setkey_string(outRec, "T_STOP", tstr);
00820         sprint_at(tstr, tearth);
00821         drms_setkey_string(outRec, "T_EARTH", tstr);
00822         drms_setkey_int(outRec, "CAR_ROT", cr);
00823         drms_setkey_double(outRec, "CARRTIME", carrtime);
00824         HeliographicLocation(trot, &i, &l, &b);
00825         drms_setkey_double(outRec, "B0_ROT", b);
00826         HeliographicLocation(tstart, &i, &l, &b);
00827         drms_setkey_double(outRec, "B0_FRST", b);
00828         HeliographicLocation(tstop, &i, &l, &b);
00829         drms_setkey_double(outRec, "B0_LAST", b);
00830         drms_setkey_double(outRec, "EARTH_B0", bearth);
00831         l=(cr-1)*360.0;
00832         drms_setkey_double(outRec, "LON_FRST", l);
00833         l=cr*360.0-360.0/len[0];
00834         drms_setkey_double(outRec, "LON_LAST", l);
00835         l=-360.0/len[0];
00836         drms_setkey_double(outRec, "LON_STEP", l);
00837         l=center;
00838         drms_setkey_double(outRec, "W_OFFSET", l);
00839         drms_setkey_string(outRec, "W_WEIGHT", "Even");
00840         i=lastidx-firstidx+1;
00841         drms_setkey_int(outRec, "MAG_NUM", i);
00842         sprint_at(tstr, imrec[firstidx].tobs);
00843         drms_setkey_string(outRec, "MAG_FRST", tstr);
00844         sprint_at(tstr, imrec[lastidx].tobs);
00845         drms_setkey_string(outRec, "MAG_LAST", tstr);
00846         sprint_at(tstr, magtime);
00847         drms_setkey_string(outRec, "MAG_ROT", tstr);
00848         drms_setkey_float(outRec, HWNWIDTH, halfWindow);
00849         drms_setkey_float(outRec, EQPOINTS, nEquivPtsReq);
00850         drms_setkey_float(outRec, "NSIGMA", nsig);
00851 //        drms_setkey_float(outRec, "NNOISE", noiseS); -- noiseS not be used so far.
00852         
00853 
00854 // save the differential rotation correction parameters 
00855         if (carrStretch > 0)
00856           {
00857             drms_setkey_int(outRec, kCARSTRCH, carrStretch);
00858             drms_setkey_float(outRec, kDIFROT_A, diffrotA);
00859             drms_setkey_float(outRec, kDIFROT_B, diffrotB);
00860             drms_setkey_float(outRec, kDIFROT_C, diffrotC);
00861           }
00862 
00863 
00864     if (sortedMagCol)
00865       {
00866         for (i = 0; i < len[0]; i++)
00867           {
00868         if (sortedMagCol[i])
00869           {
00870             free(sortedMagCol[i]);
00871           }
00872           } 
00873          free(sortedMagCol);
00874        }
00875 
00876   drms_close_records(inRD, DRMS_FREE_RECORD);
00877   drms_close_records(outRD, DRMS_INSERT_RECORD);
00878   return 0;
00879 }
00880 
00881 int CalcSynCols(int start,
00882         int end,
00883         int incr,
00884         MagCol_t **smc,
00885         float *synop,
00886         int *wt,
00887         char *ww,
00888         float *epts,
00889         float *losSynop,
00890         int *len,
00891         float center,
00892         float nEquivPtsReq, 
00893         int los,
00894         int radialFound,
00895         float sensAdj,
00896         float noiseS, 
00897         float nsig,
00898         float maxNoiseAdj, 
00899         int minOutPts,
00900         int dlog)
00901 {
00902    float midSynRow;
00903    float cosrho; /* equal cos(magnetogram latitude) * cos(magnetogram dlatitude) 
00904           *   where dlatitude is positive delta between CM and point.
00905           * approximate dlatitude with 'center' (assume for all 
00906           * magnetograms, data comes from CM - center, even though
00907           *   they will come from points surrounding this column.
00908           */
00909    float cosCenter;
00910    char *rejDescFormat;
00911    int col;
00912    int row;
00913    int i;
00914    int j;
00915    static int nrej = 0;
00916    int marked;
00917    float nEquivPts;
00918 
00919    midSynRow = (float)(len[1] - 1) / 2.0;
00920    cosCenter = cos(center * M_PI / 180);
00921    rejDescFormat = "      magnetogram: ds=%d, sn=%d, magcol=%d, magrow=%d\n";
00922    if (minOutPts < 2) minOutPts = 2;
00923 
00924    /* loop over synoptic chart column */
00925    for (col = start; incr > 0 ? col <= end : col >= end; col += incr) 
00926    {
00927       float *val = (float *)malloc(sizeof(float) * wt[col]);
00928       float *ept = (float *)malloc(sizeof(float) * wt[col]);
00929       MagCol_t *pcols = NULL;
00930       marked = 0;
00931 
00932       if (dlog)
00933       {
00934      pcols = (MagCol_t *)malloc(sizeof(MagCol_t) * wt[col]);
00935      memset(pcols, 0, sizeof(MagCol_t) * wt[col]);
00936       }
00937 
00938       if (!val || !ept) printf("MALLOC FAILED"); 
00939 
00940       /* sort MagCol_ts */
00941       SortMagCols(smc[col], wt[col]);
00942 
00943       /* print out all mags for this column */
00944 /*
00945       nEquivPts = 0.0;
00946       for (i = 0; i < wt[col]; i++) 
00947       {
00948      nEquivPts += smc[col][i].equivPts;
00949 
00950          if (i == 0) printf("it's begin\n");
00951          if (i >= 0 && i < 20) printf("col=%d, i=%d, ds=%d, dist=%f\n", col, i, smc[col][i].ds, smc[col][i].dist);
00952          if (i == wt[col]-1) printf("it's end\n");
00953 
00954      if (nEquivPts >= nEquivPtsReq && !marked) marked = 1;
00955       }
00956 */
00957 
00958       /* loop over synoptic chart row */
00959       for (row = 0; row < len[1]; ++row) 
00960       {
00961      float sinLat = (row - midSynRow) / midSynRow;
00962      float cosLat = sqrt(1 - sinLat * sinLat);
00963      float sum;
00964      int npts;
00965 
00966      if (radialFound) cosrho = cosLat * cosCenter;
00967 
00968      sum = 0.0;
00969      npts = 0;
00970      j = 0;
00971      nEquivPts = 0.0;
00972       
00973      /* Loop over synoptic chart stack above col, row, averaging the values. */
00974      /* Take the magnetograms with the best values first (smallest dist to CM) 
00975       * and take only the fewest possible until nEquivPtsReq one-minute-equivalent 
00976       * points are attained. */
00977      for (i = 0; i < wt[col] && nEquivPts < nEquivPtsReq; i++) 
00978      {
00979         float magVal = *(smc[col][i].datacol + row);
00980         if (!drms_ismissing_float(magVal)) 
00981         {
00982            nEquivPts += smc[col][i].equivPts;
00983            sum += magVal;
00984            val[j] = magVal;
00985            ept[j] = smc[col][i].equivPts;
00986 
00987            if (dlog && pcols) pcols[j] = smc[col][i];
00988 
00989            j++;
00990            npts++;
00991         }
00992      }
00993 
00994      /* Calculate summary statistics */
00995      if (npts > 1 && nsig > 0.0) {
00996         double avg = 0.0;
00997         double ssqr = 0.0;
00998         double sig;
00999         float med;
01000         int nStatPts = 0;
01001         float *statVals;
01002         float noiseLevel;
01003         float *outlier = NULL;
01004         int statsDone = 0;
01005         int sIdx = 0;
01006         float dev = 0.0;
01007         int nPtsB4Rej = npts;
01008 
01009         /* calculate threshold - below this, point values are within noise levels -
01010            we don't want to reject them under any circumstances */
01011         noiseLevel = kNOISE_EQ * noiseS * sensAdj; // noiseLevel not be used so far.
01012 
01013 //printf("noiseLevel = %f, sensAdj= %f\n", noiseLevel, sensAdj);
01014 
01015         if (radialFound) noiseLevel = noiseLevel * MIN(1 / cosrho, maxNoiseAdj);
01016 
01017         /* copy the vals - will be destructively analyzing them */
01018         statVals = (float *)malloc(sizeof(float) * npts);
01019 
01020         if (!statVals) printf("MALLOC_FAILED"); 
01021            
01022         memcpy(statVals, val, sizeof(float) * npts);
01023         nStatPts = magStats(&statVals, npts, sum, minOutPts, &avg, &med);
01024        
01025         /* Reject outliers whose values exceeds the noise threshold */
01026         for (j = 0; j < nPtsB4Rej; ++j) 
01027         {         
01028            if (fabs(val[j]) > noiseLevel)
01029            {
01030           if (!statsDone)
01031           {
01032              for (sIdx = 0; sIdx < nStatPts; sIdx++) 
01033              {
01034             ssqr += statVals[sIdx] * statVals[sIdx];
01035              }
01036 
01037              sig = sqrt((ssqr - nStatPts * avg * avg) / (nStatPts - 1));
01038              statsDone = 1;
01039           }
01040 
01041           dev = fabs(val[j] - med);
01042 
01043           if (dev > nsig * sig)
01044           {
01045              sum -= val[j];
01046              nEquivPts -= ept[j];
01047              --npts;
01048              ++nrej;
01049              if (dlog)
01050              {
01051             char rejDescBuf[128];
01052 
01053             snprintf(rejDescBuf, 
01054                  sizeof(rejDescBuf), 
01055                  rejDescFormat, 
01056                  pcols[j].ds,
01057 //               pcols[j].sn,
01058                  pcols[j].col,
01059                  row);
01060 
01061              } 
01062           }
01063            }
01064         }
01065 
01066         free(statVals);
01067      } /* npts > 1 && nsig > 0.0 */
01068 
01069     
01070       
01071      /* Calcuate the average value for each x,y in the stack */
01072      if (npts)
01073      {
01074 
01075         float synVal = sum/npts;
01076         float minVal = (SHRT_MIN + 1);  // * kOutScale;
01077         float maxVal = (SHRT_MAX - 1);  // * kOutScale;
01078 
01079         /* The desired output will be a 16 bpp image with bscale == 0.1.
01080          * If 16 bits cannot hold the float value, then cap the value at 
01081          * the largest value (in magnitude) that 16 bits can hold. */
01082         if (synVal < minVal)
01083         {
01084            synop[row * len[0] + col] = minVal;
01085         }
01086         else if (synVal > maxVal)
01087         {
01088            synop[row * len[0] + col] = maxVal;
01089         }
01090         else
01091         {
01092            synop[row * len[0] + col] = synVal;
01093         }
01094      }
01095      else
01096      {
01097         synop[row * len[0] + col] = DRMS_MISSING_FLOAT;
01098      }
01099       
01100      ww[row * len[0] + col] = npts;
01101      epts[row * len[0] + col] = nEquivPts;
01102      if (los && radialFound)
01103      {
01104         int offset = row * len[0] + col;
01105         if (drms_ismissing_float(synop[offset]))
01106         {
01107            losSynop[offset] = DRMS_MISSING_FLOAT;
01108         }
01109         else
01110         {
01111            losSynop[offset] = synop[offset] * cosLat;
01112         }
01113      }
01114       } /* row loop */
01115 
01116       free(val);
01117       val = NULL;
01118       free(ept);
01119       ept = NULL;
01120 
01121       if (pcols)
01122       {
01123      free(pcols);
01124      pcols = NULL;
01125       }
01126    } /* col loop */
01127    return 0;
01128 }
01129 
01130 int cmp(const void *a, const void *b)
01131 {
01132     float aa = *(float *)a;
01133     float bb = *(float *)b;
01134     return (aa<bb) ? -1 : (aa==bb) ? 0 : 1;
01135 }
01136 
01137 /* To simplify, assumes one outlier only */
01138 int magStats(float **val, int npts, double sum, int outThreshold, double *avg, float *med)
01139 {
01140    int actPts = npts;
01141    float *out = NULL;
01142    float *in = *val;
01143    int idx = 0;
01144    int iOut = 0;
01145    float medianInt = 0.0;
01146    float first;
01147    float last;
01148 
01149    qsort((void *)in, npts, sizeof(float), cmp);
01150 
01151    if (npts > outThreshold && npts > 1)
01152    {
01153       actPts = npts - 1;
01154       out = malloc(sizeof(float) * actPts);
01155 
01156       medianInt = in[npts/2];
01157 
01158       first = in[0];
01159       last = in[npts - 1];
01160 
01161       if (fabs(first - medianInt) <= fabs(last - medianInt))
01162       {
01163      idx = 0;
01164      sum -= last;
01165       }
01166       else
01167       {
01168      idx = 1;
01169      sum -= first;
01170       }
01171 
01172       for (iOut = 0; idx < npts && iOut < actPts; idx++, iOut++)
01173       {
01174      out[iOut] = in[idx];
01175       }
01176 
01177       free(in);
01178       *val = out;
01179       in = *val;
01180    }
01181 
01182    *avg = sum / actPts;
01183    *med = in[actPts / 2];
01184    
01185    return actPts;
01186 }
01187 
01188 int CmpMagCol(const void *a, const void *b)
01189 {
01190    MagCol_t *aa = (MagCol_t *)a;
01191    MagCol_t *bb = (MagCol_t *)b;
01192 
01193    if (aa->valid && bb->valid)
01194    {
01195       return (fabsf(aa->dist) < fabsf(bb->dist)) ? -1 : 
01196     ((fabsf(aa->dist) == fabsf(bb->dist)) ? 0 : 1);
01197    }
01198    else if (aa->valid)
01199    {
01200       return -1;
01201    }
01202    else if (bb->valid)
01203    {
01204       return 1;
01205    }
01206    else
01207    {
01208       return 0;
01209    }
01210 }
01211 
01212 void SortMagCols(MagCol_t *mc, int nelem)
01213 {
01214 
01215     MagCol_t *mc_cpy = NULL;
01216     int i;
01217     mc_cpy = (MagCol_t *)malloc(sizeof(MagCol_t) * nelem);
01218     for (i = 0; i < nelem; i++)
01219         mc_cpy[i] = mc[i];
01220     m_sort(mc, mc_cpy, nelem);
01221     free(mc_cpy);
01222 }
01223 
01224 void m_sort(MagCol_t *mc, MagCol_t *mc_cpy, int sz)
01225 {
01226     int n, i, j, mid;
01227         mid = sz/2;
01228     if (mid > 1) m_sort(mc_cpy, mc, mid);
01229     if (sz-mid > 1) m_sort(mc_cpy + mid, mc + mid, sz - mid);
01230 
01231     // merge 
01232     MagCol_t *mc_cpy2 = mc_cpy + mid;
01233     for (n=0, i=0, j=0; i<mid && j<sz-mid; n++)
01234         mc[n] = (fabs(mc_cpy[i].dist) < fabs(mc_cpy2[j].dist)) ? mc_cpy[i++] : mc_cpy2[j++];
01235     while (i<mid)
01236         mc[n++] = mc_cpy[i++];
01237     while (j<sz-mid)
01238         mc[n++] = mc_cpy2[j++];
01239 }
01240 
01241    /* qsort((void *)mc, nelem, sizeof(MagCol_t), CmpMagCol); */
01242    
01243    /* Try merge sort, since mc will be essentially two sorted lists. */
01244 //   int topNeg = -1;
01245 //   int topPos = -1;
01246 //   int i, j;
01247 //   MagCol_t *res = NULL;
01248 //   MagCol_t temp;
01249 //   int idx;
01250 
01251 /*
01252    for (i = 0; i < nelem; i++)
01253    {
01254       for (j = (i+1); j < nelem; j++)
01255       {
01256         if (fabs(mc[i].dist) > fabs(mc[j].dist))
01257         {
01258          temp = mc[i];
01259          mc[i] = mc[j];
01260          mc[j] = temp;
01261         }
01262       }
01263     }
01264     return;
01265 }
01266 
01267 
01268    for (idx = 0; idx < nelem; idx++)
01269    {
01270       if (mc[idx].dist >= 0)
01271       {
01272      topPos = idx;
01273      break;
01274       }
01275    }
01276 
01277    if (idx == nelem)
01278    {
01279       // only negative dists 
01280       topNeg = nelem - 1;
01281    }
01282    else if (topPos > 0)
01283    {
01284       // both negative and positive dists 
01285       topNeg = topPos - 1;
01286    }
01287 
01288    if (topNeg >= 0 || topPos >= 0)
01289    {
01290       res = (MagCol_t *)malloc(sizeof(MagCol_t) * nelem);
01291       memset(res, 0, sizeof(MagCol_t) * nelem);
01292    }
01293 
01294    if (res)
01295    {
01296       idx = 0;
01297       while (topNeg >=0 && topPos >= 0 && topPos < nelem)
01298       {
01299      if (fabsf(mc[topNeg].dist) < fabsf(mc[topPos].dist))
01300      {
01301         res[idx] = mc[topNeg];
01302         topNeg--;
01303      }
01304      else
01305      {
01306         res[idx] = mc[topPos];
01307         topPos++;
01308      }
01309 
01310      idx++;
01311       }
01312 
01313       while (topNeg >= 0)
01314       {
01315      res[idx] = mc[topNeg];
01316      topNeg--;
01317      idx++;
01318       }
01319 
01320       while (topPos >= 0 && topPos < nelem)
01321       {
01322      res[idx] = mc[topPos];
01323      topPos++;
01324      idx++;
01325       }
01326 
01327       memcpy(mc, res, sizeof(MagCol_t) * nelem);
01328    }
01329 }
01330 */
01331 
01332 void FreeMagColsData(int start, int end, int incr, int *n, MagCol_t **mc)
01333 {
01334    if (mc)
01335    {
01336       float *data = NULL;
01337       int cIdx;
01338       int idx;
01339 
01340       for (cIdx = start; incr > 0 ? cIdx <= end : cIdx >= end; cIdx += incr)
01341       {
01342      for (idx = 0; idx < n[cIdx]; idx++)
01343      {
01344         data = mc[cIdx][idx].datacol;
01345         if (data)
01346         {
01347            free(data);
01348            mc[cIdx][idx].datacol = NULL;
01349            mc[cIdx][idx].valid = 0;
01350         }
01351      }
01352       }
01353    }
01354 }
01355 
01356 /*
01357 double local_earth_meridian_crossing(double L, int cr)
01358 {
01359     double ta, tb, t, la, lb, l;
01360     double ephem[EPHEM_SIZE];
01361     double trot;
01362     
01363     trot = meridian_crossing(L, cr);
01364     ta = trot - 1800.0;
01365     tb = trot + 1800.0;
01366     sun_ephemeris(ta, ephem, 0.0, 0.0);
01367     la = fmod(ephem[EPHEM_L0],360.0)-L;
01368     sun_ephemeris(tb, ephem, 0.0, 0.0);
01369     lb = fmod(ephem[EPHEM_L0],360.0)-L;
01370     while (fabs(tb-ta) > 0.01) {
01371     t = ta - (tb - ta) * la / (lb - la);
01372     sun_ephemeris(t, ephem, 0.0, 0.0);
01373     l = fmod(ephem[EPHEM_L0],360.0)-L;
01374     ta = tb;
01375     la = lb;
01376     tb = t;
01377     lb = l;
01378     }
01379 
01380     return t;
01381 }
01382 
01383 double earth_B(TIME t)
01384 {
01385     double ephem[EPHEM_SIZE];
01386     sun_ephemeris(t, ephem, 0.0, 0.0);
01387     return ephem[EPHEM_B0]*180.0/M_PI;
01388 }
01389 
01390 */
01391 
01392 /*
01393 $Id: hmisynop.c,v 1.3 2014/02/10 21:24:47 yliu Exp $
01394 $Source: /home/cvsuser/cvsroot/JSOC/proj/mag/synop/apps/hmisynop.c,v $
01395 $Author: yliu $
01396 */
01397 /* $Log: hmisynop.c,v $
01398 /* Revision 1.3  2014/02/10 21:24:47  yliu
01399 /* Remove several include-lines Feb. 10, 2014
01400 /*
01401 /* Revision 1.2  2011/11/20 01:27:32  xudong
01402 /* updated fstats.h
01403 /*
01404 /* Revision 1.1  2011/04/15 01:44:09  xudong
01405 /* for yang
01406 /*
01407  * Revision 1.24  2007/10/26 17:51:39  arta
01408  * Fix bug where for loop limit was changed within loop.
01409  *
01410  * Revision 1.23  2007/10/18 16:45:58  arta
01411  * Fix bug where there was an attempt to read pixels just off the left or right edge of a remapped image when center was either 60 or -60.
01412  *
01413  * Revision 1.22  2007/09/27 17:34:25  arta
01414  * Change name of degrees per halfwindow parameter to HWNWIDTH.
01415  *
01416  * Revision 1.21  2007/09/27 17:12:28  arta
01417  * Change the W_HWIDTH attribute type to a float, to match the variable used to populate it.
01418  *
01419  * Revision 1.20  2007/09/26 16:06:41  arta
01420  * Make log more accurately show the stack of mags for each synoptic chart column.
01421  *
01422  * Revision 1.19  2007/09/26 00:28:38  arta
01423  * Use dist from CM, not dist from mid col.
01424  *
01425  * Revision 1.18  2007/08/20 16:59:41  arta
01426  * Add output to help identify data set name.  Fix sensAdj detection.
01427  *
01428  * Revision 1.17  2007/08/10 23:58:00  arta
01429  * Write out carr stretch keywords.
01430  *
01431  * Revision 1.16  2007/08/09 15:10:49  arta
01432  * Do proper check for out-of-bounds array index.
01433  *
01434  * Revision 1.15  2007/08/06 18:49:03  arta
01435  * Fix for crash in sgi4.
01436  *
01437  * Revision 1.14  2007/07/18  14:53:27  arta
01438  * Improve performance.  Use merge-sort, not quicksort, as data were essentially read in into two sorted lists.  Decrease memory footprint - do synop columen the data are ready, don't wait until the end.  Clear memory right after each column is done.
01439  *
01440  * Revision 1.13  2007/07/16 21:01:49  arta
01441  * new outlier detection.  Calculate noise at various cos(rho).
01442  *
01443  * Revision 1.12  2007/07/11 16:30:56  arta
01444  * magsynop - add line-of-sight flag (to output line-of-sight images in addition to radial ones), add force flag (to allow processing of obsolete data), fix bug in radial image detection
01445  * fixplatescale_v0 - write out soi_vers_num in BLDVER18 attribute
01446  * fdmagcal and fdradial - set pixels outside of rsun to missing (do not leave them unmodified), do not write BLDVER18 in fdmagcal.
01447  *
01448  * Revision 1.11  2007/05/14 18:44:49  arta
01449  * Ensure v2helio retains the BLDVER18 keyword for magsynop use.
01450  *
01451  * Revision 1.10  2007/05/14 16:57:39  arta
01452  * Use BLDVER18, not BLDVER15 when assessing how current data are.
01453  *
01454  * Revision 1.9  2007/05/14 05:10:19  arta
01455  * Fix sgi build.
01456  *
01457  * Revision 1.8  2007/05/08 23:39:18  arta
01458  * Check for consistent MDI CM build in magsynop.  Requires passing along keywords in v2helio.  Also, ensure that all images are either radial or line-of-sight, but not a mixture.
01459  *
01460  * Revision 1.7  2007/05/08 17:23:40  arta
01461  * Store epts file in a more compact form.
01462  *
01463  * Revision 1.6  2007/05/05 03:31:05  arta
01464  * Cap values larger than can be contained by max value.
01465  *
01466  * Revision 1.5  2007/03/26 19:09:49  arta
01467  * Track magnetogram row, col, dataseries, recordnum info to diagnose missing data.
01468  *
01469  * Revision 1.4  2007/03/16 16:35:10  arta
01470  * Modify algorithm to use the BEST magnetograms for each synoptic column.  BEST is closest to CM.  Use as many magnetograms as required to reach 20 1-minute magnetogram equivalence points.  Crop magnetograms at +/- 20 of CM.
01471  *
01472  * Revision 1.3  2007/01/26 18:25:28  rick
01473  * changed name of earth_merid... to local_earth_merid to avoid lib conflict
01474  *
01475  * Revision 1.2  2000/05/30 16:44:14  kehcheng
01476  * initial revision
01477  * */