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 xudong 1.2   // 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 xudong 1.1   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 xudong 1.1   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 xudong 1.1   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 xudong 1.2   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 xudong 1.1   
 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 xudong 1.1   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