sharp/apps/errorprop.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


int errorprop (float *bTotal, float *bAzim, float *bInc, float *ebT, float *ebA, float *ebI,
               float *ebTbA, float *ebTbI, float *ebIbA, 
               double lon, double lat, double lonc, double latc, double pAng, int nx, int ny, double x, double y,
               double *BtVariance, double *BpVariance, double *BrVariance) {

/* Inputs:
*  bTotal, bFill, bInc, bAzim: magnetic field inverted by an inversion code. They are
*                              field strength, inclination, and azimuth.
* ebT, ebF, ebI, ebA:          variance of bTotal, bInc, bAzim.
* ebTbF, ebFbI, ebFbA:         covariance of bTotal, bFill, bInc, bAzim.
* ebTbI, ebTbA, ebIbA:         covariance of bTotal, bFill, bInc, bAzim.
* lon, lat:                    heliographic coordinates of the location of map that the data
*                              is projected.
* lonc, latc:                  heliographic coordinates of the image disk center.
* pAng:                        position angle of the heliographic north pole, measured eastward
*                              from the north.
* nx, ny:                      size of the array.
* x, y:                        location of the pixel in image coordinates.
*
* Outputs:
* BrVariance, BtVariance, BpVariance: variances of these three components of magnetic
*                                     field.
*
* Note:                 Interpolation useed here is a cubic convolution interpolation.
*/
  
  // Xudong Oct 18 2011: Fill factor variances covariances are all NaNs. removed
  // NNB interpolation, just 1 point
  // Fixed definition of azi for all derivatives
  // 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
 
  static double raddeg = M_PI / 180.;
  double b, inc, azim;
  double a11, a12, a13, a21, a22, a23, a31, a32, a33;
  double dBrdBtotal, dBrdInc, dBrdAzim;
  double dBtdBtotal, dBtdInc, dBtdAzim;
  double dBpdBtotal, dBpdInc, dBpdAzim;
  double errBT, errINC, errAZ;
  double errBTINC, errBTAZ, errINCAZ;
  double BrSigma2, BtSigma2, BpSigma2;
  int ix, iy;

  if (x < 1. || x >= (float)(nx-2) || y < 1. || y >= (float)(ny-2))
  return (1);
  
  a11 = - sin(latc) * sin(pAng) * sin(lon - lonc) + cos(pAng) * cos(lon - lonc);
  a12 =  sin(latc) * cos(pAng) * sin(lon - lonc) + sin(pAng) * cos(lon - lonc);
  a13 = -cos(latc) * sin(lon - lonc);
  a21 = -sin(lat) * (sin(latc) * sin(pAng) * cos(lon - lonc) + cos(pAng) * sin(lon - lonc))
        - cos(lat) * cos(latc) * sin(pAng);
  a22 =  sin(lat) * (sin(latc) * cos(pAng) * cos(lon - lonc) - sin(pAng) * sin(lon - lonc))
        + cos(lat) * cos(latc) * cos(pAng);
  a23 = -cos(latc) * sin(lat) * cos(lon - lonc) + sin(latc) * cos(lat);
  a31 =  cos(lat) * (sin(latc) * sin(pAng) * cos(lon - lonc) + cos(pAng) * sin(lon - lonc))
        - sin(lat) * cos(latc) * sin(pAng);
  a32 = -cos(lat) * (sin(latc) * cos(pAng) * cos(lon - lonc) - sin(pAng) * sin(lon - lonc))
        + sin(lat) * cos(latc) * cos(pAng);
  a33 =  cos(lat) * cos(latc) * cos(lon - lonc) + sin(lat) * sin(latc);

  ix = (int)x;
  iy = (int)y;

  int iAll = iy * nx + ix;

  if (isnan(bTotal[iAll])) return(1);
  b = bTotal[iAll]; /* field strength from the pixels used for the interpolation */
  if (isnan(bInc[iAll])) return(1);
  inc = raddeg * bInc[iAll]; /* inclination, in deg */
  if (isnan(bAzim[iAll])) return(1);
  azim = raddeg * bAzim[iAll]; /* azimuth, in deg */
  if (isnan(ebT[iAll])) return(1);
  errBT = ebT[iAll]; /* variance of field strength, in G^2 */
  if (isnan(ebI[iAll])) return(1);
  errINC = ebI[iAll]; /* variance of inclination, in rad^2 */
  if (isnan(ebA[iAll])) return(1);
  errAZ = ebA[iAll]; /* variance of azimuth, in rad^2 */
  if (isnan(ebTbI[iAll])) return(1);
  errBTINC = ebTbI[iAll]; /* covariance of field and inclination, in G.rad */
  if (isnan(ebTbA[iAll])) return(1);
  errBTAZ = ebTbA[iAll];  /* covariance of field and azimuth, in G.rad */
  if (isnan(ebIbA[iAll])) return(1);
  errINCAZ = ebIbA[iAll]; /* covariance of inclination and azimuth, in rad^2 */


  dBpdBtotal = (- a11 * sin(inc) * sin(azim) + a12 * sin(inc) * cos(azim) + a13 * cos(inc));
  dBpdInc = b * (- a11 * cos(inc) * sin(azim) + a12 * cos(inc) * cos(azim) - a13 * sin(inc));
  dBpdAzim = b * (- a11 * sin(inc) * cos(azim) - a12 * sin(inc) * sin(azim));
  
  dBtdBtotal = (- a21 * sin(inc) * sin(azim) + a22 * sin(inc) * cos(azim) + a23 * cos(inc)) * (-1);
  dBtdInc = b * (- a21 * cos(inc) * sin(azim) + a22 * cos(inc) * cos(azim) - a23 * sin(inc)) * (-1);
  dBtdAzim = b * (- a21 * sin(inc) * cos(azim) - a22 * sin(inc) * sin(azim)) * (-1);
  
  dBrdBtotal = (- a31 * sin(inc) * sin(azim) + a32 * sin(inc) * cos(azim) + a33 * cos(inc));
  dBrdInc = b * (- a31 * cos(inc) * sin(azim) + a32 * cos(inc) * cos(azim) - a33 * sin(inc));
  dBrdAzim = b * (- a31 * sin(inc) * cos(azim) - a32 * sin(inc) * sin(azim));

  BrSigma2 = 0.0;
  BtSigma2 = 0.0;
  BpSigma2 = 0.0;

  BtSigma2 = BtSigma2 + dBtdBtotal * dBtdBtotal * errBT + 
                        dBtdInc * dBtdInc * errINC + 
                        dBtdAzim * dBtdAzim * errAZ + 
                        2.0 * dBtdBtotal * dBtdInc * errBTINC +
                        2.0 * dBtdBtotal * dBtdAzim * errBTAZ +
                        2.0 * dBtdInc * dBtdAzim * errINCAZ;
  BpSigma2 = BpSigma2 + dBpdBtotal * dBpdBtotal * errBT +
                        dBpdInc * dBpdInc * errINC + 
                        dBpdAzim * dBpdAzim * errAZ + 
                        2.0 * dBpdBtotal * dBpdInc * errBTINC +
                        2.0 * dBpdBtotal * dBpdAzim * errBTAZ +
                        2.0 * dBpdInc * dBpdAzim * errINCAZ;
  BrSigma2 = BrSigma2 + dBrdBtotal * dBrdBtotal * errBT + 
                        dBrdInc * dBrdInc * errINC + 
                        dBrdAzim * dBrdAzim * errAZ + 
                        2.0 * dBrdBtotal * dBrdInc * errBTINC +
                        2.0 * dBrdBtotal * dBrdAzim * errBTAZ +
                        2.0 * dBrdInc * dBrdAzim * errINCAZ;

  *BrVariance = BrSigma2;
  *BtVariance = BtSigma2;
  *BpVariance = BpSigma2;
  
  return 0;
}


Revision:	1.2
Committed:	Sat Feb 2 04:53:58 2019 UTC (4 years, 7 months ago) by xudong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Ver_LATEST, Ver_9-3, Ver_9-41, Ver_9-5, Ver_9-4, HEAD
Changes since 1.1:	+5 -4 lines
Log Message:	fixed signs of dBt/dXX
#	User	Rev	Content
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			* 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			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	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			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