JSOC_Documentation: proj/sharp/apps/sharp.c Source File

00001 /*
00002  *  sharp.c   
00003  *
00004  *  This module creates the pipeline Space Weather Active Region Patches (SHARPs).
00005  *  It is a hard-coded strip-down version of bmap.c.
00006  *  It takes the Mharp and Bharp series and create the following quantities:
00007  *
00008  *      Series 1: Sharp_CEA
00009  *            CEA remapped magnetogram, bitmap, continuum, doppler (same size in map coordinate, need manual spec?)
00010  *            CEA remapped vector field (Br, Bt, Bp) (same as above)
00011  *                Space weather indices based on vector cutouts (step 2)
00012  *
00013  *      Series 2: Sharp_cutout:
00014  *            cutouts of magnetogram, bitmap, continuum, doppler (HARP defined, various sizes in CCD pixels)
00015  *            cutouts of all vector data segments (same as above)
00016  *           
00017  *  Author:
00018  *      Xudong Sun; Monica Bobra
00019  *
00020  *  Version:
00021  *              v0.0 Jul 02 2012
00022  *              v0.1 Jul 23 2012
00023  *              v0.2 Sep 04 2012
00024  *              v0.3 Dec 18 2012
00025  *              v0.4 Jan 02 2013
00026  *              v0.5 Jan 23 2013
00027  *              v0.6 Aug 12 2013
00028  *              v0.7 Jan 02 2014
00029  *              v0.8 Feb 12 2014
00030  *              v0.9 Mar 04 2014
00031  *
00032  *  Notes:
00033  *      v0.0
00034  *      Mharp & Bharp must be fully specified; other input are series names only
00035  *      All input records need to match, otherwise quit
00036  *      Mapping parameters depend on keywords of each record only, not necessarily consistent for now
00037  *      Cutout doesn't work for char segments yet (drms_segment_readslice bug)
00038  *      SW indices require ephemeris info which is not passed properly as of now
00039  *      v0.1
00040  *      Fixed char I/O thanks to Art
00041  *      SW indices fixed
00042  *      Added doppler and continuum
00043  *              Added other keywords: HEADER (populated by cvs build version), DATE_B
00044  *      v0.3
00045  *      Fixed memory leakage of 0.15G per rec; denoted with "Dec 18"
00046  *      v0.4
00047  *      Took out convert_inplace(). Was causing all the images to be int
00048  *              v0.5
00049  *              Corrected ephemeris keywords, added argument mInfo for setKeys()
00050  *      v0.6
00051  *      Changes in remapping of bitmap and conf_disambig, now near neighbor without anti-aliasing
00052  *              v0.7
00053  *              Added full disk as "b"
00054  *              Global flag fullDisk is set if "b" is set
00055  *              Utilize BharpRS and BharpRec all around
00056  *              Pass mharpRec to set_keys() too in case of full disk
00057  *              Fixed Bunit (removed from copy_me_keys(), added loops for Bunits in set_keys() here)
00058  *              Error for CEA still does account for disambiguation yet
00059  *      v0.8
00060  *      Added disambig to azimuth during error propagation
00061  *      Changed usage for disambig: bit 2 (radial acute) for full disk, bit 0 for patch
00062  *      Fixed disambig cutout for patch: 0 for even, 7 for odd
00063  *      v0.9
00064  *      Fixed unit
00065  *      Check whether in PATCH of FD mode, so the error propagation uses disambiguated azimuth or not
00066  *
00067  *
00068  *  Example Calls:
00069  *      [I]   B (full disk disambiguation)
00070  *      > sharp "mharp=hmi.Mharp_720s[1832][2012.07.12_15:24]" "b=hmi_test.B_720s[2012.07.12_15:24]" "dop=hmi.V_720s[2012.07.12_15:24]" "cont=hmi.Ic_720s[2012.07.12_15:24]" "sharp_cea=su_xudong.sharp_cea_720s" "sharp_cut=su_xudong.sharp_720s"
00071  *      [II]  BHARP (patch disambiguation)
00072  *      > sharp "mharp=hmi.Mharp_720s[1832][2012.07.12_15:24]" "bharp=hmi.Bharp_720s[1832][2012.07.12_15:24]" "dop=hmi.V_720s[2012.07.12_15:24]" "cont=hmi.Ic_720s[2012.07.12_15:24]" "sharp_cea=su_xudong.sharp_cea_720s" "sharp_cut=su_xudong.sharp_720s"
00073  *
00074  *
00075  */
00076 
00077 #include <stdio.h>
00078 #include <stdlib.h>
00079 #include <time.h>
00080 #include <sys/time.h>
00081 #include <math.h>
00082 #include <string.h>
00083 #include "jsoc_main.h"
00084 #include "astro.h"
00085 #include "fstats.h"
00086 #include "cartography.c"
00087 #include "fresize.h"
00088 #include "finterpolate.h"
00089 #include "img2helioVector.c"
00090 #include "copy_me_keys.c"
00091 #include "errorprop.c"
00092 #include "sw_functions.c"
00093 
00094 //#include <mkl.h> // Comment out mkl.h, which can only run on solar3
00095 #include <mkl_blas.h>
00096 #include <mkl_service.h>
00097 #include <mkl_lapack.h>
00098 #include <mkl_vml_functions.h>
00099 #include <omp.h>
00100 
00101 #define PI              (M_PI)
00102 #define RADSINDEG       (PI/180.)
00103 #define RAD2ARCSEC      (648000./M_PI)
00104 #define SECINDAY        (86400.)
00105 #define FOURK           (4096)
00106 #define FOURK2    (16777216)
00107 
00108 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00109 
00110 // Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
00111 #define NYQVIST     (0.015)
00112 
00113 // Maximum variation of LONDTMAX-LONDTMIN
00114 #define MAXLONDIFF  (1.2e-4)
00115 
00116 // Some other things
00117 #ifndef MIN
00118 #define MIN(a,b) (((a)<(b)) ? (a) : (b))
00119 #endif
00120 #ifndef MAX
00121 #define MAX(a,b) (((a)>(b)) ? (a) : (b))
00122 #endif
00123 
00124 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
00125 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00126 
00127 #define kNotSpecified "Not Specified"
00128 
00129 // Macros for WCS transformations.  assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
00130 // and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
00131 // PIX_X and PIX_Y are CCD pixel addresses, WX and WY are arc-sec W and N on the Sun from disk center.
00132 #define PIX_X(wx,wy) ((((wx-crvalx)*cosa + (wy-crvaly)*sina)/cdelt)+crpix1)
00133 #define PIX_Y(wx,wy) ((((wy-crvaly)*cosa - (wx-crvalx)*sina)/cdelt)+crpix2)
00134 #define WX(pix_x,pix_y) (((pix_x-crpix1)*cosa - (pix_y-crpix2)*sina)*cdelt+crvalx)
00135 #define WY(pix_x,pix_y) (((pix_y-crpix2)*cosa + (pix_x-crpix1)*sina)*cdelt+crvaly)
00136 
00137 #define DISAMB_AZI      1
00138 #define XSCALE          0.03
00139 #define YSCALE          0.03
00140 #define NBIN            3
00141 #define INTERP          0
00142 #define dpath    "/home/jsoc/cvs/Development/JSOC"
00143 
00144 
00145 /* ========================================================================================================== */
00146 
00147 // Space weather keywords
00148 struct swIndex {
00149     float mean_vf;
00150     float count_mask;
00151     float absFlux;
00152     float mean_hf;
00153     float mean_gamma;
00154     float mean_derivative_btotal;
00155     float mean_derivative_bh;
00156     float mean_derivative_bz;
00157     float mean_jz;
00158     float us_i;
00159     float mean_alpha;
00160     float mean_ih;
00161     float total_us_ih;
00162     float total_abs_ih;
00163     float totaljz;
00164     float totpot;
00165     float meanpot;
00166     float area_w_shear_gt_45;
00167     float meanshear_angle;
00168     float area_w_shear_gt_45h;
00169     float meanshear_angleh;
00170     float mean_derivative_btotal_err;
00171     float mean_vf_err;
00172     float mean_gamma_err;
00173     float mean_derivative_bh_err;
00174     float mean_derivative_bz_err;
00175     float mean_jz_err;
00176     float us_i_err;
00177     float mean_alpha_err;
00178     float mean_ih_err;
00179     float total_us_ih_err;
00180     float total_abs_ih_err;
00181     float totaljz_err;
00182     float meanpot_err;
00183     float totpot_err;
00184     float meanshear_angle_err;
00185     float Rparam;
00186     float totfx;
00187     float totfy;
00188     float totfz;
00189     float totbsq;
00190     float epsx;
00191     float epsy;
00192     float epsz;
00193 };
00194 
00195 // Mapping method
00196 enum projection {
00197     carree,
00198     cassini,
00199     mercator,
00200     cyleqa,
00201     sineqa,
00202     gnomonic,
00203     postel,
00204     stereographic,
00205     orthographic,
00206     lambert
00207 };
00208 
00209 // WSC code
00210 char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
00211     "SIN", "ZEA"};
00212 
00213 // Ephemeris information
00214 struct ephemeris {
00215     double disk_lonc, disk_latc;
00216     double disk_xc, disk_yc;
00217     double rSun, asd, pa;
00218 };
00219 
00220 // Mapping information
00221 struct mapInfo {
00222     float xc, yc;       // reference point: center
00223     int nrow, ncol;     // size
00224     float xscale, yscale;   // scale
00225     int nbin;
00226     enum projection proj;   // projection method
00227     struct ephemeris ephem;     // ephemeris info
00228     float *xi_out, *zeta_out;   // coordinate on full disk image to sample at
00229 };
00230 
00231 /* ========================================================================================================== */
00232 
00233 /* Get all input data series */
00234 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, DRMS_RecordSet_t **bharpRS_ptr,
00235                char *mharpQuery, char *bharpQuery);
00236 
00237 /* Check if Mharp and Bharp match */
00238 int compareHarp(DRMS_RecordSet_t *mharpRS, DRMS_RecordSet_t *bharpRS);
00239 
00240 /* Get other data series */
00241 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS);
00242 
00243 /* Find record from record set with given T_rec */
00244 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
00245 
00246 /* Create CEA record */
00247 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
00248                     DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
00249                     DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
00250 
00251 /* Mapping single segment, wrapper */
00252 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
00253               struct mapInfo *mInfo, char *segName);
00254 
00255 /* Mapping vector magnetogram, wrapper */
00256 int mapVectorB(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
00257 
00258 /* Mapping errors of vector magnetogram, wraper */
00259 int mapVectorBErr(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
00260 
00261 /* Determine reference point coordinate and patch size according to input */
00262 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo);
00263 
00264 /* Get ephemeris information */
00265 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem);
00266 
00267 /* Compute the coordinates at which the full disk image is sampled */
00268 void findCoord(struct mapInfo *mInfo);
00269 
00270 /* Mapping function */
00271 int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
00272 
00273 /* Performing local vector transformation */
00274 void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
00275 
00276 /* Map and propogate errors */
00277 int getBErr(float *bx_err, float *by_err, float *bz_err,
00278             DRMS_Record_t *inRec, struct mapInfo *mInfo);
00279 
00280 /* Read full disk vector magnetogram */
00281 int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img);
00282 
00283 /* Read variances and covariances of vector magnetograms */
00284 int readVectorBErr(DRMS_Record_t *bharpRec,
00285                    float *bT, float *bI, float *bA,
00286                    float *errbT, float *errbI, float *errbA,
00287                    float *errbTbI, float *errbTbA, float *errbIbA);
00288 
00289 // ===================
00290 
00291 /* Create Cutout record */
00292 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
00293                     DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
00294                     DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
00295 
00296 /* Get cutout and write segment */
00297 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName);
00298 
00299 // ===================
00300 
00301 /* Compute space weather indices, no error checking for now */
00302 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo);
00303 
00304 /* Set space weather indices, no error checking for now */
00305 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
00306 
00307 /* Set all keywords, no error checking for now */
00308 // Changed Dec 30 XS
00309 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
00310 
00311 // ===================
00312 
00313 /* Nearest neighbor interpolation */
00314 float nnb (float *f, int nx, int ny, double x, double y);
00315 
00316 /* Wrapper for Jesper's rebin code */
00317 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss);
00318 
00319 /* ========================================================================================================== */
00320 
00321 /* Remap segment names */
00322 #define BR_SEG_CEA      "Br"
00323 #define BT_SEG_CEA      "Bt"
00324 #define BP_SEG_CEA      "Bp"
00325 #define BR_ERR_SEG_CEA      "Br_err"
00326 #define BT_ERR_SEG_CEA      "Bt_err"
00327 #define BP_ERR_SEG_CEA      "Bp_err"
00328 
00329 /* Cutout segment names, input identical to output */
00330 char *MharpSegs[] = {"magnetogram", "bitmap"};
00331 char *BharpSegs[] = {"inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
00332     "damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
00333     "conv_flag", // fixed
00334     "info_map", "confid_map",
00335     "inclination_err", "azimuth_err", "field_err", "vlos_err", "alpha_err",
00336     "field_inclination_err", "field_az_err", "inclin_azimuth_err",
00337     "field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
00338     "disambig", "conf_disambig"};
00339 // For stats
00340 char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
00341     "inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
00342     "damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
00343     "conv_flag", // fixed
00344     "info_map", "confid_map",
00345     "inclination_err", "azimuth_err", "field_err", "vlos_err", "alpha_err",
00346     "field_inclination_err", "field_az_err", "inclin_azimuth_err",
00347     "field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
00348     "disambig", "conf_disambig"};
00349 char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
00350     BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA, "conf_disambig"};
00351 // For Bunits, added Dec 30 XS
00352 char *CutBunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",
00353     "degree", "degree", "Mx/cm^2", "cm/s", "mA", " ",
00354     "length units", "DN/s", "DN/s", " ", " ",
00355     " ",
00356     " ", " ",
00357     "degree", "degree", "Mx/cm^2", "cm/s", " ",
00358     " ", " ", " ",
00359     " ", " ", " ",
00360     " ", " "};
00361 char *CEABunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",
00362     "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", " "};      // Mar 4 2014 XS
00363 
00364 /* ========================================================================================================== */
00365 
00366 char *module_name = "sharp";
00367 int seed;
00368 
00369 int fullDisk;       // full disk mode
00370 int amb4err;      // Use azimuth disambiguation for error propagation, default is 0 for patch and 1 for FD
00371 
00372 ModuleArgs_t module_args[] =
00373 {
00374     {ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
00375     {ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},
00376     {ARG_STRING, "b", kNotSpecified, "Input B series, if set, overrides bharp."},
00377     {ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},
00378     {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
00379     {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
00380     {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
00381     {ARG_INT,    "seed", "987654", "Seed for the random number generator."},
00382     {ARG_INT,   "f_amb4err", "0", "Force using disambiguation in error propagation"},     // Mar 4 2014 XS
00383     {ARG_END}
00384 };
00385 
00386 int DoIt(void)
00387 {
00388     int errbufstat = setvbuf(stderr, NULL, _IONBF, BUFSIZ);
00389     int outbufstat = setvbuf(stdout, NULL, _IONBF, BUFSIZ);
00390     
00391     int status = DRMS_SUCCESS;
00392     int nrecs, irec;
00393     
00394     char *mharpQuery, *bharpQuery, *bQuery;
00395     char *dopQuery, *contQuery;
00396     char *sharpCeaQuery, *sharpCutQuery;
00397     
00398     DRMS_RecordSet_t *mharpRS = NULL, *bharpRS = NULL;
00399     DRMS_RecordSet_t *dopRS = NULL, *contRS = NULL;
00400     
00401     /* Get parameters */
00402     
00403     mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
00404     bharpQuery = (char *) params_get_str(&cmdparams, "bharp");
00405     bQuery = (char *) params_get_str(&cmdparams, "b");
00406     dopQuery = (char *) params_get_str(&cmdparams, "dop");
00407     contQuery = (char *) params_get_str(&cmdparams, "cont");
00408     sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
00409     sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
00410     
00411     seed = params_get_int(&cmdparams, "seed");
00412     int f_amb4err = params_get_int(&cmdparams, "f_amb4err");
00413     
00414     /* Get input data, check everything */
00415     
00416     // Full disk mode if "b" is set
00417     if (strcmp(bQuery, kNotSpecified)) {
00418         fullDisk = 1;
00419         bharpQuery = bQuery;
00420         //        SHOW(bharpQuery); SHOW("\n");
00421         SHOW("Full disk mode\n");
00422     } else {
00423         fullDisk = 0;
00424         SHOW("Harp mode\n");
00425     }
00426     
00427     // Mar 4 2014
00428     if (f_amb4err == 0) {         // no forcing, 0 for patch and 1 for FD
00429         amb4err = fullDisk ? 1 : 0;
00430     } else {
00431         amb4err = 1;
00432     }
00433     printf("amb4err=%d\n", amb4err);
00434     
00435     // Bharp point to B if full disk
00436     if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
00437         DIE("Input harp data error.");
00438     nrecs = mharpRS->n;
00439     
00440     if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
00441         DIE("Input doppler data error.");
00442     
00443     if (getInputRS_aux(&contRS, contQuery, mharpRS))
00444         DIE("Input continuum data error.");
00445     
00446     /* Start */
00447     
00448     printf("==============\nStart. %d image(s) in total.\n", nrecs);
00449     
00450     for (irec = 0; irec < nrecs; irec++) {
00451         
00452         /* Records in work */
00453         
00454         DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;
00455         
00456         mharpRec = mharpRS->records[irec];
00457         
00458         TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
00459         
00460         if (!fullDisk) {
00461             bharpRec = bharpRS->records[irec];
00462         } else {
00463             if (getInputRec_aux(&bharpRec, bharpRS, trec)) {     // Bharp point to full disk B
00464                 printf("Fetching B failed, image #%d skipped.\n", irec);
00465                 continue;
00466             }
00467         }
00468         
00469         struct swIndex swKeys;
00470         
00471         DRMS_Record_t *dopRec = NULL, *contRec = NULL;
00472         
00473         if (getInputRec_aux(&dopRec, dopRS, trec)) {
00474             printf("Fetching Doppler failed, image #%d skipped.\n", irec);
00475             continue;
00476         }
00477         if (getInputRec_aux(&contRec, contRS, trec)) {
00478             printf("Fetching continuum failed, image #%d skipped.\n", irec);
00479             continue;
00480         }
00481         
00482         /* Create CEA record */
00483         
00484         DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
00485         if (status) {       // if failed
00486             printf("Creating CEA failed, image #%d skipped.\n", irec);
00487             continue;
00488         }
00489         
00490         if (createCeaRecord(mharpRec, bharpRec, dopRec, contRec, sharpCeaRec, &swKeys)) {       // do the work
00491             printf("Creating CEA failed, image #%d skipped.\n", irec);
00492             drms_close_record(sharpCeaRec, DRMS_FREE_RECORD);
00493             continue;
00494         }       // swKeys updated here
00495         
00496         drms_close_record(sharpCeaRec, DRMS_INSERT_RECORD);
00497         
00498         /* Create Cutout record */
00499         
00500         DRMS_Record_t *sharpCutRec = drms_create_record(drms_env, sharpCutQuery, DRMS_PERMANENT, &status);
00501         if (status) {       // if failed
00502             printf("Creating cutout failed, image #%d skipped.\n", irec);
00503             continue;
00504         }
00505         
00506         if (createCutRecord(mharpRec, bharpRec, dopRec, contRec, sharpCutRec, &swKeys)) {       // do the work
00507             printf("Creating cutout failed, image #%d skipped.\n", irec);
00508             drms_close_record(sharpCutRec, DRMS_FREE_RECORD);
00509             continue;
00510         }       // swKeys used here
00511         
00512         drms_close_record(sharpCutRec, DRMS_INSERT_RECORD);
00513         
00514         /* Done */
00515         
00516         printf("Image #%d done.\n", irec);
00517         
00518     } // irec
00519     
00520     
00521     drms_close_records(mharpRS, DRMS_FREE_RECORD);
00522     drms_close_records(bharpRS, DRMS_FREE_RECORD);
00523     drms_close_records(dopRS, DRMS_FREE_RECORD);                // Dec 18 2012
00524     drms_close_records(contRS, DRMS_FREE_RECORD);               // Dec 18 2012
00525     
00526     return 0;
00527     
00528 }   // DoIt
00529 
00530 
00531 // ===================================================================
00532 // ===================================================================
00533 // ===================================================================
00534 
00535 
00536 /*
00537  * Get input data series, including mHarp and bharp
00538  * Need all records to match, otherwise quit
00539  *
00540  */
00541 
00542 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, DRMS_RecordSet_t **bharpRS_ptr,
00543                char *mharpQuery, char *bharpQuery)
00544 {
00545     
00546     int status = 0;
00547     
00548     *mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
00549     if (status || (*mharpRS_ptr)->n == 0) return 1;
00550     
00551     if (fullDisk) {
00552         if (getInputRS_aux(bharpRS_ptr, bharpQuery, *mharpRS_ptr)) return 1;
00553     } else {
00554         *bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
00555         if (status || (*bharpRS_ptr)->n == 0) return 1;
00556         if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
00557     }
00558     
00559     return 0;
00560     
00561 }
00562 
00563 /*
00564  * Check if Mharp and Bharp match
00565  *
00566  */
00567 
00568 int compareHarp(DRMS_RecordSet_t *mharpRS, DRMS_RecordSet_t *bharpRS)
00569 {
00570     
00571     int status = 0;
00572     int nrecs = mharpRS->n;
00573     
00574     DRMS_Record_t *mharpRec_t = NULL, *bharpRec_t = NULL;       // temporary recs for utility
00575     
00576     if (bharpRS->n != nrecs) {
00577         return 1;       // return 1 if different
00578     }
00579     
00580     for (int i = 0; i < nrecs; i++) {
00581         mharpRec_t = mharpRS->records[i];
00582         bharpRec_t = bharpRS->records[i];
00583         if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
00584              drms_getkey_int(bharpRec_t, "HARPNUM", &status)) ||
00585             (drms_getkey_time(mharpRec_t, "T_REC", &status) !=
00586              drms_getkey_time(bharpRec_t, "T_REC", &status)))
00587         {
00588             return 1;
00589         }
00590     }
00591     
00592     return 0;
00593     
00594 }
00595 
00596 /*
00597  * Get other data series, check all T_REC are available
00598  *
00599  */
00600 
00601 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS)
00602 {
00603     
00604     int status = 0;
00605     
00606     *inRS_ptr = drms_open_records(drms_env, inQuery, &status);
00607     if (status || (*inRS_ptr)->n == 0) return status;
00608     
00609     // Check if all T_rec are available, need to match both ways
00610     int n = harpRS->n, n0 = (*inRS_ptr)->n;
00611     
00612     for (int i0 = 0; i0 < n0; i0++) {
00613         DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
00614         TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
00615         TIME trec = 0;
00616         for (int i = 0; i < n; i++) {
00617             DRMS_Record_t *harpRec = harpRS->records[i];
00618             trec = drms_getkey_time(harpRec, "T_REC", &status);
00619             if (fabs(trec0 - trec) < 10) break;
00620         }
00621         if (fabs(trec0 - trec) >= 10) return 1;
00622     }
00623     
00624     for (int i = 0; i < n; i++) {
00625         DRMS_Record_t *harpRec = harpRS->records[i];
00626         TIME trec = drms_getkey_time(harpRec, "T_REC", &status);
00627         TIME trec0 = 0;
00628         for (int i0 = 0; i0 < n0; i0++) {
00629             DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
00630             trec0 = drms_getkey_time(inRec, "T_REC", &status);
00631             if (fabs(trec0 - trec) < 10) break;
00632         }
00633         if (fabs(trec0 - trec) >= 10) return 1;
00634     }
00635     
00636     return 0;
00637     
00638 }
00639 
00640 /*
00641  * Find record from record set with given T_rec
00642  *
00643  */
00644 
00645 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec)
00646 {
00647     
00648     int status = 0;
00649     
00650     int n = inRS->n;
00651     for (int i = 0; i < n; i++) {
00652         *inRec_ptr = inRS->records[i];
00653         TIME trec0 = drms_getkey_time((*inRec_ptr), "T_REC", &status);
00654         if (fabs(trec0 - trec) < 10) return 0;
00655     }
00656     
00657     return 1;
00658     
00659 }
00660 
00661 
00662 
00663 
00664 /*
00665  * Create CEA record: top level subroutine
00666  * Also compute all the space weather keywords here
00667  *
00668  */
00669 
00670 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
00671                     DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
00672                     DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
00673 {
00674     
00675     int status = 0;
00676     DRMS_Segment_t *inSeg;
00677     DRMS_Array_t *inArray;
00678     
00679     struct mapInfo mInfo;
00680     mInfo.proj = (enum projection) cyleqa;      // projection method
00681     mInfo.xscale = XSCALE;
00682     mInfo.yscale = YSCALE;
00683     
00684     int ncol0, nrow0;       // oversampled map size
00685     
00686     // Get ephemeris
00687     
00688     if (getEphemeris(mharpRec, &(mInfo.ephem))) {
00689         SHOW("CEA: get ephemeris error\n");
00690         return 1;
00691     }
00692     
00693     // Find position
00694     
00695     if (findPosition(mharpRec, &mInfo)) {
00696         SHOW("CEA: find position error\n");
00697         return 1;
00698     }
00699     
00700     // ========================================
00701     // Do this for all bitmaps, Aug 12 2013 XS
00702     // ========================================
00703     
00704     mInfo.nbin = 1;         // for bitmaps. suppress anti-aliasing
00705     ncol0 = mInfo.ncol;
00706     nrow0 = mInfo.nrow;
00707     
00708     mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
00709     mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
00710     
00711     findCoord(&mInfo);      // compute it here so it could be shared by the following 4 functions
00712     
00713     if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
00714         SHOW("CEA: mapping bitmap error\n");
00715         return 1;
00716     }
00717     printf("Bitmap mapping done.\n");
00718     
00719     if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
00720         SHOW("CEA: mapping conf_disambig error\n");
00721         return 1;
00722     }
00723     printf("Conf disambig mapping done.\n");
00724     
00725     free(mInfo.xi_out);
00726     free(mInfo.zeta_out);
00727     
00728     // ========================================
00729     // Do this again for floats, Aug 12 2013 XS
00730     // ========================================
00731     // Create xi_out, zeta_out array in mInfo:
00732     // Coordinates to sample in original full disk image
00733     
00734     mInfo.nbin = NBIN;
00735     ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
00736     nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
00737     
00738     mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
00739     mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
00740     
00741     findCoord(&mInfo);      // compute it here so it could be shared by the following 4 functions
00742     
00743     // Mapping single segment: Mharp, etc.
00744     
00745     if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
00746         SHOW("CEA: mapping magnetogram error\n");
00747         return 1;
00748     }
00749     printf("Magnetogram mapping done.\n");
00750     
00751     if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
00752         SHOW("CEA: mapping dopplergram error\n");
00753         return 1;
00754     }
00755     printf("Dopplergram mapping done.\n");
00756     
00757     if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {
00758         SHOW("CEA: mapping continuum error\n");
00759         return 1;
00760     }
00761     printf("Intensitygram mapping done.\n");
00762     
00763     // Mapping vector B
00764     
00765     if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
00766         SHOW("CEA: mapping vector B error\n");
00767         return 1;
00768     }
00769     printf("Vector B mapping done.\n");
00770     
00771     // Mapping vector B errors
00772     
00773     if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) {
00774         SHOW("CEA: mapping vector B uncertainty error\n");
00775         return 1;
00776     }
00777     printf("Vector B error done.\n");
00778     
00779     // Keywords & Links
00780     
00781     drms_copykey(sharpRec, mharpRec, "T_REC");
00782     drms_copykey(sharpRec, mharpRec, "HARPNUM");
00783     
00784     if (fullDisk) {
00785         DRMS_Link_t *bLink = hcon_lookup_lower(&sharpRec->links, "B");
00786         if (bLink) drms_link_set("B", sharpRec, bharpRec);
00787     } else {
00788         DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
00789         if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
00790     }
00791     DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
00792     if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);
00793     
00794     setKeys(sharpRec, mharpRec, bharpRec, &mInfo);            // Set all other keywords
00795     drms_copykey(sharpRec, mharpRec, "QUALITY");        // copied from los records
00796     
00797     // Space weather
00798     
00799     computeSWIndex(swKeys_ptr, sharpRec, &mInfo);       // compute it!
00800     printf("Space weather indices done.\n");
00801     
00802     setSWIndex(sharpRec, swKeys_ptr);   // Set space weather indices
00803     
00804     // Stats
00805     
00806     int nCEASegs = ARRLENGTH(CEASegs);
00807     for (int iSeg = 0; iSeg < nCEASegs; iSeg++) {
00808         DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
00809         DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
00810         int stat = set_statistics(outSeg, outArray, 1);
00811         //      printf("%d => %d\n", iSeg, stat);
00812         drms_free_array(outArray);
00813     }
00814     
00815     free(mInfo.xi_out);
00816     free(mInfo.zeta_out);
00817     return 0;
00818     
00819 }
00820 
00821 
00822 /*
00823  * Mapping a single segment
00824  * Read in full disk image, utilize mapImage for mapping
00825  * then write the segment out, segName same in in/out Rec
00826  *
00827  */
00828 
00829 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
00830               struct mapInfo *mInfo, char *segName)
00831 {
00832     
00833     int status = 0;
00834     int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
00835     int dims[2] = {nx, ny};
00836     int interpOpt = INTERP;     // Aug 12 XS, default, overridden below for bitmaps and conf_disambig
00837     
00838     // Input full disk array
00839     
00840     DRMS_Segment_t *inSeg = NULL;
00841     inSeg = drms_segment_lookup(inRec, segName);
00842     if (!inSeg) return 1;
00843     
00844     DRMS_Array_t *inArray = NULL;
00845     inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00846     if (!inArray) return 1;
00847     
00848     if (!strcmp(segName, "conf_disambig") || !strcmp(segName, "bitmap")) {
00849         // Moved out so it works for FD conf_disambig as well
00850         // Jan 2 2014 XS
00851         interpOpt = 3;      // Aug 12 XS, near neighbor
00852     }
00853     
00854     float *inData;
00855     int xsz = inArray->axis[0], ysz = inArray->axis[1];
00856     if ((xsz != FOURK) || (ysz != FOURK)) {     // for bitmap, make tmp full disk
00857         float *inData0 = (float *) inArray->data;
00858         inData = (float *) (calloc(FOURK2, sizeof(float)));
00859         int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
00860         int y0 = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1;
00861         int ind_map;
00862         for (int row = 0; row < ysz; row++) {
00863             for (int col = 0; col < xsz; col++) {
00864                 ind_map = (row + y0) * FOURK + (col + x0);
00865                 inData[ind_map] = inData0[row * xsz + col];
00866             }
00867         }
00868         drms_free_array(inArray); inArray = NULL;
00869     } else {
00870         inData = (float *) inArray->data;
00871     }
00872     
00873     // Mapping
00874     
00875     float *map = (float *) (malloc(nxny * sizeof(float)));
00876     if (performSampling(map, inData, mInfo, interpOpt))     // Add interpOpt for different types, Aug 12 XS
00877     {if (inArray) drms_free_array(inArray); free(map); return 1;}
00878     
00879     // Write out
00880     
00881     DRMS_Segment_t *outSeg = NULL;
00882     outSeg = drms_segment_lookup(sharpRec, segName);
00883     if (!outSeg) return 1;
00884     
00885     //  DRMS_Type_t arrayType = outSeg->info->type;
00886     DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
00887     if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
00888     
00889     // convert to needed data type
00890     
00891     //  drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);     // Jan 02 2013
00892     
00893     outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
00894     //  outArray->parent_segment = outSeg;
00895     outArray->israw = 0;        // always compressed
00896     outArray->bzero = outSeg->bzero;
00897     outArray->bscale = outSeg->bscale;
00898     
00899     status = drms_segment_write(outSeg, outArray, 0);
00900     if (status) return 0;
00901     
00902     if (inArray) drms_free_array(inArray);
00903     if ((xsz != FOURK) || (ysz != FOURK)) free(inData);         // Dec 18 2012
00904     if (outArray) drms_free_array(outArray);
00905     return 0;
00906     
00907 }
00908 
00909 
00910 /*
00911  * Mapping vector magnetogram
00912  *
00913  */
00914 
00915 int mapVectorB(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
00916 {
00917     
00918     int status = 0;
00919     int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
00920     int dims[2] = {nx, ny};
00921     
00922     // Read in segments, filling factor assume to be 1
00923     
00924     float *bx_img = (float *) (malloc(FOURK2 * sizeof(float)));
00925     float *by_img = (float *) (malloc(FOURK2 * sizeof(float)));
00926     float *bz_img = (float *) (malloc(FOURK2 * sizeof(float)));
00927     
00928     if (readVectorB(bharpRec, bx_img, by_img, bz_img)) {
00929         printf("Read full disk image error\n");
00930         free(bx_img); free(by_img); free(bz_img);
00931         return 1;
00932     }
00933     
00934     // Mapping
00935     
00936     float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation
00937     
00938     bx_map = (float *) (malloc(nxny * sizeof(float)));
00939     if (performSampling(bx_map, bx_img, mInfo, INTERP))
00940     {free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
00941     
00942     by_map = (float *) (malloc(nxny * sizeof(float)));
00943     if (performSampling(by_map, by_img, mInfo, INTERP))
00944     {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
00945     
00946     bz_map = (float *) (malloc(nxny * sizeof(float)));
00947     if (performSampling(bz_map, bz_img, mInfo, INTERP))
00948     {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
00949     
00950     free(bx_img); free(by_img); free(bz_img);
00951     
00952     // Vector transform
00953     
00954     vectorTransform(bx_map, by_map, bz_map, mInfo);
00955     
00956     for (int i = 0; i < nxny; i++) by_map[i] *= -1;     // positive theta pointing south
00957     
00958     // Write out
00959     
00960     DRMS_Segment_t *outSeg;
00961     DRMS_Array_t *outArray;
00962     
00963     float *data_prt[3] = {bx_map, by_map, bz_map};
00964     char *segName[3] = {BP_SEG_CEA, BT_SEG_CEA, BR_SEG_CEA};
00965     
00966     for (int iSeg = 0; iSeg < 3; iSeg++) {
00967         outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
00968         outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
00969         outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
00970         //      outArray->parent_segment = outSeg;
00971         outArray->israw = 0;
00972         outArray->bzero = outSeg->bzero;
00973         outArray->bscale = outSeg->bscale;
00974         status = drms_segment_write(outSeg, outArray, 0);
00975         if (status) return 1;
00976         drms_free_array(outArray);
00977     }
00978     
00979     //
00980     
00981     return 0;
00982     
00983 }
00984 
00985 
00986 /*
00987  * Mapping vector magnetogram errors
00988  *
00989  */
00990 
00991 int mapVectorBErr(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
00992 {
00993     
00994     int status = 0;
00995     
00996     int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
00997     int dims[2] = {nx, ny};
00998     
00999     // Compute propogated errors, using nearest neighbour interpolation
01000     
01001     float *bx_err = (float *) (malloc(nxny * sizeof(float)));
01002     float *by_err = (float *) (malloc(nxny * sizeof(float)));
01003     float *bz_err = (float *) (malloc(nxny * sizeof(float)));
01004     
01005     if (getBErr(bx_err, by_err, bz_err, bharpRec, mInfo)) {
01006         free(bx_err); free(by_err); free(bz_err);
01007         return 1;
01008     }
01009     
01010     // Write out
01011     
01012     DRMS_Segment_t *outSeg;
01013     DRMS_Array_t *outArray;
01014     
01015     float *data_prt[3] = {bx_err, by_err, bz_err};
01016     char *segName[3] = {BP_ERR_SEG_CEA, BT_ERR_SEG_CEA, BR_ERR_SEG_CEA};
01017     
01018     for (int iSeg = 0; iSeg < 3; iSeg++) {
01019         outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
01020         outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
01021         outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
01022         //      outArray->parent_segment = outSeg;
01023         outArray->israw = 0;
01024         outArray->bzero = outSeg->bzero;
01025         outArray->bscale = outSeg->bscale;
01026         status = drms_segment_write(outSeg, outArray, 0);
01027         if (status) return 1;
01028         drms_free_array(outArray);
01029     }
01030     
01031     //
01032     
01033     return 0;
01034     
01035 }
01036 
01037 
01038 
01039 /*
01040  * Determine reference point coordinate and patch size according to keywords
01041  * xc, yc are the coordinate of patch center, in degrees
01042  * ncol and nrow are the final size
01043  *
01044  */
01045 
01046 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo)
01047 {
01048     
01049     int status = 0;
01050     int harpnum = drms_getkey_int(inRec, "HARPNUM", &status);
01051     TIME trec = drms_getkey_time(inRec, "T_REC", &status);
01052     float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
01053     
01054     /* Center coord */
01055     // Changed into double Jun 16 2014 XS
01056     
01057     double minlon = drms_getkey_double(inRec, "LONDTMIN", &status); if (status) return 1;       // Stonyhurst lon
01058     double maxlon = drms_getkey_double(inRec, "LONDTMAX", &status); if (status) return 1;
01059     double minlat = drms_getkey_double(inRec, "LATDTMIN", &status); if (status) return 1;
01060     double maxlat = drms_getkey_double(inRec, "LATDTMAX", &status); if (status) return 1;
01061     
01062     // A bug fixer for HARP (per M. Turmon)
01063     // When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
01064     // Also keywords such as "SIZE" will be NaN
01065     // We compute minlon & minlat then by
01066     // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
01067     
01068     //  float psize = drms_getkey_float(inRec, "SIZE", &status);
01069     //  if (psize != psize) {
01070     
01071     if (minlon != minlon || maxlon != maxlon) {     // check lons instead of SIZE
01072         TIME t0 = drms_getkey_time(inRec, "T_FRST1", &status); if (status) return 1;            // changed from T_FRST to T_FRST1, T_FRST may not exist
01073         double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
01074         char firstRecQuery[100], t0_str[100];
01075         sprint_time(t0_str, t0, "TAI", 0);
01076         snprintf(firstRecQuery, 100, "%s[%d][%s]", inRec->seriesinfo->seriesname, harpnum, t0_str);
01077         DRMS_RecordSet_t *tmpRS = drms_open_records(drms_env, firstRecQuery, &status);
01078         if (status || tmpRS->n != 1) return 1;
01079         DRMS_Record_t *tmpRec = tmpRS->records[0];
01080         double minlon0 = drms_getkey_double(tmpRec, "LONDTMIN", &status); if (status) return 1;
01081         double maxlon0 = drms_getkey_double(tmpRec, "LONDTMAX", &status); if (status) return 1;
01082         minlon = minlon0 + (trec - t0) * omega / SECINDAY;
01083         maxlon = maxlon0 + (trec - t0) * omega / SECINDAY;
01084         printf("%s, %f, %f\n", firstRecQuery, minlon, maxlon);
01085     }
01086     
01087     mInfo->xc = (maxlon + minlon) / 2. + disk_lonc;
01088     mInfo->yc = (maxlat + minlat) / 2.;
01089     
01090     /* Size */
01091     // Rounded to 1.d3 precision first. Jun 16 2014 XS
01092     // The previous fix does not work. LONDTMAX-LONDTMIN varies from frame to frame
01093     // Need to find out the maximum possible difference, MAXLONDIFF (1.2e-4)
01094     // Now, ncol = (maxlon-minlon)/xscale, if the decimal part is outside 0.5 \pm (MAXLONDIFF/xscale)
01095     // proceed as it is. else, all use floor on ncol
01096     
01097     float dpix = (MAXLONDIFF / mInfo->xscale) * 1.5;        // "danger zone"
01098     float ncol = (maxlon - minlon) / mInfo->xscale;
01099     float d_ncol = fabs(ncol - floor(ncol) - 0.5);          // distance to 0.5
01100     if (d_ncol < dpix) {
01101         mInfo->ncol = floor(ncol);
01102     } else {
01103         mInfo->ncol = round(ncol);
01104     }
01105 
01106     mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);
01107     
01108     printf("xcol=%f, ncol=%d, nrow=%d\n", ncol, mInfo->ncol, mInfo->nrow);
01109     
01110     return 0;
01111     
01112 }
01113 
01114 
01115 /*
01116  * Fetch ephemeris info from a DRMS record
01117  * No error checking for now
01118  *
01119  */
01120 
01121 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem)
01122 {
01123     
01124     int status = 0;
01125     
01126     float crota2 = drms_getkey_float(inRec, "CROTA2", &status); // rotation
01127     double sina = sin(crota2 * RADSINDEG);
01128     double cosa = cos(crota2 * RADSINDEG);
01129     
01130     ephem->pa = - crota2 * RADSINDEG;
01131     ephem->disk_latc = drms_getkey_float(inRec, "CRLT_OBS", &status) * RADSINDEG;
01132     ephem->disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status) * RADSINDEG;
01133     
01134     float crvalx = 0.0;
01135     float crvaly = 0.0;
01136     float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
01137     float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
01138     float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy
01139     ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;      // Center of disk in pixel, starting at 0
01140     ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
01141     
01142     float dSun = drms_getkey_float(inRec, "DSUN_OBS", &status);
01143     float rSun_ref = drms_getkey_float(inRec, "RSUN_REF", &status);
01144     if (status) rSun_ref = 6.96e8;
01145     
01146     ephem->asd = asin(rSun_ref/dSun);
01147     ephem->rSun = asin(rSun_ref / dSun) * RAD2ARCSEC / cdelt;
01148     
01149     return 0;
01150     
01151 }
01152 
01153 
01154 /*
01155  * Compute the coordinates to be sampled on full disk image
01156  * mInfo->xi_out & mInfo->zeta_out
01157  * This is oversampled, its size is ncol0 & nrow0 as shown below
01158  *
01159  *
01160  */
01161 
01162 void findCoord(struct mapInfo *mInfo)
01163 {
01164     
01165     int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN
01166     int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
01167     
01168     float xscale0 = mInfo->xscale / mInfo->nbin * RADSINDEG;        // oversampling resolution
01169     float yscale0 = mInfo->yscale / mInfo->nbin * RADSINDEG;        // in rad
01170     
01171     double lonc = mInfo->xc * RADSINDEG;    // in rad
01172     double latc = mInfo->yc * RADSINDEG;
01173     
01174     double disk_lonc = (mInfo->ephem).disk_lonc;
01175     double disk_latc = (mInfo->ephem).disk_latc;
01176     
01177     double rSun = (mInfo->ephem).rSun;
01178     double disk_xc = (mInfo->ephem).disk_xc / rSun;
01179     double disk_yc = (mInfo->ephem).disk_yc / rSun;
01180     double pa = (mInfo->ephem).pa;
01181     
01182     // Temp pointers
01183     
01184     float *xi_out = mInfo->xi_out;
01185     float *zeta_out = mInfo->zeta_out;
01186     
01187     // start
01188     
01189     double x, y;        // map coord
01190     double lat, lon;    // helio coord
01191     double xi, zeta;    // image coord (for one point)
01192     
01193     int ind_map;
01194     
01195     for (int row0 = 0; row0 < nrow0; row0++) {
01196         for (int col0 = 0; col0 < ncol0; col0++) {
01197             
01198             ind_map = row0 * ncol0 + col0;
01199             
01200             x = (col0 + 0.5 - ncol0/2.) * xscale0;      // in rad
01201             y = (row0 + 0.5 - nrow0/2.) * yscale0;
01202             
01203             /* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
01204              * the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
01205              * projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
01206              * is the heliographic longitude and latitude of the map center. Both are in degree.
01207              */
01208             
01209             if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
01210                 xi_out[ind_map] = -1;
01211                 zeta_out[ind_map] = -1;
01212                 continue;
01213             }
01214             
01215             /* map the grid [lon, lat] in the heliographic coordinates to [xi, zeta], a point in the
01216              * image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
01217              */
01218             
01219             if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
01220                             disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
01221                 xi_out[ind_map] = -1;
01222                 zeta_out[ind_map] = -1;
01223                 continue;
01224             }
01225             
01226             xi_out[ind_map] = xi * rSun;
01227             zeta_out[ind_map] = zeta * rSun;
01228             
01229         }
01230     }
01231     
01232 }
01233 
01234 
01235 /*
01236  * Sampling function
01237  * oversampling by nbin, then binning using a Gaussian
01238  * save results in outData, always of float type
01239  *
01240  */
01241 
01242 int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
01243 {
01244     
01245     int status = 0;
01246     int ind_map;
01247     
01248     int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN
01249     int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
01250     
01251     // Changed Aug 12 2013, XS, for bitmaps
01252     float *outData0;
01253     if (interpOpt == 3 && mInfo->nbin == 1) {
01254         outData0 = outData;
01255     } else {
01256         outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
01257     }
01258     
01259     float *xi_out = mInfo->xi_out;
01260     float *zeta_out = mInfo->zeta_out;
01261     
01262     // Interpolation
01263     
01264     struct fint_struct pars;
01265     // Aug 12 2013, passed in as argument now
01266     
01267     switch (interpOpt) {
01268         case 0:         // Wiener, 6 order, 1 constraint
01269             init_finterpolate_wiener(&pars, 6, 1, 6, 2, 1, 1, NULL, dpath);
01270             break;
01271         case 1:         // Cubic convolution
01272             init_finterpolate_cubic_conv(&pars, 1., 3.);
01273             break;
01274         case 2:         // Bilinear
01275             init_finterpolate_linear(&pars, 1.);
01276             break;
01277         case 3:         // Near neighbor
01278             break;
01279         default:
01280             return 1;
01281     }
01282     
01283     printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
01284     if (interpOpt == 3) {           // Aug 6 2013, Xudong
01285         for (int row0 = 0; row0 < nrow0; row0++) {
01286             for (int col0 = 0; col0 < ncol0; col0++) {
01287                 ind_map = row0 * ncol0 + col0;
01288                 outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
01289             }
01290         }
01291     } else {
01292         finterpolate(&pars, inData, xi_out, zeta_out, outData0,
01293                      FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
01294     }
01295     
01296     // Rebinning, smoothing
01297     
01298     if (interpOpt == 3 && mInfo->nbin == 1) {
01299         return 0;
01300     } else {
01301         frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);        // Gaussian
01302         free(outData0);     // Dec 18 2012
01303     }
01304     
01305     //
01306     
01307     return 0;
01308     
01309 }
01310 
01311 
01312 /*
01313  * Performing local vector transformation
01314  *  xyz: z refers to vertical (radial) component, x EW (phi), y NS
01315  *
01316  */
01317 
01318 void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo)
01319 {
01320     
01321     int ncol = mInfo->ncol;
01322     int nrow = mInfo->nrow;
01323     
01324     float xscale = mInfo->xscale * RADSINDEG;       // in rad
01325     float yscale = mInfo->yscale * RADSINDEG;
01326     
01327     double lonc = mInfo->xc * RADSINDEG;    // in rad
01328     double latc = mInfo->yc * RADSINDEG;
01329     
01330     double disk_lonc = (mInfo->ephem).disk_lonc;
01331     double disk_latc = (mInfo->ephem).disk_latc;
01332     
01333     double rSun = (mInfo->ephem).rSun;
01334     double disk_xc = (mInfo->ephem).disk_xc / rSun;
01335     double disk_yc = (mInfo->ephem).disk_yc / rSun;
01336     double pa = (mInfo->ephem).pa;
01337     
01338     int ind_map;
01339     double x, y;
01340     double lat, lon;    // lat / lon for current point
01341     
01342     double bx_tmp, by_tmp, bz_tmp;
01343     
01344     //
01345     
01346     for (int row = 0; row < mInfo->nrow; row++) {
01347         for (int col = 0; col < mInfo->ncol; col++) {
01348             
01349             ind_map = row * mInfo->ncol + col;
01350             
01351             x = (col + 0.5 - ncol / 2.) * xscale;
01352             y = (row + 0.5 - nrow / 2.) * yscale;
01353             
01354             if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
01355                 bx_map[ind_map] = DRMS_MISSING_FLOAT;
01356                 by_map[ind_map] = DRMS_MISSING_FLOAT;
01357                 bz_map[ind_map] = DRMS_MISSING_FLOAT;
01358                 continue;
01359             }
01360             
01361             bx_tmp = by_tmp = bz_tmp = 0;
01362             
01363             img2helioVector (bx_map[ind_map], by_map[ind_map], bz_map[ind_map],
01364                              &bx_tmp, &by_tmp, &bz_tmp,
01365                              lon, lat, disk_lonc, disk_latc, pa);
01366             
01367             bx_map[ind_map] = bx_tmp;
01368             by_map[ind_map] = by_tmp;
01369             bz_map[ind_map] = bz_tmp;
01370             
01371         }
01372     }
01373     
01374 }
01375 
01376 
01377 
01378 /*
01379  * Map and propogate vector field errors
01380  *
01381  */
01382 
01383 int getBErr(float *bx_err, float *by_err, float *bz_err,
01384             DRMS_Record_t *inRec, struct mapInfo *mInfo)
01385 {
01386     
01387     int status = 0;
01388     
01389     // Get variances and covariances, filling factor assume to be 1
01390     
01391     float *bT = (float *) (malloc(FOURK2 * sizeof(float))); // field
01392     float *bI = (float *) (malloc(FOURK2 * sizeof(float))); // inclination
01393     float *bA = (float *) (malloc(FOURK2 * sizeof(float))); // azimuth
01394     
01395     float *errbT = (float *) (malloc(FOURK2 * sizeof(float)));
01396     float *errbI = (float *) (malloc(FOURK2 * sizeof(float)));
01397     float *errbA = (float *) (malloc(FOURK2 * sizeof(float)));
01398     
01399     float *errbTbI = (float *) (malloc(FOURK2 * sizeof(float)));
01400     float *errbTbA = (float *) (malloc(FOURK2 * sizeof(float)));
01401     float *errbIbA = (float *) (malloc(FOURK2 * sizeof(float)));
01402     
01403     if (readVectorBErr(inRec,
01404                        bT, bI, bA,
01405                        errbT, errbI, errbA,
01406                        errbTbI, errbTbA, errbIbA)) {
01407         printf("Read full disk variances & covariances error\n");
01408         free(bT); free(bI); free(bA);
01409         free(errbT); free(errbI); free(errbA);
01410         free(errbTbI); free(errbTbA); free(errbIbA);
01411         return 1;
01412     }
01413     
01414     // Size
01415     
01416     int ncol = mInfo->ncol;
01417     int nrow = mInfo->nrow;
01418     
01419     float xscale = mInfo->xscale * RADSINDEG;       // in rad
01420     float yscale = mInfo->yscale * RADSINDEG;
01421     
01422     double lonc = mInfo->xc * RADSINDEG;    // in rad
01423     double latc = mInfo->yc * RADSINDEG;
01424     
01425     double disk_lonc = (mInfo->ephem).disk_lonc;
01426     double disk_latc = (mInfo->ephem).disk_latc;
01427     
01428     double rSun = (mInfo->ephem).rSun;
01429     double disk_xc = (mInfo->ephem).disk_xc / rSun;
01430     double disk_yc = (mInfo->ephem).disk_yc / rSun;
01431     double pa = (mInfo->ephem).pa;
01432     
01433     // Start
01434     
01435     double x, y;          // map coord
01436     double lat, lon;      // spherical coord
01437     double xi, zeta;      // image coord, round to full pixel
01438     
01439     int ind_map, ind_img;
01440     
01441     double bpSigma2, btSigma2, brSigma2;        // variances after prop
01442     
01443     for (int row = 0; row < mInfo->nrow; row++) {
01444         for (int col = 0; col < mInfo->ncol; col++) {
01445             
01446             ind_map = row * mInfo->ncol + col;
01447             
01448             x = (col + 0.5 - ncol / 2.) * xscale;
01449             y = (row + 0.5 - nrow / 2.) * yscale;
01450             
01451             if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
01452                 bx_err[ind_map] = DRMS_MISSING_FLOAT;
01453                 by_err[ind_map] = DRMS_MISSING_FLOAT;
01454                 bz_err[ind_map] = DRMS_MISSING_FLOAT;
01455                 continue;
01456             }
01457             
01458             if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
01459                             disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
01460                 bx_err[ind_map] = DRMS_MISSING_FLOAT;
01461                 by_err[ind_map] = DRMS_MISSING_FLOAT;
01462                 bz_err[ind_map] = DRMS_MISSING_FLOAT;       // Mar 7
01463                 continue;
01464             }
01465             
01466             xi *= rSun; xi = round(xi);
01467             zeta *= rSun; zeta = round(zeta);     // nearest neighbor
01468             
01469             ind_img = round(zeta * FOURK + xi);
01470             
01471             if (errorprop(bT, bA, bI,
01472                           errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
01473                           lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
01474                           &btSigma2, &bpSigma2, &brSigma2)) {
01475                 bx_err[ind_map] = DRMS_MISSING_FLOAT;
01476                 by_err[ind_map] = DRMS_MISSING_FLOAT;
01477                 bz_err[ind_map] = DRMS_MISSING_FLOAT;
01478                 continue;
01479             }
01480             
01481             bx_err[ind_map] = sqrt(bpSigma2);
01482             by_err[ind_map] = sqrt(btSigma2);
01483             bz_err[ind_map] = sqrt(brSigma2);
01484             
01485         }
01486     }
01487     
01488     //
01489     
01490     free(bT); free(bI); free(bA);
01491     free(errbT); free(errbI); free(errbA);
01492     free(errbTbI); free(errbTbA); free(errbIbA);
01493     return 0;
01494     
01495 }
01496 
01497 
01498 
01499 /*
01500  * Read full disk vector magnetograms
01501  * Fill factor is 1, use default disambiguity resolution
01502  *
01503  */
01504 
01505 int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img)
01506 {
01507     
01508     int status = 0;
01509     
01510     DRMS_Segment_t *inSeg;
01511     DRMS_Array_t *inArray_ambig;
01512     DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
01513     
01514     char *ambig;
01515     float *bTotal, *bAzim, *bIncl;
01516     
01517     inSeg = drms_segment_lookup(inRec, "disambig");
01518     inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);
01519     if (status) return 1;
01520     ambig = (char *)inArray_ambig->data;
01521     
01522     inSeg = drms_segment_lookup(inRec, "field");
01523     inArray_bTotal = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
01524     if (status) return 1;
01525     bTotal = (float *)inArray_bTotal->data;
01526     
01527     inSeg = drms_segment_lookup(inRec, "azimuth");
01528     inArray_bAzim = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
01529     if (status) return 1;
01530     bAzim = (float *)inArray_bAzim->data;
01531     
01532     inSeg = drms_segment_lookup(inRec, "inclination");
01533     inArray_bIncl = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
01534     if (status) return 1;
01535     bIncl = (float *)inArray_bIncl->data;
01536     
01537     // Convert CCD xyz
01538     
01539     int llx, lly;       // lower-left corner
01540     int bmx, bmy;       // bitmap size
01541     
01542     if (fullDisk) {
01543         llx = lly = 0;
01544         bmx = bmy = FOURK;
01545     } else {
01546         llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
01547         lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
01548         bmx = inArray_ambig->axis[0];
01549         bmy = inArray_ambig->axis[1];
01550     }
01551     
01552     int kx, ky, kOff;
01553     int ix = 0, jy = 0, yOff = 0, iData = 0;
01554     int xDim = FOURK, yDim = FOURK;
01555     int amb = 0;
01556     
01557     for (jy = 0; jy < yDim; jy++)
01558     {
01559         ix = 0;
01560         yOff = jy * xDim;
01561         ky = jy - lly;
01562         for (ix = 0; ix < xDim; ix++)
01563         {
01564             iData = yOff + ix;
01565             kx = ix - llx;
01566             
01567             // zero azi pointing up, zero incl pointing out from sun
01568             bx_img[iData] = - bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * sin(bAzim[iData] * RADSINDEG);
01569             by_img[iData] = bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * cos(bAzim[iData] * RADSINDEG);
01570             bz_img[iData] = bTotal[iData] * cos(bIncl[iData] * RADSINDEG);
01571             
01572             // Disambiguation
01573             
01574             if (kx < 0 || kx >= bmx || ky < 0 || ky >= bmy) {
01575                 continue;
01576             } else {
01577                 kOff = ky * bmx + kx;
01578                 //              if (ambig[kOff] % 2) {      // 180
01579                 // Feb 12 2014, use bit #2 for full disk, lowest bit for patch
01580                 if (fullDisk) { amb = (ambig[kOff] / 4) % 2; } else { amb = ambig[kOff] % 2; }
01581                 if (amb) {              // Feb 12 2014, use bit #2
01582                     bx_img[iData] *= -1.; by_img[iData] *= -1.;
01583                 }
01584             }
01585         }
01586     }
01587     
01588     // Clean up
01589     
01590     drms_free_array(inArray_ambig);
01591     drms_free_array(inArray_bTotal);
01592     drms_free_array(inArray_bAzim);
01593     drms_free_array(inArray_bIncl);
01594     
01595     return 0;
01596     
01597 }
01598 
01599 
01600 /*
01601  * Read variances and covariances of vector magnetograms
01602  *
01603  */
01604 
01605 int readVectorBErr(DRMS_Record_t *inRec,
01606                    float *bT, float *bI, float *bA,
01607                    float *errbT, float *errbI, float *errbA,
01608                    float *errbTbI, float *errbTbA, float *errbIbA)
01609 {
01610     
01611     int status = 0;
01612     
01613     float *data_ptr[9];
01614     char *segName[9] = {"field", "inclination", "azimuth",
01615         "field_err", "inclination_err", "azimuth_err",
01616         "field_inclination_err", "field_az_err", "inclin_azimuth_err"};
01617     DRMS_Segment_t *inSegs[9];
01618     DRMS_Array_t *inArrays[9];
01619     
01620     // Read full disk images
01621     // Do we need disambig? Dec 30 XS
01622     
01623     for (int iSeg = 0; iSeg < 9; iSeg++) {
01624         
01625         inSegs[iSeg] = drms_segment_lookup(inRec, segName[iSeg]);
01626         inArrays[iSeg] = drms_segment_read(inSegs[iSeg], DRMS_TYPE_FLOAT, &status);
01627         data_ptr[iSeg] = (float *) inArrays[iSeg]->data;
01628         
01629     }
01630     
01631     float *bT0 = data_ptr[0], *bI0 = data_ptr[1], *bA0 = data_ptr[2];
01632     float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];
01633     float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];
01634     
01635     // Add disambig, Feb 12 2014
01636     
01637     DRMS_Segment_t *inSeg;
01638     DRMS_Array_t *inArray_ambig;
01639     
01640     if (amb4err) {              // Mar 4 2014
01641     
01642         inSeg = drms_segment_lookup(inRec, "disambig");
01643         inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);
01644         if (status) return 1;
01645         char *ambig = (char *)inArray_ambig->data;
01646         
01647         int llx, lly;       // lower-left corner
01648         int bmx, bmy;       // bitmap size
01649         
01650         if (fullDisk) {
01651             llx = lly = 0;
01652             bmx = bmy = FOURK;
01653         } else {
01654             llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
01655             lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
01656             bmx = inArray_ambig->axis[0];
01657             bmy = inArray_ambig->axis[1];
01658         }
01659         
01660         int idx, idx_a;
01661         int amb;
01662         
01663         for (int j = 0; j < bmy; j++) {
01664             for (int i = 0; i < bmx; i++) {
01665                 idx_a = j * bmx + i;
01666                 idx = (j + lly) * FOURK + (i + llx);
01667                 // Feb 12 2014, use bit #2 for full disk, lowest bit for patch
01668                 if (fullDisk) { amb = (ambig[idx_a] / 4) % 2; } else { amb = ambig[idx_a] % 2; }
01669                 if (amb) { bA0[idx] += 180.; }
01670             }
01671         }
01672         
01673     }
01674     
01675     // Convert errors to variances, correlation coefficients to covariances
01676     
01677     for (int i = 0; i < FOURK2; i++) {
01678         
01679         if (fabs(errbI0[i]) > 180.) errbI0[i] = 180.;
01680         if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;
01681         
01682         bT[i] = bT0[i];
01683         bI[i] = bI0[i];     // in deg, coverted in errorprop
01684         bA[i] = bA0[i];
01685         
01686         errbT[i] = errbT0[i] * errbT0[i];
01687         errbI[i] = errbI0[i] * errbI0[i] * RADSINDEG * RADSINDEG;
01688         errbA[i] = errbA0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
01689         
01690         errbTbI[i] = errbTbI0[i] * errbT0[i] * errbI0[i] * RADSINDEG;
01691         errbTbA[i] = errbTbA0[i] * errbT0[i] * errbA0[i] * RADSINDEG;
01692         errbIbA[i] = errbIbA0[i] * errbI0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
01693         
01694     }
01695     
01696     //
01697     
01698     for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
01699     if (amb4err) drms_free_array(inArray_ambig);        // Feb 12; Mar 04 2014
01700     
01701     return 0;
01702     
01703 }
01704 
01705 
01706 /*
01707  * Create Cutout record: top level subroutine
01708  * Do the loops on segments and set the keywords here
01709  * Work is done in writeCutout routine below
01710  *
01711  */
01712 
01713 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
01714                     DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
01715                     DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
01716 {
01717     
01718     int status = 0;
01719     
01720     int iHarpSeg;
01721     int nMharpSegs = ARRLENGTH(MharpSegs), nBharpSegs = ARRLENGTH(BharpSegs);
01722     
01723     // Cutout Mharp
01724     
01725     for (iHarpSeg = 0; iHarpSeg < nMharpSegs; iHarpSeg++) {
01726         if (writeCutout(sharpRec, mharpRec, mharpRec, MharpSegs[iHarpSeg])) {
01727             printf("Mharp cutout fails for %s\n", MharpSegs[iHarpSeg]);
01728             break;
01729         }
01730     }
01731     if (iHarpSeg != nMharpSegs) {
01732         SHOW("Cutout: segment number unmatch\n");
01733         return 1;       // if failed
01734     }
01735     printf("Magnetogram cutout done.\n");
01736     
01737     // Cutout Doppler
01738     
01739     if (writeCutout(sharpRec, dopRec, mharpRec, "Dopplergram")) {
01740         printf("Doppler cutout failed\n");
01741         return 1;
01742     }
01743     printf("Dopplergram cutout done.\n");
01744     
01745     // Cutout Continuum
01746     
01747     if (writeCutout(sharpRec, contRec, mharpRec, "continuum")) {
01748         printf("Continuum cutout failed\n");
01749         return 1;
01750     }
01751     printf("Intensitygram cutout done.\n");
01752     
01753     // Coutout Bharp
01754     
01755     for (iHarpSeg = 0; iHarpSeg < nBharpSegs; iHarpSeg++) {
01756         if (writeCutout(sharpRec, bharpRec, mharpRec, BharpSegs[iHarpSeg])) {
01757             printf("Bharp cutout fails for %s\n", BharpSegs[iHarpSeg]);
01758             break;
01759         }
01760     }
01761     if (iHarpSeg != nBharpSegs) return 1;       // if failed
01762     printf("Vector B cutout done.\n");
01763     
01764     // Keywords & Links
01765     
01766     drms_copykey(sharpRec, mharpRec, "T_REC");
01767     drms_copykey(sharpRec, mharpRec, "HARPNUM");
01768     
01769     DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
01770     if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);
01771     DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
01772     if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
01773     
01774     setSWIndex(sharpRec, swKeys_ptr);   // Set space weather indices
01775     setKeys(sharpRec, mharpRec, bharpRec, NULL);              // Set all other keywords, NULL specifies cutout
01776     
01777     // Stats
01778     
01779     int nCutSegs = ARRLENGTH(CutSegs);
01780     for (int iSeg = 0; iSeg < nCutSegs; iSeg++) {
01781         DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
01782         DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
01783         set_statistics(outSeg, outArray, 1);
01784         drms_free_array(outArray);
01785     }
01786     
01787     return 0;
01788     
01789 }
01790 
01791 
01792 /*
01793  * Get cutout and write segment
01794  * Change DISAMB_AZI to apply disambiguation to azimuth
01795  *
01796  */
01797 
01798 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName)
01799 {
01800     
01801     int status = 0;
01802     
01803     DRMS_Segment_t *inSeg = NULL, *outSeg = NULL;
01804     DRMS_Array_t *cutoutArray = NULL;
01805     //  DRMS_Type_t arrayType;
01806     
01807     int ll[2], ur[2], nx, ny, nxny;     // lower-left and upper right coords
01808     
01809     /* Info */
01810     
01811     inSeg = drms_segment_lookup(inRec, SegName);
01812     if (!inSeg) return 1;
01813     
01814     nx = (int) drms_getkey_float(harpRec, "CRSIZE1", &status);
01815     ny = (int) drms_getkey_float(harpRec, "CRSIZE2", &status);
01816     nxny = nx * ny;
01817     ll[0] = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1; if (status) return 1;
01818     ll[1] = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1; if (status) return 1;
01819     ur[0] = ll[0] + nx - 1; if (status) return 1;
01820     ur[1] = ll[1] + ny - 1; if (status) return 1;
01821     
01822     if (inSeg->axis[0] == nx && inSeg->axis[1] == ny) {         // for bitmaps, infomaps, etc.
01823         cutoutArray = drms_segment_read(inSeg, DRMS_TYPE_DOUBLE, &status);
01824         if (status) return 1;
01825     } else if (inSeg->axis[0] == FOURK && inSeg->axis[1] == FOURK) {        // for full disk ones
01826         cutoutArray = drms_segment_readslice(inSeg, DRMS_TYPE_DOUBLE, ll, ur, &status);
01827         if (status) return 1;
01828     } else {
01829         return 1;
01830     }
01831     
01832     // Feb 12 2014, fool-proof, for patch, change everything to 0 or 7!!!
01833     // This is a fix for disambiguation before Aug 2013
01834     
01835     if (!strcmp(SegName, "disambig") && !fullDisk) {
01836         double *disamb = (double *) (cutoutArray->data);
01837         for (int i = 0; i < nxny; i++) {
01838             if (((int)disamb[i]) % 2) { disamb[i] = 7; } else { disamb[i] = 0; }
01839         }
01840     }
01841     
01842     /* Adding disambiguation resolution to cutout azimuth? */
01843     
01844 #if DISAMB_AZI
01845     int amb;
01846     if (!strcmp(SegName, "azimuth")) {
01847         DRMS_Segment_t *disambSeg = NULL;
01848         disambSeg = drms_segment_lookup(inRec, "disambig");
01849         if (!disambSeg) {drms_free_array(cutoutArray); return 1;}
01850         DRMS_Array_t *disambArray;
01851         if (fullDisk) { // Jan 2 2014 XS
01852             disambArray = drms_segment_readslice(disambSeg, DRMS_TYPE_CHAR, ll, ur, &status);
01853             if (status) return 1;
01854         } else {
01855             if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
01856                 disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
01857                 if (status) {drms_free_array(cutoutArray); return 1;}
01858             } else {
01859                 drms_free_array(cutoutArray);
01860                 return 1;
01861             }
01862         }
01863         double *azimuth = (double *) cutoutArray->data;
01864         char *disamb = (char *) disambArray->data;
01865         for (int n = 0; n < nxny; n++) {
01866             //          if (disamb[n] % 2) azimuth[n] += 180.;      // Nov 12 2013 Fixed!!!
01867             // Feb 12 2014, use bit #2 for full disk, lowest bit for patch
01868             if (fullDisk) { amb = (disamb[n] / 4) % 2; } else { amb = disamb[n] % 2; }
01869             if (amb) azimuth[n] += 180.;
01870         }
01871         drms_free_array(disambArray);
01872     }
01873 #endif
01874     
01875     /* Write out */
01876     
01877     outSeg = drms_segment_lookup(outRec, SegName);
01878     if (!outSeg) return 1;
01879     //  drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);  // Jan 02 2013
01880     outSeg->axis[0] = cutoutArray->axis[0];
01881     outSeg->axis[1] = cutoutArray->axis[1];
01882     //  cutoutArray->parent_segment = outSeg;
01883     cutoutArray->israw = 0;     // always compressed
01884     cutoutArray->bzero = outSeg->bzero;
01885     cutoutArray->bscale = outSeg->bscale;       // Same as inArray's
01886     status = drms_segment_write(outSeg, cutoutArray, 0);
01887     drms_free_array(cutoutArray);
01888     if (status) return 1;
01889     
01890     return 0;
01891     
01892 }
01893 
01894 
01895 /*
01896  * Compute space weather indices, no error checking for now
01897  * Based on M. Bobra's swharp_vectorB.c
01898  * No error checking for now
01899  *
01900  */
01901 
01902 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo)
01903 {
01904     
01905     int status = 0;
01906     int nx = mInfo->ncol, ny = mInfo->nrow;
01907     int nxny = nx * ny;
01908     int dims[2] = {nx, ny};
01909     
01910     // Get bx, by, bz, mask
01911     
01912     // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
01913     DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap");
01914     DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
01915     int *bitmask = (int *) bitmaskArray->data;      // get the previously made mask array
01916     
01917     //Use conf_disambig map as a threshold on spaceweather quantities
01918     DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
01919     DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
01920     int *mask = (int *) maskArray->data;        // get the previously made mask array
01921     
01922     DRMS_Segment_t *bxSeg = drms_segment_lookup(inRec, BP_SEG_CEA);
01923     DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
01924     float *bx = (float *) bxArray->data;        // bx
01925     
01926     DRMS_Segment_t *bySeg = drms_segment_lookup(inRec, BT_SEG_CEA);
01927     DRMS_Array_t *byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status);
01928     float *by = (float *) byArray->data;        // by
01929     for (int i = 0; i < nxny; i++) by[i] *= -1;
01930     
01931     DRMS_Segment_t *bzSeg = drms_segment_lookup(inRec, BR_SEG_CEA);
01932     DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
01933     float *bz = (float *) bzArray->data;        // bz
01934     
01935        //Use magnetogram map to compute R
01936        DRMS_Segment_t *losSeg = drms_segment_lookup(inRec, "magnetogram");
01937        DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
01938        float *los = (float *) losArray->data;          // los
01939     
01940     DRMS_Segment_t *bz_errSeg = drms_segment_lookup(inRec, BR_ERR_SEG_CEA);
01941     DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
01942     float *bz_err = (float *) bz_errArray->data;        // bz_err
01943     
01944     DRMS_Segment_t *by_errSeg = drms_segment_lookup(inRec, BT_ERR_SEG_CEA);
01945     DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
01946     float *by_err = (float *) by_errArray->data;        // by_err
01947     //for (int i = 0; i < nxny; i++) by_err[i] *= -1;
01948     
01949     DRMS_Segment_t *bx_errSeg = drms_segment_lookup(inRec, BP_ERR_SEG_CEA);
01950     DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
01951     float *bx_err = (float *) bx_errArray->data;        // bx_err
01952     
01953     // Get emphemeris
01954     float  cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
01955     float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
01956     double rsun_ref    = drms_getkey_double(inRec, "RSUN_REF", &status);
01957     double rsun_obs    = drms_getkey_double(inRec, "RSUN_OBS", &status);
01958     float imcrpix1     = drms_getkey_float(inRec, "IMCRPIX1", &status);
01959     float imcrpix2     = drms_getkey_float(inRec, "IMCRPIX2", &status);
01960     float crpix1       = drms_getkey_float(inRec, "CRPIX1", &status);
01961     float crpix2       = drms_getkey_float(inRec, "CRPIX2", &status);
01962     
01963         // convert cdelt1_orig from degrees to arcsec
01964         float cdelt1       = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
01965 
01966     //if (nx1 > floor((nx-1)/scale + 1) )
01967     //  DIE("X-dimension of output array in fsample() is too large.");
01968     //if (ny1 > floor((ny-1)/scale + 1) )
01969     //  DIE("Y-dimension of output array in fsample() is too large.");
01970     
01971     // Temp arrays
01972     float *bh      = (float *) (malloc(nxny * sizeof(float)));
01973     float *bt      = (float *) (malloc(nxny * sizeof(float)));
01974     float *jz      = (float *) (malloc(nxny * sizeof(float)));
01975     float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
01976     float *bpx     = (float *) (malloc(nxny * sizeof(float)));
01977     float *bpy     = (float *) (malloc(nxny * sizeof(float)));
01978     float *bpz     = (float *) (malloc(nxny * sizeof(float)));
01979     float *derx    = (float *) (malloc(nxny * sizeof(float)));
01980     float *dery    = (float *) (malloc(nxny * sizeof(float)));
01981     float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
01982     float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
01983     float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
01984     float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
01985     float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
01986     float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
01987     float *bt_err  = (float *) (malloc(nxny * sizeof(float)));
01988     float *bh_err  = (float *) (malloc(nxny * sizeof(float)));
01989         float *jz_err  = (float *) (malloc(nxny * sizeof(float)));
01990         float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
01991         float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
01992         float *jz_rms_err  = (float *) (malloc(nxny * sizeof(float)));
01993     float *err_term1   = (float *) (calloc(nxny, sizeof(float)));
01994     float *err_term2   = (float *) (calloc(nxny, sizeof(float)));
01995     float *err_termA   = (float *) (calloc(nxny, sizeof(float)));
01996     float *err_termB   = (float *) (calloc(nxny, sizeof(float)));
01997     float *err_termAt  = (float *) (calloc(nxny, sizeof(float)));
01998     float *err_termBt  = (float *) (calloc(nxny, sizeof(float)));
01999     float *err_termAh  = (float *) (calloc(nxny, sizeof(float)));
02000     float *err_termBh  = (float *) (calloc(nxny, sizeof(float)));
02001      
02002         // define some values for the R calculation
02003         int scale = round(2.0/cdelt1);
02004         int nx1 = nx/scale;
02005         int ny1 = ny/scale;
02006         int nxp = nx1+40;
02007         int nyp = ny1+40;
02008         float *rim     = (float *)malloc(nx1*ny1*sizeof(float));
02009         float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));
02010         float *p1n0    = (float *)malloc(nx1*ny1*sizeof(float));
02011         float *p1p     = (float *)malloc(nx1*ny1*sizeof(float));
02012         float *p1n     = (float *)malloc(nx1*ny1*sizeof(float));
02013         float *p1      = (float *)malloc(nx1*ny1*sizeof(float));
02014         float *pmap    = (float *)malloc(nxp*nyp*sizeof(float));
02015         float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));
02016         float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));
02017 
02018         // define some arrays for the lorentz force calculation
02019         float *fx = (float *) (malloc(nxny * sizeof(float)));
02020         float *fy = (float *) (malloc(nxny * sizeof(float)));
02021         float *fz = (float *) (malloc(nxny * sizeof(float)));
02022     
02023 
02024     //spaceweather quantities computed
02025     if (computeAbsFlux(bz_err, bz , dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),  &(swKeys_ptr->mean_vf_err),
02026                            &(swKeys_ptr->count_mask), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
02027         {
02028         swKeys_ptr->absFlux = DRMS_MISSING_FLOAT;       // If fail, fill in NaN
02029         swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
02030                 swKeys_ptr->mean_vf_err = DRMS_MISSING_FLOAT;
02031                 swKeys_ptr->count_mask  = DRMS_MISSING_INT;
02032     }
02033     
02034     for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
02035     greenpot(bpx, bpy, bpz, nx, ny);
02036     
02037     computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
02038     
02039     if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), &(swKeys_ptr->mean_gamma_err),mask, bitmask))
02040     {
02041                 swKeys_ptr->mean_gamma     =  DRMS_MISSING_FLOAT;
02042                 swKeys_ptr->mean_gamma_err =  DRMS_MISSING_FLOAT;
02043         }
02044     
02045     computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);
02046     
02047     if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, 
02048                                     dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err), err_termAt, err_termBt))
02049         {
02050         swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
02051         swKeys_ptr->mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
02052         }
02053     
02054     if (computeBhderivative(bh, bh_err, dims, &(swKeys_ptr->mean_derivative_bh), 
02055                                 &(swKeys_ptr->mean_derivative_bh_err), mask, bitmask, derx_bh, dery_bh, err_termAh, err_termBh))
02056         {
02057         swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
02058                 swKeys_ptr->mean_derivative_bh_err = DRMS_MISSING_FLOAT;
02059     }
02060     
02061     if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err), 
02062                                 mask, bitmask, derx_bz, dery_bz, err_termA, err_termB))
02063         {
02064         swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
02065                 swKeys_ptr->mean_derivative_bz_err = DRMS_MISSING_FLOAT;
02066         }
02067     
02068     computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
02069                   derx, dery, err_term1, err_term2);
02070     
02071     
02072         if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &(swKeys_ptr->mean_jz),
02073                        &(swKeys_ptr->mean_jz_err), &(swKeys_ptr->us_i), &(swKeys_ptr->us_i_err), mask, bitmask, cdelt1,
02074                        rsun_ref, rsun_obs, derx, dery))
02075         {
02076                 swKeys_ptr->mean_jz            = DRMS_MISSING_FLOAT;
02077         swKeys_ptr->us_i               = DRMS_MISSING_FLOAT;
02078                 swKeys_ptr->mean_jz_err        = DRMS_MISSING_FLOAT;
02079                 swKeys_ptr->us_i_err           = DRMS_MISSING_FLOAT;
02080     }
02081     
02082     if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &(swKeys_ptr->mean_alpha), &(swKeys_ptr->mean_alpha_err), 
02083                          mask, bitmask, cdelt1, rsun_ref, rsun_obs))
02084         {
02085         swKeys_ptr->mean_alpha         = DRMS_MISSING_FLOAT;
02086                 swKeys_ptr->mean_alpha_err     = DRMS_MISSING_FLOAT;
02087         }
02088     
02089     if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &(swKeys_ptr->mean_ih), &(swKeys_ptr->mean_ih_err), 
02090                             &(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
02091                             &(swKeys_ptr->total_us_ih_err), &(swKeys_ptr->total_abs_ih_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
02092         {
02093         swKeys_ptr->mean_ih            = DRMS_MISSING_FLOAT;
02094         swKeys_ptr->total_us_ih        = DRMS_MISSING_FLOAT;
02095         swKeys_ptr->total_abs_ih       = DRMS_MISSING_FLOAT;
02096                 swKeys_ptr->mean_ih_err        = DRMS_MISSING_FLOAT;
02097                 swKeys_ptr->total_us_ih_err    = DRMS_MISSING_FLOAT;
02098                 swKeys_ptr->total_abs_ih_err   = DRMS_MISSING_FLOAT;
02099     }
02100     
02101     if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &(swKeys_ptr->totaljz), &(swKeys_ptr->totaljz_err),
02102                      mask, bitmask, cdelt1, rsun_ref, rsun_obs))
02103         {  
02104         swKeys_ptr->totaljz            = DRMS_MISSING_FLOAT;
02105                 swKeys_ptr->totaljz_err        = DRMS_MISSING_FLOAT;
02106     }
02107     
02108     if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims,
02109                   &(swKeys_ptr->meanpot), &(swKeys_ptr->meanpot_err), &(swKeys_ptr->totpot), &(swKeys_ptr->totpot_err),
02110                   mask, bitmask, cdelt1, rsun_ref, rsun_obs))
02111         {
02112         swKeys_ptr->meanpot            = DRMS_MISSING_FLOAT; // If fail, fill in NaN
02113         swKeys_ptr->totpot             = DRMS_MISSING_FLOAT;
02114                 swKeys_ptr->meanpot_err        = DRMS_MISSING_FLOAT;
02115                 swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;
02116     }
02117     
02118     
02119     if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,
02120                   &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),
02121                   mask, bitmask)) 
02122     {
02123         swKeys_ptr->meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN
02124         swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
02125                 swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;
02126     }
02127     
02128     if (computeR(bz_err, los , dims, &(swKeys_ptr->Rparam), cdelt1, rim, p1p0, p1n0,
02129                      p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
02130         {
02131         swKeys_ptr->Rparam = DRMS_MISSING_FLOAT;        // If fail, fill in NaN
02132         }
02133 
02134     
02135     if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &(swKeys_ptr->totfx), &(swKeys_ptr->totfy), &(swKeys_ptr->totfz), &(swKeys_ptr->totbsq),
02136            &(swKeys_ptr->epsx), &(swKeys_ptr->epsy), &(swKeys_ptr->epsz), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
02137         {  
02138         swKeys_ptr->totfx             = DRMS_MISSING_FLOAT;
02139                 swKeys_ptr->totfy             = DRMS_MISSING_FLOAT;
02140         swKeys_ptr->totfz             = DRMS_MISSING_FLOAT;
02141         swKeys_ptr->totbsq            = DRMS_MISSING_FLOAT;
02142                 swKeys_ptr->epsx              = DRMS_MISSING_FLOAT;
02143         swKeys_ptr->epsy              = DRMS_MISSING_FLOAT;
02144                 swKeys_ptr->epsz              = DRMS_MISSING_FLOAT;
02145 
02146     }
02147     
02148     // Clean up the arrays
02149     
02150     drms_free_array(bitmaskArray);      // Dec 18 2012 Xudong
02151     drms_free_array(maskArray);
02152     drms_free_array(bxArray);
02153     drms_free_array(byArray);
02154     drms_free_array(bzArray);
02155         drms_free_array(losArray);              // Mar 7
02156         drms_free_array(bx_errArray);
02157     drms_free_array(by_errArray);
02158     drms_free_array(bz_errArray);
02159     
02160     free(bh); free(bt); free(jz); free(jz_smooth);
02161     free(bpx); free(bpy); free(bpz);
02162     free(derx); free(dery);
02163     free(derx_bt); free(dery_bt);
02164     free(derx_bz); free(dery_bz);
02165     free(derx_bh); free(dery_bh);
02166     free(bt_err); free(bh_err);  free(jz_err);
02167         free(jz_err_squared); free(jz_rms_err);
02168         free(jz_err_squared_smooth);
02169 
02170         // free the arrays that are related to the numerical derivatives
02171         free(err_term2);
02172         free(err_term1);
02173         free(err_termB);
02174         free(err_termA);
02175         free(err_termBt);
02176         free(err_termAt);
02177         free(err_termBh);
02178         free(err_termAh);
02179 
02180         // free the arrays that are related to the r calculation     
02181         free(rim);
02182         free(p1p0);
02183         free(p1n0);
02184         free(p1p);
02185         free(p1n);
02186         free(p1);
02187         free(pmap);
02188         free(p1pad);
02189         free(pmapn);
02190 
02191         // free the arrays that are related to the lorentz calculation
02192         free(fx); free(fy); free(fz);
02193 }
02194 
02195 /*
02196  * Set space weather indices, no error checking for now
02197  *
02198  */
02199 
02200 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
02201 {
02202     drms_setkey_float(outRec, "USFLUX",  swKeys_ptr->mean_vf);
02203     drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
02204     drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
02205     drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
02206     drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
02207     drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
02208     drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
02209     drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
02210     drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih);
02211     drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih);
02212     drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
02213     drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
02214     drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
02215     drms_setkey_float(outRec, "TOTPOT",  swKeys_ptr->totpot);
02216     drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
02217     drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
02218     drms_setkey_float(outRec, "CMASK",   swKeys_ptr->count_mask);
02219     drms_setkey_float(outRec, "ERRBT",   swKeys_ptr->mean_derivative_btotal_err);
02220     drms_setkey_float(outRec, "ERRVF",   swKeys_ptr->mean_vf_err);
02221     drms_setkey_float(outRec, "ERRGAM",  swKeys_ptr->mean_gamma_err);
02222     drms_setkey_float(outRec, "ERRBH",   swKeys_ptr->mean_derivative_bh_err);
02223     drms_setkey_float(outRec, "ERRBZ",   swKeys_ptr->mean_derivative_bz_err);
02224     drms_setkey_float(outRec, "ERRJZ",   swKeys_ptr->mean_jz_err);
02225     drms_setkey_float(outRec, "ERRUSI",  swKeys_ptr->us_i_err);
02226     drms_setkey_float(outRec, "ERRALP",  swKeys_ptr->mean_alpha_err);
02227     drms_setkey_float(outRec, "ERRMIH",  swKeys_ptr->mean_ih_err);
02228     drms_setkey_float(outRec, "ERRTUI",  swKeys_ptr->total_us_ih_err);
02229     drms_setkey_float(outRec, "ERRTAI",  swKeys_ptr->total_abs_ih_err);
02230     drms_setkey_float(outRec, "ERRJHT",  swKeys_ptr->totaljz_err);
02231     drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err);
02232     drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err);
02233     drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err);
02234     drms_setkey_float(outRec, "R_VALUE", swKeys_ptr->Rparam);
02235     drms_setkey_float(outRec, "TOTFX",   swKeys_ptr->totfx);
02236     drms_setkey_float(outRec, "TOTFY",   swKeys_ptr->totfy);
02237     drms_setkey_float(outRec, "TOTFZ",   swKeys_ptr->totfz);
02238     drms_setkey_float(outRec, "TOTBSQ",  swKeys_ptr->totbsq);
02239     drms_setkey_float(outRec, "EPSX",    swKeys_ptr->epsx);
02240     drms_setkey_float(outRec, "EPSY",    swKeys_ptr->epsy);
02241     drms_setkey_float(outRec, "EPSZ",    swKeys_ptr->epsz);
02242 };
02243 
02244 /*
02245  * Set all keywords, no error checking for now
02246  *
02247  */
02248 
02249 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
02250 {
02251     
02252     copy_me_keys(bharpRec, outRec);
02253     copy_patch_keys(mharpRec, outRec);      // Dec 30
02254     copy_geo_keys(mharpRec, outRec);        // Dec 30
02255     copy_ambig_keys(bharpRec, outRec);
02256     
02257     int status = 0;
02258     
02259     // Change a few geometry keywords for CEA & cutout records
02260     if (mInfo != NULL) {        // CEA
02261         
02262         drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
02263         drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
02264         
02265         drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
02266         drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
02267         drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
02268         drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
02269         drms_setkey_string(outRec, "CUNIT1", "degree");
02270         drms_setkey_string(outRec, "CUNIT2", "degree");
02271         
02272         char key[64];
02273         snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
02274         drms_setkey_string(outRec, "CTYPE1", key);
02275         snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
02276         drms_setkey_string(outRec, "CTYPE2", key);
02277         drms_setkey_float(outRec, "CROTA2", 0.0);
02278         
02279         // Jan 2 2014 XS
02280         int nSeg = ARRLENGTH(CEASegs);
02281         for (int iSeg = 0; iSeg < nSeg; iSeg++) {
02282             DRMS_Segment_t *outSeg = NULL;
02283             outSeg = drms_segment_lookup(outRec, CEASegs[iSeg]);
02284             if (!outSeg) continue;
02285             // Set Bunit
02286             char bunit_xxx[20];
02287             sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
02288             //printf("%s, %s\n", bunit_xxx, CEABunits[iSeg]);
02289             drms_setkey_string(outRec, bunit_xxx, CEABunits[iSeg]);
02290         }
02291         
02292     } else {        // Cutout
02293         
02294         float disk_xc, disk_yc;
02295         if (fullDisk) {
02296             disk_xc = drms_getkey_float(bharpRec, "CRPIX1", &status);
02297             disk_yc = drms_getkey_float(bharpRec, "CRPIX2", &status);
02298         } else {
02299             disk_xc = drms_getkey_float(mharpRec, "IMCRPIX1", &status);
02300             disk_yc = drms_getkey_float(mharpRec, "IMCRPIX2", &status);
02301         }
02302         float x_ll = drms_getkey_float(mharpRec, "CRPIX1", &status);
02303         float y_ll = drms_getkey_float(mharpRec, "CRPIX2", &status);
02304         // Defined as disk center's pixel address wrt lower-left of cutout
02305         drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
02306         drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
02307         // Always 0.
02308         drms_setkey_float(outRec, "CRVAL1", 0);
02309         drms_setkey_float(outRec, "CRVAL2", 0);
02310         
02311         // Jan 2 2014 XS
02312         int nSeg = ARRLENGTH(CutSegs);
02313         for (int iSeg = 0; iSeg < nSeg; iSeg++) {
02314             DRMS_Segment_t *outSeg = NULL;
02315             outSeg = drms_segment_lookup(outRec, CutSegs[iSeg]);
02316             if (!outSeg) continue;
02317             // Set Bunit
02318             char bunit_xxx[20];
02319             sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
02320             //printf("%s, %s\n", bunit_xxx, CutBunits[iSeg]);
02321             drms_setkey_string(outRec, bunit_xxx, CutBunits[iSeg]);
02322         }
02323         
02324         
02325     }
02326     
02327     // Mar 19 XS
02328     if (fullDisk) {
02329         drms_setkey_int(outRec, "AMBPATCH", 0);
02330         drms_setkey_int(outRec, "AMBWEAK", 2);
02331     } else {
02332         drms_setkey_int(outRec, "AMBPATCH", 1);
02333     }
02334     
02335     TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
02336     tnow = (double)time(NULL);
02337     tnow += UNIX_epoch;
02338     
02339     val = drms_getkey_time(bharpRec, "DATE", &status);
02340     drms_setkey_time(outRec, "DATE_B", val);
02341     drms_setkey_time(outRec, "DATE", tnow);
02342     
02343     // set cvs commit version into keyword HEADER
02344     char *cvsinfo  = strdup("$Id: sharp.c,v 1.38 2015/03/18 00:28:26 xudong Exp $");
02345     char *cvsinfo2 = sw_functions_version();
02346     char cvsinfoall[2048];
02347     strcat(cvsinfoall,cvsinfo);
02348     strcat(cvsinfoall,"\n");
02349     strcat(cvsinfoall,cvsinfo2);
02350     status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
02351     
02352 };
02353 
02354 
02355 /* ############# Nearest neighbour interpolation ############### */
02356 
02357 float nnb (float *f, int nx, int ny, double x, double y)
02358 {
02359     
02360     if (x <= -0.5 || y <= -0.5 || x > nx - 0.5 || y > ny - 0.5)
02361         return DRMS_MISSING_FLOAT;
02362     int ilow = floor (x);
02363     int jlow = floor (y);
02364     int i = ((x - ilow) > 0.5) ? ilow + 1 : ilow;
02365     int j = ((y - jlow) > 0.5) ? jlow + 1 : jlow;
02366     return f[j * nx + i];
02367     
02368 }
02369 
02370 
02371 /* ################## Wrapper for Jesper's rebin code ################## */
02372 
02373 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss)
02374 {
02375     
02376     struct fresize_struct fresizes;
02377     int nxout, nyout, xoff, yoff;
02378     int nlead = nx;
02379     
02380     nxout = nx / nbin; nyout = ny / nbin;
02381     if (gauss && nbin != 1)
02382         init_fresize_gaussian(&fresizes, (nbin / 2), (nbin / 2 * 2), nbin);     // for nbin=3, sigma=1, half truncate width=2
02383     else
02384         init_fresize_bin(&fresizes, nbin);
02385     xoff = nbin / 2 + nbin / 2;
02386     yoff = nbin / 2 + nbin / 2;
02387     fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT);
02388     
02389 }