(file) Return to swharp_vectorB.c CVS log (file) (dir) Up to [Development] / JSOC / proj / sharp / apps

Diff for /JSOC/proj/sharp/apps/swharp_vectorB.c between version 1.1 and 1.2

version 1.1, 2012/04/04 22:56:31 version 1.2, 2012/06/22 19:40:27
Line 2 
Line 2 
  *  swharp_vectorB.  *  swharp_vectorB.
  *  *
  *  Created by Xudong Sun on 8/22/11.  *  Created by Xudong Sun on 8/22/11.
  *  Modified to -- include ALL spaceweather keywords by Monica Bobra 25 Aug 2011   *
  *              -- include potential field calculation   *  Modified by Monica Bobra to
  *              -- run only on los data 5 March 2011   *              -- 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  *  Bz arrays
  *  Write out abs(B) as data segment and a few keywords as SW index  *  Write out abs(B) as data segment and a few keywords as SW index
  *  *
Line 12 
Line 14 
  *  First use the bmap module to create the in= and mask= parameters.  *  First use the bmap module to create the in= and mask= parameters.
  *  *
  *  then run this module:  *  then run this module:
  *  swharp_vectorB  "in=su_mbobra.test_mmap_me[][]" /   *  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.test_mmap_bitmap_me[][]" "out=su_mbobra.swharp_test_v2" "dzvalue=0.001"   *  "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 <jsoc_main.h>
 #include <stdio.h> #include <stdio.h>
 #include <stdlib.h> #include <stdlib.h>
 #include <math.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" #include "sharp_functions.c"
  
 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);} #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
 #define SHOW(msg) {printf("%s", msg); fflush(stdout);} #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
 #define IN_FILES        3   /* Number of input files */ #define IN_FILES        3   /* Number of input files */
 #define PI (3.141592653589793)  #define PI (3.141592653589793)        /* Ratio of circumference to diameter of a circle*/
   #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
  
 /* declaring all the functions */ /* 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);  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 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 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 readFits(char *filename, float **image, int *dims);
Line 58  ModuleArgs_t module_args[] =
Line 71  ModuleArgs_t module_args[] =
     {ARG_END}     {ARG_END}
 }; };
  
   
 int DoIt(void) int DoIt(void)
 { {
  
Line 70  int DoIt(void)
Line 82  int DoIt(void)
         DRMS_Record_t *inRec, *outRec, *maskRec;         DRMS_Record_t *inRec, *outRec, *maskRec;
         DRMS_Segment_t *inSegBx, *inSegBy, *inSegBz, *outSeg, *maskSeg;         DRMS_Segment_t *inSegBx, *inSegBy, *inSegBz, *outSeg, *maskSeg;
         DRMS_Array_t *inArrayBx, *inArrayBy, *inArrayBz, *outArray, *maskArray;         DRMS_Array_t *inArrayBx, *inArrayBy, *inArrayBz, *outArray, *maskArray;
         float *bx, *by, *bz, *outData, *bh, *bt, *jz, *bpx, *bpy, *bpz;          float *inverseMu, *mu, *bx, *by, *bz, *outData, *bh, *bt, *jz, *bpx, *bpy, *bpz;
         int *mask;         int *mask;
         int dims[2], nxny, nx, ny;              // dimensions;  NAXIS1 = dims[0] which is the number of columns.         int dims[2], nxny, nx, ny;              // dimensions;  NAXIS1 = dims[0] which is the number of columns.
         float mean_vf = 0.0;         float mean_vf = 0.0;
Line 143  int DoIt(void)
Line 155  int DoIt(void)
                 maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);                 maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
                 if (status) DIE("No mask data file found. \n");                 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;                 bx   = (float *)inArrayBx->data;
                 by   = (float *)inArrayBy->data;                 by   = (float *)inArrayBy->data;
                 bz   = (float *)inArrayBz->data;                 bz   = (float *)inArrayBz->data;
Line 179  int DoIt(void)
Line 210  int DoIt(void)
                 /*===========================================*/                 /*===========================================*/
                 /* Malloc some arrays     */                 /* Malloc some arrays     */
  
                   inverseMu   = (float *)malloc(nx*ny*sizeof(float));
                   mu   = (float *)malloc(nx*ny*sizeof(float));
                 bh   = (float *)malloc(nx*ny*sizeof(float));                 bh   = (float *)malloc(nx*ny*sizeof(float));
                 bt   = (float *)malloc(nx*ny*sizeof(float));                 bt   = (float *)malloc(nx*ny*sizeof(float));
                 jz   = (float *)malloc(nx*ny*sizeof(float));                 jz   = (float *)malloc(nx*ny*sizeof(float));
Line 189  int DoIt(void)
Line 222  int DoIt(void)
                 /*===========================================*/                 /*===========================================*/
                 /* SW Keyword computation */                 /* SW Keyword computation */
  
                 if (computeAbsFlux(bz, dims, &absFlux, &mean_vf, mask))                  computeMu(bz, dims, mu, inverseMu);
   
                   if (computeAbsFlux(bz, dims, &absFlux, &mean_vf, mask, inverseMu))
                 {                 {
                    absFlux = 0.0 / 0.0;         // If fail, fill in NaN                    absFlux = 0.0 / 0.0;         // If fail, fill in NaN
                    mean_vf = 0.0 / 0.0;                    mean_vf = 0.0 / 0.0;

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
Karen Tian
 Powered by
ViewCVS 0.9.4