sharp/apps/swharp_vectorB.c

/* 
 *  swharp_vectorB. 
 *
 *  Created by Xudong Sun on 8/22/11.
 *
 *  Modified by Monica Bobra to 
 *              -- include ALL spaceweather keywords in Leka and Barnes (2003, I and II)
 *              -- include potential field calculation from Keiji Hayashi
 *              -- run on los and vector data 15 april 2012 via sharp_functions.c
 *  Bz arrays
 *  Write out abs(B) as data segment and a few keywords as SW index
 *
 *  Use:
 *  First use the bmap module to create the in= and mask= parameters.
 *
 *  then run this module:
 *  swharp_vectorB "in=su_mbobra.bmap_fd10[401][2011.03.09_00:00:00_TAI-2011.03.10_03:00:00_TAI]" /
 *  "mask=su_mbobra.bitmap_fd10[401][2011.03.09_00:00:00_TAI-2011.03 * .10_03:00:00_TAI]" /
 *  "out=su_mbobra.sharp_fd10" "dzvalue=0.001" 
 */

#include <jsoc_main.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
float cdelt1_orig;
float cdelt1;
double rsun_ref;
double dsun_obs;
double rsun_obs;
float imcrpix1;
float imcrpix2;
float crpix1;
float crpix2;

#include "sharp_functions.c"

#define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
#define SHOW(msg) {printf("%s", msg); fflush(stdout);}
#define IN_FILES        3             /* Number of input files */
#define PI (3.141592653589793)        /* Ratio of circumference to diameter of a circle*/
#define MUNAUGHT (0.0000012566370614) /* magnetic constant */

/* declaring all the functions */
int computeMu(float *bz, int *dims, float *mu, float *inverseMu);
int computeAbsFlux(float *bz, int *dims, float *absFlux, float *mean_vf_ptr, int *mask, float *inverseMu);
int computeBh(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_hf_ptr, int *mask);
int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_gamma_ptr, int *mask);
int readFits(char *filename, float **image, int *dims);
int writeFits(char *filename, float *image, int *dims);
int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask);
int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask);
int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask);
int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask);
int computeJz(float *by, float *bx, int *dims, float *jz, float *mean_jz_ptr, float *us_i_ptr, int *mask);
int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask);
int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr, int *mask);
int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr, int *mask);
void greenpot(float *bx, float *by, float *bz, int nnx, int nny); 
int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims, float *meanpotptr, float *totpotptr, int *mask);
int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, float *meanshear_angleptr, float *area_w_shear_gt_45ptr, float *meanshear_anglehptr, float *area_w_shear_gt_45hptr, int *mask);

char *module_name = "swharp_vectorB";   /* Module name */

ModuleArgs_t module_args[] =
{
    {ARG_STRING, "in", NULL, "Input vec mag recordset."},
    {ARG_STRING, "mask", NULL, "Input bitmap recordset."},
    {ARG_STRING, "out", NULL, "Output series."},
    {ARG_FLOAT,  "dzvalue", NULL, "Monopole depth."},
    {ARG_END}
};

int DoIt(void)
{
        
    int status = DRMS_SUCCESS;
        
        char *inQuery, *outQuery;       // input series query string
        char *maskQuery;                // mask series query string
        DRMS_RecordSet_t *inRecSet, *outRecSet, *maskRecSet;
        DRMS_Record_t *inRec, *outRec, *maskRec;
        DRMS_Segment_t *inSegBx, *inSegBy, *inSegBz, *outSeg, *maskSeg;
        DRMS_Array_t *inArrayBx, *inArrayBy, *inArrayBz, *outArray, *maskArray; 
        float *inverseMu, *mu, *bx, *by, *bz, *outData, *bh, *bt, *jz, *bpx, *bpy, *bpz;
        int *mask;
        int dims[2], nxny, nx, ny;      // dimensions;  NAXIS1 = dims[0] which is the number of columns.
        float mean_vf = 0.0; 
        float absFlux = 0.0;
        float mean_hf = 0.0;
        float mean_gamma = 0.0;
        float mean_derivative_btotal = 0.0;
        float mean_derivative_bh = 0.0; 
        float mean_derivative_bz = 0.0;
        float mean_jz = 0.0;
        float us_i = 0.0;
        float mean_alpha = 0.0;
        float mean_ih = 0.0;
        float total_us_ih = 0.0;
        float total_abs_ih = 0.0;
        float totaljz = 0.0;
        float totpot  =0.0;
        float meanpot = 0.0;
        float area_w_shear_gt_45 = 0.0;
        float meanshear_angle = 0.0;
        float area_w_shear_gt_45h = 0.0;
        float meanshear_angleh = 0.0;
        int nrecs, irec, i;

        /* Input */
        
        inQuery = (char *) params_get_str(&cmdparams, "in");
        inRecSet = drms_open_records(drms_env, inQuery, &status);
        if (status || inRecSet->n == 0) DIE("No input data found");
        nrecs = inRecSet->n;
        
        /* Mask */
        
        maskQuery = (char *) params_get_str(&cmdparams, "mask");
        maskRecSet = drms_open_records(drms_env, maskQuery, &status);
        if (status || maskRecSet->n == 0) DIE("No mask data found");
        if (maskRecSet->n != nrecs) DIE("Mask and Input series do not have a 1:1 match"); 

        /* Output */
        
        outQuery = (char *) params_get_str(&cmdparams, "out");
        outRecSet = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
        if (status) DIE("Output recordset not created");
        
        /* Do this for each record */
        
        for (irec = 0; irec < nrecs; irec++)
        {
                /* Input record and data */
                
                inRec = inRecSet->records[irec];
                printf("Input Record #%d of #%d\n", irec+1, nrecs); fflush(stdout);
                
                maskRec = maskRecSet->records[irec];
                printf("Mask Record #%d of #%d\n", irec+1, nrecs); fflush(stdout);
 
                inSegBx = drms_segment_lookupnum(inRec, 0); /* Assume this is Bx equivalent */
                inSegBy = drms_segment_lookupnum(inRec, 1); /* Assume this is By equivalent */
                inSegBz = drms_segment_lookupnum(inRec, 2); /* Assume this is Bz equivalent */

                maskSeg = drms_segment_lookupnum(maskRec, 0); /* This is the bitmap */

                inArrayBx = drms_segment_read(inSegBx, DRMS_TYPE_FLOAT, &status);
                if (status) DIE("No Bx data file found. \n");
                inArrayBy = drms_segment_read(inSegBy, DRMS_TYPE_FLOAT, &status);
                if (status) DIE("No By data file found. \n");
                inArrayBz = drms_segment_read(inSegBz, DRMS_TYPE_FLOAT, &status);
                if (status) DIE("No Bz data file found. \n");
                
                maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
                if (status) DIE("No mask data file found. \n");

                cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
                dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
                rsun_ref    = drms_getkey_double(inRec, "RSUN_REF", &status);
                rsun_obs    = drms_getkey_double(inRec, "RSUN_OBS", &status);
                imcrpix1    = drms_getkey_float(inRec, "IMCRPIX1", &status);
                imcrpix2    = drms_getkey_float(inRec, "IMCRPIX2", &status);
                crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);
                crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);

                
                cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)
                
                printf("cdelt1_orig=%f\n",cdelt1_orig);
                printf("cdelt1=%f\n",cdelt1);
                printf("rsun_obs/rsun_ref=%f\n",rsun_obs/rsun_ref);
                printf("rsun_ref/rsun_obs=%f\n",rsun_ref/rsun_obs);
                printf("test1=%f\n",((1/cdelt1_orig)*(rsun_obs/rsun_ref)*(1000000.)));
                printf("test2=%f\n",((cdelt1_orig)*(PI/180)*(rsun_ref)*(1/1000000.)));

                bx   = (float *)inArrayBx->data;
                by   = (float *)inArrayBy->data;
                bz   = (float *)inArrayBz->data;
                mask = (int *)maskArray->data;

                nx = dims[0] = inArrayBz->axis[0];
                ny = dims[1] = inArrayBz->axis[1];
                nxny = dims[0] * dims[1];
                if (maskArray->axis[0] != nx || maskArray->axis[1] != ny) DIE("Mask and Input series are not of the same size"); 

                /* This is to modify the data for each PROJECTION method */
                int flag;
                char* value1;

                value1 = drms_getkey_string(inRec, "PROJECTION", &status);
                flag = strcmp("LambertCylindrical",value1);
                if (flag == 0)
                {
                      int i, j;   
                      for (j = 0; j < ny; j++) 
                      {
                         for (i = 0; i < nx; i++) 
                         {
                            by[j * nx + i] = - by[j * nx + i];
                         }      
                      }
                }
                
                /* Output data */

                outRec = outRecSet->records[irec];
                drms_setlink_static(outRec, "SRCLINK",  inRec->recnum);

                /*===========================================*/    
                /* Malloc some arrays     */

                inverseMu   = (float *)malloc(nx*ny*sizeof(float));
                mu   = (float *)malloc(nx*ny*sizeof(float));
                bh   = (float *)malloc(nx*ny*sizeof(float));
                bt   = (float *)malloc(nx*ny*sizeof(float));
                jz   = (float *)malloc(nx*ny*sizeof(float));
                bpx  = (float *)malloc(nx*ny*sizeof(float));
                bpy  = (float *)malloc(nx*ny*sizeof(float));
                bpz  = (float *)malloc(nx*ny*sizeof(float));

                /*===========================================*/   
                /* SW Keyword computation */
                
                computeMu(bz, dims, mu, inverseMu);

                if (computeAbsFlux(bz, dims, &absFlux, &mean_vf, mask, inverseMu)) 
                {
                   absFlux = 0.0 / 0.0;         // If fail, fill in NaN
                   mean_vf = 0.0 / 0.0;
                }
                drms_setkey_float(outRec, "USFLUX", mean_vf);

                for (i=0 ;i<nxny; i++){bpz[i]=bz[i];}
                greenpot(bpx, bpy, bpz, nx, ny);
                        
                computeBh(bx, by, bz, bh, dims, &mean_hf, mask);        
               
                if (computeGamma(bx, by, bz, bh, dims, &mean_gamma, mask)) mean_gamma = 0.0 / 0.0;      
                drms_setkey_float(outRec, "MEANGAM", mean_gamma);

                computeB_total(bx, by, bz, bt, dims, mask);     
                
                if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask)) mean_derivative_btotal = 0.0 / 0.0;
                drms_setkey_float(outRec, "MEANGBT", mean_derivative_btotal);

                if (computeBhderivative(bh, dims, &mean_derivative_bh, mask)) mean_derivative_bh = 0.0 / 0.0;
                drms_setkey_float(outRec, "MEANGBH", mean_derivative_bh);
                
                if (computeBzderivative(bz, dims, &mean_derivative_bz, mask)) mean_derivative_bz = 0.0 / 0.0; // If fail, fill in NaN
                drms_setkey_float(outRec, "MEANGBZ", mean_derivative_bz);

                if(computeJz(bx, by, dims, jz, &mean_jz, &us_i, mask))
                {  
                   mean_jz = 0.0 / 0.0;
                   us_i = 0.0 /0.0; 
                } 
                drms_setkey_float(outRec, "MEANJZD", mean_jz);
                drms_setkey_float(outRec, "TOTUSJZ", us_i);

                if (computeAlpha(bz, dims, jz, &mean_alpha, mask)) mean_alpha = 0.0 / 0.0; 
                drms_setkey_float(outRec, "MEANALP", mean_alpha);

                if (computeHelicity(bz, dims, jz, &mean_ih, &total_us_ih, &total_abs_ih, mask)) 
                {  
                   mean_ih     = 0.0/0.0; 
                   total_us_ih = 0.0/0.0;
                   total_abs_ih= 0.0/0.0;
                } 
                drms_setkey_float(outRec, "MEANJZH", mean_ih);
                drms_setkey_float(outRec, "TOTUSJH", total_us_ih);
                drms_setkey_float(outRec, "ABSNJZH", total_abs_ih);

                if (computeSumAbsPerPolarity(bz, jz, dims, &totaljz, mask)) totaljz = 0.0 / 0.0;
                drms_setkey_float(outRec, "SAVNCPP", totaljz);

                if (computeFreeEnergy(bx, by, bpx, bpy, dims, &meanpot, &totpot, mask)) 
                {
                   meanpot = 0.0 / 0.0; // If fail, fill in NaN
                   totpot = 0.0 / 0.0;
                }
                drms_setkey_float(outRec, "MEANPOT", meanpot);
                drms_setkey_float(outRec, "TOTPOT", totpot);

                if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims, &meanshear_angle, &area_w_shear_gt_45, &meanshear_angleh, &area_w_shear_gt_45h, mask)) 
                {
                   meanshear_angle = 0.0 / 0.0; // If fail, fill in NaN
                   area_w_shear_gt_45 = 0.0/0.0;
                   meanshear_angleh = 0.0 / 0.0; // If fail, fill in NaN
                   area_w_shear_gt_45h = 0.0/0.0;
                }
                printf("meanshear_angle=%f, area_w_shear_gt_45=%f, meanshear_angleh=%f, area_w_shear_gt_45h=%f\n",meanshear_angle,area_w_shear_gt_45,meanshear_angleh,area_w_shear_gt_45h);
                drms_setkey_float(outRec, "MEANSHR", meanshear_angle);
                drms_setkey_float(outRec, "SHRGT45", area_w_shear_gt_45);
                //drms_setkey_float(outRec, "MEANSHRH", meanshear_angleh);
                //drms_setkey_float(outRec, "SHRGT45H", area_w_shear_gt_45h);

 
                /*===========================================*/  
                /* Set non-SW keywords */
                
                drms_copykey(outRec, inRec, "T_REC");
                drms_copykey(outRec, inRec, "HARPNUM");
                
                /* Clean up */
                drms_free_array(inArrayBz);
                drms_free_array(maskArray);

    }

        drms_close_records(inRecSet, DRMS_FREE_RECORD);
        drms_close_records(outRecSet, DRMS_INSERT_RECORD);
        
    return 0;

}
Revision:	1.2
Committed:	Fri Jun 22 19:40:27 2012 UTC (11 years, 3 months ago) by mbobra
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Ver_6-4, Ver_9-1, Ver_LATEST, Ver_9-3, Ver_9-41, Ver_9-2, Ver_8-8, Ver_8-2, Ver_8-3, Ver_8-0, Ver_8-1, Ver_8-6, Ver_8-7, Ver_8-4, Ver_8-5, Ver_7-1, Ver_7-0, Ver_9-5, Ver_9-4, Ver_8-10, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Changes since 1.1:	+53 -18 lines
Log Message:	The threshold for the mask has been changed to greater than 2 (old threshold was greater than 1). Modified Files: swharp_vectorB.c
#	User	Rev	Content
1	mbobra	1.2	/*
2			* swharp_vectorB.
3	mbobra	1.1	*
4			* Created by Xudong Sun on 8/22/11.
5	mbobra	1.2	*
6			* Modified by Monica Bobra to
7			* -- include ALL spaceweather keywords in Leka and Barnes (2003, I and II)
8			* -- include potential field calculation from Keiji Hayashi
9			* -- run on los and vector data 15 april 2012 via sharp_functions.c
10	mbobra	1.1	* Bz arrays
11			* Write out abs(B) as data segment and a few keywords as SW index
12			*
13			* Use:
14			* First use the bmap module to create the in= and mask= parameters.
15			*
16			* then run this module:
17	mbobra	1.2	* swharp_vectorB "in=su_mbobra.bmap_fd10[401][2011.03.09_00:00:00_TAI-2011.03.10_03:00:00_TAI]" /
18			* "mask=su_mbobra.bitmap_fd10[401][2011.03.09_00:00:00_TAI-2011.03 * .10_03:00:00_TAI]" /
19			* "out=su_mbobra.sharp_fd10" "dzvalue=0.001"
20	mbobra	1.1	*/
21
22			#include <jsoc_main.h>
23			#include <stdio.h>
24			#include <stdlib.h>
25			#include <math.h>
26	mbobra	1.2	float cdelt1_orig;
27			float cdelt1;
28			double rsun_ref;
29			double dsun_obs;
30			double rsun_obs;
31			float imcrpix1;
32			float imcrpix2;
33			float crpix1;
34			float crpix2;
35
36	mbobra	1.1	#include "sharp_functions.c"
37
38			#define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
39			#define SHOW(msg) {printf("%s", msg); fflush(stdout);}
40	mbobra	1.2	#define IN_FILES 3 /* Number of input files */
41			#define PI (3.141592653589793) /* Ratio of circumference to diameter of a circle*/
42			#define MUNAUGHT (0.0000012566370614) /* magnetic constant */
43	mbobra	1.1
44			/* declaring all the functions */
45	mbobra	1.2	int computeMu(float bz, int dims, float mu, float inverseMu);
46			int computeAbsFlux(float bz, int dims, float absFlux, float mean_vf_ptr, int mask, float inverseMu);
47	mbobra	1.1	int computeBh(float bx, float by, float bz, float bh, int dims, float mean_hf_ptr, int *mask);
48			int computeGamma(float bx, float by, float bz, float bh, int dims, float mean_gamma_ptr, int *mask);
49			int readFits(char filename, float image, int dims);
50			int writeFits(char filename, float image, int *dims);
51			int computeB_total(float bx, float by, float bz, float bt, int dims, int mask);
52			int computeBtotalderivative(float bt, int dims, float mean_derivative_btotal_ptr, int mask);
53			int computeBhderivative(float bh, int dims, float mean_derivative_bh_ptr, int mask);
54			int computeBzderivative(float bz, int dims, float mean_derivative_bz_ptr, int mask);
55			int computeJz(float by, float bx, int dims, float jz, float mean_jz_ptr, float us_i_ptr, int *mask);
56			int computeAlpha(float bz, int dims, float jz, float mean_alpha_ptr, int *mask);
57			int computeHelicity(float bz, int dims, float jz, float mean_ih_ptr, float total_us_ih_ptr, float total_abs_ih_ptr, int *mask);
58			int computeSumAbsPerPolarity(float bz, float jz, int dims, float totaljzptr, int *mask);
59			void greenpot(float bx, float by, float *bz, int nnx, int nny);
60			int computeFreeEnergy(float bx, float by, float bpx, float bpy, int dims, float meanpotptr, float totpotptr, int mask);
61			int computeShearAngle(float bx, float by, float bz, float bpx, float bpy, float bpz, int dims, float meanshear_angleptr, float area_w_shear_gt_45ptr, float meanshear_anglehptr, float area_w_shear_gt_45hptr, int mask);
62
63			char module_name = "swharp_vectorB"; / Module name */
64
65			ModuleArgs_t module_args[] =
66			{
67			{ARG_STRING, "in", NULL, "Input vec mag recordset."},
68			{ARG_STRING, "mask", NULL, "Input bitmap recordset."},
69			{ARG_STRING, "out", NULL, "Output series."},
70			{ARG_FLOAT, "dzvalue", NULL, "Monopole depth."},
71			{ARG_END}
72			};
73
74			int DoIt(void)
75			{
76
77			int status = DRMS_SUCCESS;
78
79	mbobra	1.2	char inQuery, outQuery; // input series query string
80	mbobra	1.1	char *maskQuery; // mask series query string
81			DRMS_RecordSet_t inRecSet, outRecSet, *maskRecSet;
82			DRMS_Record_t inRec, outRec, *maskRec;
83			DRMS_Segment_t inSegBx, inSegBy, inSegBz, outSeg, *maskSeg;
84			DRMS_Array_t inArrayBx, inArrayBy, inArrayBz, outArray, *maskArray;
85	mbobra	1.2	float inverseMu, mu, bx, by, bz, outData, bh, bt, jz, bpx, bpy, bpz;
86	mbobra	1.1	int *mask;
87	mbobra	1.2	int dims[2], nxny, nx, ny; // dimensions; NAXIS1 = dims[0] which is the number of columns.
88	mbobra	1.1	float mean_vf = 0.0;
89			float absFlux = 0.0;
90			float mean_hf = 0.0;
91			float mean_gamma = 0.0;
92			float mean_derivative_btotal = 0.0;
93			float mean_derivative_bh = 0.0;
94			float mean_derivative_bz = 0.0;
95			float mean_jz = 0.0;
96			float us_i = 0.0;
97			float mean_alpha = 0.0;
98			float mean_ih = 0.0;
99			float total_us_ih = 0.0;
100			float total_abs_ih = 0.0;
101			float totaljz = 0.0;
102			float totpot =0.0;
103			float meanpot = 0.0;
104			float area_w_shear_gt_45 = 0.0;
105			float meanshear_angle = 0.0;
106			float area_w_shear_gt_45h = 0.0;
107			float meanshear_angleh = 0.0;
108			int nrecs, irec, i;
109
110			/* Input */
111
112			inQuery = (char *) params_get_str(&cmdparams, "in");
113			inRecSet = drms_open_records(drms_env, inQuery, &status);
114			if (status \|\| inRecSet->n == 0) DIE("No input data found");
115			nrecs = inRecSet->n;
116
117			/* Mask */
118
119			maskQuery = (char *) params_get_str(&cmdparams, "mask");
120			maskRecSet = drms_open_records(drms_env, maskQuery, &status);
121			if (status \|\| maskRecSet->n == 0) DIE("No mask data found");
122			if (maskRecSet->n != nrecs) DIE("Mask and Input series do not have a 1:1 match");
123
124			/* Output */
125
126			outQuery = (char *) params_get_str(&cmdparams, "out");
127			outRecSet = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
128			if (status) DIE("Output recordset not created");
129
130			/* Do this for each record */
131
132			for (irec = 0; irec < nrecs; irec++)
133			{
134			/* Input record and data */
135
136			inRec = inRecSet->records[irec];
137			printf("Input Record #%d of #%d\n", irec+1, nrecs); fflush(stdout);
138
139			maskRec = maskRecSet->records[irec];
140			printf("Mask Record #%d of #%d\n", irec+1, nrecs); fflush(stdout);
141
142			inSegBx = drms_segment_lookupnum(inRec, 0); /* Assume this is Bx equivalent */
143			inSegBy = drms_segment_lookupnum(inRec, 1); /* Assume this is By equivalent */
144			inSegBz = drms_segment_lookupnum(inRec, 2); /* Assume this is Bz equivalent */
145
146			maskSeg = drms_segment_lookupnum(maskRec, 0); /* This is the bitmap */
147
148			inArrayBx = drms_segment_read(inSegBx, DRMS_TYPE_FLOAT, &status);
149			if (status) DIE("No Bx data file found. \n");
150			inArrayBy = drms_segment_read(inSegBy, DRMS_TYPE_FLOAT, &status);
151			if (status) DIE("No By data file found. \n");
152			inArrayBz = drms_segment_read(inSegBz, DRMS_TYPE_FLOAT, &status);
153			if (status) DIE("No Bz data file found. \n");
154
155			maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
156			if (status) DIE("No mask data file found. \n");
157
158	mbobra	1.2	cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status);
159			dsun_obs = drms_getkey_float(inRec, "DSUN_OBS", &status);
160			rsun_ref = drms_getkey_double(inRec, "RSUN_REF", &status);
161			rsun_obs = drms_getkey_double(inRec, "RSUN_OBS", &status);
162			imcrpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
163			imcrpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
164			crpix1 = drms_getkey_float(inRec, "CRPIX1", &status);
165			crpix2 = drms_getkey_float(inRec, "CRPIX2", &status);
166
167
168			cdelt1=( (rsun_refcdelt1_origPI/180.) / (dsun_obs) )(180./PI)(3600.); //convert cdelt1 from degrees to arcsec (approximately)
169
170			printf("cdelt1_orig=%f\n",cdelt1_orig);
171			printf("cdelt1=%f\n",cdelt1);
172			printf("rsun_obs/rsun_ref=%f\n",rsun_obs/rsun_ref);
173			printf("rsun_ref/rsun_obs=%f\n",rsun_ref/rsun_obs);
174			printf("test1=%f\n",((1/cdelt1_orig)(rsun_obs/rsun_ref)(1000000.)));
175			printf("test2=%f\n",((cdelt1_orig)(PI/180)(rsun_ref)*(1/1000000.)));
176
177	mbobra	1.1	bx = (float *)inArrayBx->data;
178			by = (float *)inArrayBy->data;
179			bz = (float *)inArrayBz->data;
180			mask = (int *)maskArray->data;
181
182			nx = dims[0] = inArrayBz->axis[0];
183			ny = dims[1] = inArrayBz->axis[1];
184			nxny = dims[0] * dims[1];
185			if (maskArray->axis[0] != nx \|\| maskArray->axis[1] != ny) DIE("Mask and Input series are not of the same size");
186
187			/* This is to modify the data for each PROJECTION method */
188			int flag;
189			char* value1;
190
191			value1 = drms_getkey_string(inRec, "PROJECTION", &status);
192			flag = strcmp("LambertCylindrical",value1);
193			if (flag == 0)
194			{
195			int i, j;
196			for (j = 0; j < ny; j++)
197			{
198			for (i = 0; i < nx; i++)
199			{
200			by[j * nx + i] = - by[j * nx + i];
201			}
202			}
203			}
204
205			/* Output data */
206
207			outRec = outRecSet->records[irec];
208			drms_setlink_static(outRec, "SRCLINK", inRec->recnum);
209
210			/===========================================/
211			/* Malloc some arrays */
212
213	mbobra	1.2	inverseMu = (float )malloc(nxny*sizeof(float));
214			mu = (float )malloc(nxny*sizeof(float));
215	mbobra	1.1	bh = (float )malloc(nxny*sizeof(float));
216			bt = (float )malloc(nxny*sizeof(float));
217			jz = (float )malloc(nxny*sizeof(float));
218			bpx = (float )malloc(nxny*sizeof(float));
219			bpy = (float )malloc(nxny*sizeof(float));
220			bpz = (float )malloc(nxny*sizeof(float));
221
222			/===========================================/
223			/* SW Keyword computation */
224
225	mbobra	1.2	computeMu(bz, dims, mu, inverseMu);
226
227			if (computeAbsFlux(bz, dims, &absFlux, &mean_vf, mask, inverseMu))
228	mbobra	1.1	{
229			absFlux = 0.0 / 0.0; // If fail, fill in NaN
230			mean_vf = 0.0 / 0.0;
231			}
232			drms_setkey_float(outRec, "USFLUX", mean_vf);
233
234			for (i=0 ;i<nxny; i++){bpz[i]=bz[i];}
235			greenpot(bpx, bpy, bpz, nx, ny);
236
237			computeBh(bx, by, bz, bh, dims, &mean_hf, mask);
238
239			if (computeGamma(bx, by, bz, bh, dims, &mean_gamma, mask)) mean_gamma = 0.0 / 0.0;
240			drms_setkey_float(outRec, "MEANGAM", mean_gamma);
241
242			computeB_total(bx, by, bz, bt, dims, mask);
243
244			if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask)) mean_derivative_btotal = 0.0 / 0.0;
245			drms_setkey_float(outRec, "MEANGBT", mean_derivative_btotal);
246
247			if (computeBhderivative(bh, dims, &mean_derivative_bh, mask)) mean_derivative_bh = 0.0 / 0.0;
248			drms_setkey_float(outRec, "MEANGBH", mean_derivative_bh);
249
250			if (computeBzderivative(bz, dims, &mean_derivative_bz, mask)) mean_derivative_bz = 0.0 / 0.0; // If fail, fill in NaN
251			drms_setkey_float(outRec, "MEANGBZ", mean_derivative_bz);
252
253			if(computeJz(bx, by, dims, jz, &mean_jz, &us_i, mask))
254			{
255			mean_jz = 0.0 / 0.0;
256			us_i = 0.0 /0.0;
257			}
258			drms_setkey_float(outRec, "MEANJZD", mean_jz);
259			drms_setkey_float(outRec, "TOTUSJZ", us_i);
260
261			if (computeAlpha(bz, dims, jz, &mean_alpha, mask)) mean_alpha = 0.0 / 0.0;
262			drms_setkey_float(outRec, "MEANALP", mean_alpha);
263
264			if (computeHelicity(bz, dims, jz, &mean_ih, &total_us_ih, &total_abs_ih, mask))
265			{
266	mbobra	1.2	mean_ih = 0.0/0.0;
267	mbobra	1.1	total_us_ih = 0.0/0.0;
268			total_abs_ih= 0.0/0.0;
269			}
270			drms_setkey_float(outRec, "MEANJZH", mean_ih);
271			drms_setkey_float(outRec, "TOTUSJH", total_us_ih);
272			drms_setkey_float(outRec, "ABSNJZH", total_abs_ih);
273
274			if (computeSumAbsPerPolarity(bz, jz, dims, &totaljz, mask)) totaljz = 0.0 / 0.0;
275			drms_setkey_float(outRec, "SAVNCPP", totaljz);
276
277			if (computeFreeEnergy(bx, by, bpx, bpy, dims, &meanpot, &totpot, mask))
278			{
279			meanpot = 0.0 / 0.0; // If fail, fill in NaN
280			totpot = 0.0 / 0.0;
281			}
282			drms_setkey_float(outRec, "MEANPOT", meanpot);
283			drms_setkey_float(outRec, "TOTPOT", totpot);
284
285			if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims, &meanshear_angle, &area_w_shear_gt_45, &meanshear_angleh, &area_w_shear_gt_45h, mask))
286			{
287			meanshear_angle = 0.0 / 0.0; // If fail, fill in NaN
288			area_w_shear_gt_45 = 0.0/0.0;
289			meanshear_angleh = 0.0 / 0.0; // If fail, fill in NaN
290			area_w_shear_gt_45h = 0.0/0.0;
291			}
292			printf("meanshear_angle=%f, area_w_shear_gt_45=%f, meanshear_angleh=%f, area_w_shear_gt_45h=%f\n",meanshear_angle,area_w_shear_gt_45,meanshear_angleh,area_w_shear_gt_45h);
293			drms_setkey_float(outRec, "MEANSHR", meanshear_angle);
294			drms_setkey_float(outRec, "SHRGT45", area_w_shear_gt_45);
295			//drms_setkey_float(outRec, "MEANSHRH", meanshear_angleh);
296			//drms_setkey_float(outRec, "SHRGT45H", area_w_shear_gt_45h);
297
298
299			/===========================================/
300			/* Set non-SW keywords */
301
302			drms_copykey(outRec, inRec, "T_REC");
303			drms_copykey(outRec, inRec, "HARPNUM");
304
305			/* Clean up */
306			drms_free_array(inArrayBz);
307			drms_free_array(maskArray);
308
309			}
310
311			drms_close_records(inRecSet, DRMS_FREE_RECORD);
312			drms_close_records(outRecSet, DRMS_INSERT_RECORD);
313
314			return 0;
315
316			}