JSOC_Documentation: proj/sharp/apps/update_sharp

00001 /*                                                                      
00002  *      MODULE NAME: update_sharp_keys.c
00003  *
00004  *  DESCRIPTION: This module recalculates SHARP keywords. This is accomplished by 
00005  *      cloning a record, recalculating keywords of choice (i.e., user input), 
00006  *      and pointing to the same segments as the old record. Associated error keys
00007  *      are computed by default. 
00008  * 
00009  *      This module accounts for versioning by appending to the keyword CODEVER7 and HISTORY:
00010  *      CODEVER7 will contain multiple lines: the production build of sharp.c, the 
00011  *      production build of include file sw_functions.c, and the production build
00012  *      of update_sharp_keys.c. HISTORY will include a human-readable sentence about which
00013  *      keywords were updated. 
00014  * 
00015  *      N.B.: This module calculates the keyword specified in keylist and then does a 
00016  *      drms_copykeys() to copy over all the remaining keywords. However, if [1] a keyword
00017  *      does not exist in the .jsd files of the input series, and [2] said keyword X is not
00018  *      specified in keylist, then [3] drms_copykeys() will copy DRMS_MISSING_VALUE into 
00019  *      keyword X for all the output series.
00020  *
00021  *      This module does not produce segments.
00022  *
00023  *      INPUTS     : -- DRMS SHARP series
00024  *                   -- DRMS SHARP CEA series
00025  *                   -- HARPNUM
00026  *                   -- comma separated list of keywords to recalculate
00027  *                   -- DEBUG flag (use like this debug=debug)
00028  *
00029  *  AUTHOR     : Monica Bobra
00030  *
00031  *  Version    :   v0.0 Jun 14 2013
00032  *                     v0.1     Feb 11 2014 Added different input and output series
00033  *
00034  *  EXAMPLE    :
00035  *      update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //
00036  *      HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT debug=debug
00037  *
00038  */
00039 
00040 /* Include files */
00041 #include <stdio.h>
00042 #include <stdlib.h>
00043 #include <time.h>
00044 #include <sys/time.h>
00045 #include <math.h>
00046 #include <string.h>
00047 #include "jsoc_main.h"
00048 #include "astro.h"
00049 #include "fstats.h"
00050 #include "cartography.c"
00051 #include "fresize.h"
00052 #include "finterpolate.h"
00053 #include "img2helioVector.c"
00054 #include "copy_me_keys.c"
00055 #include "errorprop.c"
00056 #include "sw_functions.c"
00057 #include <mkl_blas.h>
00058 #include <mkl_service.h>
00059 #include <mkl_lapack.h>
00060 #include <mkl_vml_functions.h>
00061 #include <omp.h>
00062  
00063 /* Define values */                     
00064 #define PI              (M_PI)
00065 #define RADSINDEG       (PI/180.)
00066 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00067 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
00068 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00069 
00070 char *module_name = "update_sharp_keys";  /* Module name */
00071 char *version_id  = "2014 Feb 12";        /* Version number */
00072 
00073 ModuleArgs_t module_args[] =
00074 {
00075     {ARG_STRING, "sharpseriesin",      NULL, "Input Sharp dataseries"},
00076     {ARG_STRING, "sharpceaseriesin",   NULL, "Input Sharp CEA dataseries"},
00077     {ARG_STRING, "sharpseriesout",     NULL, "Output Sharp dataseries"},
00078     {ARG_STRING, "sharpceaseriesout",  NULL, "Output Sharp CEA dataseries"},
00079     {ARG_INT,    "HARPNUM",            NULL, "HARP number"},
00080     {ARG_STRING, "keylist",            NULL, "comma separated list of keywords to update"},
00081     {ARG_STRING, "debug",              NULL, "debug mode functionality. use like this: debug=debug"},
00082     {ARG_END}
00083 };
00084 
00085 
00086 int DoIt(void)                    
00087 {
00088         int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
00089         int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
00090 
00091     int status = DRMS_SUCCESS;
00092     int nrecs, irec;
00093 
00094         /* Keywords */
00095     float mean_vf; 
00096     float absFlux;
00097         float count_mask;
00098     float mean_hf;
00099     float mean_gamma;
00100     float mean_derivative_btotal;
00101     float mean_derivative_bh; 
00102     float mean_derivative_bz;
00103     float mean_jz;
00104     float us_i;
00105     float mean_alpha;
00106     float mean_ih;
00107     float total_us_ih;
00108     float total_abs_ih;
00109     float totaljz;
00110     float totpot;
00111     float meanpot;
00112     float area_w_shear_gt_45;
00113     float meanshear_angle;
00114     float area_w_shear_gt_45h;
00115     float meanshear_angleh; 
00116         float mean_derivative_btotal_err;
00117         float mean_vf_err;
00118         float mean_gamma_err;
00119         float mean_derivative_bh_err;
00120         float mean_derivative_bz_err;
00121         float mean_jz_err;
00122         float us_i_err; 
00123         float mean_alpha_err;
00124         float mean_ih_err;
00125         float total_us_ih_err;
00126         float total_abs_ih_err;
00127         float totaljz_err;
00128         float meanpot_err;
00129         float totpot_err;
00130         float Rparam;
00131         float meanshear_angle_err;
00132         float totfx;
00133         float totfy;
00134         float totfz;
00135         float totbsq;
00136         float epsx;
00137         float epsy;
00138         float epsz;
00139 
00140     char *sharpseriesin = (char *) params_get_str(&cmdparams, "sharpseriesin");
00141     char *sharpceaseriesin = (char *) params_get_str(&cmdparams, "sharpceaseriesin");
00142     char *sharpseriesout = (char *) params_get_str(&cmdparams, "sharpseriesout");
00143     char *sharpceaseriesout = (char *) params_get_str(&cmdparams, "sharpceaseriesout");
00144 
00145         int sameflag = strcmp(sharpseriesin, sharpseriesout) == 0;
00146         int testflag = strcmp(sharpceaseriesin, sharpceaseriesout) == 0;
00147         if (testflag != sameflag)
00148           DIE("Either both outputs must be the same as their inputs, or both be different.");
00149 
00150     int harpnum = params_get_int(&cmdparams, "HARPNUM");
00151 
00152     char sharpquery[256];
00153     char sharpceaquery[256];
00154     sprintf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
00155     printf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
00156 
00157     sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
00158     sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
00159 
00160     DRMS_RecordSet_t *sharpinrecset= drms_open_records(drms_env, sharpquery, &status);
00161     if (status != DRMS_SUCCESS)
00162         DIE("Problem opening sharp recordset.");
00163     nrecs=sharpinrecset->n;
00164     if (nrecs == 0)
00165         DIE("Sharp recordset contains no records.");
00166         printf("nrecs=%d\n",nrecs);
00167 
00168     DRMS_RecordSet_t *sharpceainrecset = drms_open_records(drms_env, sharpceaquery, &status);
00169     if (status != DRMS_SUCCESS)
00170         DIE("Problem opening sharp cea recordset.");
00171     if (sharpceainrecset->n != nrecs)
00172         DIE("Sharp cea recordset differs from sharp recordset.");
00173 
00174         DRMS_RecordSet_t *sharpoutrecset, *sharpceaoutrecset;
00175         if (sameflag)
00176         {
00177     sharpoutrecset = drms_clone_records(sharpinrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
00178     if (status != DRMS_SUCCESS)
00179         DIE("Problem cloning sharp records.");
00180     sharpceaoutrecset = drms_clone_records(sharpceainrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
00181     if (status != DRMS_SUCCESS)
00182         DIE("Problem cloning sharp cea records.");
00183         }
00184         else
00185         {
00186     sharpoutrecset = drms_create_records(drms_env, nrecs, sharpseriesout, DRMS_PERMANENT, &status);
00187     if (status != DRMS_SUCCESS)
00188         DIE("Problem creating sharp records.");
00189     sharpceaoutrecset = drms_create_records(drms_env, nrecs, sharpceaseriesout, DRMS_PERMANENT, &status);
00190     if (status != DRMS_SUCCESS)
00191         DIE("Problem creating sharp cea records.");
00192         }
00193 
00194 
00195     char *keylist = (char *) params_get_str(&cmdparams, "keylist");
00196     char *debug   = (char *) params_get_str(&cmdparams, "debug");
00197 
00198         // Flags to indicate which keyword will be recalculated
00199     int meanvfflag  = (strstr(keylist,"USFLUX")  != NULL);  // generalize so that lowercase is acceptable
00200     int totpotflag  = (strstr(keylist,"TOTPOT")  != NULL);  
00201         int meangamflag = (strstr(keylist,"MEANGAM") != NULL);
00202         int meangbtflag = (strstr(keylist,"MEANGBT") != NULL);
00203         int meangbzflag = (strstr(keylist,"MEANGBZ") != NULL);
00204         int meangbhflag = (strstr(keylist,"MEANGBH") != NULL);
00205         int meanjzdflag = (strstr(keylist,"MEANJZD") != NULL);
00206         int totusjzflag = (strstr(keylist,"TOTUSJZ") != NULL);
00207         int meanalpflag = (strstr(keylist,"MEANALP") != NULL);
00208         int meanjzhflag = (strstr(keylist,"MEANJZH") != NULL);
00209         int totusjhflag = (strstr(keylist,"TOTUSJH") != NULL);
00210         int absnjzhflag = (strstr(keylist,"ABSNJZH") != NULL);
00211         int savncppflag = (strstr(keylist,"SAVNCPP") != NULL);
00212         int meanpotflag = (strstr(keylist,"MEANPOT") != NULL);
00213         int meanshrflag = (strstr(keylist,"MEANSHR") != NULL);
00214         int shrgt45flag = (strstr(keylist,"SHRGT45") != NULL);
00215         int r_valueflag = (strstr(keylist,"R_VALUE") != NULL);
00216         int totbsqflag  = (strstr(keylist,"TOTBSQ")  != NULL);
00217         int totfxflag   = (strstr(keylist,"TOTFX")   != NULL);
00218         int totfyflag   = (strstr(keylist,"TOTFY")   != NULL);
00219         int totfzflag   = (strstr(keylist,"TOTFZ")   != NULL);
00220         int epsxflag    = (strstr(keylist,"EPSX")    != NULL);
00221         int epsyflag    = (strstr(keylist,"EPSY")    != NULL);
00222         int epszflag    = (strstr(keylist,"EPSZ")    != NULL);
00223         int debugflag   = (strstr(debug,"debug")     != NULL);
00224 
00225     DRMS_Record_t *sharpinrec = sharpinrecset->records[0];
00226     DRMS_Record_t *sharpceainrec = sharpceainrecset->records[0];
00227     DRMS_Segment_t *inseg = drms_segment_lookup(sharpceainrec, "Br");
00228     int nx = inseg->axis[0];
00229     int ny = inseg->axis[1];
00230     int nxny = nx * ny;
00231     int dims[2] = {nx, ny};
00232 
00233     // Temp arrays  
00234     float *bh      = (float *) (malloc(nxny * sizeof(float)));
00235     float *bt      = (float *) (malloc(nxny * sizeof(float)));
00236     float *jz      = (float *) (malloc(nxny * sizeof(float)));
00237     float *bpx     = (float *) (malloc(nxny * sizeof(float)));
00238     float *bpy     = (float *) (malloc(nxny * sizeof(float)));
00239     float *bpz     = (float *) (malloc(nxny * sizeof(float)));
00240     float *derx    = (float *) (malloc(nxny * sizeof(float)));
00241     float *dery    = (float *) (malloc(nxny * sizeof(float)));
00242     float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
00243     float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
00244     float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
00245     float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
00246     float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
00247     float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
00248     float *bt_err  = (float *) (malloc(nxny * sizeof(float)));
00249     float *bh_err  = (float *) (malloc(nxny * sizeof(float)));
00250         float *jz_err  = (float *) (malloc(nxny * sizeof(float)));
00251         float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
00252         float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));
00253         float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
00254     float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
00255     float *err_term1   = (float *) (calloc(nxny, sizeof(float)));
00256     float *err_term2   = (float *) (calloc(nxny, sizeof(float)));
00257         float *fx = (float *) (malloc(nxny * sizeof(float)));
00258         float *fy = (float *) (malloc(nxny * sizeof(float)));
00259         float *fz = (float *) (malloc(nxny * sizeof(float)));
00260 
00261     for (irec=0;irec<nrecs;irec++)
00262     {
00263 
00264         //DRMS_Record_t *sharpinrec = sharpinrecset->records[irec];
00265     DRMS_Record_t *sharpceainrec = sharpceainrecset->records[irec];
00266     DRMS_Segment_t *inseg = drms_segment_lookup(sharpceainrec, "Br");
00267  
00268         // ephemeris variables
00269     float  cdelt1_orig, cdelt1, dsun_obs, imcrpix1, imcrpix2, crpix1, crpix2;
00270     double rsun_ref, rsun_obs;
00271 
00272         // outrecords 
00273     DRMS_Record_t *sharpoutrec, *sharpceaoutrec, *testrec;
00274 
00275         // prepare to set CODEVER7 (CVS Version of the SHARP module)
00276     char *cvsinfo0;
00277     char *history0;
00278     char *cvsinfo1 = strdup("$Id: update_sharp_keys.c,v 1.15 2015/03/11 22:18:05 mbobra Exp $");
00279     char *cvsinfo2 = sw_functions_version();
00280     char *cvsinfoall = (char *)malloc(2048);
00281         char historyofthemodule[2048];
00282         cvsinfo0    = drms_getkey_string(sharpinrec, "CODEVER7", &status);
00283         strcpy(cvsinfoall,cvsinfo0);
00284         strcat(cvsinfoall,"\n");
00285         strcat(cvsinfoall,cvsinfo2);
00286 
00287         // prepare to set HISTORY (History of the data)
00288     char timebuf[1024];
00289     float UNIX_epoch = -220924792.000;  
00290     sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
00291         sprintf(historyofthemodule,"The following keywords were re-computed on %s: %s.",timebuf,keylist); 
00292     //printf("historyofthemodule=%s\n",historyofthemodule);
00293     history0    = drms_getkey_string(sharpinrec, "HISTORY", &status);
00294     strcat(historyofthemodule,"\n");
00295         strcat(historyofthemodule,history0);
00296  
00297         // grab one record such that we can malloc the R-parameter calculation arrays outside of the loop
00298         // assume scale = round(2./cdelt1) does not change for any given HARP
00299         testrec  = sharpceainrecset->records[nrecs/2];
00300         dsun_obs = drms_getkey_float(testrec, "DSUN_OBS",   &status); 
00301         rsun_ref = drms_getkey_double(testrec, "RSUN_REF", &status);
00302         cdelt1_orig = drms_getkey_float(testrec, "CDELT1",   &status);
00303         cdelt1      = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
00304         int scale = round(2.0/cdelt1);
00305         int nx1 = nx/scale;
00306         int ny1 = ny/scale;
00307         int nxp = nx1+40;
00308         int nyp = ny1+40;
00309 
00310         // check to make sure that the dimensions passed into fsample() are adequate
00311     if (nx1 > floor((nx-1)/scale + 1) )
00312         DIE("X-dimension of output array in fsample() is too large.");
00313     if (ny1 > floor((ny-1)/scale + 1) )
00314             DIE("Y-dimension of output array in fsample() is too large.");
00315  
00316         // malloc some arrays for the R parameter calculation
00317         float *rim     = (float *)malloc(nx1*ny1*sizeof(float));
00318         float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));
00319         float *p1n0    = (float *)malloc(nx1*ny1*sizeof(float));
00320         float *p1p     = (float *)malloc(nx1*ny1*sizeof(float));
00321         float *p1n     = (float *)malloc(nx1*ny1*sizeof(float));
00322         float *p1      = (float *)malloc(nx1*ny1*sizeof(float));
00323         float *pmap    = (float *)malloc(nxp*nyp*sizeof(float));
00324         float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));
00325         float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));
00326     
00327        // Get emphemeris
00328        sharpinrec    = sharpinrecset->records[irec];
00329        sharpceainrec = sharpceainrecset->records[irec];
00330        cdelt1_orig   = drms_getkey_float(sharpceainrec, "CDELT1",   &status);
00331        dsun_obs      = drms_getkey_float(sharpinrec, "DSUN_OBS",   &status);
00332        rsun_ref      = drms_getkey_double(sharpinrec, "RSUN_REF", &status);
00333        rsun_obs      = drms_getkey_double(sharpinrec, "RSUN_OBS", &status);
00334        imcrpix1      = drms_getkey_float(sharpinrec, "IMCRPIX1", &status);
00335        imcrpix2      = drms_getkey_float(sharpinrec, "IMCRPIX2", &status);
00336        crpix1        = drms_getkey_float(sharpinrec, "CRPIX1", &status);
00337        crpix2        = drms_getkey_float(sharpinrec, "CRPIX2", &status);
00338 
00339        sharpoutrec = sharpoutrecset->records[irec];
00340        sharpceaoutrec = sharpceaoutrecset->records[irec];
00341            if (!sameflag)
00342            {
00343              drms_copykeys(sharpoutrec, sharpinrec, 0, kDRMS_KeyClass_Explicit);
00344              drms_copykeys(sharpceaoutrec, sharpceainrec, 0, kDRMS_KeyClass_Explicit);
00345            }
00346 
00347            TIME trec, tnow; /* 1970.01.01_00:00:00_UTC */
00348            tnow = (double)time(NULL);
00349            tnow += UNIX_epoch;
00350 
00351            if (debugflag)
00352            {
00353               printf("i=%d\n",irec);
00354            }
00355 
00356        // set CODEVER7 and HISTORY and DATE
00357        drms_setkey_string(sharpoutrec, "CODEVER7", cvsinfoall);
00358            drms_setkey_string(sharpoutrec,"HISTORY",historyofthemodule);           
00359        drms_setkey_string(sharpceaoutrec, "CODEVER7", cvsinfoall);
00360            drms_setkey_string(sharpceaoutrec,"HISTORY",historyofthemodule);
00361            drms_setkey_time(sharpoutrec,"DATE",tnow);  
00362            drms_setkey_time(sharpceaoutrec,"DATE",tnow);
00363 
00364        // Get bx, by, bz, mask  
00365        // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
00366        DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");
00367        DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
00368        int *bitmask = (int *) bitmaskArray->data;    // get the previously made mask array  
00369 
00370        //Use conf_disambig map as a threshold on spaceweather quantities 
00371        DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceainrec, "conf_disambig");       
00372        DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);    
00373            if (status != DRMS_SUCCESS)
00374               DIE("No CONF_DISAMBIG segment.");
00375        int *mask = (int *) maskArray->data;     // get the previously made mask array   
00376 
00377        //Use magnetogram map to compute R
00378        DRMS_Segment_t *losSeg = drms_segment_lookup(sharpceainrec, "magnetogram");
00379        DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
00380        float *los = (float *) losArray->data;       // los
00381 
00382        DRMS_Segment_t *bxSeg = drms_segment_lookup(sharpceainrec, "Bp");
00383        DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
00384        float *bx = (float *) bxArray->data;     // bx
00385     
00386        DRMS_Segment_t *bySeg = drms_segment_lookup(sharpceainrec, "Bt");
00387        DRMS_Array_t *byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status);
00388        float *by = (float *) byArray->data;     // by
00389        for (int i = 0; i < nxny; i++) by[i] *= -1;
00390     
00391        DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceainrec, "Br");
00392        DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
00393        float *bz = (float *) bzArray->data;     // bz
00394 
00395        DRMS_Segment_t *bz_errSeg = drms_segment_lookup(sharpceainrec, "Br_err");
00396        DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
00397        float *bz_err = (float *) bz_errArray->data;     // bz_err
00398 
00399        DRMS_Segment_t *by_errSeg = drms_segment_lookup(sharpceainrec, "Bt_err");
00400        DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
00401        float *by_err = (float *) by_errArray->data;     // by_err
00402 
00403        DRMS_Segment_t *bx_errSeg = drms_segment_lookup(sharpceainrec, "Bp_err");
00404        DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
00405        float *bx_err = (float *) bx_errArray->data;     // bx_err
00406 
00407            /***** USFLUX, Example Function 1 *************************************/ 
00408        if (meanvfflag)
00409        {
00410               // Compute unsigned flux 
00411               if (computeAbsFlux(bz_err, bz , dims, &absFlux, &mean_vf,  &mean_vf_err, 
00412                                  &count_mask, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
00413               {
00414              mean_vf           = DRMS_MISSING_FLOAT;
00415                  mean_vf_err       = DRMS_MISSING_FLOAT;
00416                  count_mask        = DRMS_MISSING_INT;
00417           }
00418 
00419            drms_setkey_float(sharpoutrec, "USFLUX",  mean_vf);
00420            drms_setkey_float(sharpoutrec, "CMASK",   count_mask);   
00421            drms_setkey_float(sharpoutrec, "ERRVF",   mean_vf_err);
00422 
00423            drms_setkey_float(sharpceaoutrec, "USFLUX",  mean_vf);
00424            drms_setkey_float(sharpceaoutrec, "CMASK",   count_mask);    
00425            drms_setkey_float(sharpceaoutrec, "ERRVF",   mean_vf_err);
00426            }
00427 
00428             /***** RPARAM, Example Function 14 *************************************/ 
00429        if (r_valueflag)
00430        {
00431               if (computeR(bz_err, los , dims, &Rparam, cdelt1, rim, p1p0, p1n0,
00432                            p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
00433 
00434               {
00435          Rparam            = DRMS_MISSING_FLOAT;
00436           }
00437 
00438            drms_setkey_float(sharpoutrec, "R_VALUE",  Rparam);
00439            drms_setkey_float(sharpceaoutrec, "R_VALUE", Rparam);
00440 
00441            }
00442   
00443            /***** MEANPOT and TOTPOT, Example Function 13  (Requires Keiji's code) **/
00444        
00445        if (totpotflag || meanpotflag)
00446        {
00447               // First compute potential field  
00448               for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
00449           greenpot(bpx, bpy, bpz, nx, ny);
00450  
00451               // Then compute energy   
00452           if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims, 
00453                                     &meanpot, &meanpot_err, &totpot, &totpot_err, 
00454                         mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 
00455               {
00456              meanpot           = DRMS_MISSING_FLOAT; 
00457              totpot            = DRMS_MISSING_FLOAT;
00458                  meanpot_err       = DRMS_MISSING_FLOAT;
00459                  totpot_err        = DRMS_MISSING_FLOAT;
00460           }
00461        
00462               // Set sharp keys          
00463               drms_setkey_float(sharpoutrec, "MEANPOT",  meanpot); 
00464               drms_setkey_float(sharpoutrec, "TOTPOT",   totpot);   
00465               drms_setkey_float(sharpoutrec, "ERRMPOT",  meanpot_err); 
00466               drms_setkey_float(sharpoutrec, "ERRTPOT",  totpot_err);  
00467 
00468               // Set sharp cea keys
00469               drms_setkey_float(sharpceaoutrec, "MEANPOT",  meanpot); 
00470               drms_setkey_float(sharpceaoutrec, "TOTPOT",   totpot);    
00471               drms_setkey_float(sharpceaoutrec, "ERRMPOT",  meanpot_err); 
00472               drms_setkey_float(sharpceaoutrec, "ERRTPOT",  totpot_err);
00473            }
00474        
00475            /***** MEANGAM, Example Function 3 ************************************/
00476    
00477            if (meangamflag)
00478            {
00479               // First compute horizontal field
00480               computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask);
00481 
00482               // Then compute inclination angle, gamma
00483               if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &mean_gamma, &mean_gamma_err,mask, bitmask))
00484           { 
00485                  mean_gamma                 =  DRMS_MISSING_FLOAT;
00486                  mean_gamma_err             =  DRMS_MISSING_FLOAT;
00487               }
00488 
00489               // Set sharp keys    
00490               drms_setkey_float(sharpoutrec, "MEANGAM",  mean_gamma); 
00491               drms_setkey_float(sharpoutrec, "ERRGAM",   mean_gamma_err);   
00492 
00493               // Set sharp cea keys
00494               drms_setkey_float(sharpceaoutrec, "MEANGAM", mean_gamma); 
00495               drms_setkey_float(sharpceaoutrec, "ERRGAM",  mean_gamma_err);     
00496            }
00497 
00498            /***** MEANGBT, Example Function 5 (Requires Function 4) *************/
00499 
00500            if (meangbtflag)
00501            {
00502               // First, compute Bt 
00503               computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);
00504 
00505               // Then take the derivative of Bt
00506           if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask, bitmask, derx_bt, dery_bt, bt_err, 
00507                                           &mean_derivative_btotal_err, err_term1, err_term2))
00508               {
00509              mean_derivative_btotal     = DRMS_MISSING_FLOAT;
00510              mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
00511               }
00512 
00513               // Set sharp keys    
00514               drms_setkey_float(sharpoutrec, "MEANGBT",  mean_derivative_btotal); 
00515               drms_setkey_float(sharpoutrec, "ERRBT",  mean_derivative_btotal_err); 
00516     
00517               // Set sharp cea keys  
00518               drms_setkey_float(sharpceaoutrec, "MEANGBT", mean_derivative_btotal); 
00519               drms_setkey_float(sharpceaoutrec, "ERRBT", mean_derivative_btotal_err);   
00520            }
00521    
00522            /***** MEANGBH, Example Function 6 (Requires Function 2) *************/
00523 
00524            if (meangbhflag)
00525            {
00526               // First, compute Bh
00527           computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask);
00528 
00529               // Then take the derivative of Bh
00530           if (computeBhderivative(bh, bh_err, dims, &mean_derivative_bh, &mean_derivative_bh_err, mask, bitmask, derx_bh, dery_bh, err_term1, err_term2))
00531               {
00532              mean_derivative_bh       = DRMS_MISSING_FLOAT;
00533                  mean_derivative_bh_err   = DRMS_MISSING_FLOAT;
00534           }
00535 
00536               // Set sharp keys   
00537               drms_setkey_float(sharpoutrec, "MEANGBH", mean_derivative_bh); 
00538               drms_setkey_float(sharpoutrec, "ERRBH", mean_derivative_bh_err); 
00539     
00540               // Set sharp cea keys  
00541               drms_setkey_float(sharpceaoutrec, "MEANGBH", mean_derivative_bh); 
00542               drms_setkey_float(sharpceaoutrec, "ERRBH", mean_derivative_bh_err);   
00543            }
00544 
00545            /***** MEANGBZ, Example Function 7 ************************************/
00546 
00547            if (meangbzflag)
00548            {
00549               // Compute Bz derivative
00550               if (computeBzderivative(bz, bz_err, dims, &mean_derivative_bz, &mean_derivative_bz_err, mask, bitmask, derx_bz, dery_bz, err_term1, err_term2))
00551               {
00552              mean_derivative_bz     = DRMS_MISSING_FLOAT; 
00553                  mean_derivative_bz_err = DRMS_MISSING_FLOAT; 
00554               }
00555     
00556               // Set sharp keys   
00557               drms_setkey_float(sharpoutrec, "MEANGBZ",  mean_derivative_bz); 
00558               drms_setkey_float(sharpoutrec, "ERRBZ",  mean_derivative_bz_err); 
00559     
00560               // Set sharp cea keys   
00561               drms_setkey_float(sharpceaoutrec, "MEANGBZ", mean_derivative_bz); 
00562               drms_setkey_float(sharpceaoutrec, "ERRBZ", mean_derivative_bz_err);   
00563            }
00564 
00565            /***** MEANJZD and TOTUSJZ, Example Function 9 (Requires Function 8) ***/
00566 
00567            if (meanjzdflag || totusjzflag)
00568            {
00569               // First, compute Jz on all pixels
00570           computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
00571                      derx, dery, err_term1, err_term2);
00572    
00573               // Then take the sums of the appropriate pixels
00574               if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &mean_jz,
00575                        &mean_jz_err, &us_i, &us_i_err, mask, bitmask, cdelt1,
00576                        rsun_ref, rsun_obs, derx, dery))
00577 
00578 
00579               {
00580                  mean_jz                = DRMS_MISSING_FLOAT;
00581              us_i                   = DRMS_MISSING_FLOAT;
00582                  mean_jz_err            = DRMS_MISSING_FLOAT;
00583                  us_i_err               = DRMS_MISSING_FLOAT;
00584           }
00585 
00586               // Set sharp keys   
00587               drms_setkey_float(sharpoutrec, "MEANJZD",  mean_jz); 
00588               drms_setkey_float(sharpoutrec, "TOTUSJZ",  us_i);     
00589               drms_setkey_float(sharpoutrec, "ERRJZ",  mean_jz_err); 
00590               drms_setkey_float(sharpoutrec, "ERRUSI",  us_i_err);  
00591 
00592               // Set sharp cea keys   
00593               drms_setkey_float(sharpceaoutrec, "MEANJZD", mean_jz); 
00594               drms_setkey_float(sharpceaoutrec, "TOTUSJZ", us_i);   
00595               drms_setkey_float(sharpceaoutrec, "ERRJZ", mean_jz_err); 
00596               drms_setkey_float(sharpceaoutrec, "ERRUSI", us_i_err);    
00597            }
00598 
00599            /***** MEANALP, Example Function 10 (Requires Function 8)*********/
00600 
00601            if (meanalpflag)
00602            {
00603               // First, compute Jz on all pixels
00604           computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
00605                         derx, dery, err_term1, err_term2);
00606 
00607               // Then compute alpha quantities on select pixels
00608           if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &mean_alpha, &mean_alpha_err, 
00609                                mask, bitmask, cdelt1, rsun_ref, rsun_obs))
00610               {
00611              mean_alpha             = DRMS_MISSING_FLOAT;
00612                  mean_alpha_err         = DRMS_MISSING_FLOAT;
00613               }
00614 
00615               // Set sharp keys     
00616               drms_setkey_float(sharpoutrec, "MEANALP", mean_alpha); 
00617               drms_setkey_float(sharpoutrec, "ERRALP", mean_alpha_err);     
00618 
00619               // Set sharp cea keys   
00620               drms_setkey_float(sharpceaoutrec, "MEANALP", mean_alpha); 
00621               drms_setkey_float(sharpceaoutrec, "ERRALP", mean_alpha_err);  
00622 
00623            }
00624 
00625            /***** MEANJZH, TOTUSJH, ABSNJZH, Example Function 11 (Requires Function 8) ***/
00626 
00627            if (meanjzhflag || totusjhflag || absnjzhflag)
00628            {
00629               // First, compute Jz on all pixels
00630               computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
00631                         derx, dery, err_term1, err_term2);
00632 
00633               // Then compute helicity quantities on select pixels
00634            if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &mean_ih, &mean_ih_err, &total_us_ih, 
00635                                &total_abs_ih, &total_us_ih_err, &total_abs_ih_err, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
00636               {  
00637              mean_ih                = DRMS_MISSING_FLOAT; 
00638              total_us_ih            = DRMS_MISSING_FLOAT;
00639              total_abs_ih           = DRMS_MISSING_FLOAT;
00640                  mean_ih_err            = DRMS_MISSING_FLOAT;
00641                  total_us_ih_err        = DRMS_MISSING_FLOAT;
00642                  total_abs_ih_err       = DRMS_MISSING_FLOAT;
00643           }
00644 
00645               // Set sharp keys 
00646               drms_setkey_float(sharpoutrec, "MEANJZH",  mean_ih); 
00647               drms_setkey_float(sharpoutrec, "TOTUSJH",  total_us_ih);  
00648               drms_setkey_float(sharpoutrec, "ABSNJZH",  total_abs_ih); 
00649               drms_setkey_float(sharpoutrec, "ERRMIH",  mean_ih_err);   
00650               drms_setkey_float(sharpoutrec, "ERRTUI",  total_us_ih_err); 
00651               drms_setkey_float(sharpoutrec, "ERRTAI",  total_abs_ih_err);  
00652 
00653               // Set sharp cea keys 
00654               drms_setkey_float(sharpceaoutrec, "MEANJZH", mean_ih); 
00655               drms_setkey_float(sharpceaoutrec, "TOTUSJH", total_us_ih);    
00656               drms_setkey_float(sharpceaoutrec, "ABSNJZH", total_abs_ih); 
00657               drms_setkey_float(sharpceaoutrec, "ERRMIH", mean_ih_err);     
00658               drms_setkey_float(sharpceaoutrec, "ERRTUI", total_us_ih_err); 
00659               drms_setkey_float(sharpceaoutrec, "ERRTAI", total_abs_ih_err);    
00660            }
00661 
00662            /***** SAVNCPP, Example Function 12 (Requires Function 8) *******************/
00663    
00664            if (savncppflag)
00665            {
00666               // First, compute Jz on all pixels
00667           computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
00668                         derx, dery, err_term1, err_term2);
00669 
00670               // Then compute sums of Jz on select pixels
00671           if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &totaljz, &totaljz_err, 
00672                            mask, bitmask, cdelt1, rsun_ref, rsun_obs))
00673               {
00674              totaljz                = DRMS_MISSING_FLOAT;
00675                  totaljz_err            = DRMS_MISSING_FLOAT;
00676           }
00677 
00678               // Set sharp keys     
00679               drms_setkey_float(sharpoutrec, "SAVNCPP",  totaljz); 
00680               drms_setkey_float(sharpoutrec, "ERRJHT",  totaljz_err);   
00681 
00682               // Set sharp cea keys 
00683               drms_setkey_float(sharpceaoutrec, "SAVNCPP", totaljz); 
00684               drms_setkey_float(sharpceaoutrec, "ERRJHT", totaljz_err);     
00685            }
00686 
00687            /***** MEANSHR and SHRGT45, Example Function 14 (Requires Keiji's code) **********/
00688 
00689            if (meanshrflag || shrgt45flag)
00690            {
00691               // First compute potential field
00692           for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
00693           greenpot(bpx, bpy, bpz, nx, ny); 
00694 
00695               // Then compute shear angles
00696           if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims, 
00697                     &meanshear_angle, &meanshear_angle_err, &area_w_shear_gt_45,  
00698                     mask, bitmask)) 
00699               {
00700               meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN
00701           area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
00702                   meanshear_angle_err= DRMS_MISSING_FLOAT;
00703           }
00704 
00705               // Set sharp keys     
00706               drms_setkey_float(sharpoutrec, "MEANSHR",  meanshear_angle); 
00707               drms_setkey_float(sharpoutrec, "SHRGT45",  area_w_shear_gt_45);   
00708               drms_setkey_float(sharpoutrec, "ERRMSHA",  meanshear_angle_err);  
00709 
00710               // Set sharp cea keys 
00711               drms_setkey_float(sharpceaoutrec, "MEANSHR", meanshear_angle); 
00712               drms_setkey_float(sharpceaoutrec, "SHRGT45", area_w_shear_gt_45);     
00713               drms_setkey_float(sharpceaoutrec, "ERRMSHA", meanshear_angle_err);      
00714            }
00715 
00716           /***** TOTFX, TOTFY, TOTFX, TOTBSQ, EPSX, EPSY, EPSZ , Example Function 16 (Lorentz forces) **********/
00717 
00718            if (totfxflag || totfyflag || totfzflag || totbsqflag || epsxflag || epsyflag || epszflag)
00719            {
00720 
00721                 // Compute Lorentz forces
00722             if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &totfx, &totfy, &totfz, &totbsq,
00723                     &epsx, &epsy, &epsz, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
00724                 {  
00725            totfx             = DRMS_MISSING_FLOAT;
00726                    totfy             = DRMS_MISSING_FLOAT;
00727            totfz             = DRMS_MISSING_FLOAT;
00728            totbsq            = DRMS_MISSING_FLOAT;
00729                    epsx              = DRMS_MISSING_FLOAT;
00730            epsy              = DRMS_MISSING_FLOAT;
00731                    epsz              = DRMS_MISSING_FLOAT;
00732             }
00733 
00734               // Set sharp keys     
00735               drms_setkey_float(sharpoutrec, "TOTFX", totfx); 
00736               drms_setkey_float(sharpoutrec, "TOTFY", totfy);   
00737               drms_setkey_float(sharpoutrec, "TOTFZ", totfz);   
00738               drms_setkey_float(sharpoutrec, "TOTBSQ",totbsq); 
00739               drms_setkey_float(sharpoutrec, "EPSX",  epsx);    
00740               drms_setkey_float(sharpoutrec, "EPSY",  epsy);    
00741               drms_setkey_float(sharpoutrec, "EPSZ",  epsz);    
00742 
00743               // Set sharp cea keys 
00744               drms_setkey_float(sharpceaoutrec, "TOTFX", totfx); 
00745               drms_setkey_float(sharpceaoutrec, "TOTFY", totfy);    
00746               drms_setkey_float(sharpceaoutrec, "TOTFZ", totfz);    
00747               drms_setkey_float(sharpceaoutrec, "TOTBSQ",totbsq); 
00748               drms_setkey_float(sharpceaoutrec, "EPSX",  epsx);     
00749               drms_setkey_float(sharpceaoutrec, "EPSY",  epsy);     
00750               drms_setkey_float(sharpceaoutrec, "EPSZ",  epsz);     
00751               
00752            }
00753 
00754            /******************************* END FLAGS **********************************/
00755 
00756     drms_free_array(bitmaskArray);      
00757     drms_free_array(maskArray);
00758     drms_free_array(bxArray);           
00759     drms_free_array(byArray);
00760     drms_free_array(bzArray);
00761     drms_free_array(bx_errArray);           
00762     drms_free_array(by_errArray);
00763     drms_free_array(bz_errArray);
00764     drms_free_array(losArray);
00765     free(cvsinfoall);
00766 
00767         free(rim);
00768         free(p1p0);
00769         free(p1n0);
00770         free(p1p);
00771         free(p1n);
00772         free(p1);
00773         free(pmap);
00774         free(p1pad); 
00775         free(pmapn);
00776 
00777     } //endfor
00778 
00779         free(fx); free(fy); free(fz);
00780     free(bh); free(bt); free(jz); 
00781     free(bpx); free(bpy); free(bpz); 
00782     free(derx); free(dery);  
00783     free(derx_bt); free(dery_bt);   
00784     free(derx_bz); free(dery_bz);    
00785     free(derx_bh); free(dery_bh); 
00786     free(bt_err); free(bh_err);  free(jz_err); 
00787         free(jz_err_squared); free(jz_rms_err);
00788         free(jz_err_squared_smooth);
00789         free(jz_smooth);
00790         free(err_term2); 
00791         free(err_term1); 
00792     drms_close_records(sharpinrecset, DRMS_FREE_RECORD);
00793     drms_close_records(sharpceainrecset, DRMS_FREE_RECORD);
00794 
00795     drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD);
00796     drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD);
00797 
00798     printf("Success=%d\n",sameflag);
00799 
00800 return 0;    
00801 
00802 }   // DoIt
proj/sharp/apps/update_sharp_keys.c
