version 1.1, 2012/08/23 08:38:33
|
version 1.14, 2013/05/21 17:06:32
|
|
|
* Version: | * Version: |
* v0.0 Jul 02 2012 | * v0.0 Jul 02 2012 |
* v0.1 Jul 23 2012 | * v0.1 Jul 23 2012 |
|
* v0.2 Sep 04 2012 |
|
* v0.3 Dec 18 2012 |
|
* v0.4 Jan 02 2013 |
* | * |
* Notes: | * Notes: |
* v0.0 | * v0.0 |
|
|
* Fixed char I/O thanks to Art | * Fixed char I/O thanks to Art |
* SW indices fixed | * SW indices fixed |
* Added doppler and continuum | * Added doppler and continuum |
|
* Added other keywords: HEADER (populated by cvs build version), DATE_B |
|
* v0.3 |
|
* Fixed memory leakage of 0.15G per rec; denoted with "Dec 18" |
|
* v0.4 |
|
* Took out convert_inplace(). Was causing all the images to be int |
|
|
* | * |
* Example: | * Example: |
* sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \ | * sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \ |
|
|
#include "errorprop.c" | #include "errorprop.c" |
#include "sw_functions.c" | #include "sw_functions.c" |
| |
|
//#include <mkl.h> // Comment out mkl.h, which can only run on solar3 |
#include <mkl_blas.h> | #include <mkl_blas.h> |
#include <mkl_service.h> | #include <mkl_service.h> |
#include <mkl_lapack.h> | #include <mkl_lapack.h> |
|
|
// Space weather keywords | // Space weather keywords |
struct swIndex { | struct swIndex { |
float mean_vf; | float mean_vf; |
|
float count_mask; |
float absFlux; | float absFlux; |
float mean_hf; | float mean_hf; |
float mean_gamma; | float mean_gamma; |
Line 131 struct swIndex { |
|
Line 142 struct swIndex { |
|
float meanshear_angle; | float meanshear_angle; |
float area_w_shear_gt_45h; | float area_w_shear_gt_45h; |
float meanshear_angleh; | float meanshear_angleh; |
|
float mean_derivative_btotal_err; |
|
float mean_vf_err; |
|
float mean_gamma_err; |
|
float mean_derivative_bh_err; |
|
float mean_derivative_bz_err; |
|
float mean_jz_err; |
|
float us_i_err; |
|
float mean_alpha_err; |
|
float mean_ih_err; |
|
float total_us_ih_err; |
|
float total_abs_ih_err; |
|
float totaljz_err; |
|
float meanpot_err; |
|
float totpot_err; |
|
float meanshear_angle_err; |
}; | }; |
| |
// Mapping method | // Mapping method |
Line 147 enum projection { |
|
Line 173 enum projection { |
|
lambert | lambert |
}; | }; |
| |
// Ephemeris |
// Ephemeris information |
struct ephemeris { | struct ephemeris { |
double disk_lonc, disk_latc; | double disk_lonc, disk_latc; |
double disk_xc, disk_yc; | double disk_xc, disk_yc; |
Line 181 int getInputRS_aux(DRMS_RecordSet_t **in |
|
Line 207 int getInputRS_aux(DRMS_RecordSet_t **in |
|
/* Find record from record set with given T_rec */ | /* Find record from record set with given T_rec */ |
int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec); | int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec); |
| |
// =================== |
|
|
|
/* Create CEA record */ | /* Create CEA record */ |
int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, | int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, |
DRMS_Record_t *dopRec, DRMS_Record_t *contRec, | DRMS_Record_t *dopRec, DRMS_Record_t *contRec, |
Line 289 char *CEASegs[] = {"magnetogram", "bitma |
|
Line 313 char *CEASegs[] = {"magnetogram", "bitma |
|
BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA}; | BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA}; |
/* ========================================================================================================== */ | /* ========================================================================================================== */ |
| |
|
|
|
|
char *module_name = "sharp"; | char *module_name = "sharp"; |
char *version_id = "2012 Jul 02"; /* Version number */ |
char *version_id = "2012 Dec 18"; /* Version number */ |
|
int seed; |
| |
ModuleArgs_t module_args[] = | ModuleArgs_t module_args[] = |
{ | { |
Line 302 ModuleArgs_t module_args[] = |
|
Line 325 ModuleArgs_t module_args[] = |
|
{ARG_STRING, "cont", kNotSpecified, "Input Continuum series."}, | {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."}, |
{ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."}, | {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."}, |
{ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."}, | {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."}, |
|
{ARG_INT, "seed", "987654", "Seed for the random number generator."}, |
{ARG_END} | {ARG_END} |
}; | }; |
| |
int DoIt(void) | int DoIt(void) |
{ | { |
|
int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ); |
|
int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ); |
| |
int status = DRMS_SUCCESS; | int status = DRMS_SUCCESS; |
int nrecs, irec; | int nrecs, irec; |
|
|
contQuery = (char *) params_get_str(&cmdparams, "cont"); | contQuery = (char *) params_get_str(&cmdparams, "cont"); |
sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea"); | sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea"); |
sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut"); | sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut"); |
|
sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut"); |
|
|
|
seed = params_get_int(&cmdparams, "seed"); |
| |
/* Get input data, check everything */ | /* Get input data, check everything */ |
| |
|
|
| |
} // irec | } // irec |
| |
|
|
drms_close_records(mharpRS, DRMS_FREE_RECORD); | drms_close_records(mharpRS, DRMS_FREE_RECORD); |
drms_close_records(bharpRS, DRMS_FREE_RECORD); | drms_close_records(bharpRS, DRMS_FREE_RECORD); |
|
drms_close_records(dopRS, DRMS_FREE_RECORD); // Dec 18 2012 |
|
drms_close_records(contRS, DRMS_FREE_RECORD); // Dec 18 2012 |
| |
return 0; | return 0; |
| |
Line 564 int createCeaRecord(DRMS_Record_t *mharp |
|
Line 596 int createCeaRecord(DRMS_Record_t *mharp |
|
| |
// Get ephemeris | // Get ephemeris |
| |
if (getEphemeris(mharpRec, &(mInfo.ephem))) return 1; |
if (getEphemeris(mharpRec, &(mInfo.ephem))) { |
|
SHOW("CEA: get ephemeris error\n"); |
|
return 1; |
|
} |
| |
// Find position | // Find position |
| |
if (findPosition(mharpRec, &mInfo)) return 1; |
if (findPosition(mharpRec, &mInfo)) { |
|
SHOW("CEA: find position error\n"); |
|
return 1; |
|
} |
| |
// Create xi_out, zeta_out array in mInfo: | // Create xi_out, zeta_out array in mInfo: |
// Coordinates to sample in original full disk image | // Coordinates to sample in original full disk image |
Line 584 int createCeaRecord(DRMS_Record_t *mharp |
|
Line 622 int createCeaRecord(DRMS_Record_t *mharp |
|
| |
// Mapping single segment: Mharp, etc. | // Mapping single segment: Mharp, etc. |
| |
if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) return 1; |
if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) { |
if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) return 1; |
SHOW("CEA: mapping magnetogram error\n"); |
|
return 1; |
|
} |
|
if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) { |
|
SHOW("CEA: mapping bitmap error\n"); |
|
return 1; |
|
} |
printf("Magnetogram mapping done.\n"); | printf("Magnetogram mapping done.\n"); |
|
if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) { |
if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) return 1; |
SHOW("CEA: mapping dopplergram error\n"); |
|
return 1; |
|
} |
printf("Dopplergram mapping done.\n"); | printf("Dopplergram mapping done.\n"); |
| |
if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) return 1; |
if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) { |
|
SHOW("CEA: mapping continuum error\n"); |
|
return 1; |
|
} |
printf("Intensitygram mapping done.\n"); | printf("Intensitygram mapping done.\n"); |
| |
//if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) return 1; |
if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) { |
//printf("Conf disambig mapping done.\n"); |
SHOW("CEA: mapping conf_disambig error\n"); |
|
return 1; |
if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "disambig")) return 1; |
} |
printf("disambig mapping done.\n"); |
printf("Conf disambig mapping done.\n"); |
| |
// Mapping vector B | // Mapping vector B |
| |
if (mapVectorB(sharpRec, bharpRec, &mInfo)) return 1; |
if (mapVectorB(sharpRec, bharpRec, &mInfo)) { |
|
SHOW("CEA: mapping vector B error\n"); |
|
return 1; |
|
} |
printf("Vector B mapping done.\n"); | printf("Vector B mapping done.\n"); |
| |
// Mapping vector B errors | // Mapping vector B errors |
| |
if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) return 1; |
if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) { |
|
SHOW("CEA: mapping vector B uncertainty error\n"); |
|
return 1; |
|
} |
printf("Vector B error done.\n"); | printf("Vector B error done.\n"); |
| |
// Keywords & Links | // Keywords & Links |
Line 621 int createCeaRecord(DRMS_Record_t *mharp |
|
Line 676 int createCeaRecord(DRMS_Record_t *mharp |
|
if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec); | if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec); |
| |
setKeys(sharpRec, bharpRec); // Set all other keywords | setKeys(sharpRec, bharpRec); // Set all other keywords |
|
drms_copykey(sharpRec, mharpRec, "QUALITY"); // copied from los records |
| |
// Space weather | // Space weather |
| |
Line 703 int mapScaler(DRMS_Record_t *sharpRec, D |
|
Line 759 int mapScaler(DRMS_Record_t *sharpRec, D |
|
outSeg = drms_segment_lookup(sharpRec, segName); | outSeg = drms_segment_lookup(sharpRec, segName); |
if (!outSeg) return 1; | if (!outSeg) return 1; |
| |
DRMS_Type_t arrayType = outSeg->info->type; |
// DRMS_Type_t arrayType = outSeg->info->type; |
DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status); | DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status); |
if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;} | if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;} |
| |
// convert to needed data type | // convert to needed data type |
| |
drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray); |
// drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray); // Jan 02 2013 |
| |
outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1]; | outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1]; |
outArray->parent_segment = outSeg; |
// outArray->parent_segment = outSeg; |
outArray->israw = 0; // always compressed | outArray->israw = 0; // always compressed |
outArray->bzero = outSeg->bzero; | outArray->bzero = outSeg->bzero; |
outArray->bscale = outSeg->bscale; | outArray->bscale = outSeg->bscale; |
Line 721 int mapScaler(DRMS_Record_t *sharpRec, D |
|
Line 777 int mapScaler(DRMS_Record_t *sharpRec, D |
|
if (status) return 0; | if (status) return 0; |
| |
if (inArray) drms_free_array(inArray); | if (inArray) drms_free_array(inArray); |
|
if ((xsz != FOURK) || (ysz != FOURK)) free(inData); // Dec 18 2012 |
if (outArray) drms_free_array(outArray); | if (outArray) drms_free_array(outArray); |
return 0; | return 0; |
| |
Line 787 int mapVectorB(DRMS_Record_t *sharpRec, |
|
Line 844 int mapVectorB(DRMS_Record_t *sharpRec, |
|
outSeg = drms_segment_lookup(sharpRec, segName[iSeg]); | outSeg = drms_segment_lookup(sharpRec, segName[iSeg]); |
outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status); | outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status); |
outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1]; | outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1]; |
outArray->parent_segment = outSeg; |
// outArray->parent_segment = outSeg; |
outArray->israw = 0; | outArray->israw = 0; |
outArray->bzero = outSeg->bzero; | outArray->bzero = outSeg->bzero; |
outArray->bscale = outSeg->bscale; | outArray->bscale = outSeg->bscale; |
Line 839 int mapVectorBErr(DRMS_Record_t *sharpRe |
|
Line 896 int mapVectorBErr(DRMS_Record_t *sharpRe |
|
outSeg = drms_segment_lookup(sharpRec, segName[iSeg]); | outSeg = drms_segment_lookup(sharpRec, segName[iSeg]); |
outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status); | outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status); |
outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1]; | outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1]; |
outArray->parent_segment = outSeg; |
// outArray->parent_segment = outSeg; |
outArray->israw = 0; | outArray->israw = 0; |
outArray->bzero = outSeg->bzero; | outArray->bzero = outSeg->bzero; |
outArray->bscale = outSeg->bscale; | outArray->bscale = outSeg->bscale; |
Line 884 int findPosition(DRMS_Record_t *inRec, s |
|
Line 941 int findPosition(DRMS_Record_t *inRec, s |
|
// We compute minlon & minlat then by | // We compute minlon & minlat then by |
// LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT | // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT |
| |
float psize = drms_getkey_float(inRec, "SIZE", &status); |
// float psize = drms_getkey_float(inRec, "SIZE", &status); |
if (psize != psize) { |
// if (psize != psize) { |
TIME t0 = drms_getkey_time(inRec, "T_FRST", &status); if (status) return 1; |
|
|
if (minlon != minlon || maxlon != maxlon) { // check lons instead of SIZE |
|
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 |
double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1; | double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1; |
char firstRecQuery[100], t0_str[100]; | char firstRecQuery[100], t0_str[100]; |
sprint_time(t0_str, t0, "TAI", 0); | sprint_time(t0_str, t0, "TAI", 0); |
Line 1080 int performSampling(float *outData, floa |
|
Line 1139 int performSampling(float *outData, floa |
|
// Rebinning, smoothing | // Rebinning, smoothing |
| |
frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian | frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian |
|
free(outData0); // Dec 18 2012 |
| |
// | // |
| |
Line 1457 int createCutRecord(DRMS_Record_t *mharp |
|
Line 1517 int createCutRecord(DRMS_Record_t *mharp |
|
break; | break; |
} | } |
} | } |
if (iHarpSeg != nMharpSegs) return 1; // if failed |
if (iHarpSeg != nMharpSegs) { |
|
SHOW("Cutout: segment number unmatch\n"); |
|
return 1; // if failed |
|
} |
printf("Magnetogram cutout done.\n"); | printf("Magnetogram cutout done.\n"); |
| |
// Cutout Doppler | // Cutout Doppler |
Line 1559 int writeCutout(DRMS_Record_t *outRec, D |
|
Line 1622 int writeCutout(DRMS_Record_t *outRec, D |
|
| |
#if DISAMB_AZI | #if DISAMB_AZI |
if (!strcmp(SegName, "azimuth")) { | if (!strcmp(SegName, "azimuth")) { |
DRMS_Segment_t *disambSeg = drms_segment_lookup(inRec, "disambig"); |
DRMS_Segment_t *disambSeg = NULL; |
|
disambSeg = drms_segment_lookup(inRec, "disambig"); |
if (!disambSeg) {drms_free_array(cutoutArray); return 1;} | if (!disambSeg) {drms_free_array(cutoutArray); return 1;} |
DRMS_Array_t *disambArray; | DRMS_Array_t *disambArray; |
if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) { | if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) { |
disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status); | disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status); |
if (status) {drms_free_array(cutoutArray); return 1;} | if (status) {drms_free_array(cutoutArray); return 1;} |
} else { | } else { |
|
drms_free_array(cutoutArray); |
return 1; | return 1; |
} | } |
double *azimuth = (double *) cutoutArray->data; | double *azimuth = (double *) cutoutArray->data; |
Line 1581 int writeCutout(DRMS_Record_t *outRec, D |
|
Line 1646 int writeCutout(DRMS_Record_t *outRec, D |
|
| |
outSeg = drms_segment_lookup(outRec, SegName); | outSeg = drms_segment_lookup(outRec, SegName); |
if (!outSeg) return 1; | if (!outSeg) return 1; |
drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray); |
// drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray); // Jan 02 2013 |
outSeg->axis[0] = cutoutArray->axis[0]; | outSeg->axis[0] = cutoutArray->axis[0]; |
outSeg->axis[1] = cutoutArray->axis[1]; | outSeg->axis[1] = cutoutArray->axis[1]; |
cutoutArray->parent_segment = outSeg; |
// cutoutArray->parent_segment = outSeg; |
cutoutArray->israw = 0; // always compressed | cutoutArray->israw = 0; // always compressed |
cutoutArray->bzero = outSeg->bzero; | cutoutArray->bzero = outSeg->bzero; |
cutoutArray->bscale = outSeg->bscale; // Same as inArray's | cutoutArray->bscale = outSeg->bscale; // Same as inArray's |
Line 1611 void computeSWIndex(struct swIndex *swKe |
|
Line 1676 void computeSWIndex(struct swIndex *swKe |
|
int nx = mInfo->ncol, ny = mInfo->nrow; | int nx = mInfo->ncol, ny = mInfo->nrow; |
int nxny = nx * ny; | int nxny = nx * ny; |
int dims[2] = {nx, ny}; | int dims[2] = {nx, ny}; |
|
|
// Get bx, by, bz, mask | // Get bx, by, bz, mask |
| |
// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities | // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities |
//DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "bitmap"); |
DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap"); |
//DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status); |
DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status); |
//int *mask = (int *) maskArray->data; // get the previously made mask array |
int *bitmask = (int *) bitmaskArray->data; // get the previously made mask array |
| |
//Use conf_disambig map as a threshold on spaceweather quantities | //Use conf_disambig map as a threshold on spaceweather quantities |
//DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig"); |
DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig"); |
DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "disambig"); |
|
DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status); | DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status); |
int *mask = (int *) maskArray->data; // get the previously made mask array | int *mask = (int *) maskArray->data; // get the previously made mask array |
| |
Line 1637 void computeSWIndex(struct swIndex *swKe |
|
Line 1702 void computeSWIndex(struct swIndex *swKe |
|
DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status); | DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status); |
float *bz = (float *) bzArray->data; // bz | float *bz = (float *) bzArray->data; // bz |
| |
|
DRMS_Segment_t *bz_errSeg = drms_segment_lookup(inRec, BR_ERR_SEG_CEA); |
|
DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status); |
|
float *bz_err = (float *) bz_errArray->data; // bz_err |
|
|
|
DRMS_Segment_t *by_errSeg = drms_segment_lookup(inRec, BT_ERR_SEG_CEA); |
|
DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status); |
|
float *by_err = (float *) by_errArray->data; // by_err |
|
//for (int i = 0; i < nxny; i++) by_err[i] *= -1; |
|
|
|
DRMS_Segment_t *bx_errSeg = drms_segment_lookup(inRec, BP_ERR_SEG_CEA); |
|
DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status); |
|
float *bx_err = (float *) bx_errArray->data; // bx_err |
|
|
// Get emphemeris | // Get emphemeris |
| |
//float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status); | //float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status); |
Line 1649 void computeSWIndex(struct swIndex *swKe |
|
Line 1727 void computeSWIndex(struct swIndex *swKe |
|
float crpix1 = drms_getkey_float(inRec, "CRPIX1", &status); | float crpix1 = drms_getkey_float(inRec, "CRPIX1", &status); |
float crpix2 = drms_getkey_float(inRec, "CRPIX2", &status); | float crpix2 = drms_getkey_float(inRec, "CRPIX2", &status); |
| |
//float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately) |
|
|
|
printf("cdelt1=%f\n",cdelt1); |
|
printf("rsun_ref=%f\n",rsun_ref); |
|
printf("rsun_obs=%f\n",rsun_obs); |
|
printf("dsun_obs=%f\n",dsun_obs); |
|
|
|
// Temp arrays | // Temp arrays |
| |
float *bh = (float *) (malloc(nxny * sizeof(float))); | float *bh = (float *) (malloc(nxny * sizeof(float))); |
float *bt = (float *) (malloc(nxny * sizeof(float))); | float *bt = (float *) (malloc(nxny * sizeof(float))); |
float *jz = (float *) (malloc(nxny * sizeof(float))); | float *jz = (float *) (malloc(nxny * sizeof(float))); |
|
float *jz_smooth = (float *) (malloc(nxny * sizeof(float))); |
float *bpx = (float *) (malloc(nxny * sizeof(float))); | float *bpx = (float *) (malloc(nxny * sizeof(float))); |
float *bpy = (float *) (malloc(nxny * sizeof(float))); | float *bpy = (float *) (malloc(nxny * sizeof(float))); |
float *bpz = (float *) (malloc(nxny * sizeof(float))); | float *bpz = (float *) (malloc(nxny * sizeof(float))); |
Line 1672 void computeSWIndex(struct swIndex *swKe |
|
Line 1744 void computeSWIndex(struct swIndex *swKe |
|
float *dery_bh = (float *) (malloc(nxny * sizeof(float))); | float *dery_bh = (float *) (malloc(nxny * sizeof(float))); |
float *derx_bz = (float *) (malloc(nxny * sizeof(float))); | float *derx_bz = (float *) (malloc(nxny * sizeof(float))); |
float *dery_bz = (float *) (malloc(nxny * sizeof(float))); | float *dery_bz = (float *) (malloc(nxny * sizeof(float))); |
|
float *bt_err = (float *) (malloc(nxny * sizeof(float))); |
|
float *bh_err = (float *) (malloc(nxny * sizeof(float))); |
|
float *jz_err = (float *) (malloc(nxny * sizeof(float))); |
|
float *jz_err_squared = (float *) (malloc(nxny * sizeof(float))); |
|
float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float))); |
|
float *jz_rms_err = (float *) (malloc(nxny * sizeof(float))); |
|
//spaceweather quantities computed |
| |
// Compute |
|
| |
if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf), |
if (computeAbsFlux(bz_err, bz , dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf), &(swKeys_ptr->mean_vf_err), |
mask, cdelt1, rsun_ref, rsun_obs)){ |
&(swKeys_ptr->count_mask), mask, bitmask, cdelt1, rsun_ref, rsun_obs)) |
|
{ |
swKeys_ptr->absFlux = DRMS_MISSING_FLOAT; // If fail, fill in NaN | swKeys_ptr->absFlux = DRMS_MISSING_FLOAT; // If fail, fill in NaN |
swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT; | swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->mean_vf_err = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->count_mask = DRMS_MISSING_INT; |
} | } |
| |
for (int i = 0; i < nxny; i++) bpz[i] = bz[i]; | for (int i = 0; i < nxny; i++) bpz[i] = bz[i]; |
greenpot(bpx, bpy, bpz, nx, ny); | greenpot(bpx, bpy, bpz, nx, ny); |
| |
computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask); |
computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask); |
| |
if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask)) |
if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), &(swKeys_ptr->mean_gamma_err),mask, bitmask)) |
|
{ |
swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT; | swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->mean_gamma_err = DRMS_MISSING_FLOAT; |
|
} |
| |
computeB_total(bx, by, bz, bt, dims, mask); |
computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask); |
| |
if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, derx_bt, dery_bt)) |
if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err))) |
|
{ |
swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT; | swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->mean_derivative_btotal_err = DRMS_MISSING_FLOAT; |
|
} |
| |
if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, derx_bh, dery_bh)) |
if (computeBhderivative(bh, bh_err, dims, &(swKeys_ptr->mean_derivative_bh), &(swKeys_ptr->mean_derivative_bh_err), mask, bitmask, derx_bh, dery_bh)) |
|
{ |
swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT; | swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->mean_derivative_bh_err = DRMS_MISSING_FLOAT; |
|
} |
| |
if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, derx_bz, dery_bz)) |
if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err), mask, bitmask, derx_bz, dery_bz)) |
|
{ |
swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN | swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN |
|
swKeys_ptr->mean_derivative_bz_err = DRMS_MISSING_FLOAT; |
|
} |
| |
|
computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, |
|
derx, dery); |
| |
| |
if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, |
if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &(swKeys_ptr->mean_jz), |
cdelt1, rsun_ref, rsun_obs, derx, dery)) { |
&(swKeys_ptr->mean_jz_err), &(swKeys_ptr->us_i), &(swKeys_ptr->us_i_err), mask, bitmask, cdelt1, |
|
rsun_ref, rsun_obs, derx, dery)) |
|
{ |
swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT; | swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT; |
swKeys_ptr->us_i = DRMS_MISSING_FLOAT; | swKeys_ptr->us_i = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->mean_jz_err = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->us_i_err = DRMS_MISSING_FLOAT; |
} | } |
| |
printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz); |
if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &(swKeys_ptr->mean_alpha), &(swKeys_ptr->mean_alpha_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs)) |
|
{ |
if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, cdelt1, rsun_ref, rsun_obs)) |
|
swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT; | swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->mean_alpha_err = DRMS_MISSING_FLOAT; |
|
} |
| |
if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih), |
if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &(swKeys_ptr->mean_ih), &(swKeys_ptr->mean_ih_err), &(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih), |
&(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih), |
&(swKeys_ptr->total_us_ih_err), &(swKeys_ptr->total_abs_ih_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs)) |
mask, cdelt1, rsun_ref, rsun_obs)) { |
{ |
swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT; | swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT; |
swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT; | swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT; |
swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT; | swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->mean_ih_err = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->total_us_ih_err = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->total_abs_ih_err = DRMS_MISSING_FLOAT; |
} | } |
| |
if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz), |
if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &(swKeys_ptr->totaljz), &(swKeys_ptr->totaljz_err), |
mask, cdelt1, rsun_ref, rsun_obs)) |
mask, bitmask, cdelt1, rsun_ref, rsun_obs)) |
|
{ |
swKeys_ptr->totaljz = DRMS_MISSING_FLOAT; | swKeys_ptr->totaljz = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->totaljz_err = DRMS_MISSING_FLOAT; |
|
} |
| |
|
if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims, |
if (computeFreeEnergy(bx, by, bpx, bpy, dims, |
&(swKeys_ptr->meanpot), &(swKeys_ptr->meanpot_err), &(swKeys_ptr->totpot), &(swKeys_ptr->totpot_err), |
&(swKeys_ptr->meanpot), &(swKeys_ptr->totpot), |
mask, bitmask, cdelt1, rsun_ref, rsun_obs)) |
mask, cdelt1, rsun_ref, rsun_obs)) { |
{ |
swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN | swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN |
swKeys_ptr->totpot = DRMS_MISSING_FLOAT; | swKeys_ptr->totpot = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->meanpot_err = DRMS_MISSING_FLOAT; |
|
swKeys_ptr->totpot_err = DRMS_MISSING_FLOAT; |
} | } |
| |
if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims, |
if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims, |
&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45), |
&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45), |
&(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h), |
mask, bitmask)) { |
mask)) { |
|
swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN | swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN |
swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT; | swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT; |
swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN |
swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT; |
swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT; |
|
} | } |
| |
// Clean up |
// Clean up the arrays |
| |
|
drms_free_array(bitmaskArray); // Dec 18 2012 Xudong |
drms_free_array(maskArray); | drms_free_array(maskArray); |
drms_free_array(bxArray); | drms_free_array(bxArray); |
drms_free_array(byArray); | drms_free_array(byArray); |
drms_free_array(bzArray); | drms_free_array(bzArray); |
| |
free(bh); free(bt); free(jz); |
free(bh); free(bt); free(jz); free(jz_smooth); |
free(bpx); free(bpy); free(bpz); | free(bpx); free(bpy); free(bpz); |
free(derx); free(dery); | free(derx); free(dery); |
free(derx_bt); free(dery_bt); | free(derx_bt); free(dery_bt); |
free(derx_bz); free(dery_bz); | free(derx_bz); free(dery_bz); |
free(derx_bh); free(dery_bh); | free(derx_bh); free(dery_bh); |
|
free(bt_err); free(bh_err); free(jz_err); |
|
free(jz_err_squared); free(jz_rms_err); |
|
free(jz_err_squared_smooth); |
} | } |
| |
|
|
/* | /* |
* Set space weather indices, no error checking for now | * Set space weather indices, no error checking for now |
* | * |
Line 1783 void setSWIndex(DRMS_Record_t *outRec, s |
|
Line 1891 void setSWIndex(DRMS_Record_t *outRec, s |
|
drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot); | drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot); |
drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle); | drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle); |
drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45); | drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45); |
|
drms_setkey_float(outRec, "CMASK", swKeys_ptr->count_mask); |
|
drms_setkey_float(outRec, "ERRBT", swKeys_ptr->mean_derivative_btotal_err); |
|
drms_setkey_float(outRec, "ERRVF", swKeys_ptr->mean_vf_err); |
|
drms_setkey_float(outRec, "ERRGAM", swKeys_ptr->mean_gamma_err); |
|
drms_setkey_float(outRec, "ERRBH", swKeys_ptr->mean_derivative_bh_err); |
|
drms_setkey_float(outRec, "ERRBZ", swKeys_ptr->mean_derivative_bz_err); |
|
drms_setkey_float(outRec, "ERRJZ", swKeys_ptr->mean_jz_err); |
|
drms_setkey_float(outRec, "ERRUSI", swKeys_ptr->us_i_err); |
|
drms_setkey_float(outRec, "ERRALP", swKeys_ptr->mean_alpha_err); |
|
drms_setkey_float(outRec, "ERRMIH", swKeys_ptr->mean_ih_err); |
|
drms_setkey_float(outRec, "ERRTUI", swKeys_ptr->total_us_ih_err); |
|
drms_setkey_float(outRec, "ERRTAI", swKeys_ptr->total_abs_ih_err); |
|
drms_setkey_float(outRec, "ERRJHT", swKeys_ptr->totaljz_err); |
|
drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err); |
|
drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err); |
|
drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err); |
}; | }; |
| |
|
|
/* | /* |
* Set all keywords, no error checking for now | * Set all keywords, no error checking for now |
* | * |
Line 1797 void setKeys(DRMS_Record_t *outRec, DRMS |
|
Line 1920 void setKeys(DRMS_Record_t *outRec, DRMS |
|
copy_patch_keys(inRec, outRec); | copy_patch_keys(inRec, outRec); |
copy_geo_keys(inRec, outRec); | copy_geo_keys(inRec, outRec); |
copy_ambig_keys(inRec, outRec); | copy_ambig_keys(inRec, outRec); |
|
|
|
char timebuf[1024]; |
|
float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */ |
|
double val; |
|
int status = DRMS_SUCCESS; |
|
|
|
val = drms_getkey_double(inRec, "DATE",&status); |
|
drms_setkey_double(outRec, "DATE_B", val); |
|
sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0); |
|
drms_setkey_string(outRec, "DATE", timebuf); |
|
|
|
// set cvs commit version into keyword HEADER |
|
char *cvsinfo = strdup("$Id$"); |
|
// status = drms_setkey_string(outRec, "HEADER", cvsinfo); |
|
status = drms_setkey_string(outRec, "CODEVER7", cvsinfo); |
|
|
}; | }; |
| |
// | // |