1 xudong 1.1 #include <math.h>
 2            
 3            int img2helioVector (double bxImg, double byImg, double bzImg, double *bxHelio, 
 4                                 double *byHelio, double *bzHelio, double lon, double lat,
 5                                 double lonc, double latc, double pAng) {
 6            /*
 7            *     perform tramsformation of a vector from image location (lon, lat) to
 8            *     heliographic center. The formula is from Hagyard (1987), and further
 9            *     developed by Gary & Hagyard (1990).
10            *
11            * Arguments:
12            *
13            *    bxImg, byImg, bzImg: three components of vector magnetic field on image
14            *                         coordinates.
15            *    lon, lat:            heliographic coordinates of the location where the vector field
16            *                         measured. They are in radians.
17            *    lonc, latc:          heliographic coordinates of the image disk center. They are in
18            *                         radians.
19            *    pAng:                position angle of the heliographic north pole, measured eastward
20            *                         from the north. It's in radians.   
21            */
22 xudong 1.1 
23                static double raddeg = M_PI / 180.;
24                double a11, a12, a13, a21, a22, a23, a31, a32, a33;
25            
26                a11 = -sin(latc) * sin(pAng) * sin(lon - lonc) + cos(pAng) * cos(lon - lonc);
27                a12 =  sin(latc) * cos(pAng) * sin(lon - lonc) + sin(pAng) * cos(lon - lonc);
28                a13 = -cos(latc) * sin(lon - lonc);
29                a21 = -sin(lat) * (sin(latc) * sin(pAng) * cos(lon - lonc) + cos(pAng) * sin(lon - lonc)) 
30                      - cos(lat) * cos(latc) * sin(pAng);
31                a22 =  sin(lat) * (sin(latc) * cos(pAng) * cos(lon - lonc) - sin(pAng) * sin(lon - lonc))
32                      + cos(lat) * cos(latc) * cos(pAng);
33                a23 = -cos(latc) * sin(lat) * cos(lon - lonc) + sin(latc) * cos(lat);
34                a31 =  cos(lat) * (sin(latc) * sin(pAng) * cos(lon - lonc) + cos(pAng) * sin(lon - lonc))
35                      - sin(lat) * cos(latc) * sin(pAng);
36                a32 = -cos(lat) * (sin(latc) * cos(pAng) * cos(lon - lonc) - sin(pAng) * sin(lon - lonc))
37                      + sin(lat) * cos(latc) * cos(pAng);
38                a33 =  cos(lat) * cos(latc) * cos(lon - lonc) + sin(lat) * sin(latc);
39            
40                *bxHelio = a11 * bxImg + a12 * byImg + a13 * bzImg;
41                *byHelio = a21 * bxImg + a22 * byImg + a23 * bzImg;
42                *bzHelio = a31 * bxImg + a32 * byImg + a33 * bzImg;
43 xudong 1.1 
44            return 0;
45            }