1 | /*!\file Xy2llx.cpp
|
---|
2 | */
|
---|
3 |
|
---|
4 | #include "./Xy2llx.h"
|
---|
5 | #include "../../shared/shared.h"
|
---|
6 | #include "../../include/include.h"
|
---|
7 | #include "../../toolkits/toolkits.h"
|
---|
8 | #include "../../EnumDefinitions/EnumDefinitions.h"
|
---|
9 |
|
---|
10 | int Xy2llx(double* lat, double* lon,
|
---|
11 | double* x, double* y, int ncoord,
|
---|
12 | int sgn){
|
---|
13 |
|
---|
14 | /* This is a cpp conversion of the following:
|
---|
15 | %XY2LL - converts xy to lat long
|
---|
16 | %
|
---|
17 | % Converts Polar Stereographic (X,Y) coordinates for the polar regions to
|
---|
18 | % latitude and longitude Stereographic (X,Y) coordinates for the polar
|
---|
19 | % regions.
|
---|
20 | % Author: Michael P. Schodlok, December 2003 (map2xy.m)
|
---|
21 | %
|
---|
22 | % Usage:
|
---|
23 | % [lat,lon] = xy2ll(x,y,sgn);
|
---|
24 | % [lat,lon] = xy2ll(x,y,sgn,central_meridian,standard_parallel);
|
---|
25 | %
|
---|
26 | % - sgn = Sign of latitude +1 : north latitude (default is mer=45 lat=70)
|
---|
27 | % -1 : south latitude (default is mer=0 lat=71)
|
---|
28 | */
|
---|
29 |
|
---|
30 | double delta,slat;
|
---|
31 | bool flag=true;
|
---|
32 |
|
---|
33 | /* Get central_meridian and standard_parallel depending on hemisphere */
|
---|
34 | if (sgn == 1) {
|
---|
35 | delta = 45;
|
---|
36 | slat = 70;
|
---|
37 | _printf_(flag,"Warning: expecting coordinates in polar stereographic (Std Latitude: 70N Meridian: 45).\n");
|
---|
38 | }
|
---|
39 | else if (sgn == -1) {
|
---|
40 | delta = 0;
|
---|
41 | slat = 71;
|
---|
42 | _printf_(flag,"Warning: expecting coordinates in polar stereographic (Std Latitude: 71S Meridian: 0).\n");
|
---|
43 | }
|
---|
44 | else
|
---|
45 | _error_("Sign should be either +1 or -1.\n");
|
---|
46 |
|
---|
47 | return(Xy2llx(lat,lon,
|
---|
48 | x,y,ncoord,
|
---|
49 | sgn,delta,slat));
|
---|
50 | }
|
---|
51 |
|
---|
52 | int Xy2llx(double* lat, double* lon,
|
---|
53 | double* x, double* y, int ncoord,
|
---|
54 | int sgn, double central_meridian, double standard_parallel){
|
---|
55 |
|
---|
56 | /* This is a cpp conversion of the following:
|
---|
57 | %XY2LL - converts xy to lat long
|
---|
58 | %
|
---|
59 | % Converts Polar Stereographic (X,Y) coordinates for the polar regions to
|
---|
60 | % latitude and longitude Stereographic (X,Y) coordinates for the polar
|
---|
61 | % regions.
|
---|
62 | % Author: Michael P. Schodlok, December 2003 (map2xy.m)
|
---|
63 | %
|
---|
64 | % Usage:
|
---|
65 | % [lat,lon] = xy2ll(x,y,sgn);
|
---|
66 | % [lat,lon] = xy2ll(x,y,sgn,central_meridian,standard_parallel);
|
---|
67 | %
|
---|
68 | % - sgn = Sign of latitude +1 : north latitude (default is mer=45 lat=70)
|
---|
69 | % -1 : south latitude (default is mer=0 lat=71)
|
---|
70 | */
|
---|
71 |
|
---|
72 | int i,iret=0;
|
---|
73 | double delta,slat;
|
---|
74 | double cde,re,ex2,ex,pi;
|
---|
75 | double sl,rho,cm,T,chi;
|
---|
76 |
|
---|
77 | if ((sgn != 1) && (sgn != -1))
|
---|
78 | _error_("Sign should be either +1 or -1.\n");
|
---|
79 |
|
---|
80 | delta = central_meridian;
|
---|
81 | slat = standard_parallel;
|
---|
82 |
|
---|
83 | /* Conversion constant from degrees to radians */
|
---|
84 | cde = 57.29577951;
|
---|
85 | /* Radius of the earth in meters */
|
---|
86 | re = 6378.273*pow(10.,3.);
|
---|
87 | /* Eccentricity of the Hughes ellipsoid squared */
|
---|
88 | ex2 = .006693883;
|
---|
89 | /* Eccentricity of the Hughes ellipsoid */
|
---|
90 | ex = sqrt(ex2);
|
---|
91 | /* pi is not defined in cpp */
|
---|
92 | pi = 4.*atan(1.);
|
---|
93 |
|
---|
94 | /* loop over all the coordinate pairs */
|
---|
95 |
|
---|
96 | for (i=0; i<ncoord; i++) {
|
---|
97 | sl = slat*pi/180.;
|
---|
98 | cm = cos(sl) / sqrt(1.0 - ex2 * (pow(sin(sl),2.)));
|
---|
99 | rho = sqrt(pow(x[i],2.) + pow(y[i],2.));
|
---|
100 | T = tan((pi / 4.0) - (sl / 2.0)) / pow(((1.0 - ex * sin(sl)) / (1.0 + ex * sin(sl))),(ex / 2.0));
|
---|
101 |
|
---|
102 | if (fabs(slat-90.) < 1.e-5)
|
---|
103 | T = rho * sqrt(pow((1. + ex),(1. + ex)) * pow((1. - ex),(1. - ex))) / 2. / re;
|
---|
104 | else
|
---|
105 | T = rho * T / (re * cm);
|
---|
106 |
|
---|
107 | chi = (pi / 2.0) - 2.0 * atan(T);
|
---|
108 | lat[i] = chi + ((ex2 / 2.0) + (5.0 * pow(ex2,2.0) / 24.0) + (pow(ex2,3.0) / 12.0)) *
|
---|
109 | sin(2.0 * chi) + ((7.0 * pow(ex2,2.0) / 48.0) + (29.0 * pow(ex2,3.0) / 240.0)) *
|
---|
110 | sin(4.0 * chi) + (7.0 * pow(ex2,3.0) / 120.0) * sin(6.0 * chi) ;
|
---|
111 |
|
---|
112 | lat[i] = (double)sgn * lat[i];
|
---|
113 | lon[i] = atan2((double)sgn * x[i],-(double)sgn * y[i]);
|
---|
114 | lon[i] = (double)sgn * lon[i];
|
---|
115 |
|
---|
116 | if (rho <= 0.1) {
|
---|
117 | lat[i] = 90. * (double)sgn;
|
---|
118 | lon[i] = 0.0;
|
---|
119 | iret=1;
|
---|
120 | }
|
---|
121 |
|
---|
122 | lon[i] = lon[i] * 180. / pi;
|
---|
123 | lat[i] = lat[i] * 180. / pi;
|
---|
124 | lon[i] = lon[i] - delta;
|
---|
125 | }
|
---|
126 |
|
---|
127 | return(iret);
|
---|
128 | }
|
---|
129 |
|
---|