| 1 |
#include <stdio.h> |
| 2 |
#include <stdlib.h> |
| 3 |
#include <math.h> |
| 4 |
|
| 5 |
|
| 6 |
int errorprop (float *bTotal, float *bAzim, float *bInc, float *ebT, float *ebA, float *ebI, |
| 7 |
float *ebTbA, float *ebTbI, float *ebIbA, |
| 8 |
double lon, double lat, double lonc, double latc, double pAng, int nx, int ny, double x, double y, |
| 9 |
double *BtVariance, double *BpVariance, double *BrVariance) { |
| 10 |
|
| 11 |
/* Inputs: |
| 12 |
* bTotal, bFill, bInc, bAzim: magnetic field inverted by an inversion code. They are |
| 13 |
* field strength, inclination, and azimuth. |
| 14 |
* ebT, ebF, ebI, ebA: variance of bTotal, bInc, bAzim. |
| 15 |
* ebTbF, ebFbI, ebFbA: covariance of bTotal, bFill, bInc, bAzim. |
| 16 |
* ebTbI, ebTbA, ebIbA: covariance of bTotal, bFill, bInc, bAzim. |
| 17 |
* lon, lat: heliographic coordinates of the location of map that the data |
| 18 |
* is projected. |
| 19 |
* lonc, latc: heliographic coordinates of the image disk center. |
| 20 |
* pAng: position angle of the heliographic north pole, measured eastward |
| 21 |
* from the north. |
| 22 |
* nx, ny: size of the array. |
| 23 |
* x, y: location of the pixel in image coordinates. |
| 24 |
* |
| 25 |
* Outputs: |
| 26 |
* BrVariance, BtVariance, BpVariance: variances of these three components of magnetic |
| 27 |
* field. |
| 28 |
* |
| 29 |
* Note: Interpolation useed here is a cubic convolution interpolation. |
| 30 |
*/ |
| 31 |
|
| 32 |
// Xudong Oct 18 2011: Fill factor variances covariances are all NaNs. removed |
| 33 |
// NNB interpolation, just 1 point |
| 34 |
// Fixed definition of azi for all derivatives |
| 35 |
// Feb 1 2019, dBt/dXX is actually dBy/dXX, so the sign is wrong; however they are used in pairs (squared) so final result was okay |
| 36 |
|
| 37 |
static double raddeg = M_PI / 180.; |
| 38 |
double b, inc, azim; |
| 39 |
double a11, a12, a13, a21, a22, a23, a31, a32, a33; |
| 40 |
double dBrdBtotal, dBrdInc, dBrdAzim; |
| 41 |
double dBtdBtotal, dBtdInc, dBtdAzim; |
| 42 |
double dBpdBtotal, dBpdInc, dBpdAzim; |
| 43 |
double errBT, errINC, errAZ; |
| 44 |
double errBTINC, errBTAZ, errINCAZ; |
| 45 |
double BrSigma2, BtSigma2, BpSigma2; |
| 46 |
int ix, iy; |
| 47 |
|
| 48 |
if (x < 1. || x >= (float)(nx-2) || y < 1. || y >= (float)(ny-2)) |
| 49 |
return (1); |
| 50 |
|
| 51 |
a11 = - sin(latc) * sin(pAng) * sin(lon - lonc) + cos(pAng) * cos(lon - lonc); |
| 52 |
a12 = sin(latc) * cos(pAng) * sin(lon - lonc) + sin(pAng) * cos(lon - lonc); |
| 53 |
a13 = -cos(latc) * sin(lon - lonc); |
| 54 |
a21 = -sin(lat) * (sin(latc) * sin(pAng) * cos(lon - lonc) + cos(pAng) * sin(lon - lonc)) |
| 55 |
- cos(lat) * cos(latc) * sin(pAng); |
| 56 |
a22 = sin(lat) * (sin(latc) * cos(pAng) * cos(lon - lonc) - sin(pAng) * sin(lon - lonc)) |
| 57 |
+ cos(lat) * cos(latc) * cos(pAng); |
| 58 |
a23 = -cos(latc) * sin(lat) * cos(lon - lonc) + sin(latc) * cos(lat); |
| 59 |
a31 = cos(lat) * (sin(latc) * sin(pAng) * cos(lon - lonc) + cos(pAng) * sin(lon - lonc)) |
| 60 |
- sin(lat) * cos(latc) * sin(pAng); |
| 61 |
a32 = -cos(lat) * (sin(latc) * cos(pAng) * cos(lon - lonc) - sin(pAng) * sin(lon - lonc)) |
| 62 |
+ sin(lat) * cos(latc) * cos(pAng); |
| 63 |
a33 = cos(lat) * cos(latc) * cos(lon - lonc) + sin(lat) * sin(latc); |
| 64 |
|
| 65 |
ix = (int)x; |
| 66 |
iy = (int)y; |
| 67 |
|
| 68 |
int iAll = iy * nx + ix; |
| 69 |
|
| 70 |
if (isnan(bTotal[iAll])) return(1); |
| 71 |
b = bTotal[iAll]; /* field strength from the pixels used for the interpolation */ |
| 72 |
if (isnan(bInc[iAll])) return(1); |
| 73 |
inc = raddeg * bInc[iAll]; /* inclination, in deg */ |
| 74 |
if (isnan(bAzim[iAll])) return(1); |
| 75 |
azim = raddeg * bAzim[iAll]; /* azimuth, in deg */ |
| 76 |
if (isnan(ebT[iAll])) return(1); |
| 77 |
errBT = ebT[iAll]; /* variance of field strength, in G^2 */ |
| 78 |
if (isnan(ebI[iAll])) return(1); |
| 79 |
errINC = ebI[iAll]; /* variance of inclination, in rad^2 */ |
| 80 |
if (isnan(ebA[iAll])) return(1); |
| 81 |
errAZ = ebA[iAll]; /* variance of azimuth, in rad^2 */ |
| 82 |
if (isnan(ebTbI[iAll])) return(1); |
| 83 |
errBTINC = ebTbI[iAll]; /* covariance of field and inclination, in G.rad */ |
| 84 |
if (isnan(ebTbA[iAll])) return(1); |
| 85 |
errBTAZ = ebTbA[iAll]; /* covariance of field and azimuth, in G.rad */ |
| 86 |
if (isnan(ebIbA[iAll])) return(1); |
| 87 |
errINCAZ = ebIbA[iAll]; /* covariance of inclination and azimuth, in rad^2 */ |
| 88 |
|
| 89 |
|
| 90 |
dBpdBtotal = (- a11 * sin(inc) * sin(azim) + a12 * sin(inc) * cos(azim) + a13 * cos(inc)); |
| 91 |
dBpdInc = b * (- a11 * cos(inc) * sin(azim) + a12 * cos(inc) * cos(azim) - a13 * sin(inc)); |
| 92 |
dBpdAzim = b * (- a11 * sin(inc) * cos(azim) - a12 * sin(inc) * sin(azim)); |
| 93 |
|
| 94 |
dBtdBtotal = (- a21 * sin(inc) * sin(azim) + a22 * sin(inc) * cos(azim) + a23 * cos(inc)) * (-1); |
| 95 |
dBtdInc = b * (- a21 * cos(inc) * sin(azim) + a22 * cos(inc) * cos(azim) - a23 * sin(inc)) * (-1); |
| 96 |
dBtdAzim = b * (- a21 * sin(inc) * cos(azim) - a22 * sin(inc) * sin(azim)) * (-1); |
| 97 |
|
| 98 |
dBrdBtotal = (- a31 * sin(inc) * sin(azim) + a32 * sin(inc) * cos(azim) + a33 * cos(inc)); |
| 99 |
dBrdInc = b * (- a31 * cos(inc) * sin(azim) + a32 * cos(inc) * cos(azim) - a33 * sin(inc)); |
| 100 |
dBrdAzim = b * (- a31 * sin(inc) * cos(azim) - a32 * sin(inc) * sin(azim)); |
| 101 |
|
| 102 |
BrSigma2 = 0.0; |
| 103 |
BtSigma2 = 0.0; |
| 104 |
BpSigma2 = 0.0; |
| 105 |
|
| 106 |
BtSigma2 = BtSigma2 + dBtdBtotal * dBtdBtotal * errBT + |
| 107 |
dBtdInc * dBtdInc * errINC + |
| 108 |
dBtdAzim * dBtdAzim * errAZ + |
| 109 |
2.0 * dBtdBtotal * dBtdInc * errBTINC + |
| 110 |
2.0 * dBtdBtotal * dBtdAzim * errBTAZ + |
| 111 |
2.0 * dBtdInc * dBtdAzim * errINCAZ; |
| 112 |
BpSigma2 = BpSigma2 + dBpdBtotal * dBpdBtotal * errBT + |
| 113 |
dBpdInc * dBpdInc * errINC + |
| 114 |
dBpdAzim * dBpdAzim * errAZ + |
| 115 |
2.0 * dBpdBtotal * dBpdInc * errBTINC + |
| 116 |
2.0 * dBpdBtotal * dBpdAzim * errBTAZ + |
| 117 |
2.0 * dBpdInc * dBpdAzim * errINCAZ; |
| 118 |
BrSigma2 = BrSigma2 + dBrdBtotal * dBrdBtotal * errBT + |
| 119 |
dBrdInc * dBrdInc * errINC + |
| 120 |
dBrdAzim * dBrdAzim * errAZ + |
| 121 |
2.0 * dBrdBtotal * dBrdInc * errBTINC + |
| 122 |
2.0 * dBrdBtotal * dBrdAzim * errBTAZ + |
| 123 |
2.0 * dBrdInc * dBrdAzim * errINCAZ; |
| 124 |
|
| 125 |
*BrVariance = BrSigma2; |
| 126 |
*BtVariance = BtSigma2; |
| 127 |
*BpVariance = BpSigma2; |
| 128 |
|
| 129 |
return 0; |
| 130 |
} |
| 131 |
|
| 132 |
|
| 133 |
|
