(file) Return to update_sharp_keys.c CVS log (file) (dir) Up to [Development] / JSOC / proj / sharp / apps

Diff for /JSOC/proj/sharp/apps/update_sharp_keys.c between version 1.8 and 1.19

version 1.8, 2014/06/05 21:27:32 version 1.19, 2021/05/26 19:25:17
Line 1 
Line 1 
 /* /*
  *      MODULE NAME: update_sharp_keys.c  *      MODULE NAME: update_sharp_keys.c
  *  *
  *      DESCRIPTION: This module recalculates SHARP keywords. This is accomplished by   *          DESCRIPTION: This module recalculates SHARP keywords.
  *      cloning a record, recalculating keywords of choice (i.e., user input),   *      This is accomplished by cloning a record, recalculating keywords of choice
  *      and pointing to the same segments as the old record. Associated error keys   *      (i.e., user input), and pointing to the same segments as the old record.
  *      are computed by default.   *      Associated error keys are computed by default.
  *  *
  *      This module accounts for versioning by appending to the keyword CODEVER7 and HISTORY:  *      This module accounts for versioning by appending to the keyword CODEVER7 and HISTORY:
  *      CODEVER7 will contain multiple lines: the production build of sharp.c, the  *      CODEVER7 will contain multiple lines: the production build of sharp.c, the
Line 12 
Line 12 
  *      of update_sharp_keys.c. HISTORY will include a human-readable sentence about which  *      of update_sharp_keys.c. HISTORY will include a human-readable sentence about which
  *      keywords were updated.  *      keywords were updated.
  *  *
    *      N.B.: This module calculates the keyword specified in keylist and then does a
    *      drms_copykeys() to copy over all the remaining keywords. However, if [1] a keyword
    *      does not exist in the .jsd files of the input series, and [2] said keyword X is not
    *      specified in keylist, then [3] drms_copykeys() will copy DRMS_MISSING_VALUE into
    *      keyword X for all the output series.
    *
  *      This module does not produce segments.  *      This module does not produce segments.
  *  *
  *      INPUTS     : -- DRMS SHARP series  *      INPUTS     : -- DRMS SHARP series
  *                   -- DRMS SHARP CEA series  *                   -- DRMS SHARP CEA series
  *                   -- HARPNUM  *                   -- HARPNUM
  *                   -- comma separated list of keywords to recalculate  *                   -- comma separated list of keywords to recalculate
    *                   -- DEBUG flag
  *  *
  *      AUTHOR     : Monica Bobra  *      AUTHOR     : Monica Bobra
  *  *
  *      Version    :   v0.0     Jun 14 2013  
  *                     v0.1     Feb 11 2014 Added different input and output series  
  *  
  *      EXAMPLE    :  *      EXAMPLE    :
  *      update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //  *      update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //
  *      HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT   *      HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT debug=debug
  *  *
  */  */
  
Line 59 
Line 63 
 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0])) #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);} #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
 #define SHOW(msg) {printf("%s", msg); fflush(stdout);} #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
   #define DEBUG(true)
  
 char *module_name = "update_sharp_keys";  /* Module name */ char *module_name = "update_sharp_keys";  /* Module name */
 char *version_id  = "2014 Feb 12";        /* Version number */  char *version_id  = "2021 May 25";        /* Version number */
  
 ModuleArgs_t module_args[] = ModuleArgs_t module_args[] =
 { {
Line 70  ModuleArgs_t module_args[] =
Line 75  ModuleArgs_t module_args[] =
         {ARG_STRING, "sharpseriesout",     NULL, "Output Sharp dataseries"},         {ARG_STRING, "sharpseriesout",     NULL, "Output Sharp dataseries"},
         {ARG_STRING, "sharpceaseriesout",  NULL, "Output Sharp CEA dataseries"},         {ARG_STRING, "sharpceaseriesout",  NULL, "Output Sharp CEA dataseries"},
         {ARG_INT,    "HARPNUM",         NULL, "HARP number"},         {ARG_INT,    "HARPNUM",         NULL, "HARP number"},
         {ARG_STRING, "keylist",         NULL, "comma separated list of keywords to update"},          {ARG_STRING, "keylist",            NULL, "comma separated list of keywords to update (for LOS keywords enter MEANGBZLOS or USFLUXLOS)"},
           {ARG_STRING, "debug",              NULL, "debug mode functionality. use like this: debug=debug"},
         {ARG_END}         {ARG_END}
 }; };
  
Line 85  int DoIt(void)
Line 91  int DoIt(void)
  
         /* Keywords */         /* Keywords */
         float mean_vf;         float mean_vf;
       float mean_vf_los;
         float absFlux;         float absFlux;
       float absFlux_los;
         float count_mask;         float count_mask;
       float count_mask_los;
         float mean_hf;         float mean_hf;
         float mean_gamma;         float mean_gamma;
         float mean_derivative_btotal;         float mean_derivative_btotal;
         float mean_derivative_bh;         float mean_derivative_bh;
         float mean_derivative_bz;         float mean_derivative_bz;
           float mean_derivative_los;
         float mean_jz;         float mean_jz;
         float us_i;         float us_i;
         float mean_alpha;         float mean_alpha;
Line 143  int DoIt(void)
Line 153  int DoIt(void)
  
         char sharpquery[256];         char sharpquery[256];
         char sharpceaquery[256];         char sharpceaquery[256];
         printf("sameflag = %d\n", sameflag);  
         printf("testflag = %d\n", testflag);  
         sprintf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);         sprintf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
         printf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);         printf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
  
Line 185  int DoIt(void)
Line 193  int DoIt(void)
                 DIE("Problem creating sharp cea records.");                 DIE("Problem creating sharp cea records.");
         }         }
  
   
         char *keylist = (char *) params_get_str(&cmdparams, "keylist");         char *keylist = (char *) params_get_str(&cmdparams, "keylist");
           char *debug   = (char *) params_get_str(&cmdparams, "debug");
  
         // Flags to indicate which keyword will be recalculated         // Flags to indicate which keyword will be recalculated
         int meanvfflag  = (strstr(keylist,"USFLUX")  != NULL);  // generalize so that lowercase is acceptable      int meanvflosflag = (strstr(keylist,"USFLUXL")  != NULL);
       int meanglosflag  = (strstr(keylist,"MEANGBL") != NULL);
       int meanvfflag    = (strstr(keylist,"USFLUX")  != NULL);
         int totpotflag  = (strstr(keylist,"TOTPOT")  != NULL);         int totpotflag  = (strstr(keylist,"TOTPOT")  != NULL);
         int meangamflag = (strstr(keylist,"MEANGAM") != NULL);         int meangamflag = (strstr(keylist,"MEANGAM") != NULL);
         int meangbtflag = (strstr(keylist,"MEANGBT") != NULL);         int meangbtflag = (strstr(keylist,"MEANGBT") != NULL);
Line 213  int DoIt(void)
Line 223  int DoIt(void)
         int epsxflag    = (strstr(keylist,"EPSX")    != NULL);         int epsxflag    = (strstr(keylist,"EPSX")    != NULL);
         int epsyflag    = (strstr(keylist,"EPSY")    != NULL);         int epsyflag    = (strstr(keylist,"EPSY")    != NULL);
         int epszflag    = (strstr(keylist,"EPSZ")    != NULL);         int epszflag    = (strstr(keylist,"EPSZ")    != NULL);
       int debugflag     = (strstr(debug,"debug")     != NULL);
  
         DRMS_Record_t *sharpinrec = sharpinrecset->records[0];          for (irec=0;irec<nrecs;irec++)
         DRMS_Record_t *sharpceainrec = sharpceainrecset->records[0];          {
           DRMS_Record_t *sharpinrec = sharpinrecset->records[irec];
               DRMS_Record_t *sharpceainrec = sharpceainrecset->records[irec];
         DRMS_Segment_t *inseg = drms_segment_lookup(sharpceainrec, "Br");         DRMS_Segment_t *inseg = drms_segment_lookup(sharpceainrec, "Br");
         int nx = inseg->axis[0];         int nx = inseg->axis[0];
         int ny = inseg->axis[1];         int ny = inseg->axis[1];
Line 244  int DoIt(void)
Line 257  int DoIt(void)
         float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));         float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));
         float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));         float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
         float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));         float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
               float *err_term1 = (float *) (calloc(nxny, sizeof(float)));
               float *err_term2 = (float *) (calloc(nxny, sizeof(float)));
           float *fx        = (float *) (malloc(nxny * sizeof(float)));
           float *fy        = (float *) (malloc(nxny * sizeof(float)));
           float *fz        = (float *) (malloc(nxny * sizeof(float)));
           float *derx_los  = (float *) (malloc(nxny * sizeof(float)));
               float *dery_los  = (float *) (malloc(nxny * sizeof(float)));
  
         // ephemeris variables         // ephemeris variables
         float  cdelt1_orig, cdelt1, dsun_obs, imcrpix1, imcrpix2, crpix1, crpix2;         float  cdelt1_orig, cdelt1, dsun_obs, imcrpix1, imcrpix2, crpix1, crpix2;
Line 269  int DoIt(void)
Line 289  int DoIt(void)
         float UNIX_epoch = -220924792.000;         float UNIX_epoch = -220924792.000;
         sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);         sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
         sprintf(historyofthemodule,"The following keywords were re-computed on %s: %s.",timebuf,keylist);         sprintf(historyofthemodule,"The following keywords were re-computed on %s: %s.",timebuf,keylist);
         printf("historyofthemodule=%s\n",historyofthemodule);              //printf("historyofthemodule=%s\n",historyofthemodule);
         history0    = drms_getkey_string(sharpinrec, "HISTORY", &status);         history0    = drms_getkey_string(sharpinrec, "HISTORY", &status);
         strcat(historyofthemodule,"\n");         strcat(historyofthemodule,"\n");
         strcat(historyofthemodule,history0);         strcat(historyofthemodule,history0);
Line 280  int DoIt(void)
Line 300  int DoIt(void)
         dsun_obs = drms_getkey_float(testrec, "DSUN_OBS",   &status);         dsun_obs = drms_getkey_float(testrec, "DSUN_OBS",   &status);
         rsun_ref = drms_getkey_double(testrec, "RSUN_REF", &status);         rsun_ref = drms_getkey_double(testrec, "RSUN_REF", &status);
         cdelt1_orig = drms_getkey_float(testrec, "CDELT1",   &status);         cdelt1_orig = drms_getkey_float(testrec, "CDELT1",   &status);
         cdelt1      = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);  
         int scale = round(2.0/cdelt1);         int scale = round(2.0/cdelt1);
         int nx1 = nx/scale;         int nx1 = nx/scale;
         int ny1 = ny/scale;         int ny1 = ny/scale;
         int nxp = nx1+40;         int nxp = nx1+40;
         int nyp = ny1+40;         int nyp = ny1+40;
  
           // check to make sure that the dimensions passed into fsample() are adequate
         if (nx1 > floor((nx-1)/scale + 1) )         if (nx1 > floor((nx-1)/scale + 1) )
                 DIE("X-dimension of output array in fsample() is too large.");                 DIE("X-dimension of output array in fsample() is too large.");
         if (ny1 > floor((ny-1)/scale + 1) )         if (ny1 > floor((ny-1)/scale + 1) )
                 DIE("Y-dimension of output array in fsample() is too large.");                 DIE("Y-dimension of output array in fsample() is too large.");
   
         // malloc some arrays for the R parameter calculation         // malloc some arrays for the R parameter calculation
         float *rim     = (float *)malloc(nx1*ny1*sizeof(float));         float *rim     = (float *)malloc(nx1*ny1*sizeof(float));
         float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));         float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));
Line 303  int DoIt(void)
Line 322  int DoIt(void)
         float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));         float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));
         float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));         float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));
  
         // define some arrays for the lorentz force calculation  
         float *fx = (float *) (malloc(nxny * sizeof(float)));  
         float *fy = (float *) (malloc(nxny * sizeof(float)));  
         float *fz = (float *) (malloc(nxny * sizeof(float)));  
   
   
         for (irec=0;irec<nrecs;irec++)  
         {  
            // Get emphemeris            // Get emphemeris
            sharpinrec    = sharpinrecset->records[irec];            sharpinrec    = sharpinrecset->records[irec];
            sharpceainrec = sharpceainrecset->records[irec];            sharpceainrec = sharpceainrecset->records[irec];
Line 322  int DoIt(void)
Line 333  int DoIt(void)
            imcrpix2      = drms_getkey_float(sharpinrec, "IMCRPIX2", &status);            imcrpix2      = drms_getkey_float(sharpinrec, "IMCRPIX2", &status);
            crpix1        = drms_getkey_float(sharpinrec, "CRPIX1", &status);            crpix1        = drms_getkey_float(sharpinrec, "CRPIX1", &status);
            crpix2        = drms_getkey_float(sharpinrec, "CRPIX2", &status);            crpix2        = drms_getkey_float(sharpinrec, "CRPIX2", &status);
           cdelt1      = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
  
            sharpoutrec = sharpoutrecset->records[irec];            sharpoutrec = sharpoutrecset->records[irec];
            sharpceaoutrec = sharpceaoutrecset->records[irec];            sharpceaoutrec = sharpceaoutrecset->records[irec];
Line 331  int DoIt(void)
Line 343  int DoIt(void)
              drms_copykeys(sharpceaoutrec, sharpceainrec, 0, kDRMS_KeyClass_Explicit);              drms_copykeys(sharpceaoutrec, sharpceainrec, 0, kDRMS_KeyClass_Explicit);
            }            }
  
            TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */             TIME trec, tnow; /* 1970.01.01_00:00:00_UTC */
            tnow = (double)time(NULL);            tnow = (double)time(NULL);
            tnow += UNIX_epoch;            tnow += UNIX_epoch;
  
              if (debugflag)
              {
                 printf("i=%d\n",irec);
              }
   
            // set CODEVER7 and HISTORY and DATE            // set CODEVER7 and HISTORY and DATE
            drms_setkey_string(sharpoutrec, "CODEVER7", cvsinfoall);            drms_setkey_string(sharpoutrec, "CODEVER7", cvsinfoall);
            drms_setkey_string(sharpoutrec,"HISTORY",historyofthemodule);            drms_setkey_string(sharpoutrec,"HISTORY",historyofthemodule);
Line 343  int DoIt(void)
Line 360  int DoIt(void)
            drms_setkey_time(sharpoutrec,"DATE",tnow);            drms_setkey_time(sharpoutrec,"DATE",tnow);
            drms_setkey_time(sharpceaoutrec,"DATE",tnow);            drms_setkey_time(sharpceaoutrec,"DATE",tnow);
  
            // Get bx, by, bz, mask              // Get all the segments if computing vector magnetic field keywords
            // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities  
            DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");            DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");
            DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);            DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
            int *bitmask = (int *) bitmaskArray->data;    // get the previously made mask array            int *bitmask = (int *) bitmaskArray->data;    // get the previously made mask array
  
            //Use conf_disambig map as a threshold on spaceweather quantities  
            DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceainrec, "conf_disambig");            DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceainrec, "conf_disambig");
            DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);            DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
            if (status != DRMS_SUCCESS)            if (status != DRMS_SUCCESS)
               DIE("No CONF_DISAMBIG segment.");               DIE("No CONF_DISAMBIG segment.");
            int *mask = (int *) maskArray->data;         // get the previously made mask array            int *mask = (int *) maskArray->data;         // get the previously made mask array
  
            //Use magnetogram map to compute R  
            DRMS_Segment_t *losSeg = drms_segment_lookup(sharpceainrec, "magnetogram");            DRMS_Segment_t *losSeg = drms_segment_lookup(sharpceainrec, "magnetogram");
            DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);            DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
               if (status != DRMS_SUCCESS)
               DIE("No magnetogram segment.");
            float *los = (float *) losArray->data;               // los            float *los = (float *) losArray->data;               // los
               //printf("Line 387");
  
            DRMS_Segment_t *bxSeg = drms_segment_lookup(sharpceainrec, "Bp");            DRMS_Segment_t *bxSeg = drms_segment_lookup(sharpceainrec, "Bp");
            DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);            DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
               if (status != DRMS_SUCCESS)
               DIE("No Bp segment.");
            float *bx = (float *) bxArray->data;         // bx            float *bx = (float *) bxArray->data;         // bx
               //printf("Line 393");
  
            DRMS_Segment_t *bySeg = drms_segment_lookup(sharpceainrec, "Bt");            DRMS_Segment_t *bySeg = drms_segment_lookup(sharpceainrec, "Bt");
            DRMS_Array_t *byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status);            DRMS_Array_t *byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status);
               if (status != DRMS_SUCCESS)
               DIE("No Bt segment.");
            float *by = (float *) byArray->data;         // by            float *by = (float *) byArray->data;         // by
            for (int i = 0; i < nxny; i++) by[i] *= -1;            for (int i = 0; i < nxny; i++) by[i] *= -1;
               //printf("Line 399");
  
            DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceainrec, "Br");            DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceainrec, "Br");
               //printf("Line 408\n");
               if (bzSeg == NULL)
               DIE("somekind of pointer problem.");
   
            DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);            DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
               //printf("status=%d\n",status);
               if (status != DRMS_SUCCESS)
               DIE("No Br segment.");
            float *bz = (float *) bzArray->data;         // bz            float *bz = (float *) bzArray->data;         // bz
               //printf("Line 414");
  
            DRMS_Segment_t *bz_errSeg = drms_segment_lookup(sharpceainrec, "Br_err");            DRMS_Segment_t *bz_errSeg = drms_segment_lookup(sharpceainrec, "Br_err");
            DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);            DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
               if (status != DRMS_SUCCESS)
               DIE("No Br_err segment.");
            float *bz_err = (float *) bz_errArray->data;         // bz_err            float *bz_err = (float *) bz_errArray->data;         // bz_err
               //printf("Line 408");
  
            DRMS_Segment_t *by_errSeg = drms_segment_lookup(sharpceainrec, "Bt_err");            DRMS_Segment_t *by_errSeg = drms_segment_lookup(sharpceainrec, "Bt_err");
            DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);            DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
               if (status != DRMS_SUCCESS)
               DIE("No Bt_err segment.");
            float *by_err = (float *) by_errArray->data;         // by_err            float *by_err = (float *) by_errArray->data;         // by_err
               //printf("Line 412");
  
            DRMS_Segment_t *bx_errSeg = drms_segment_lookup(sharpceainrec, "Bp_err");            DRMS_Segment_t *bx_errSeg = drms_segment_lookup(sharpceainrec, "Bp_err");
            DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);            DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
               if (status != DRMS_SUCCESS)
               DIE("No Bp_err segment.");
            float *bx_err = (float *) bx_errArray->data;         // bx_err            float *bx_err = (float *) bx_errArray->data;         // bx_err
  
            /***** USFLUX, Example Function 1 *************************************/            /***** USFLUX, Example Function 1 *************************************/
Line 419  int DoIt(void)
Line 459  int DoIt(void)
  
            drms_setkey_float(sharpoutrec, "R_VALUE",  Rparam);            drms_setkey_float(sharpoutrec, "R_VALUE",  Rparam);
            drms_setkey_float(sharpceaoutrec, "R_VALUE", Rparam);            drms_setkey_float(sharpceaoutrec, "R_VALUE", Rparam);
           }
  
           /***** MEANGBZ, Example Function 7 ************************************/
           if (meangbzflag)
           {
               // Compute Bz derivative
               if (computeBzderivative(bz, bz_err, dims, &mean_derivative_bz, &mean_derivative_bz_err, mask, bitmask, derx_bz, dery_bz, err_term1, err_term2))
               {
                   mean_derivative_bz     = DRMS_MISSING_FLOAT;
                   mean_derivative_bz_err = DRMS_MISSING_FLOAT;
            }            }
  
            /***** MEANPOT and TOTPOT, Example Function 13  (Requires Keiji's code) **/              // Set sharp keys
               drms_setkey_float(sharpoutrec, "MEANGBZ",  mean_derivative_bz);
               drms_setkey_float(sharpoutrec, "ERRBZ",  mean_derivative_bz_err);
   
               // Set sharp cea keys
               drms_setkey_float(sharpceaoutrec, "MEANGBZ", mean_derivative_bz);
               drms_setkey_float(sharpceaoutrec, "ERRBZ", mean_derivative_bz_err);
           }
  
           /***** MEANPOT and TOTPOT, Example Function 13  (Requires Keiji's code) **/
            if (totpotflag || meanpotflag)            if (totpotflag || meanpotflag)
            {            {
               // First compute potential field               // First compute potential field
Line 455  int DoIt(void)
Line 512  int DoIt(void)
            }            }
  
            /***** MEANGAM, Example Function 3 ************************************/            /***** MEANGAM, Example Function 3 ************************************/
   
            if (meangamflag)            if (meangamflag)
            {            {
               // First compute horizontal field               // First compute horizontal field
Line 478  int DoIt(void)
Line 534  int DoIt(void)
            }            }
  
            /***** MEANGBT, Example Function 5 (Requires Function 4) *************/            /***** MEANGBT, Example Function 5 (Requires Function 4) *************/
   
            if (meangbtflag)            if (meangbtflag)
            {            {
               // First, compute Bt               // First, compute Bt
Line 486  int DoIt(void)
Line 541  int DoIt(void)
  
               // Then take the derivative of Bt               // Then take the derivative of Bt
               if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask, bitmask, derx_bt, dery_bt, bt_err,               if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask, bitmask, derx_bt, dery_bt, bt_err,
                                           &mean_derivative_btotal_err))                                          &mean_derivative_btotal_err, err_term1, err_term2))
               {               {
                  mean_derivative_btotal     = DRMS_MISSING_FLOAT;                  mean_derivative_btotal     = DRMS_MISSING_FLOAT;
                  mean_derivative_btotal_err = DRMS_MISSING_FLOAT;                  mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
Line 502  int DoIt(void)
Line 557  int DoIt(void)
            }            }
  
            /***** MEANGBH, Example Function 6 (Requires Function 2) *************/            /***** MEANGBH, Example Function 6 (Requires Function 2) *************/
   
            if (meangbhflag)            if (meangbhflag)
            {            {
               // First, compute Bh               // First, compute Bh
               computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask);               computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask);
  
               // Then take the derivative of Bh               // Then take the derivative of Bh
               if (computeBhderivative(bh, bh_err, dims, &mean_derivative_bh, &mean_derivative_bh_err, mask, bitmask, derx_bh, dery_bh))                  if (computeBhderivative(bh, bh_err, dims, &mean_derivative_bh, &mean_derivative_bh_err, mask, bitmask, derx_bh, dery_bh, err_term1, err_term2))
               {               {
                  mean_derivative_bh       = DRMS_MISSING_FLOAT;                  mean_derivative_bh       = DRMS_MISSING_FLOAT;
                  mean_derivative_bh_err   = DRMS_MISSING_FLOAT;                  mean_derivative_bh_err   = DRMS_MISSING_FLOAT;
Line 524  int DoIt(void)
Line 578  int DoIt(void)
               drms_setkey_float(sharpceaoutrec, "ERRBH", mean_derivative_bh_err);               drms_setkey_float(sharpceaoutrec, "ERRBH", mean_derivative_bh_err);
            }            }
  
            /***** MEANGBZ, Example Function 7 ************************************/  
   
            if (meangbzflag)  
            {  
               // Compute Bz derivative  
               if (computeBzderivative(bz, bz_err, dims, &mean_derivative_bz, &mean_derivative_bz_err, mask, bitmask, derx_bz, dery_bz))  
               {  
                  mean_derivative_bz     = DRMS_MISSING_FLOAT;  
                  mean_derivative_bz_err = DRMS_MISSING_FLOAT;  
               }  
   
               // Set sharp keys  
               drms_setkey_float(sharpoutrec, "MEANGBZ",  mean_derivative_bz);  
               drms_setkey_float(sharpoutrec, "ERRBZ",  mean_derivative_bz_err);  
   
               // Set sharp cea keys  
               drms_setkey_float(sharpceaoutrec, "MEANGBZ", mean_derivative_bz);  
               drms_setkey_float(sharpceaoutrec, "ERRBZ", mean_derivative_bz_err);  
            }  
   
            /***** MEANJZD and TOTUSJZ, Example Function 9 (Requires Function 8) ***/            /***** MEANJZD and TOTUSJZ, Example Function 9 (Requires Function 8) ***/
   
            if (meanjzdflag || totusjzflag)            if (meanjzdflag || totusjzflag)
            {            {
               // First, compute Jz on all pixels               // First, compute Jz on all pixels
               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
                      derx, dery);                        derx, dery, err_term1, err_term2);
  
               // Then take the sums of the appropriate pixels               // Then take the sums of the appropriate pixels
               if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &mean_jz,               if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &mean_jz,
                        &mean_jz_err, &us_i, &us_i_err, mask, bitmask, cdelt1,                        &mean_jz_err, &us_i, &us_i_err, mask, bitmask, cdelt1,
                        rsun_ref, rsun_obs, derx, dery))                        rsun_ref, rsun_obs, derx, dery))
  
   
               {               {
                  mean_jz                = DRMS_MISSING_FLOAT;                  mean_jz                = DRMS_MISSING_FLOAT;
                  us_i                   = DRMS_MISSING_FLOAT;                  us_i                   = DRMS_MISSING_FLOAT;
Line 579  int DoIt(void)
Line 611  int DoIt(void)
            }            }
  
            /***** MEANALP, Example Function 10 (Requires Function 8)*********/            /***** MEANALP, Example Function 10 (Requires Function 8)*********/
   
            if (meanalpflag)            if (meanalpflag)
            {            {
               // First, compute Jz on all pixels               // First, compute Jz on all pixels
               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
                         derx, dery);                        derx, dery, err_term1, err_term2);
  
               // Then compute alpha quantities on select pixels               // Then compute alpha quantities on select pixels
               if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &mean_alpha, &mean_alpha_err,               if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &mean_alpha, &mean_alpha_err,
Line 605  int DoIt(void)
Line 636  int DoIt(void)
            }            }
  
            /***** MEANJZH, TOTUSJH, ABSNJZH, Example Function 11 (Requires Function 8) ***/            /***** MEANJZH, TOTUSJH, ABSNJZH, Example Function 11 (Requires Function 8) ***/
   
            if (meanjzhflag || totusjhflag || absnjzhflag)            if (meanjzhflag || totusjhflag || absnjzhflag)
            {            {
               // First, compute Jz on all pixels               // First, compute Jz on all pixels
               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
                         derx, dery);                        derx, dery, err_term1, err_term2);
  
               // Then compute helicity quantities on select pixels               // Then compute helicity quantities on select pixels
            if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &mean_ih, &mean_ih_err, &total_us_ih,            if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &mean_ih, &mean_ih_err, &total_us_ih,
Line 642  int DoIt(void)
Line 672  int DoIt(void)
            }            }
  
            /***** SAVNCPP, Example Function 12 (Requires Function 8) *******************/            /***** SAVNCPP, Example Function 12 (Requires Function 8) *******************/
   
            if (savncppflag)            if (savncppflag)
            {            {
               // First, compute Jz on all pixels               // First, compute Jz on all pixels
               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
                         derx, dery);                        derx, dery, err_term1, err_term2);
  
               // Then compute sums of Jz on select pixels               // Then compute sums of Jz on select pixels
               if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &totaljz, &totaljz_err,               if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &totaljz, &totaljz_err,
Line 667  int DoIt(void)
Line 696  int DoIt(void)
            }            }
  
            /***** MEANSHR and SHRGT45, Example Function 14 (Requires Keiji's code) **********/            /***** MEANSHR and SHRGT45, Example Function 14 (Requires Keiji's code) **********/
   
            if (meanshrflag || shrgt45flag)            if (meanshrflag || shrgt45flag)
            {            {
               // First compute potential field               // First compute potential field
Line 696  int DoIt(void)
Line 724  int DoIt(void)
            }            }
  
           /***** TOTFX, TOTFY, TOTFX, TOTBSQ, EPSX, EPSY, EPSZ , Example Function 16 (Lorentz forces) **********/           /***** TOTFX, TOTFY, TOTFX, TOTBSQ, EPSX, EPSY, EPSZ , Example Function 16 (Lorentz forces) **********/
   
            if (totfxflag || totfyflag || totfzflag || totbsqflag || epsxflag || epsyflag || epszflag)            if (totfxflag || totfyflag || totfzflag || totbsqflag || epsxflag || epsyflag || epszflag)
            {            {
   
                 // Compute Lorentz forces                 // Compute Lorentz forces
                 if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &totfx, &totfy, &totfz, &totbsq,                 if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &totfx, &totfy, &totfz, &totbsq,
                     &epsx, &epsy, &epsz, mask, bitmask, cdelt1, rsun_ref, rsun_obs))                     &epsx, &epsy, &epsz, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
Line 733  int DoIt(void)
Line 759  int DoIt(void)
  
            }            }
  
            /******************************* END FLAGS **********************************/          /***** USFLUX, Example Function 17 *************************************/
               if (meanvflosflag)
               {
                  //printf("772");
                  if (computeAbsFlux_los(los, dims, &absFlux_los, &mean_vf_los,
                              &count_mask_los, bitmask, cdelt1, rsun_ref, rsun_obs))
                  {
                   absFlux_los = DRMS_MISSING_FLOAT;               // If fail, fill in NaN
                   mean_vf_los = DRMS_MISSING_FLOAT;
                   count_mask_los = DRMS_MISSING_INT;
                  }
               //printf("line757\n");
               //printf("mean_vf_los=%f\n",mean_vf_los);
               //printf("cdelt1=%f\n",cdelt1);
               //printf("rsun_ref=%f\n",rsun_ref);
               //printf("rsun_obs=%f\n",rsun_obs);
               //printf("Key #%d done.\n", irec);
               drms_setkey_float(sharpoutrec, "USFLUXL",  mean_vf_los);
               drms_setkey_float(sharpoutrec, "CMASKL",   count_mask_los);
   
               drms_setkey_float(sharpceaoutrec, "USFLUXL",  mean_vf_los);
               drms_setkey_float(sharpceaoutrec, "CMASKL",   count_mask_los);
           }
   
           /***** MEANGBZ, Example Function 18 ************************************/
           if (meanglosflag)
           {
           //printf("line 794\n");
               // Compute Bz derivative
               if (computeLOSderivative(los, dims, &mean_derivative_los, bitmask, derx_los, dery_los))
               {
                   //printf("line 799\n");
                   mean_derivative_los     = DRMS_MISSING_FLOAT;
               }
   
               // Set sharp keys
               //printf("mean_derivative_los=%f\n",mean_derivative_los);
               drms_setkey_float(sharpoutrec, "MEANGBL",  mean_derivative_los);
  
               // Set sharp cea keys
               drms_setkey_float(sharpceaoutrec, "MEANGBL", mean_derivative_los);
           }
   
           /******************************* END FLAGS **********************************/
           //printf("line809\n");
            drms_free_array(bitmaskArray);            drms_free_array(bitmaskArray);
            drms_free_array(maskArray);            drms_free_array(maskArray);
            drms_free_array(bxArray);            drms_free_array(bxArray);
Line 744  int DoIt(void)
Line 813  int DoIt(void)
            drms_free_array(by_errArray);            drms_free_array(by_errArray);
            drms_free_array(bz_errArray);            drms_free_array(bz_errArray);
            drms_free_array(losArray);            drms_free_array(losArray);
         } //endfor  
   
         drms_close_records(sharpinrecset, DRMS_FREE_RECORD);  
         drms_close_records(sharpceainrecset, DRMS_FREE_RECORD);  
   
         drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD);  
         drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD);  
   
         //used for select calculations  
         free(bh); free(bt); free(jz);  
         free(bpx); free(bpy); free(bpz);  
         free(derx); free(dery);  
         free(derx_bt); free(dery_bt);  
         free(derx_bz); free(dery_bz);  
         free(derx_bh); free(dery_bh);  
         free(bt_err); free(bh_err);  free(jz_err);  
         free(jz_err_squared); free(jz_rms_err);  
         free(cvsinfoall);         free(cvsinfoall);
   
         free(rim);         free(rim);
         free(p1p0);         free(p1p0);
         free(p1n0);         free(p1n0);
Line 772  int DoIt(void)
Line 823  int DoIt(void)
         free(pmap);         free(pmap);
         free(p1pad);         free(p1pad);
         free(pmapn);         free(pmapn);
   
         // free the arrays that are related to the lorentz calculation  
         free(fx); free(fy); free(fz);         free(fx); free(fy); free(fz);
       free(bh); free(bt); free(jz);
       free(bpx); free(bpy); free(bpz);
       free(derx); free(dery);
       free(derx_bt); free(dery_bt);
       free(derx_los); free(dery_los);
       free(derx_bz); free(dery_bz);
       free(derx_bh); free(dery_bh);
       free(bt_err); free(bh_err);  free(jz_err);
       free(jz_err_squared); free(jz_rms_err);
       free(jz_err_squared_smooth);
       free(jz_smooth);
       free(err_term2);
       free(err_term1);
   
           } //endfor
   
       // Close all the records
           drms_close_records(sharpinrecset, DRMS_FREE_RECORD);
           drms_close_records(sharpceainrecset, DRMS_FREE_RECORD);
   
           drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD);
           drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD);
   
           printf("Success=%d\n",sameflag);
  
 return 0; return 0;
   
 }       // DoIt }       // DoIt


Legend:
Removed from v.1.8  
changed lines
  Added in v.1.19

Karen Tian
Powered by
ViewCVS 0.9.4