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
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 #include <math.h>
00092 #include "TransverseMercator.h"
00093 #include "MapProjection5Parameters.h"
00094 #include "MapProjectionCoordinates.h"
00095 #include "GeodeticCoordinates.h"
00096 #include "CoordinateConversionException.h"
00097 #include "ErrorMessages.h"
00098 #include "WarningMessages.h"
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109 using namespace MSP::CCS;
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 #define PI 3.14159265358979323e0
00123 #define PI_OVER_2 (PI/2.0e0)
00124 #define MAX_LAT ((PI * 89.99)/180.0)
00125 #define MAX_DELTA_LONG ((PI * 90)/180.0)
00126 #define MIN_SCALE_FACTOR 0.3
00127 #define MAX_SCALE_FACTOR 3.0
00128
00129
00130
00131
00132
00133
00134 TransverseMercator::TransverseMercator(
00135 double ellipsoidSemiMajorAxis,
00136 double ellipsoidFlattening,
00137 double centralMeridian,
00138 double latitudeOfTrueScale,
00139 double falseEasting,
00140 double falseNorthing,
00141 double scaleFactor ) :
00142 CoordinateSystem(),
00143 TranMerc_es( 0.0066943799901413800 ),
00144 TranMerc_ebs( 0.0067394967565869 ),
00145 TranMerc_Origin_Long( 0.0 ),
00146 TranMerc_Origin_Lat( 0.0 ),
00147 TranMerc_False_Easting( 0.0 ),
00148 TranMerc_False_Northing( 0.0 ),
00149 TranMerc_Scale_Factor( 1.0 ),
00150 TranMerc_ap( 6367449.1458008 ),
00151 TranMerc_bp( 16038.508696861 ),
00152 TranMerc_cp( 16.832613334334 ),
00153 TranMerc_dp( 0.021984404273757 ),
00154 TranMerc_ep( 3.1148371319283e-005 ),
00155 TranMerc_Delta_Easting( 40000000.0 ),
00156 TranMerc_Delta_Northing( 40000000.0 )
00157 {
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 double tn;
00176 double tn2;
00177 double tn3;
00178 double tn4;
00179 double tn5;
00180 double TranMerc_b;
00181 double inv_f = 1 / ellipsoidFlattening;
00182 char errorStatus[500];
00183 errorStatus[0] = '\0';
00184
00185 if (ellipsoidSemiMajorAxis <= 0.0)
00186 {
00187 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00188 }
00189 if ((inv_f < 250) || (inv_f > 350))
00190 {
00191 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00192 }
00193 if ((latitudeOfTrueScale < -PI_OVER_2) || (latitudeOfTrueScale > PI_OVER_2))
00194 {
00195 strcat( errorStatus, ErrorMessages::originLatitude );
00196 }
00197 if ((centralMeridian < -PI) || (centralMeridian > (2*PI)))
00198 {
00199 strcat( errorStatus, ErrorMessages::centralMeridian );
00200 }
00201 if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
00202 {
00203 strcat( errorStatus, ErrorMessages::scaleFactor );
00204 }
00205
00206 if( strlen( errorStatus ) > 0)
00207 throw CoordinateConversionException( errorStatus );
00208
00209 semiMajorAxis = ellipsoidSemiMajorAxis;
00210 flattening = ellipsoidFlattening;
00211
00212 TranMerc_Origin_Lat = latitudeOfTrueScale;
00213 if (centralMeridian > PI)
00214 centralMeridian -= (2*PI);
00215 TranMerc_Origin_Long = centralMeridian;
00216 TranMerc_False_Northing = falseNorthing;
00217 TranMerc_False_Easting = falseEasting;
00218 TranMerc_Scale_Factor = scaleFactor;
00219
00220
00221 TranMerc_es = 2 * flattening - flattening * flattening;
00222
00223 TranMerc_ebs = (1 / (1 - TranMerc_es)) - 1;
00224
00225 TranMerc_b = semiMajorAxis * (1 - flattening);
00226
00227 tn = (semiMajorAxis - TranMerc_b) / (semiMajorAxis + TranMerc_b);
00228 tn2 = tn * tn;
00229 tn3 = tn2 * tn;
00230 tn4 = tn3 * tn;
00231 tn5 = tn4 * tn;
00232
00233 TranMerc_ap = semiMajorAxis *
00234 (1.e0 - tn + 5.e0 * (tn2 - tn3)/4.e0 + 81.e0 * (tn4 - tn5)/64.e0 );
00235 TranMerc_bp = 3.e0 * semiMajorAxis *
00236 (tn - tn2 + 7.e0 * (tn3 - tn4) /8.e0 + 55.e0 * tn5/64.e0 ) / 2.e0;
00237 TranMerc_cp = 15.e0 * semiMajorAxis *
00238 (tn2 - tn3 + 3.e0 * (tn4 - tn5 )/4.e0) /16.0;
00239 TranMerc_dp = 35.e0 * semiMajorAxis *
00240 (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
00241 TranMerc_ep = 315.e0 * semiMajorAxis * (tn4 - tn5) / 512.e0;
00242
00243 MapProjectionCoordinates* deltaCoordinates = 0;
00244 try
00245 {
00246 GeodeticCoordinates gcDeltaNorthing(
00247 CoordinateType::geodetic,
00248 MAX_DELTA_LONG + TranMerc_Origin_Long, MAX_LAT );
00249 deltaCoordinates = convertFromGeodetic( &gcDeltaNorthing );
00250 TranMerc_Delta_Northing = deltaCoordinates->northing();
00251 TranMerc_Delta_Northing ++;
00252
00253 delete deltaCoordinates;
00254 deltaCoordinates = 0;
00255
00256 GeodeticCoordinates gcDeltaEasting(
00257 CoordinateType::geodetic, MAX_DELTA_LONG + TranMerc_Origin_Long, 0 );
00258 deltaCoordinates = convertFromGeodetic( &gcDeltaEasting );
00259 TranMerc_Delta_Easting = deltaCoordinates->easting();
00260 TranMerc_Delta_Easting++;
00261
00262 delete deltaCoordinates;
00263 }
00264 catch(CoordinateConversionException e)
00265 {
00266 delete deltaCoordinates;
00267
00268 TranMerc_Delta_Easting = 40000000.0;
00269 TranMerc_Delta_Northing = 40000000.0;
00270 }
00271 }
00272
00273
00274 TransverseMercator::TransverseMercator( const TransverseMercator &tm )
00275 {
00276 semiMajorAxis = tm.semiMajorAxis;
00277 flattening = tm.flattening;
00278 TranMerc_es = tm.TranMerc_es;
00279 TranMerc_ebs = tm.TranMerc_ebs;
00280 TranMerc_Origin_Long = tm.TranMerc_Origin_Long;
00281 TranMerc_Origin_Lat = tm.TranMerc_Origin_Lat;
00282 TranMerc_False_Easting = tm.TranMerc_False_Easting;
00283 TranMerc_False_Northing = tm.TranMerc_False_Northing;
00284 TranMerc_Scale_Factor = tm.TranMerc_Scale_Factor;
00285 TranMerc_ap = tm.TranMerc_ap;
00286 TranMerc_bp = tm.TranMerc_bp;
00287 TranMerc_cp = tm.TranMerc_cp;
00288 TranMerc_dp = tm.TranMerc_dp;
00289 TranMerc_ep = tm.TranMerc_ep;
00290 TranMerc_Delta_Easting = tm.TranMerc_Delta_Easting;
00291 TranMerc_Delta_Northing = tm.TranMerc_Delta_Northing;
00292 }
00293
00294
00295 TransverseMercator::~TransverseMercator()
00296 {
00297 }
00298
00299
00300 TransverseMercator& TransverseMercator::operator=( const TransverseMercator &tm )
00301 {
00302 if( this != &tm )
00303 {
00304 semiMajorAxis = tm.semiMajorAxis;
00305 flattening = tm.flattening;
00306 TranMerc_es = tm.TranMerc_es;
00307 TranMerc_ebs = tm.TranMerc_ebs;
00308 TranMerc_Origin_Long = tm.TranMerc_Origin_Long;
00309 TranMerc_Origin_Lat = tm.TranMerc_Origin_Lat;
00310 TranMerc_False_Easting = tm.TranMerc_False_Easting;
00311 TranMerc_False_Northing = tm.TranMerc_False_Northing;
00312 TranMerc_Scale_Factor = tm.TranMerc_Scale_Factor;
00313 TranMerc_ap = tm.TranMerc_ap;
00314 TranMerc_bp = tm.TranMerc_bp;
00315 TranMerc_cp = tm.TranMerc_cp;
00316 TranMerc_dp = tm.TranMerc_dp;
00317 TranMerc_ep = tm.TranMerc_ep;
00318 TranMerc_Delta_Easting = tm.TranMerc_Delta_Easting;
00319 TranMerc_Delta_Northing = tm.TranMerc_Delta_Northing;
00320 }
00321
00322 return *this;
00323 }
00324
00325
00326 MapProjection5Parameters* TransverseMercator::getParameters() const
00327 {
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343 return new MapProjection5Parameters(
00344 CoordinateType::transverseMercator,
00345 TranMerc_Origin_Long, TranMerc_Origin_Lat, TranMerc_Scale_Factor,
00346 TranMerc_False_Easting, TranMerc_False_Northing );
00347 }
00348
00349
00350 MSP::CCS::MapProjectionCoordinates* TransverseMercator::convertFromGeodetic(
00351 MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00352 {
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366 double c;
00367 double c2;
00368 double c3;
00369 double c5;
00370 double c7;
00371 double dlam;
00372 double eta;
00373 double eta2;
00374 double eta3;
00375 double eta4;
00376 double s;
00377 double sn;
00378 double t;
00379 double tan2;
00380 double tan3;
00381 double tan4;
00382 double tan5;
00383 double tan6;
00384 double t1;
00385 double t2;
00386 double t3;
00387 double t4;
00388 double t5;
00389 double t6;
00390 double t7;
00391 double t8;
00392 double t9;
00393 double tmd;
00394 double tmdo;
00395 double temp_Origin;
00396 double temp_Long;
00397 char errorStatus[256];
00398 errorStatus[0] = '\0';
00399
00400 double longitude = geodeticCoordinates->longitude();
00401 double latitude = geodeticCoordinates->latitude();
00402
00403 if ((latitude < -MAX_LAT) || (latitude > MAX_LAT))
00404 {
00405 strcat( errorStatus, ErrorMessages::latitude );
00406 }
00407 if (longitude > PI)
00408 longitude -= (2 * PI);
00409 if ((longitude < (TranMerc_Origin_Long - MAX_DELTA_LONG))
00410 || (longitude > (TranMerc_Origin_Long + MAX_DELTA_LONG)))
00411 {
00412 if (longitude < 0)
00413 temp_Long = longitude + 2 * PI;
00414 else
00415 temp_Long = longitude;
00416 if (TranMerc_Origin_Long < 0)
00417 temp_Origin = TranMerc_Origin_Long + 2 * PI;
00418 else
00419 temp_Origin = TranMerc_Origin_Long;
00420 if ((temp_Long < (temp_Origin - MAX_DELTA_LONG))
00421 || (temp_Long > (temp_Origin + MAX_DELTA_LONG)))
00422 strcat( errorStatus, ErrorMessages::longitude );
00423 }
00424
00425 if( strlen( errorStatus ) > 0)
00426 throw CoordinateConversionException( errorStatus );
00427
00428
00429
00430
00431 dlam = longitude - TranMerc_Origin_Long;
00432
00433 if (fabs(dlam) > (9.0 * PI / 180))
00434 {
00435 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00436 }
00437
00438 if (dlam > PI)
00439 dlam -= (2 * PI);
00440 if (dlam < -PI)
00441 dlam += (2 * PI);
00442 if (fabs(dlam) < 2.e-10)
00443 dlam = 0.0;
00444
00445 s = sin(latitude);
00446 c = cos(latitude);
00447 c2 = c * c;
00448 c3 = c2 * c;
00449 c5 = c3 * c2;
00450 c7 = c5 * c2;
00451 t = tan (latitude);
00452 tan2 = t * t;
00453 tan3 = tan2 * t;
00454 tan4 = tan3 * t;
00455 tan5 = tan4 * t;
00456 tan6 = tan5 * t;
00457 eta = TranMerc_ebs * c2;
00458 eta2 = eta * eta;
00459 eta3 = eta2 * eta;
00460 eta4 = eta3 * eta;
00461
00462
00463 sn = sphsn(s);
00464
00465
00466 tmd = sphtmd(latitude, s, c);
00467
00468
00469 tmdo = sphtmd(
00470 TranMerc_Origin_Lat, sin(TranMerc_Origin_Lat), cos(TranMerc_Origin_Lat) );
00471
00472
00473 t1 = (tmd - tmdo) * TranMerc_Scale_Factor;
00474 t2 = sn * s * c * TranMerc_Scale_Factor/ 2.e0;
00475 t3 = sn * s * c3 * TranMerc_Scale_Factor *
00476 (5.e0 - tan2 + 9.e0 * eta + 4.e0 * eta2) /24.e0;
00477
00478 t4 = sn * s * c5 * TranMerc_Scale_Factor *
00479 (61.e0 - 58.e0 * tan2
00480 + tan4 + 270.e0 * eta - 330.e0 * tan2 * eta + 445.e0 * eta2
00481 + 324.e0 * eta3 -680.e0 * tan2 * eta2 + 88.e0 * eta4
00482 -600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4) / 720.e0;
00483
00484 t5 = sn * s * c7 * TranMerc_Scale_Factor *
00485 (1385.e0 - 3111.e0 * tan2 + 543.e0 * tan4 - tan6) / 40320.e0;
00486
00487 double dlamSq = dlam * dlam;
00488
00489 double northing = TranMerc_False_Northing + t1 +
00490 (((t5 * dlamSq + t4) * dlamSq + t3) * dlamSq + t2) * dlamSq;
00491
00492
00493 t6 = sn * c * TranMerc_Scale_Factor;
00494 t7 = sn * c3 * TranMerc_Scale_Factor * (1.e0 - tan2 + eta ) /6.e0;
00495 t8 = sn * c5 * TranMerc_Scale_Factor * (
00496 5.e0 - 18.e0 * tan2 + tan4 + 14.e0 * eta - 58.e0 * tan2 * eta +
00497 13.e0 * eta2 + 4.e0 * eta3 - 64.e0 * tan2 * eta2 - 24.e0 * tan2 * eta3 )
00498 / 120.e0;
00499 t9 = sn * c7 * TranMerc_Scale_Factor *
00500 ( 61.e0 - 479.e0 * tan2 + 179.e0 * tan4 - tan6 ) /5040.e0;
00501
00502 double easting = TranMerc_False_Easting +
00503 (((t9 * dlamSq + t8) * dlamSq + t7) * dlamSq + t6) * dlam;
00504
00505 return new MapProjectionCoordinates(
00506 CoordinateType::transverseMercator, errorStatus, easting, northing );
00507 }
00508
00509
00510 MSP::CCS::GeodeticCoordinates* TransverseMercator::convertToGeodetic(
00511 MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00512 {
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526 double c;
00527 double de;
00528 double dlam;
00529 double eta;
00530 double eta2;
00531 double eta3;
00532 double eta4;
00533 double ftphi;
00534 int i;
00535 double s;
00536 double sn;
00537 double sr;
00538 double t;
00539 double tan2;
00540 double tan4;
00541 double t10;
00542 double t11;
00543 double t12;
00544 double t13;
00545 double t14;
00546 double t15;
00547 double t16;
00548 double t17;
00549 double tmd;
00550 double tmdo;
00551 char errorStatus[256];
00552 errorStatus[0] = '\0';
00553
00554 double easting = mapProjectionCoordinates->easting();
00555 double northing = mapProjectionCoordinates->northing();
00556
00557 if ((easting < (TranMerc_False_Easting - TranMerc_Delta_Easting))
00558 ||(easting > (TranMerc_False_Easting + TranMerc_Delta_Easting)))
00559 {
00560 strcat( errorStatus, ErrorMessages::easting );
00561 }
00562
00563 if ((northing < (TranMerc_False_Northing - TranMerc_Delta_Northing))
00564 || (northing > (TranMerc_False_Northing + TranMerc_Delta_Northing)))
00565 {
00566 strcat( errorStatus, ErrorMessages::northing );
00567 }
00568
00569 if( strlen( errorStatus ) > 0)
00570 throw CoordinateConversionException( errorStatus );
00571
00572
00573 tmdo = sphtmd(
00574 TranMerc_Origin_Lat, sin(TranMerc_Origin_Lat), cos(TranMerc_Origin_Lat) );
00575
00576
00577 tmd = tmdo + (northing - TranMerc_False_Northing) / TranMerc_Scale_Factor;
00578
00579
00580 sr = sphsr(0.e0);
00581 ftphi = tmd/sr;
00582
00583 for (i = 0; i < 5 ; i++)
00584 {
00585 s = sin(ftphi);
00586 c = cos(ftphi);
00587 t10 = sphtmd(ftphi, s, c);
00588 sr = sphsr( s);
00589 ftphi = ftphi + (tmd - t10) / sr;
00590 }
00591
00592
00593 s = sin(ftphi);
00594 c = cos(ftphi);
00595
00596
00597 sr = sphsr(s);
00598
00599
00600 sn = sphsn(s);
00601
00602
00603 t = s/c;
00604 tan2 = t * t;
00605 tan4 = tan2 * tan2;
00606 eta = TranMerc_ebs * c * c;
00607 eta2 = eta * eta;
00608 eta3 = eta2 * eta;
00609 eta4 = eta3 * eta;
00610
00611 double sn2 = sn * sn;
00612 double sn3 = sn2 * sn;
00613 double sn5 = sn2 * sn3;
00614 double sn7 = sn2 * sn5;
00615
00616 de = easting - TranMerc_False_Easting;
00617 if (fabs(de) < 0.0001)
00618 de = 0.0;
00619
00620
00621 double TranMerc_Scale_Factor2 =
00622 TranMerc_Scale_Factor * TranMerc_Scale_Factor;
00623 double TranMerc_Scale_Factor3 =
00624 TranMerc_Scale_Factor2 * TranMerc_Scale_Factor;
00625 double TranMerc_Scale_Factor4 =
00626 TranMerc_Scale_Factor2 * TranMerc_Scale_Factor2;
00627 double TranMerc_Scale_Factor5 =
00628 TranMerc_Scale_Factor3 * TranMerc_Scale_Factor2;
00629 double TranMerc_Scale_Factor6 =
00630 TranMerc_Scale_Factor4 * TranMerc_Scale_Factor2;
00631 double TranMerc_Scale_Factor7 =
00632 TranMerc_Scale_Factor6 * TranMerc_Scale_Factor;
00633 double TranMerc_Scale_Factor8 =
00634 TranMerc_Scale_Factor6 * TranMerc_Scale_Factor2;
00635
00636 t10 = t / (2.e0 * sr * sn * TranMerc_Scale_Factor2);
00637
00638 t11 = t * (5.e0 + 3.e0 * tan2 + eta - 4.e0 * eta2 - 9.e0 * tan2 * eta) /
00639 (24.e0 * sr * sn3 * TranMerc_Scale_Factor4);
00640
00641 t12 = t * (61.e0 + 90.e0 * tan2 + 46.e0 * eta + 45.E0 * tan4
00642 - 252.e0 * tan2 * eta - 3.e0 * eta2 + 100.e0
00643 * eta3 - 66.e0 * tan2 * eta2 - 90.e0 * tan4
00644 * eta + 88.e0 * eta4 + 225.e0 * tan4 * eta2
00645 + 84.e0 * tan2* eta3 - 192.e0 * tan2 * eta4)
00646 / ( 720.e0 * sr * sn5 * TranMerc_Scale_Factor6 );
00647
00648 double t6 = t * t * t * t * t * t;
00649 t13 = t * ( 1385.e0 + 3633.e0 * tan2 + 4095.e0 * tan4 + 1575.e0 * t6 ) /
00650 (40320.e0 * sr * sn7 * TranMerc_Scale_Factor8);
00651
00652 double de2 = de * de;
00653 double latitude = ftphi +
00654 (((t13 * de2 - t12) * de2 + t11) * de2 - t10) *de2;
00655
00656
00657 t14 = 1.e0 / (sn * c * TranMerc_Scale_Factor);
00658
00659 t15 = (1.e0 + 2.e0 * tan2 + eta) /
00660 (6.e0 * sn3 * c * TranMerc_Scale_Factor3);
00661
00662 t16 = (5.e0 + 6.e0 * eta + 28.e0 * tan2 - 3.e0 * eta2
00663 + 8.e0 * tan2 * eta + 24.e0 * tan4 - 4.e0
00664 * eta3 + 4.e0 * tan2 * eta2 + 24.e0 * tan2 * eta3) /
00665 (120.e0 * sn5 * c * TranMerc_Scale_Factor5);
00666
00667 t17 = (61.e0 + 662.e0 * tan2 + 1320.e0 * tan4 + 720.e0 * t6 ) /
00668 (5040.e0 * sn7 * c * TranMerc_Scale_Factor7);
00669
00670
00671 dlam = (((-t17 * de2 + t16) * de2 - t15) * de2 + t14) *de;
00672
00673
00674 double longitude = TranMerc_Origin_Long + dlam;
00675
00676 if(fabs(latitude) > (90.0 * PI / 180.0))
00677 throw CoordinateConversionException( ErrorMessages::northing );
00678
00679 if((longitude) > (PI))
00680 {
00681 longitude -= (2 * PI);
00682 if(fabs(longitude) > PI)
00683 throw CoordinateConversionException( ErrorMessages::easting );
00684 }
00685 else if((longitude) < (-PI))
00686 {
00687 longitude += (2 * PI);
00688 if(fabs(longitude) > PI)
00689 throw CoordinateConversionException( ErrorMessages::easting );
00690 }
00691
00692 if (fabs(dlam) > (9.0 * PI / 180))
00693 {
00694
00695
00696
00697 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00698 }
00699
00700 return new GeodeticCoordinates(
00701 CoordinateType::geodetic, errorStatus, longitude, latitude );
00702 }
00703
00704
00705
00706
00707
00708
00709
00710 double TransverseMercator::sphtmd(
00711 double latitude,
00712 double sinLat,
00713 double cosLat )
00714 {
00715
00716
00717 double sin2Lat = 2.0 * sinLat * cosLat;
00718 double cos2Lat = 2.0 * cosLat * cosLat - 1.0;
00719 double sin2LatCubed = sin2Lat * sin2Lat * sin2Lat;
00720 double sin4Lat = 2.0 * sin2Lat * cos2Lat;
00721 double sin6Lat = 3.0 * sin2Lat - 4.0 * sin2LatCubed;
00722 double sin8Lat = cos2Lat * ( 4.0 * sin2Lat - 8.0 * sin2LatCubed);
00723
00724 return TranMerc_ap * latitude
00725 - TranMerc_bp * sin2Lat + TranMerc_cp * sin4Lat
00726 - TranMerc_dp * sin6Lat + TranMerc_ep * sin8Lat;
00727
00728
00729
00730
00731 }
00732
00733
00734 double TransverseMercator::sphsn( double sinLat )
00735 {
00736
00737 return semiMajorAxis / sqrt( 1.e0 - TranMerc_es * sinLat * sinLat);
00738 }
00739
00740
00741 double TransverseMercator::sphsr( double sinLat )
00742 {
00743
00744 double denom = sqrt(1.e0 - TranMerc_es * sinLat * sinLat);
00745 return semiMajorAxis * (1.e0 - TranMerc_es) / denom / denom / denom;
00746 }
00747
00748
00749
00750