sharp/apps/swharp_losB.c

/*
 *  swharp_losB.
 *
 *  Created by Xudong Sun on 8/22/11.
 *  Modified to -- include ALL spaceweather keywords by Monica Bobra 25 Aug 2011
 *              -- include potential field calculation 
 *              -- run only on los data 5 March 2011
 *  Bz arrays
 *  Write out abs(B) as data segment and a few keywords as SW index
 *
 *  Use:
 *  First use the mmap module to create the in= and mask= parameters:
 *  /home/xudong/cvs/JSOC/bin/linux_x86_64/mmap "in=hmi.M_720s[2012.02.01_03:48:00_TAI]" "harp=su_turmon.Mharpv9_720s[1][2012.02.01_03:48:00_TAI]" "out=su_mbobra.test_mmap" "map=Postel" -a -e -v
 *  /home/xudong/cvs/JSOC/bin/linux_x86_64/mmap "harp=su_turmon.Mharpv9_720s[1][2012.02.01_03:48:00_TAI]" -z  "out=su_mbobra.test_mmap_bitmap" "map=Postel" "segment=bitmap" -a -v
 *
 *  then run this module:
 *  swharp_losB  "in=su_mbobra.test_mmap[][2011.10.06_23:59:60_TAI/1d]" /
 *  "mask=su_mbobra.test_mmap_bitmap[][2011.10.06_23:59:60_TAI/1d]" "out=su_mbobra.swharp_test_v1" "dzvalue=0.001"
 */


#include <jsoc_main.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#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)
#define CMPERPIX (0.504277*696000000.0*100.)/(943.) /* cm/pixel = (CDELT1*RSUN_REF*100./RSUN_OBS) */

/* CMSQUARED = CMX*CMY
   CMX = CDELT1*RSUN_REF*(#of pixels in the x-direction)*100/RSUN_OBS
   CMY = CDELT2*RSUN_REF*(#of pixels in the y-direction)*100/RSUN_OBS */


/* declaring all the functions */

int computeAbsFlux(float *bz, int *dims, float *absFlux, float *mean_vf_ptr, int *mask);
int computeBh(float *bpx, float *bpy, float *bz, float *bh, int *dims, float *mean_hf_ptr, int *mask);
int computeGamma(float *bpx, float *bpy, 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 *bpx, float *bpy, 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);
void greenpot(float *bx, float *by, float *bz, int nnx, int nny); 

char *module_name = "swharp_losB";      /* 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 *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);
                
                inSegBz = drms_segment_lookupnum(inRec, 0); /* Assume this is Bz equivalent */

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

                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");

                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     */

                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 */
                
                if (computeAbsFlux(bz, dims, &absFlux, &mean_vf, mask)) 
                {
                   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(bpx, bpy, bz, bh, dims, &mean_hf, mask);      
               
                if (computeGamma(bpx, bpy, bz, bh, dims, &mean_gamma, mask)) mean_gamma = 0.0 / 0.0;    
                drms_setkey_float(outRec, "MEANGAM", mean_gamma);

                computeB_total(bpx, bpy, 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);

                /*===========================================*/  
                /* 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.1
Committed:	Wed Apr 4 22:55:46 2012 UTC (11 years, 5 months ago) by mbobra
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Ver_6-2, Ver_6-3, 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
Log Message:	module to calculate the LOS B SHARP data
#	Content
1	/*
2	* swharp_losB.
3	*
4	* Created by Xudong Sun on 8/22/11.
5	* Modified to -- include ALL spaceweather keywords by Monica Bobra 25 Aug 2011
6	* -- include potential field calculation
7	* -- run only on los data 5 March 2011
8	* Bz arrays
9	* Write out abs(B) as data segment and a few keywords as SW index
10	*
11	* Use:
12	* First use the mmap module to create the in= and mask= parameters:
13	* /home/xudong/cvs/JSOC/bin/linux_x86_64/mmap "in=hmi.M_720s[2012.02.01_03:48:00_TAI]" "harp=su_turmon.Mharpv9_720s[1][2012.02.01_03:48:00_TAI]" "out=su_mbobra.test_mmap" "map=Postel" -a -e -v
14	* /home/xudong/cvs/JSOC/bin/linux_x86_64/mmap "harp=su_turmon.Mharpv9_720s[1][2012.02.01_03:48:00_TAI]" -z "out=su_mbobra.test_mmap_bitmap" "map=Postel" "segment=bitmap" -a -v
15	*
16	* then run this module:
17	* swharp_losB "in=su_mbobra.test_mmap[][2011.10.06_23:59:60_TAI/1d]" /
18	* "mask=su_mbobra.test_mmap_bitmap[][2011.10.06_23:59:60_TAI/1d]" "out=su_mbobra.swharp_test_v1" "dzvalue=0.001"
19	*/
20
21
22	#include <jsoc_main.h>
23	#include <stdio.h>
24	#include <stdlib.h>
25	#include <math.h>
26	#include "sharp_functions.c"
27
28	#define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
29	#define SHOW(msg) {printf("%s", msg); fflush(stdout);}
30	#define IN_FILES 3 /* Number of input files */
31	#define PI (3.141592653589793)
32	#define CMPERPIX (0.504277696000000.0100.)/(943.) /* cm/pixel = (CDELT1RSUN_REF100./RSUN_OBS) */
33
34	/* CMSQUARED = CMX*CMY
35	CMX = CDELT1RSUN_REF(#of pixels in the x-direction)*100/RSUN_OBS
36	CMY = CDELT2RSUN_REF(#of pixels in the y-direction)100/RSUN_OBS /
37
38
39	/* declaring all the functions */
40
41	int computeAbsFlux(float bz, int dims, float absFlux, float mean_vf_ptr, int *mask);
42	int computeBh(float bpx, float bpy, float bz, float bh, int dims, float mean_hf_ptr, int *mask);
43	int computeGamma(float bpx, float bpy, float bz, float bh, int dims, float mean_gamma_ptr, int *mask);
44	int readFits(char filename, float image, int dims);
45	int writeFits(char filename, float image, int *dims);
46	int computeB_total(float bpx, float bpy, float bz, float bt, int dims, int mask);
47	int computeBtotalderivative(float bt, int dims, float mean_derivative_btotal_ptr, int mask);
48	int computeBhderivative(float bh, int dims, float mean_derivative_bh_ptr, int mask);
49	int computeBzderivative(float bz, int dims, float mean_derivative_bz_ptr, int mask);
50	void greenpot(float bx, float by, float *bz, int nnx, int nny);
51
52	char module_name = "swharp_losB"; / Module name */
53
54	ModuleArgs_t module_args[] =
55	{
56	{ARG_STRING, "in", NULL, "Input vec mag recordset."},
57	{ARG_STRING, "mask", NULL, "Input bitmap recordset."},
58	{ARG_STRING, "out", NULL, "Output series."},
59	{ARG_FLOAT, "dzvalue", NULL, "Monopole depth."},
60	{ARG_END}
61	};
62
63
64	int DoIt(void)
65	{
66
67	int status = DRMS_SUCCESS;
68
69	char inQuery, outQuery; // input series query string
70	char *maskQuery; // mask series query string
71	DRMS_RecordSet_t inRecSet, outRecSet, *maskRecSet;
72	DRMS_Record_t inRec, outRec, *maskRec;
73	DRMS_Segment_t inSegBx, inSegBy, inSegBz, outSeg, *maskSeg;
74	DRMS_Array_t inArrayBx, inArrayBy, inArrayBz, outArray, *maskArray;
75	float bx, by, bz, outData, bh, bt, jz, bpx, bpy, bpz;
76	int *mask;
77	int dims[2], nxny, nx, ny; // dimensions; NAXIS1 = dims[0] which is the number of columns.
78	float mean_vf = 0.0;
79	float absFlux = 0.0;
80	float mean_hf = 0.0;
81	float mean_gamma = 0.0;
82	float mean_derivative_btotal = 0.0;
83	float mean_derivative_bh = 0.0;
84	float mean_derivative_bz = 0.0;
85	float mean_jz = 0.0;
86	float us_i = 0.0;
87	float mean_alpha = 0.0;
88	float mean_ih = 0.0;
89	float total_us_ih = 0.0;
90	float total_abs_ih = 0.0;
91	float totaljz = 0.0;
92	float totpot =0.0;
93	float meanpot = 0.0;
94	float area_w_shear_gt_45 = 0.0;
95	float meanshear_angle = 0.0;
96	float area_w_shear_gt_45h = 0.0;
97	float meanshear_angleh = 0.0;
98	int nrecs, irec, i;
99
100	/* Input */
101
102	inQuery = (char *) params_get_str(&cmdparams, "in");
103	inRecSet = drms_open_records(drms_env, inQuery, &status);
104	if (status \|\| inRecSet->n == 0) DIE("No input data found");
105	nrecs = inRecSet->n;
106
107	/* Mask */
108
109	maskQuery = (char *) params_get_str(&cmdparams, "mask");
110	maskRecSet = drms_open_records(drms_env, maskQuery, &status);
111	if (status \|\| maskRecSet->n == 0) DIE("No mask data found");
112	if (maskRecSet->n != nrecs) DIE("Mask and Input series do not have a 1:1 match");
113
114	/* Output */
115
116	outQuery = (char *) params_get_str(&cmdparams, "out");
117	outRecSet = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
118	if (status) DIE("Output recordset not created");
119
120	/* Do this for each record */
121
122	for (irec = 0; irec < nrecs; irec++)
123	{
124
125	/* Input record and data */
126
127	inRec = inRecSet->records[irec];
128	printf("Input Record #%d of #%d\n", irec+1, nrecs); fflush(stdout);
129
130	maskRec = maskRecSet->records[irec];
131	printf("Mask Record #%d of #%d\n", irec+1, nrecs); fflush(stdout);
132
133	inSegBz = drms_segment_lookupnum(inRec, 0); /* Assume this is Bz equivalent */
134
135	maskSeg = drms_segment_lookupnum(maskRec, 0); /* This is the bitmap */
136
137	inArrayBz = drms_segment_read(inSegBz, DRMS_TYPE_FLOAT, &status);
138	if (status) DIE("No Bz data file found. \n");
139
140	maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
141	if (status) DIE("No mask data file found. \n");
142
143	bz = (float *)inArrayBz->data;
144	mask = (int *)maskArray->data;
145
146	nx = dims[0] = inArrayBz->axis[0];
147	ny = dims[1] = inArrayBz->axis[1];
148	nxny = dims[0] * dims[1];
149	if (maskArray->axis[0] != nx \|\| maskArray->axis[1] != ny) DIE("Mask and Input series are not of the same size");
150
151	/* This is to modify the data for each PROJECTION method */
152	int flag;
153	char* value1;
154	value1 = drms_getkey_string(inRec, "PROJECTION", &status);
155	flag = strcmp("LambertCylindrical",value1);
156	if (flag == 0)
157	{
158	int i, j;
159	for (j = 0; j < ny; j++)
160	{
161	for (i = 0; i < nx; i++)
162	{
163	by[j * nx + i] = - by[j * nx + i];
164	}
165	}
166	}
167
168	/* Output data */
169
170	outRec = outRecSet->records[irec];
171	drms_setlink_static(outRec, "SRCLINK", inRec->recnum);
172
173	/===========================================/
174	/* Malloc some arrays */
175
176	bh = (float )malloc(nxny*sizeof(float));
177	bt = (float )malloc(nxny*sizeof(float));
178	jz = (float )malloc(nxny*sizeof(float));
179	bpx = (float )malloc(nxny*sizeof(float));
180	bpy = (float )malloc(nxny*sizeof(float));
181	bpz = (float )malloc(nxny*sizeof(float));
182
183	/===========================================/
184	/* SW Keyword computation */
185
186	if (computeAbsFlux(bz, dims, &absFlux, &mean_vf, mask))
187	{
188	absFlux = 0.0 / 0.0; // If fail, fill in NaN
189	mean_vf = 0.0 / 0.0;
190	}
191	drms_setkey_float(outRec, "USFLUX", mean_vf);
192
193	for (i=0 ;i<nxny; i++){bpz[i]=bz[i];}
194	greenpot(bpx, bpy, bpz, nx, ny);
195
196	computeBh(bpx, bpy, bz, bh, dims, &mean_hf, mask);
197
198	if (computeGamma(bpx, bpy, bz, bh, dims, &mean_gamma, mask)) mean_gamma = 0.0 / 0.0;
199	drms_setkey_float(outRec, "MEANGAM", mean_gamma);
200
201	computeB_total(bpx, bpy, bz, bt, dims, mask);
202
203	if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask)) mean_derivative_btotal = 0.0 / 0.0;
204	drms_setkey_float(outRec, "MEANGBT", mean_derivative_btotal);
205
206	if (computeBhderivative(bh, dims, &mean_derivative_bh, mask)) mean_derivative_bh = 0.0 / 0.0;
207	drms_setkey_float(outRec, "MEANGBH", mean_derivative_bh);
208
209	if (computeBzderivative(bz, dims, &mean_derivative_bz, mask)) mean_derivative_bz = 0.0 / 0.0; // If fail, fill in NaN
210	drms_setkey_float(outRec, "MEANGBZ", mean_derivative_bz);
211
212	/===========================================/
213	/* Set non-SW keywords */
214
215	drms_copykey(outRec, inRec, "T_REC");
216	drms_copykey(outRec, inRec, "HARPNUM");
217
218	/* Clean up */
219	drms_free_array(inArrayBz);
220	drms_free_array(maskArray);
221
222	}
223
224	drms_close_records(inRecSet, DRMS_FREE_RECORD);
225	drms_close_records(outRecSet, DRMS_INSERT_RECORD);
226
227	return 0;
228
229	}