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

### Diff for /JSOC/proj/sharp/apps/sw_functions.c between version 1.36 and 1.37

version 1.36, 2015/10/30 18:53:02 version 1.37, 2015/10/30 20:21:28
 Line 1
 Line 1

/*=========================================== /*===========================================

The following 14 functions calculate the following spaceweather indices:  The following 14 functions calculate the following spaceweather indices:
 Line 196  int computeGamma(float *bz_err, float *b
 Line 197  int computeGamma(float *bz_err, float *b
//printf("MEANGAM=%f\n",*mean_gamma_ptr);     //printf("MEANGAM=%f\n",*mean_gamma_ptr);
//printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);     //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
printf("bh[200,200],bh_err[200,200]=%f,%f\n",bh[200,200],bh_err[200,200]);

return 0;     return 0;
} }

 Line 249  int computeB_total(float *bx_err, float
 Line 247  int computeB_total(float *bx_err, float
/*===========================================*/ /*===========================================*/
/* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */ /* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */

int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr, float *err_term1, float *err_term2)  int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr, float *err_termAt, float *err_termBt)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 269  int computeBtotalderivative(float *bt, i
 Line 267  int computeBtotalderivative(float *bt, i
for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
{         {
derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;            derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
err_term1[j * nx + i] = (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) ;             err_termAt[j * nx + i] = (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) ;
}         }
}     }

 Line 279  int computeBtotalderivative(float *bt, i
 Line 277  int computeBtotalderivative(float *bt, i
for (j = 1; j <= ny-2; j++)         for (j = 1; j <= ny-2; j++)
{         {
dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;            dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
err_term2[j * nx + i] = (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) ;             err_termBt[j * nx + i] = (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) ;
}         }
}     }

 Line 326  int computeBtotalderivative(float *bt, i
 Line 324  int computeBtotalderivative(float *bt, i
if isnan(derx_bt[j * nx + i]) continue;             if isnan(derx_bt[j * nx + i]) continue;
if isnan(dery_bt[j * nx + i]) continue;             if isnan(dery_bt[j * nx + i]) continue;
sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */             sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */
err += err_term2[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+              err += err_termBt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+
err_term1[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;                     err_termAt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;
}         }
}     }
 Line 344  int computeBtotalderivative(float *bt, i
 Line 342  int computeBtotalderivative(float *bt, i
/*===========================================*/ /*===========================================*/
/* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */ /* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */

int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh, float *err_term1, float *err_term2)  int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh, float *err_termAh, float *err_termBh)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 364  int computeBhderivative(float *bh, float
 Line 362  int computeBhderivative(float *bh, float
for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
{         {
derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;            derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
err_term1[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));             err_termAh[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));
}         }
}     }

 Line 374  int computeBhderivative(float *bh, float
 Line 372  int computeBhderivative(float *bh, float
for (j = 1; j <= ny-2; j++)        for (j = 1; j <= ny-2; j++)
{        {
dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;           dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
err_term2[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));            err_termBh[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));
}        }
}     }

 Line 421  int computeBhderivative(float *bh, float
 Line 419  int computeBhderivative(float *bh, float
if isnan(derx_bh[j * nx + i]) continue;             if isnan(derx_bh[j * nx + i]) continue;
if isnan(dery_bh[j * nx + i]) continue;             if isnan(dery_bh[j * nx + i]) continue;
sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */             sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */
err += err_term2[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+              err += err_termBh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+
err_term1[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;                     err_termAh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;
}         }
}     }
 Line 438  int computeBhderivative(float *bh, float
 Line 436  int computeBhderivative(float *bh, float
/*===========================================*/ /*===========================================*/
/* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */ /* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */

int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz, float *err_term1, float *err_term2)  int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz, float *err_termA, float *err_termB)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 458  int computeBzderivative(float *bz, float
 Line 456  int computeBzderivative(float *bz, float
for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
{         {
derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;            derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
err_term1[j * nx + i] = (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)]));             err_termA[j * nx + i] = (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)]));
}         }
}     }

 Line 468  int computeBzderivative(float *bz, float
 Line 466  int computeBzderivative(float *bz, float
for (j = 1; j <= ny-2; j++)         for (j = 1; j <= ny-2; j++)
{         {
dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;            dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
err_term2[j * nx + i] = (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i]));             err_termB[j * nx + i] = (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i]));
}         }
}     }

 Line 515  int computeBzderivative(float *bz, float
 Line 513  int computeBzderivative(float *bz, float
if isnan(derx_bz[j * nx + i]) continue;             if isnan(derx_bz[j * nx + i]) continue;
if isnan(dery_bz[j * nx + i]) continue;             if isnan(dery_bz[j * nx + i]) continue;
sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i] ); /* Units of Gauss */             sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i] ); /* Units of Gauss */
err += err_term2[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) +              err += err_termB[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) +
err_term1[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;                     err_termA[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;
}         }
}     }
 Line 992  int computeShearAngle(float *bx_err, flo
 Line 990  int computeShearAngle(float *bx_err, flo
dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);             dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);
magnitude_potential   = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));             magnitude_potential   = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));
magnitude_vector      = sqrt( (bx[j * nx + i]*bx[j * nx + i])   + (by[j * nx + i]*by[j * nx + i])   + (bz[j * nx + i]*bz[j * nx + i]) );             magnitude_vector      = sqrt( (bx[j * nx + i]*bx[j * nx + i])   + (by[j * nx + i]*by[j * nx + i])   + (bz[j * nx + i]*bz[j * nx + i]) );

//printf("i,j=%d,%d\n",i,j);
//printf("dotproduct=%f\n",dotproduct);             //printf("dotproduct=%f\n",dotproduct);
//printf("magnitude_potential=%f\n",magnitude_potential);             //printf("magnitude_potential=%f\n",magnitude_potential);
//printf("magnitude_vector=%f\n",magnitude_vector);             //printf("magnitude_vector=%f\n",magnitude_vector);
//printf("dotproduct/(magnitude_potential*magnitude_vector)=%f\n",dotproduct/(magnitude_potential*magnitude_vector));
//printf("acos(dotproduct/(magnitude_potential*magnitude_vector))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector)));
//printf("acos(dotproduct/(magnitude_potential*magnitude_vector)*(180./PI))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI));

shear_angle            = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);             shear_angle            = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
sumsum                  += shear_angle;             sumsum                  += shear_angle;
 Line 1322  void greenpot(float *bx, float *by, floa
 Line 1315  void greenpot(float *bx, float *by, floa
i2e = inx + iwindow;             i2e = inx + iwindow;
if (i2s <   0){i2s =   0;}             if (i2s <   0){i2s =   0;}
if (i2e > nnx){i2e = nnx;}             if (i2e > nnx){i2e = nnx;}
//            for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
for (j2=j2s;j2<2;j2++){for (i2=i2s;i2<2;i2++)              for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
{             {
float val1 = bztmp[nnx*j2 + i2];                 float val1 = bztmp[nnx*j2 + i2];
float rr, r1, r2;
//        r1 = (float)(i2 - inx);
//        r2 = (float)(j2 - iny);
//        rr = r1*r1 + r2*r2;
//        if (rr < rwindow)
//        {
int   di, dj;                 int   di, dj;
di = abs(i2 - inx);                 di = abs(i2 - inx);
dj = abs(j2 - iny);                 dj = abs(j2 - iny);
sum = sum + val1 * rdist[nnx * dj + di] * dz;                 sum = sum + val1 * rdist[nnx * dj + di] * dz;
printf("sum,val1,rdist[nnx * dj + di]=%f,%f,%f\n",sum,val1,rdist[nnx * dj + di]);                  //        }
} }             } }
pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.             pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
}         }
} } // end of OpenMP parallelism     } } // end of OpenMP parallelism

printf("bztmp[0]=%f\n",bztmp[0]);
printf("bztmp[91746]=%f\n",bztmp[91746]);
printf("pfpot[0]=%f\n",pfpot[0]);
printf("pfpot[91746]=%f\n",pfpot[91746]);

// #pragma omp parallel for private(inx)     // #pragma omp parallel for private(inx)
for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)     for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
{     {
 Line 1348  void greenpot(float *bx, float *by, floa
 Line 1342  void greenpot(float *bx, float *by, floa
by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;         by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
} } // end of OpenMP parallelism     } } // end of OpenMP parallelism

printf("bx[0]=%f\n",bx[0]);
printf("by[0]=%f\n",by[0]);
printf("bx[91746]=%f\n",bx[91746]);
printf("by[91746]=%f\n",by[91746]);
free(rdist);     free(rdist);
free(pfpot);     free(pfpot);
free(bztmp);     free(bztmp);
 Line 1365  char *sw_functions_version() // Returns
 Line 1355  char *sw_functions_version() // Returns
return strdup("\$Id");     return strdup("\$Id");
} }

/* ---------------- end of this file ----------------*/ /* ---------------- end of this file ----------------*/

Legend:
 Removed from v.1.36 changed lines Added in v.1.37