1 mbobra 1.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 mbobra 1.1 #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.504277*696000000.0*100.)/(943.) /* cm/pixel = (CDELT1*RSUN_REF*100./RSUN_OBS) */
 33            
 34            /* CMSQUARED = CMX*CMY
 35               CMX = CDELT1*RSUN_REF*(#of pixels in the x-direction)*100/RSUN_OBS
 36               CMY = CDELT2*RSUN_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 mbobra 1.1 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 mbobra 1.1 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 mbobra 1.1         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 mbobra 1.1 	
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 mbobra 1.1                 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 mbobra 1.1 		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 mbobra 1.1 
170            		outRec = outRecSet->records[irec];
171                            drms_setlink_static(outRec, "SRCLINK",  inRec->recnum);
172            
173                            /*===========================================*/    
174                            /* Malloc some arrays     */
175            
176                            bh   = (float *)malloc(nx*ny*sizeof(float));
177                            bt   = (float *)malloc(nx*ny*sizeof(float));
178                            jz   = (float *)malloc(nx*ny*sizeof(float));
179                            bpx  = (float *)malloc(nx*ny*sizeof(float));
180                            bpy  = (float *)malloc(nx*ny*sizeof(float));
181                            bpz  = (float *)malloc(nx*ny*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 mbobra 1.1                 }
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 mbobra 1.1 
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            }