source: issm/trunk/src/c/modules/Xy2llx/Xy2llx.cpp@ 8698

Last change on this file since 8698 was 8698, checked in by Mathieu Morlighem, 14 years ago

some cosmetics

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