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

Diff for /JSOC/proj/sharp/apps/sw_functions.c between version 1.31 and 1.33

version 1.31, 2014/06/05 21:27:04 version 1.33, 2015/01/23 01:18:13
Line 563  int computeBzderivative(float *bz, float
Line 563  int computeBzderivative(float *bz, float
 //              float *noiseby, float *noisebz) //              float *noiseby, float *noisebz)
  
 int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared, int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
               int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)                int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *err_term1, float *err_term2)
  
  
 { {
Line 584  int computeJz(float *bx_err, float *by_e
Line 584  int computeJz(float *bx_err, float *by_e
         {         {
             if isnan(by[j * nx + i]) continue;             if isnan(by[j * nx + i]) continue;
             derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;             derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
               err_term1[j * nx + i] = (by_err[j * nx + i+1])*(by_err[j * nx + i+1]) + (by_err[j * nx + i-1])*(by_err[j * nx + i-1]);
         }         }
     }     }
  
Line 593  int computeJz(float *bx_err, float *by_e
Line 594  int computeJz(float *bx_err, float *by_e
         {         {
             if isnan(bx[j * nx + i]) continue;             if isnan(bx[j * nx + i]) continue;
             dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;             dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
               err_term2[j * nx + i] = (bx_err[(j+1) * nx + i])*(bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i])*(bx_err[(j-1) * nx + i]);
         }         }
     }     }
  
Line 602  int computeJz(float *bx_err, float *by_e
Line 604  int computeJz(float *bx_err, float *by_e
     {     {
         if isnan(by[j * nx + i]) continue;         if isnan(by[j * nx + i]) continue;
         derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;         derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
           err_term1[j * nx + i] = ( (3*by_err[j * nx + i])*(3*by_err[j * nx + i]) + (4*by_err[j * nx + (i+1)])*(4*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])*(by_err[j * nx + (i+2)]) );
     }     }
  
     i=nx-1;     i=nx-1;
Line 609  int computeJz(float *bx_err, float *by_e
Line 612  int computeJz(float *bx_err, float *by_e
     {     {
         if isnan(by[j * nx + i]) continue;         if isnan(by[j * nx + i]) continue;
         derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;         derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
           err_term1[j * nx + i] = ( (3*by_err[j * nx + i])*(3*by_err[j * nx + i]) + (4*by_err[j * nx + (i+1)])*(4*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])*(by_err[j * nx + (i+2)]) );
     }     }
  
     j=0;     j=0;
Line 616  int computeJz(float *bx_err, float *by_e
Line 620  int computeJz(float *bx_err, float *by_e
     {     {
         if isnan(bx[j * nx + i]) continue;         if isnan(bx[j * nx + i]) continue;
         dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;         dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
           err_term2[j * nx + i] = ( (3*bx_err[j*nx + i])*(3*bx_err[j*nx + i]) + (4*bx_err[(j+1) * nx + i])*(4*bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])*(bx_err[(j+2) * nx + i]) );
     }     }
  
     j=ny-1;     j=ny-1;
Line 623  int computeJz(float *bx_err, float *by_e
Line 628  int computeJz(float *bx_err, float *by_e
     {     {
         if isnan(bx[j * nx + i]) continue;         if isnan(bx[j * nx + i]) continue;
         dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;         dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
           err_term2[j * nx + i] = ( (3*bx_err[j*nx + i])*(3*bx_err[j*nx + i]) + (4*bx_err[(j+1) * nx + i])*(4*bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])*(bx_err[(j+2) * nx + i]) );
   
     }     }
  
     for (i = 1; i <= nx-2; i++)  
       for (i = 0; i <= nx-1; i++)
     {     {
         for (j = 1; j <= ny-2; j++)          for (j = 0; j <= ny-1; j++)
         {         {
             // calculate jz at all points             // calculate jz at all points
   
             jz[j * nx + i]            = (derx[j * nx + i]-dery[j * nx + i]);       // jz is in units of Gauss/pix             jz[j * nx + i]            = (derx[j * nx + i]-dery[j * nx + i]);       // jz is in units of Gauss/pix
             jz_err[j * nx + i]        = 0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +              jz_err[j * nx + i]        = 0.5*sqrt( err_term1[j * nx + i] + err_term2[j * nx + i] ) ;
                                                  (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;  
             jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);             jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);
             count_mask++;             count_mask++;
   
         }         }
     }     }
         return 0;         return 0;
Line 832  int computeHelicity(float *jz_err, float
Line 837  int computeHelicity(float *jz_err, float
 /* Example function 12:  Sum of Absolute Value per polarity  */ /* Example function 12:  Sum of Absolute Value per polarity  */
  
 //  The Sum of the Absolute Value per polarity is defined as the following: //  The Sum of the Absolute Value per polarity is defined as the following:
 //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.  //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per arcsecond.
 //  The units of jz are in G/pix. In this case, we would have the following: //  The units of jz are in G/pix. In this case, we would have the following:
 //  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF), //  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
 //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS) //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
Line 869  int computeSumAbsPerPolarity(float *jz_e
Line 874  int computeSumAbsPerPolarity(float *jz_e
         }         }
     }     }
  
         *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */      *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are Amperes per arcsecond */
     *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));     *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
     //printf("SAVNCPP=%g\n",*totaljzptr);     //printf("SAVNCPP=%g\n",*totaljzptr);
     //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);     //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
Line 1160  int computeR(float *bz_err, float *los,
Line 1165  int computeR(float *bz_err, float *los,
         for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
         {         {
             index = j * nx1 + i;             index = j * nx1 + i;
               if isnan(pmapn[index]) continue;
               if isnan(rim[index]) continue;
             sum += pmapn[index]*abs(rim[index]);             sum += pmapn[index]*abs(rim[index]);
         }         }
     }     }
Line 1169  int computeR(float *bz_err, float *los,
Line 1176  int computeR(float *bz_err, float *los,
     else     else
         *Rparam = log10(sum);         *Rparam = log10(sum);
  
       //printf("R_VALUE=%f\n",*Rparam);
   
     free_fresize(&fresboxcar);     free_fresize(&fresboxcar);
     free_fresize(&fresgauss);     free_fresize(&fresgauss);
  
Line 1201  int computeLorentz(float *bx, float *by
Line 1210  int computeLorentz(float *bx, float *by
     double k_h = -1.0 * area / (4. * PI) / 1.0e20;     double k_h = -1.0 * area / (4. * PI) / 1.0e20;
     double k_z = area / (8. * PI) / 1.0e20;     double k_z = area / (8. * PI) / 1.0e20;
  
     /* Multiplier */  
     float vectorMulti[] = {1.,-1.,1.};  
   
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
  
     for (int i = 0; i < nxny; i++)     for (int i = 0; i < nxny; i++)
Line 1212  int computeLorentz(float *bx, float *by
Line 1218  int computeLorentz(float *bx, float *by
        if isnan(bx[i]) continue;        if isnan(bx[i]) continue;
        if isnan(by[i]) continue;        if isnan(by[i]) continue;
        if isnan(bz[i]) continue;        if isnan(bz[i]) continue;
        fx[i]  = (bx[i] * vectorMulti[0]) * (bz[i] * vectorMulti[2]) * k_h;         fx[i]  = bx[i] * bz[i] * k_h;
        fy[i]  = (by[i] * vectorMulti[1]) * (bz[i] * vectorMulti[2]) * k_h;         fy[i]  = by[i] * bz[i] * k_h;
        fz[i]  = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z;        fz[i]  = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z;
        bsq    = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i];        bsq    = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i];
        totfx  += fx[i]; totfy += fy[i]; totfz += fz[i];        totfx  += fx[i]; totfy += fy[i]; totfz += fz[i];
Line 1228  int computeLorentz(float *bx, float *by
Line 1234  int computeLorentz(float *bx, float *by
     *epsy_ptr   = (totfy / k_h) / totbsq;     *epsy_ptr   = (totfy / k_h) / totbsq;
     *epsz_ptr   = (totfz / k_z) / totbsq;     *epsz_ptr   = (totfz / k_z) / totbsq;
  
       //printf("TOTBSQ=%f\n",*totbsq_ptr);
   
     return 0;     return 0;
  
 } }

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