| 513 |
|
{ |
| 514 |
|
for (j = 0; j <= ny-1; j++) |
| 515 |
|
{ |
| 516 |
< |
if isnan(by[j * nx + i]) continue; |
| 516 |
> |
//if isnan(by[j * nx + i]) continue; |
| 517 |
|
derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5; |
| 518 |
|
} |
| 519 |
|
} |
| 522 |
|
{ |
| 523 |
|
for (j = 1; j <= ny-2; j++) |
| 524 |
|
{ |
| 525 |
< |
if isnan(bx[j * nx + i]) continue; |
| 525 |
> |
//if isnan(bx[j * nx + i]) continue; |
| 526 |
|
dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5; |
| 527 |
|
} |
| 528 |
|
} |
| 550 |
|
} |
| 551 |
|
|
| 552 |
|
j=ny-1; |
| 553 |
< |
for (i = 0; i <= nx-1; i++) |
| 553 |
> |
for (i = 0; i <= nx-1; i++) |
| 554 |
|
{ |
| 555 |
|
if isnan(bx[j * nx + i]) continue; |
| 556 |
|
dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5; |
| 562 |
|
for (j = 0; j <= ny-1; j++) |
| 563 |
|
{ |
| 564 |
|
// calculate jz at all points |
| 565 |
< |
jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix |
| 566 |
< |
|
| 565 |
> |
jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix |
| 566 |
> |
|
| 567 |
> |
//printf("jz[j * nx + i]=%f,i=%d,j=%d\n,",jz[j * nx + i],i,j); // tmp.txt |
| 568 |
> |
//printf("i=%d,j=%d,jz[j * nx + i]=%f\n,",i,j,jz[j * nx + i]); // tmp1.txt |
| 569 |
> |
|
| 570 |
> |
|
| 571 |
|
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]) + |
| 572 |
|
(by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ; |
| 573 |
|
jz_err_squared[j * nx + i]=(jz_err[j * nx + i]*jz_err[j * nx + i]); |
| 581 |
|
/*===========================================*/ |
| 582 |
|
|
| 583 |
|
|
| 584 |
< |
/* Example function 9: Compute quantities on smoothed Jz array */ |
| 585 |
< |
|
| 582 |
< |
// All of the subsequent functions, including this one, use a smoothed Jz array. The smoothing is performed by Jesper's |
| 583 |
< |
// fresize routines. These routines are located at /cvs/JSOC/proj/libs/interpolate. A Gaussian with a FWHM of 4 pixels |
| 584 |
< |
// and truncation width of 12 pixels is used to smooth the array; however, a quick analysis shows that the mean values |
| 585 |
< |
// of qualities like Jz and helicity do not change much as a result of smoothing. The smoothed array will, of course, |
| 586 |
< |
// give a lower total Jz as the stron field pixels have been smoothed out to neighboring weaker field pixels. |
| 584 |
> |
/* Example function 9: Compute quantities on Jz array */ |
| 585 |
> |
// Compute mean and total current on Jz array. |
| 586 |
|
|
| 587 |
|
int computeJzsmooth(float *bx, float *by, int *dims, float *jz, float *jz_smooth, float *jz_err, float *jz_rms_err, float *jz_err_squared_smooth, |
| 588 |
|
float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask, |
| 603 |
|
{ |
| 604 |
|
for (j = 0; j <= ny-1; j++) |
| 605 |
|
{ |
| 607 |
– |
//printf("%f ",jz_smooth[j * nx + i]); |
| 606 |
|
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue; |
| 607 |
|
if isnan(derx[j * nx + i]) continue; |
| 608 |
|
if isnan(dery[j * nx + i]) continue; |
| 611 |
– |
//if isnan(jz_smooth[j * nx + i]) continue; |
| 612 |
– |
//curl += (jz_smooth[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */ |
| 613 |
– |
//us_i += fabs(jz_smooth[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A */ |
| 614 |
– |
//jz_rms_err[j * nx + i] = sqrt(jz_err_squared_smooth[j * nx + i]); |
| 615 |
– |
//err += (jz_rms_err[j * nx + i]*jz_rms_err[j * nx + i]); |
| 609 |
|
if isnan(jz[j * nx + i]) continue; |
| 610 |
|
curl += (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */ |
| 611 |
|
us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A */ |
| 612 |
|
err += (jz_err[j * nx + i]*jz_err[j * nx + i]); |
| 613 |
|
count_mask++; |
| 614 |
|
} |
| 622 |
– |
//printf("\n"); |
| 615 |
|
} |
| 616 |
|
|
| 617 |
|
/* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */ |
| 820 |
|
/*===========================================*/ |
| 821 |
|
/* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */ |
| 822 |
|
// The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV |
| 823 |
< |
// automatically yields erg per cubic centimeter for an input B in Gauss. |
| 823 |
> |
// automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8*PI can come out of the integral; thus, |
| 824 |
> |
// the integral is over B^2 dV and the 8*PI is divided at the end. |
| 825 |
|
// |
| 826 |
|
// Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert |
| 827 |
|
// ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm: |
| 851 |
|
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue; |
| 852 |
|
if isnan(bx[j * nx + i]) continue; |
| 853 |
|
if isnan(by[j * nx + i]) continue; |
| 854 |
< |
sum += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) ) / 8.*PI; |
| 854 |
> |
sum += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0); |
| 855 |
|
err += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]); |
| 863 |
– |
//err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i]; |
| 856 |
|
count_mask++; |
| 857 |
|
} |
| 858 |
|
} |
| 859 |
|
|
| 860 |
< |
*meanpotptr = (sum) / (count_mask); /* Units are ergs per cubic centimeter */ |
| 860 |
> |
*meanpotptr = (sum/(8.*PI)) / (count_mask); /* Units are ergs per cubic centimeter */ |
| 861 |
|
*meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask) |
| 870 |
– |
//*mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask) |
| 862 |
|
|
| 863 |
|
/* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */ |
| 864 |
< |
*totpotptr = sum*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/8.*PI)) ; |
| 865 |
< |
*totpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/8.*PI)); |
| 864 |
> |
*totpotptr = (sum)/(8.*PI); |
| 865 |
> |
*totpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI))); |
| 866 |
|
|
| 867 |
|
printf("MEANPOT=%g\n",*meanpotptr); |
| 868 |
|
printf("MEANPOT_err=%g\n",*meanpot_err_ptr); |
