(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.1 and 1.19

version 1.1, 2013/07/04 02:20:46 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  
  *  
  *      EXAMPLE    :  *      EXAMPLE    :
  *      update_sharp_keys sharpseries=hmi.sharp_720s sharpceaseries=hmi.sharp_cea_720s //   *      update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //
  *      HARPNUM=1 keylist=USFLUX,TOTPOT   *      HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT debug=debug
  *  *
  */  */
  
Line 58 
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  = "2013 Jun 14";        /* Version number */  char *version_id  = "2021 May 25";        /* Version number */
  
 ModuleArgs_t module_args[] = ModuleArgs_t module_args[] =
 { {
         {ARG_STRING, "sharpseries",     NULL, "Input/output Sharp dataseries"},          {ARG_STRING, "sharpseriesin",      NULL, "Input Sharp dataseries"},
         {ARG_STRING, "sharpceaseries",  NULL, "Input/output Sharp CEA dataseries"},          {ARG_STRING, "sharpceaseriesin",   NULL, "Input Sharp CEA dataseries"},
           {ARG_STRING, "sharpseriesout",     NULL, "Output Sharp 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 82  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 116  int DoIt(void)
Line 129  int DoIt(void)
         float totaljz_err;         float totaljz_err;
         float meanpot_err;         float meanpot_err;
         float totpot_err;         float totpot_err;
       float Rparam;
         float meanshear_angle_err;         float meanshear_angle_err;
         char *sharpseries = (char *) params_get_str(&cmdparams, "sharpseries");      float totfx;
         char *sharpceaseries = (char *) params_get_str(&cmdparams, "sharpceaseries");      float totfy;
       float totfz;
       float totbsq;
       float epsx;
       float epsy;
       float epsz;
   
           char *sharpseriesin = (char *) params_get_str(&cmdparams, "sharpseriesin");
           char *sharpceaseriesin = (char *) params_get_str(&cmdparams, "sharpceaseriesin");
           char *sharpseriesout = (char *) params_get_str(&cmdparams, "sharpseriesout");
           char *sharpceaseriesout = (char *) params_get_str(&cmdparams, "sharpceaseriesout");
   
       int sameflag = strcmp(sharpseriesin, sharpseriesout) == 0;
       int testflag = strcmp(sharpceaseriesin, sharpceaseriesout) == 0;
       if (testflag != sameflag)
       DIE("Either both outputs must be the same as their inputs, or both be different.");
  
         int harpnum = params_get_int(&cmdparams, "HARPNUM");         int harpnum = params_get_int(&cmdparams, "HARPNUM");
   
         char sharpquery[256];         char sharpquery[256];
         char sharpceaquery[256];         char sharpceaquery[256];
           sprintf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
           printf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
  
         sprintf(sharpceaquery, "%s[%d]\n", sharpceaseries,harpnum);          sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
         printf(sharpceaquery, "%s[%d]\n", sharpceaseries,harpnum);          sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
   
         sprintf(sharpquery, "%s[%d]\n", sharpseries,harpnum);  
         printf(sharpquery, "%s[%d]\n", sharpseries,harpnum);  
  
         DRMS_RecordSet_t *sharpinrecset = drms_open_records(drms_env, sharpquery, &status);         DRMS_RecordSet_t *sharpinrecset = drms_open_records(drms_env, sharpquery, &status);
         if (status != DRMS_SUCCESS)         if (status != DRMS_SUCCESS)
                 DIE("Problem opening sharp recordset.");                 DIE("Problem opening sharp recordset.");
         if (sharpinrecset->n == 0)          nrecs=sharpinrecset->n;
           if (nrecs == 0)
                 DIE("Sharp recordset contains no records.");                 DIE("Sharp recordset contains no records.");
         DRMS_RecordSet_t *sharpoutrecset = drms_clone_records(sharpinrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);          printf("nrecs=%d\n",nrecs);
         if (status != DRMS_SUCCESS)  
                 DIE("Problem cloning sharp records.");  
  
         DRMS_RecordSet_t *sharpceainrecset = drms_open_records(drms_env, sharpceaquery, &status);         DRMS_RecordSet_t *sharpceainrecset = drms_open_records(drms_env, sharpceaquery, &status);
         if (status != DRMS_SUCCESS)         if (status != DRMS_SUCCESS)
                 DIE("Problem opening sharp cea recordset.");                 DIE("Problem opening sharp cea recordset.");
         if (sharpceainrecset->n == 0)          if (sharpceainrecset->n != nrecs)
                 DIE("Sharp cea recordset contains no records.");                  DIE("Sharp cea recordset differs from sharp recordset.");
         DRMS_RecordSet_t *sharpceaoutrecset = drms_clone_records(sharpceainrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);  
           DRMS_RecordSet_t *sharpoutrecset, *sharpceaoutrecset;
           if (sameflag)
           {
           sharpoutrecset = drms_clone_records(sharpinrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
           if (status != DRMS_SUCCESS)
                   DIE("Problem cloning sharp records.");
           sharpceaoutrecset = drms_clone_records(sharpceainrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
         if (status != DRMS_SUCCESS)         if (status != DRMS_SUCCESS)
                 DIE("Problem cloning sharp cea records.");                 DIE("Problem cloning sharp cea records.");
           }
           else
           {
           sharpoutrecset = drms_create_records(drms_env, nrecs, sharpseriesout, DRMS_PERMANENT, &status);
           if (status != DRMS_SUCCESS)
                   DIE("Problem creating sharp records.");
           sharpceaoutrecset = drms_create_records(drms_env, nrecs, sharpceaseriesout, DRMS_PERMANENT, &status);
           if (status != DRMS_SUCCESS)
                   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 167  int DoIt(void)
Line 215  int DoIt(void)
         int meanpotflag = (strstr(keylist,"MEANPOT") != NULL);         int meanpotflag = (strstr(keylist,"MEANPOT") != NULL);
         int meanshrflag = (strstr(keylist,"MEANSHR") != NULL);         int meanshrflag = (strstr(keylist,"MEANSHR") != NULL);
         int shrgt45flag = (strstr(keylist,"SHRGT45") != NULL);         int shrgt45flag = (strstr(keylist,"SHRGT45") != NULL);
         DRMS_Record_t *sharpinrec = sharpinrecset->records[0];      int r_valueflag   = (strstr(keylist,"R_VALUE") != NULL);
         DRMS_Record_t *sharpceainrec = sharpceainrecset->records[0];      int totbsqflag    = (strstr(keylist,"TOTBSQ")  != NULL);
       int totfxflag     = (strstr(keylist,"TOTFX")   != NULL);
       int totfyflag     = (strstr(keylist,"TOTFY")   != NULL);
       int totfzflag     = (strstr(keylist,"TOTFZ")   != NULL);
       int epsxflag      = (strstr(keylist,"EPSX")    != NULL);
       int epsyflag      = (strstr(keylist,"EPSY")    != NULL);
       int epszflag      = (strstr(keylist,"EPSZ")    != NULL);
       int debugflag     = (strstr(debug,"debug")     != NULL);
   
           for (irec=0;irec<nrecs;irec++)
           {
           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 197  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;
         double rsun_ref, rsun_obs;         double rsun_ref, rsun_obs;
  
         // outrecords         // outrecords
         DRMS_Record_t *sharpoutrec, *sharpceaoutrec;              DRMS_Record_t *sharpoutrec, *sharpceaoutrec, *testrec;
         nrecs=sharpoutrecset->n;  
         printf("nrecs=%d\n",nrecs);  
  
         // prepare to set CODEVER7 (CVS Version of the SHARP module)         // prepare to set CODEVER7 (CVS Version of the SHARP module)
         char *cvsinfo0;         char *cvsinfo0;
Line 224  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);
  
         // no longer need inrecs          // grab one record such that we can malloc the R-parameter calculation arrays outside of the loop
         drms_close_records(sharpinrecset, DRMS_FREE_RECORD);          // assume scale = round(2./cdelt1) does not change for any given HARP
         drms_close_records(sharpceainrecset, DRMS_FREE_RECORD);          testrec  = sharpceainrecset->records[nrecs/2];
           dsun_obs = drms_getkey_float(testrec, "DSUN_OBS",   &status);
           rsun_ref = drms_getkey_double(testrec, "RSUN_REF", &status);
           cdelt1_orig = drms_getkey_float(testrec, "CDELT1",   &status);
           int scale = round(2.0/cdelt1);
           int nx1 = nx/scale;
           int ny1 = ny/scale;
           int nxp = nx1+40;
           int nyp = ny1+40;
   
           // check to make sure that the dimensions passed into fsample() are adequate
               if (nx1 > floor((nx-1)/scale + 1) )
                   DIE("X-dimension of output array in fsample() is too large.");
               if (ny1 > floor((ny-1)/scale + 1) )
                   DIE("Y-dimension of output array in fsample() is too large.");
           // malloc some arrays for the R parameter calculation
           float *rim     = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1n0    = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1p     = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1n     = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1      = (float *)malloc(nx1*ny1*sizeof(float));
           float *pmap    = (float *)malloc(nxp*nyp*sizeof(float));
           float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));
           float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));
  
         for (irec=0;irec<nrecs;irec++)  
         {  
            // Get emphemeris            // Get emphemeris
           sharpinrec    = sharpinrecset->records[irec];
           sharpceainrec = sharpceainrecset->records[irec];
           cdelt1_orig   = drms_getkey_float(sharpceainrec, "CDELT1",   &status);
           dsun_obs      = drms_getkey_float(sharpinrec, "DSUN_OBS",   &status);
           rsun_ref      = drms_getkey_double(sharpinrec, "RSUN_REF", &status);
           rsun_obs      = drms_getkey_double(sharpinrec, "RSUN_OBS", &status);
           imcrpix1      = drms_getkey_float(sharpinrec, "IMCRPIX1", &status);
           imcrpix2      = drms_getkey_float(sharpinrec, "IMCRPIX2", &status);
           crpix1        = drms_getkey_float(sharpinrec, "CRPIX1", &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];
            cdelt1_orig = drms_getkey_float(sharpceaoutrec, "CDELT1",   &status);             if (!sameflag)
            dsun_obs    = drms_getkey_float(sharpoutrec, "DSUN_OBS",   &status);             {
            rsun_ref    = drms_getkey_double(sharpoutrec, "RSUN_REF", &status);               drms_copykeys(sharpoutrec, sharpinrec, 0, kDRMS_KeyClass_Explicit);
            rsun_obs    = drms_getkey_double(sharpoutrec, "RSUN_OBS", &status);               drms_copykeys(sharpceaoutrec, sharpceainrec, 0, kDRMS_KeyClass_Explicit);
            imcrpix1    = drms_getkey_float(sharpoutrec, "IMCRPIX1", &status);             }
            imcrpix2    = drms_getkey_float(sharpoutrec, "IMCRPIX2", &status);  
            crpix1      = drms_getkey_float(sharpoutrec, "CRPIX1", &status);             TIME trec, tnow; /* 1970.01.01_00:00:00_UTC */
            crpix2      = drms_getkey_float(sharpoutrec, "CRPIX2", &status);             tnow = (double)time(NULL);
            cdelt1      = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);             tnow += UNIX_epoch;
  
            // set CODEVER7 and HISTORY             if (debugflag)
              {
                 printf("i=%d\n",irec);
              }
   
               // 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);
            drms_setkey_string(sharpceaoutrec, "CODEVER7", cvsinfoall);            drms_setkey_string(sharpceaoutrec, "CODEVER7", cvsinfoall);
            drms_setkey_string(sharpceaoutrec,"HISTORY",historyofthemodule);            drms_setkey_string(sharpceaoutrec,"HISTORY",historyofthemodule);
               drms_setkey_time(sharpoutrec,"DATE",tnow);
               drms_setkey_time(sharpceaoutrec,"DATE",tnow);
   
               // Get all the segments if computing vector magnetic field keywords
  
            // Get bx, by, bz, mask              DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");
            // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities  
            DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceaoutrec, "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(sharpceaoutrec, "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
  
            DRMS_Segment_t *bxSeg = drms_segment_lookup(sharpceaoutrec, "Bp");              DRMS_Segment_t *losSeg = drms_segment_lookup(sharpceainrec, "magnetogram");
               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
               //printf("Line 387");
   
               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(sharpceaoutrec, "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");
               //printf("Line 408\n");
               if (bzSeg == NULL)
               DIE("somekind of pointer problem.");
  
            DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceaoutrec, "Br");  
            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(sharpceaoutrec, "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(sharpceaoutrec, "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(sharpceaoutrec, "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 ************************************sdfdsf*/          /***** USFLUX, Example Function 1 *************************************/
            if (meanvfflag)            if (meanvfflag)
            {            {
               // Compute unsigned flux               // Compute unsigned flux
Line 313  int DoIt(void)
Line 447  int DoIt(void)
            drms_setkey_float(sharpceaoutrec, "ERRVF",   mean_vf_err);            drms_setkey_float(sharpceaoutrec, "ERRVF",   mean_vf_err);
            }            }
  
            /***** MEANPOT and TOTPOT, Example Function 13  (Requires Keiji's code) **/          /***** RPARAM, Example Function 14 *************************************/
               if (r_valueflag)
               {
               if (computeR(bz_err, los , dims, &Rparam, cdelt1, rim, p1p0, p1n0,
                           p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
  
               {
                           Rparam = DRMS_MISSING_FLOAT;
                   }
   
               drms_setkey_float(sharpoutrec, "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;
               }
   
               // 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 346  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 369  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 377  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 393  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 415  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 470  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 496  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 533  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 558  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 586  int DoIt(void)
Line 723  int DoIt(void)
               drms_setkey_float(sharpceaoutrec, "ERRMSHA", meanshear_angle_err);               drms_setkey_float(sharpceaoutrec, "ERRMSHA", meanshear_angle_err);
            }            }
  
            /******************************* END FLAGS **********************************/          /***** TOTFX, TOTFY, TOTFX, TOTBSQ, EPSX, EPSY, EPSZ , Example Function 16 (Lorentz forces) **********/
           if (totfxflag || totfyflag || totfzflag || totbsqflag || epsxflag || epsyflag || epszflag)
           {
               // Compute Lorentz forces
                   if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &totfx, &totfy, &totfz, &totbsq,
                                  &epsx, &epsy, &epsz, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
               {
                   totfx = DRMS_MISSING_FLOAT;
                   totfy = DRMS_MISSING_FLOAT;
                   totfz = DRMS_MISSING_FLOAT;
                   totbsq = DRMS_MISSING_FLOAT;
                   epsx = DRMS_MISSING_FLOAT;
                   epsy = DRMS_MISSING_FLOAT;
                   epsz = DRMS_MISSING_FLOAT;
               }
  
               // Set sharp keys
               drms_setkey_float(sharpoutrec, "TOTFX", totfx);
               drms_setkey_float(sharpoutrec, "TOTFY", totfy);
               drms_setkey_float(sharpoutrec, "TOTFZ", totfz);
               drms_setkey_float(sharpoutrec, "TOTBSQ",totbsq);
               drms_setkey_float(sharpoutrec, "EPSX",  epsx);
               drms_setkey_float(sharpoutrec, "EPSY",  epsy);
               drms_setkey_float(sharpoutrec, "EPSZ",  epsz);
   
               // Set sharp cea keys
               drms_setkey_float(sharpceaoutrec, "TOTFX", totfx);
               drms_setkey_float(sharpceaoutrec, "TOTFY", totfy);
               drms_setkey_float(sharpceaoutrec, "TOTFZ", totfz);
               drms_setkey_float(sharpceaoutrec, "TOTBSQ",totbsq);
               drms_setkey_float(sharpceaoutrec, "EPSX",  epsx);
               drms_setkey_float(sharpceaoutrec, "EPSY",  epsy);
               drms_setkey_float(sharpceaoutrec, "EPSZ",  epsz);
   
           }
   
           /***** 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 596  int DoIt(void)
Line 812  int DoIt(void)
            drms_free_array(bx_errArray);            drms_free_array(bx_errArray);
            drms_free_array(by_errArray);            drms_free_array(by_errArray);
            drms_free_array(bz_errArray);            drms_free_array(bz_errArray);
         } //endfor          drms_free_array(losArray);
           free(cvsinfoall);
         drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD);      free(rim);
         drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD);      free(p1p0);
       free(p1n0);
         //used for select calculations      free(p1p);
       free(p1n);
       free(p1);
       free(pmap);
       free(p1pad);
       free(pmapn);
       free(fx); free(fy); free(fz);
         free(bh); free(bt); free(jz);         free(bh); free(bt); free(jz);
         free(bpx); free(bpy); free(bpz);         free(bpx); free(bpy); free(bpz);
         free(derx); free(dery);         free(derx); free(dery);
         free(derx_bt); free(dery_bt);         free(derx_bt); free(dery_bt);
       free(derx_los); free(dery_los);
         free(derx_bz); free(dery_bz);         free(derx_bz); free(dery_bz);
         free(derx_bh); free(dery_bh);         free(derx_bh); free(dery_bh);
         free(bt_err); free(bh_err);  free(jz_err);         free(bt_err); free(bh_err);  free(jz_err);
         free(jz_err_squared); free(jz_rms_err);         free(jz_err_squared); free(jz_rms_err);
         free(cvsinfoall);      free(jz_err_squared_smooth);
       free(jz_smooth);
       free(err_term2);
       free(err_term1);
  
 return 0;          } //endfor
 }       // DoIt  
  
       // 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;
   
   }       // DoIt


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

Karen Tian
Powered by
ViewCVS 0.9.4