Changeset 8698
- Timestamp:
- 06/23/11 11:01:20 (14 years ago)
- Location:
- issm/trunk/src/c/modules
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/modules/Ll2xyx/Ll2xyx.cpp
r8689 r8698 3 3 4 4 #include "./Ll2xyx.h" 5 #include "../../include/include.h" 5 6 #include "../../shared/shared.h" 6 #include "../../include/include.h" 7 #include "../../toolkits/toolkits.h" 8 #include "../../EnumDefinitions/EnumDefinitions.h" 7 #include <math.h> 9 8 10 int Ll2xyx(double* x, double* y, 11 double* lat, double* lon, int ncoord, 12 int sgn){ 13 9 int Ll2xyx(double* x, double* y, double* lat, double* lon, int ncoord, int sgn){ 14 10 /* This is a cpp conversion of the following: 15 11 %LL2XY - converts lat long to polar stereographic … … 26 22 % -1 : south latitude (default is mer=0 lat=71) 27 23 */ 28 29 24 double delta,slat; 30 25 bool flag=true; 31 26 32 27 /* Get central_meridian and standard_parallel depending on hemisphere */ 33 if 34 delta 28 if (sgn == 1) { 29 delta= 45; 35 30 slat = 70; 36 31 _printf_(flag,"Info: creating coordinates in polar stereographic (Std Latitude: 70N Meridian: 45).\n"); 37 32 } 38 33 else if (sgn == -1) { 39 delta 34 delta= 0; 40 35 slat = 71; 41 36 _printf_(flag,"Info: creating coordinates in polar stereographic (Std Latitude: 71S Meridian: 0).\n"); 42 37 } 43 else 44 _error_("Sign should be either +1 or -1.\n"); 38 else _error_("Sign should be either +1 or -1.\n"); 45 39 46 return(Ll2xyx(x,y, 47 lat,lon,ncoord, 48 sgn,delta,slat)); 40 return(Ll2xyx(x,y, lat,lon,ncoord, sgn,delta,slat)); 49 41 } 50 42 51 int Ll2xyx(double* x, double* y, 52 double* lat, double* lon, int ncoord, 53 int sgn, double central_meridian, double standard_parallel){ 54 43 int Ll2xyx(double* x, double* y, double* lat, double* lon, int ncoord, int sgn, double central_meridian, double standard_parallel){ 55 44 /* This is a cpp conversion of the following: 56 45 %LL2XY - converts lat long to polar stereographic … … 74 63 double T,rho,sl,tc,mc; 75 64 65 <<<<<<< .mine 66 if((sgn!=1) && (sgn!=-1)) _error_("Sign should be either +1 or -1"); 67 ======= 76 68 if ((sgn != 1) && (sgn != -1)) 77 69 _error_("Sign should be either +1 or -1.\n"); 70 >>>>>>> .r8690 78 71 79 72 delta = central_meridian; … … 81 74 82 75 /* Conversion constant from degrees to radians */ 83 cde 76 cde = 57.29577951; 84 77 /* Radius of the earth in meters */ 85 re 78 re = 6378.273*pow(10.,3.); 86 79 /* Eccentricity of the Hughes ellipsoid squared */ 87 ex2 =.006693883;80 ex2 = 0.006693883; 88 81 /* Eccentricity of the Hughes ellipsoid */ 89 ex = sqrt(ex2); 90 /* pi is not defined in cpp */ 91 pi = 4.*atan(1.); 82 ex = sqrt(ex2); 92 83 93 84 /* loop over all the coordinate pairs */ 94 95 for (i=0; i<ncoord; i++) { 96 latitude = fabs(lat[i]) * pi/180.; 97 longitude = (lon[i] + delta) * pi/180.; 85 for(i=0; i<ncoord; i++){ 86 latitude = fabs(lat[i]) * PI/180.; 87 longitude = (lon[i] + delta) * PI/180.; 98 88 99 89 /* compute X and Y in grid coordinates. */ 100 T = tan( pi/4.-latitude/2.) / pow(((1.-ex*sin(latitude))/(1.+ex*sin(latitude))),(ex/2.));90 T = tan(PI/4.-latitude/2.) / pow(((1.-ex*sin(latitude))/(1.+ex*sin(latitude))),(ex/2.)); 101 91 102 if ((90. - slat) < 92 if ((90. - slat) < 1.e-5) 103 93 rho = 2.*re*T/sqrt(pow((1.+ex),(1.+ex))*pow((1.-ex),(1.-ex))); 104 94 else { 105 sl = slat* pi/180.;106 tc = tan( pi/4.-sl/2.)/pow(((1.-ex*sin(sl))/(1.+ex*sin(sl))),(ex/2.));95 sl = slat*PI/180.; 96 tc = tan(PI/4.-sl/2.)/pow(((1.-ex*sin(sl))/(1.+ex*sin(sl))),(ex/2.)); 107 97 mc = cos(sl)/sqrt(1.0-ex2*(pow(sin(sl),2.))); 108 98 rho = re*mc*T/tc; 109 99 } 110 100 111 y[i] = -rho * (double)sgn * cos((double)sgn*longitude);112 x[i] = rho * (double)sgn * sin((double)sgn*longitude);101 y[i]= -rho*(double)sgn*cos(sgn*longitude); 102 x[i]= rho*(double)sgn*sin(sgn*longitude); 113 103 114 if (latitude >= pi / 2.){104 if (latitude>= PI/2.){ 115 105 x[i] = 0.0; 116 106 y[i] = 0.0; … … 118 108 } 119 109 } 120 121 110 return(iret); 122 111 } 123 -
issm/trunk/src/c/modules/Ll2xyx/Ll2xyx.h
r8673 r8698 7 7 8 8 /* local prototypes: */ 9 int Ll2xyx(double* x, double* y, 10 double* lat, double* lon, int ncoord, 11 int sgn); 12 13 int Ll2xyx(double* x, double* y, 14 double* lat, double* lon, int ncoord, 15 int sgn, double central_meridian, double standard_parallel); 9 int Ll2xyx(double* x, double* y, double* lat, double* lon, int ncoord, int sgn); 10 int Ll2xyx(double* x, double* y, double* lat, double* lon, int ncoord, int sgn, double central_meridian, double standard_parallel); 16 11 17 12 #endif /* _LL2XYX_H */ 18 -
issm/trunk/src/c/modules/Xy2llx/Xy2llx.cpp
r8689 r8698 3 3 4 4 #include "./Xy2llx.h" 5 #include "../../include/include.h" 5 6 #include "../../shared/shared.h" 6 #include "../../include/include.h" 7 #include "../../toolkits/toolkits.h" 8 #include "../../EnumDefinitions/EnumDefinitions.h" 7 #include <math.h> 9 8 10 int Xy2llx(double* lat, double* lon, 11 double* x, double* y, int ncoord, 12 int sgn){ 13 9 int Xy2llx(double* lat, double* lon, double* x, double* y, int ncoord, int sgn){ 14 10 /* This is a cpp conversion of the following: 15 11 %XY2LL - converts xy to lat long … … 32 28 33 29 /* Get central_meridian and standard_parallel depending on hemisphere */ 34 if (sgn ==1) {35 delta 30 if (sgn == 1) { 31 delta= 45; 36 32 slat = 70; 37 33 _printf_(flag,"Warning: expecting coordinates in polar stereographic (Std Latitude: 70N Meridian: 45).\n"); 38 34 } 39 35 else if (sgn == -1) { 40 delta 36 delta= 0; 41 37 slat = 71; 42 38 _printf_(flag,"Warning: expecting coordinates in polar stereographic (Std Latitude: 71S Meridian: 0).\n"); 43 39 } 44 else 45 _error_("Sign should be either +1 or -1.\n"); 40 else _error_("Sign should be either +1 or -1.\n"); 46 41 47 return(Xy2llx(lat,lon, 48 x,y,ncoord, 49 sgn,delta,slat)); 42 return(Xy2llx(lat,lon,x,y,ncoord,sgn,delta,slat)); 50 43 } 51 44 52 int Xy2llx(double* lat, double* lon, 53 double* x, double* y, int ncoord, 54 int sgn, double central_meridian, double standard_parallel){ 55 45 int Xy2llx(double* lat, double* lon, double* x, double* y, int ncoord, int sgn, double central_meridian, double standard_parallel){ 56 46 /* This is a cpp conversion of the following: 57 47 %XY2LL - converts xy to lat long … … 75 65 double sl,rho,cm,T,chi; 76 66 77 if ((sgn != 1) && (sgn != -1)) 78 _error_("Sign should be either +1 or -1.\n"); 67 if((sgn!=1) && (sgn!=-1)) _error_("Sign should be either +1 or -1.\n"); 79 68 80 69 delta = central_meridian; … … 86 75 re = 6378.273*pow(10.,3.); 87 76 /* Eccentricity of the Hughes ellipsoid squared */ 88 ex2 =.006693883;77 ex2 = 0.006693883; 89 78 /* Eccentricity of the Hughes ellipsoid */ 90 ex = sqrt(ex2); 91 /* pi is not defined in cpp */ 92 pi = 4.*atan(1.); 79 ex = sqrt(ex2); 93 80 94 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)); 95 87 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)); 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); 101 92 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); 93 chi = (PI / 2.0) - 2.0 * atan(T); 108 94 lat[i] = chi + ((ex2 / 2.0) + (5.0 * pow(ex2,2.0) / 24.0) + (pow(ex2,3.0) / 12.0)) * 109 95 sin(2.0 * chi) + ((7.0 * pow(ex2,2.0) / 48.0) + (29.0 * pow(ex2,3.0) / 240.0)) * … … 114 100 lon[i] = (double)sgn * lon[i]; 115 101 116 if (rho <= 0.1){102 if(rho <= 0.1){ 117 103 lat[i] = 90. * (double)sgn; 118 104 lon[i] = 0.0; … … 127 113 return(iret); 128 114 } 129 -
issm/trunk/src/c/modules/Xy2llx/Xy2llx.h
r8673 r8698 7 7 8 8 /* local prototypes: */ 9 int Xy2llx(double* lat, double* lon, 10 double* x, double* y, int ncoord, 11 int sgn); 12 13 int Xy2llx(double* lat, double* lon, 14 double* x, double* y, int ncoord, 15 int sgn, double central_meridian, double standard_parallel); 9 int Xy2llx(double* lat, double* lon, double* x, double* y, int ncoord, int sgn); 10 int Xy2llx(double* lat, double* lon, double* x, double* y, int ncoord, int sgn, double central_meridian, double standard_parallel); 16 11 17 12 #endif /* _XY2LLX_H */ 18
Note:
See TracChangeset
for help on using the changeset viewer.