1 xudong 1.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 xudong 1.1 * 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            
 36              static double raddeg = M_PI / 180.;
 37              double b, inc, azim;
 38              double a11, a12, a13, a21, a22, a23, a31, a32, a33;
 39              double dBrdBtotal, dBrdInc, dBrdAzim;
 40              double dBtdBtotal, dBtdInc, dBtdAzim;
 41              double dBpdBtotal, dBpdInc, dBpdAzim;
 42              double errBT, errINC, errAZ;
 43 xudong 1.1   double errBTINC, errBTAZ, errINCAZ;
 44              double BrSigma2, BtSigma2, BpSigma2;
 45              int ix, iy;
 46            
 47              if (x < 1. || x >= (float)(nx-2) || y < 1. || y >= (float)(ny-2))
 48              return (1);
 49              
 50              a11 = - sin(latc) * sin(pAng) * sin(lon - lonc) + cos(pAng) * cos(lon - lonc);
 51              a12 =  sin(latc) * cos(pAng) * sin(lon - lonc) + sin(pAng) * cos(lon - lonc);
 52              a13 = -cos(latc) * sin(lon - lonc);
 53              a21 = -sin(lat) * (sin(latc) * sin(pAng) * cos(lon - lonc) + cos(pAng) * sin(lon - lonc))
 54                    - cos(lat) * cos(latc) * sin(pAng);
 55              a22 =  sin(lat) * (sin(latc) * cos(pAng) * cos(lon - lonc) - sin(pAng) * sin(lon - lonc))
 56                    + cos(lat) * cos(latc) * cos(pAng);
 57              a23 = -cos(latc) * sin(lat) * cos(lon - lonc) + sin(latc) * cos(lat);
 58              a31 =  cos(lat) * (sin(latc) * sin(pAng) * cos(lon - lonc) + cos(pAng) * sin(lon - lonc))
 59                    - sin(lat) * cos(latc) * sin(pAng);
 60              a32 = -cos(lat) * (sin(latc) * cos(pAng) * cos(lon - lonc) - sin(pAng) * sin(lon - lonc))
 61                    + sin(lat) * cos(latc) * cos(pAng);
 62              a33 =  cos(lat) * cos(latc) * cos(lon - lonc) + sin(lat) * sin(latc);
 63            
 64 xudong 1.1   ix = (int)x;
 65              iy = (int)y;
 66            
 67              int iAll = iy * nx + ix;
 68            
 69              if (isnan(bTotal[iAll])) return(1);
 70              b = bTotal[iAll]; /* field strength from the pixels used for the interpolation */
 71              if (isnan(bInc[iAll])) return(1);
 72              inc = raddeg * bInc[iAll]; /* inclination, in deg */
 73              if (isnan(bAzim[iAll])) return(1);
 74              azim = raddeg * bAzim[iAll]; /* azimuth, in deg */
 75              if (isnan(ebT[iAll])) return(1);
 76              errBT = ebT[iAll]; /* variance of field strength, in G^2 */
 77              if (isnan(ebI[iAll])) return(1);
 78              errINC = ebI[iAll]; /* variance of inclination, in rad^2 */
 79              if (isnan(ebA[iAll])) return(1);
 80              errAZ = ebA[iAll]; /* variance of azimuth, in rad^2 */
 81              if (isnan(ebTbI[iAll])) return(1);
 82              errBTINC = ebTbI[iAll]; /* covariance of field and inclination, in G.rad */
 83              if (isnan(ebTbA[iAll])) return(1);
 84              errBTAZ = ebTbA[iAll];  /* covariance of field and azimuth, in G.rad */
 85 xudong 1.1   if (isnan(ebIbA[iAll])) return(1);
 86              errINCAZ = ebIbA[iAll]; /* covariance of inclination and azimuth, in rad^2 */
 87            
 88            
 89              dBpdBtotal = (- a11 * sin(inc) * sin(azim) + a12 * sin(inc) * cos(azim) + a13 * cos(inc));
 90              dBpdInc = b * (- a11 * cos(inc) * sin(azim) + a12 * cos(inc) * cos(azim) - a13 * sin(inc));
 91              dBpdAzim = b * (- a11 * sin(inc) * cos(azim) - a12 * sin(inc) * sin(azim));
 92              
 93              dBtdBtotal = (- a21 * sin(inc) * sin(azim) + a22 * sin(inc) * cos(azim) + a23 * cos(inc));
 94              dBtdInc = b * (- a21 * cos(inc) * sin(azim) + a22 * cos(inc) * cos(azim) - a23 * sin(inc));
 95              dBtdAzim = b * (- a21 * sin(inc) * cos(azim) - a22 * sin(inc) * sin(azim));
 96              
 97              dBrdBtotal = (- a31 * sin(inc) * sin(azim) + a32 * sin(inc) * cos(azim) + a33 * cos(inc));
 98              dBrdInc = b * (- a31 * cos(inc) * sin(azim) + a32 * cos(inc) * cos(azim) - a33 * sin(inc));
 99              dBrdAzim = b * (- a31 * sin(inc) * cos(azim) - a32 * sin(inc) * sin(azim));
100            
101              BrSigma2 = 0.0;
102              BtSigma2 = 0.0;
103              BpSigma2 = 0.0;
104            
105              BtSigma2 = BtSigma2 + dBtdBtotal * dBtdBtotal * errBT + 
106 xudong 1.1                         dBtdInc * dBtdInc * errINC + 
107                                    dBtdAzim * dBtdAzim * errAZ + 
108                                    2.0 * dBtdBtotal * dBtdInc * errBTINC +
109                                    2.0 * dBtdBtotal * dBtdAzim * errBTAZ +
110                                    2.0 * dBtdInc * dBtdAzim * errINCAZ;
111              BpSigma2 = BpSigma2 + dBpdBtotal * dBpdBtotal * errBT +
112                                    dBpdInc * dBpdInc * errINC + 
113                                    dBpdAzim * dBpdAzim * errAZ + 
114                                    2.0 * dBpdBtotal * dBpdInc * errBTINC +
115                                    2.0 * dBpdBtotal * dBpdAzim * errBTAZ +
116                                    2.0 * dBpdInc * dBpdAzim * errINCAZ;
117              BrSigma2 = BrSigma2 + dBrdBtotal * dBrdBtotal * errBT + 
118                                    dBrdInc * dBrdInc * errINC + 
119                                    dBrdAzim * dBrdAzim * errAZ + 
120                                    2.0 * dBrdBtotal * dBrdInc * errBTINC +
121                                    2.0 * dBrdBtotal * dBrdAzim * errBTAZ +
122                                    2.0 * dBrdInc * dBrdAzim * errINCAZ;
123            
124              *BrVariance = BrSigma2;
125              *BtVariance = BtSigma2;
126              *BpVariance = BpSigma2;
127 xudong 1.1   
128              return 0;
129            }
130            
131            
132