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