00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080 #include <math.h>
00081 #include "Sinusoidal.h"
00082 #include "MapProjection3Parameters.h"
00083 #include "MapProjectionCoordinates.h"
00084 #include "GeodeticCoordinates.h"
00085 #include "CoordinateConversionException.h"
00086 #include "ErrorMessages.h"
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 using namespace MSP::CCS;
00099
00100
00101
00102
00103
00104
00105
00106 const double PI = 3.14159265358979323e0;
00107 const double PI_OVER_2 = (PI / 2.0);
00108 const double TWO_PI = (2.0 * PI);
00109
00110
00111 double sinuCoeffTimesSine( double coeff, double x, double latit )
00112 {
00113 return coeff * sin(x * latit);
00114 }
00115
00116
00117 double floatEq( double x, double v, double epsilon )
00118 {
00119 return ((v - epsilon) < x) && (x < (v + epsilon));
00120 }
00121
00122
00123
00124
00125
00126
00127
00128 Sinusoidal::Sinusoidal( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
00129 CoordinateSystem(),
00130 es2( 0.0066943799901413800 ),
00131 es4( 4.4814723452405e-005 ),
00132 es6( 3.0000678794350e-007 ),
00133 c0( .99832429845280 ),
00134 c1( .0025146070605187 ),
00135 c2( 2.6390465943377e-006 ),
00136 c3( 3.4180460865959e-009 ),
00137 a0( .0025188265843907 ),
00138 a1( 3.7009490356205e-006 ),
00139 a2( 7.4478137675038e-009 ),
00140 a3( 1.7035993238596e-011 ),
00141 Sinu_Origin_Long( 0.0 ),
00142 Sinu_False_Easting( 0.0 ),
00143 Sinu_False_Northing( 0.0 ),
00144 Sinu_Max_Easting( 20037509.0 ),
00145 Sinu_Min_Easting( -20037509.0 ),
00146 Sinu_Delta_Northing( 10001966.0 )
00147 {
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165 double j;
00166 double One_MINUS_es2, Sqrt_One_MINUS_es2, e1, e2, e3, e4;
00167 double inv_f = 1 / ellipsoidFlattening;
00168 char errorStatus[500] = "";
00169
00170 if (ellipsoidSemiMajorAxis <= 0.0)
00171 {
00172 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00173 }
00174 if ((inv_f < 250) || (inv_f > 350))
00175 {
00176 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00177 }
00178 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00179 {
00180 strcat( errorStatus, ErrorMessages::centralMeridian );
00181 }
00182
00183 if( strlen( errorStatus ) > 0)
00184 throw CoordinateConversionException( errorStatus );
00185
00186 semiMajorAxis = ellipsoidSemiMajorAxis;
00187 flattening = ellipsoidFlattening;
00188
00189 es2 = 2 * flattening - flattening * flattening;
00190 es4 = es2 * es2;
00191 es6 = es4 * es2;
00192 j = 45.0 * es6 / 1024.0;
00193 c0 = 1.0 - es2 / 4.0 - 3.0 * es4 / 64.0 - 5.0 * es6 / 256.0;
00194 c1 = 3.0 * es2 / 8.0 + 3.0 * es4 / 32.0 + j;
00195 c2 = 15.0 * es4 / 256.0 + j;
00196 c3 = 35.0 * es6 / 3072.0;
00197 One_MINUS_es2 = 1.0 - es2;
00198 Sqrt_One_MINUS_es2 = sqrt(One_MINUS_es2);
00199 e1 = (1.0 - Sqrt_One_MINUS_es2) / (1.0 + Sqrt_One_MINUS_es2);
00200 e2 = e1 * e1;
00201 e3 = e2 * e1;
00202 e4 = e3 * e1;
00203 a0 = 3.0 * e1 / 2.0 - 27.0 * e3 / 32.0 ;
00204 a1 = 21.0 * e2 / 16.0 - 55.0 * e4 / 32.0;
00205 a2 = 151.0 * e3 / 96.0;
00206 a3 = 1097.0 * e4 / 512.0;
00207 if (centralMeridian > PI)
00208 centralMeridian -= TWO_PI;
00209 Sinu_Origin_Long = centralMeridian;
00210 Sinu_False_Northing = falseNorthing;
00211 Sinu_False_Easting = falseEasting;
00212
00213 if (Sinu_Origin_Long > 0)
00214 {
00215 Sinu_Max_Easting = 19926189.0;
00216 Sinu_Min_Easting = -20037509.0;
00217 }
00218 else if (Sinu_Origin_Long < 0)
00219 {
00220 Sinu_Max_Easting = 20037509.0;
00221 Sinu_Min_Easting = -19926189.0;
00222 }
00223 else
00224 {
00225 Sinu_Max_Easting = 20037509.0;
00226 Sinu_Min_Easting = -20037509.0;
00227 }
00228 }
00229
00230
00231 Sinusoidal::Sinusoidal( const Sinusoidal &s )
00232 {
00233 semiMajorAxis = s.semiMajorAxis;
00234 flattening = s.flattening;
00235 es2 = s.es2;
00236 es4 = s.es4;
00237 es6 = s.es6;
00238 c0 = s.c0;
00239 c1 = s.c1;
00240 c2 = s.c2;
00241 c3 = s.c3;
00242 a0 = s.a0;
00243 a1 = s.a1;
00244 a2 = s.a2;
00245 a3 = s.a3;
00246 Sinu_Origin_Long = s.Sinu_Origin_Long;
00247 Sinu_False_Easting = s.Sinu_False_Easting;
00248 Sinu_False_Northing = s.Sinu_False_Northing;
00249 Sinu_Max_Easting = s.Sinu_Max_Easting;
00250 Sinu_Min_Easting = s.Sinu_Min_Easting;
00251 Sinu_Delta_Northing = s.Sinu_Delta_Northing;
00252 }
00253
00254
00255 Sinusoidal::~Sinusoidal()
00256 {
00257 }
00258
00259
00260 Sinusoidal& Sinusoidal::operator=( const Sinusoidal &s )
00261 {
00262 if( this != &s )
00263 {
00264 semiMajorAxis = s.semiMajorAxis;
00265 flattening = s.flattening;
00266 es2 = s.es2;
00267 es4 = s.es4;
00268 es6 = s.es6;
00269 c0 = s.c0;
00270 c1 = s.c1;
00271 c2 = s.c2;
00272 c3 = s.c3;
00273 a0 = s.a0;
00274 a1 = s.a1;
00275 a2 = s.a2;
00276 a3 = s.a3;
00277 Sinu_Origin_Long = s.Sinu_Origin_Long;
00278 Sinu_False_Easting = s.Sinu_False_Easting;
00279 Sinu_False_Northing = s.Sinu_False_Northing;
00280 Sinu_Max_Easting = s.Sinu_Max_Easting;
00281 Sinu_Min_Easting = s.Sinu_Min_Easting;
00282 Sinu_Delta_Northing = s.Sinu_Delta_Northing;
00283 }
00284
00285 return *this;
00286 }
00287
00288
00289 MapProjection3Parameters* Sinusoidal::getParameters() const
00290 {
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305 return new MapProjection3Parameters( CoordinateType::sinusoidal, Sinu_Origin_Long, Sinu_False_Easting, Sinu_False_Northing );
00306 }
00307
00308
00309 MSP::CCS::MapProjectionCoordinates* Sinusoidal::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00310 {
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324 double sin2lat, sin4lat, sin6lat;
00325 double dlam;
00326 double mm;
00327 double MM;
00328 char errorStatus[50] = "";
00329
00330 double longitude = geodeticCoordinates->longitude();
00331 double latitude = geodeticCoordinates->latitude();
00332 double slat = sin(latitude);
00333
00334 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00335 {
00336 strcat( errorStatus, ErrorMessages::latitude );
00337 }
00338 if ((longitude < -PI) || (longitude > TWO_PI))
00339 {
00340 strcat( errorStatus, ErrorMessages::longitude );
00341 }
00342
00343 if( strlen( errorStatus ) > 0)
00344 throw CoordinateConversionException( errorStatus );
00345
00346 dlam = longitude - Sinu_Origin_Long;
00347 if (dlam > PI)
00348 {
00349 dlam -= TWO_PI;
00350 }
00351 if (dlam < -PI)
00352 {
00353 dlam += TWO_PI;
00354 }
00355 mm = sqrt(1.0 - es2 * slat * slat);
00356
00357 sin2lat = sinuCoeffTimesSine(c1, 2.0, latitude);
00358 sin4lat = sinuCoeffTimesSine(c2, 4.0, latitude);
00359 sin6lat = sinuCoeffTimesSine(c3, 6.0, latitude);
00360 MM = semiMajorAxis * (c0 * latitude - sin2lat + sin4lat - sin6lat);
00361
00362 double easting = semiMajorAxis * dlam * cos(latitude) / mm + Sinu_False_Easting;
00363 double northing = MM + Sinu_False_Northing;
00364
00365 return new MapProjectionCoordinates( CoordinateType::sinusoidal, easting, northing );
00366 }
00367
00368
00369 MSP::CCS::GeodeticCoordinates* Sinusoidal::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00370 {
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384 double dx;
00385 double dy;
00386 double mu;
00387 double sin2mu, sin4mu, sin6mu, sin8mu;
00388 double sin_lat;
00389 double longitude, latitude;
00390 char errorStatus[50] = "";
00391
00392 double easting = mapProjectionCoordinates->easting();
00393 double northing = mapProjectionCoordinates->northing();
00394
00395 if ((easting < (Sinu_False_Easting + Sinu_Min_Easting))
00396 || (easting > (Sinu_False_Easting + Sinu_Max_Easting)))
00397 {
00398 strcat( errorStatus, ErrorMessages::easting );
00399 }
00400 if ((northing < (Sinu_False_Northing - Sinu_Delta_Northing))
00401 || (northing > (Sinu_False_Northing + Sinu_Delta_Northing)))
00402 {
00403 strcat( errorStatus, ErrorMessages::northing );
00404 }
00405
00406 if( strlen( errorStatus ) > 0)
00407 throw CoordinateConversionException( errorStatus );
00408
00409 dy = northing - Sinu_False_Northing;
00410 dx = easting - Sinu_False_Easting;
00411
00412 mu = dy / (semiMajorAxis * c0);
00413 sin2mu = sinuCoeffTimesSine(a0, 2.0, mu);
00414 sin4mu = sinuCoeffTimesSine(a1, 4.0, mu);
00415 sin6mu = sinuCoeffTimesSine(a2, 6.0, mu);
00416 sin8mu = sinuCoeffTimesSine(a3, 8.0, mu);
00417 latitude = mu + sin2mu + sin4mu + sin6mu + sin8mu;
00418
00419 if (latitude > PI_OVER_2)
00420 latitude = PI_OVER_2;
00421 else if (latitude < -PI_OVER_2)
00422 latitude = -PI_OVER_2;
00423
00424 if (floatEq(fabs(latitude),PI_OVER_2,1.0e-8))
00425 longitude = Sinu_Origin_Long;
00426 else
00427 {
00428 sin_lat = sin(latitude);
00429 longitude = Sinu_Origin_Long + dx * sqrt(1.0 - es2 *
00430 sin_lat * sin_lat) / (semiMajorAxis * cos(latitude));
00431
00432
00433 if (longitude > PI)
00434 longitude -= TWO_PI;
00435 if (longitude < -PI)
00436 longitude += TWO_PI;
00437
00438 if (longitude > PI)
00439 longitude = PI;
00440 else if (longitude < -PI)
00441 longitude = -PI;
00442 }
00443
00444 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00445 }
00446
00447
00448
00449