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 "LambertConformalConic.h"
00093 #include "MapProjection5Parameters.h"
00094 #include "MapProjection6Parameters.h"
00095 #include "MapProjectionCoordinates.h"
00096 #include "GeodeticCoordinates.h"
00097 #include "CoordinateConversionException.h"
00098 #include "ErrorMessages.h"
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 using namespace MSP::CCS;
00111
00112
00113
00114
00115
00116
00117
00118 const double PI = 3.14159265358979323e0;
00119 const double PI_OVER_2 = (PI / 2.0);
00120 const double PI_OVER_4 = (PI / 4.0);
00121 const double PI_OVER_180 = (PI / 180.0);
00122 const double MAX_LAT = (( PI * 89.99972222222222) / 180.0);
00123 const double TWO_PI = (2.0 * PI);
00124 const double MIN_SCALE_FACTOR = 1.0e-9;
00125 const double ONE_SECOND = 4.89e-6;
00126
00127
00128
00129
00130
00131
00132
00133 LambertConformalConic::LambertConformalConic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing, double scaleFactor ) :
00134 CoordinateSystem(),
00135 coordinateType( CoordinateType::lambertConformalConic1Parallel ),
00136 es( 0.08181919084262188000 ),
00137 es_OVER_2( .040909595421311 ),
00138 Lambert_1_n( 0.70710678118655 ),
00139 Lambert_1_rho0( 6388838.2901212 ),
00140 Lambert_1_rho_olat( 6388838.2901211 ),
00141 Lambert_1_t0( 0.41618115138974 ),
00142 Lambert_Origin_Long( 0.0 ),
00143 Lambert_Origin_Latitude( (45.0 * PI / 180.0) ),
00144 Lambert_False_Easting( 0.0 ),
00145 Lambert_False_Northing( 0.0 ),
00146 Lambert_Scale_Factor( 1.0 ),
00147 Lambert_2_Origin_Lat( (45 * PI / 180) ),
00148 Lambert_2_Std_Parallel_1( (40 * PI / 180) ),
00149 Lambert_2_Std_Parallel_2( (50 * PI / 180) ),
00150 Lambert_Delta_Easting( 400000000.0 ),
00151 Lambert_Delta_Northing( 400000000.0 )
00152 {
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 double es2;
00170 double inv_f = 1 / ellipsoidFlattening;
00171 char errorStatus[500] = "";
00172
00173 if (ellipsoidSemiMajorAxis <= 0.0)
00174 {
00175 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00176 }
00177 if ((inv_f < 250) || (inv_f > 350))
00178 {
00179 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00180 }
00181 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00182 {
00183 strcat( errorStatus, ErrorMessages::centralMeridian );
00184 }
00185
00186 if( strlen( errorStatus ) > 0)
00187 throw CoordinateConversionException( errorStatus );
00188
00189 semiMajorAxis = ellipsoidSemiMajorAxis;
00190 flattening = ellipsoidFlattening;
00191
00192 if (centralMeridian > PI)
00193 centralMeridian -= TWO_PI;
00194 Lambert_Origin_Long = centralMeridian;
00195 Lambert_False_Easting = falseEasting;
00196
00197 es2 = 2.0 * flattening - flattening * flattening;
00198 es = sqrt(es2);
00199 es_OVER_2 = es / 2.0;
00200
00201 CommonParameters* parameters = setCommonLambert1StandardParallelParameters(originLatitude, falseNorthing, scaleFactor);
00202
00203 Lambert_1_n = parameters->_lambertN;
00204 Lambert_1_rho0 = parameters->_lambertRho0;
00205 Lambert_1_rho_olat = parameters->_lambertRhoOlat;
00206 Lambert_1_t0 = parameters->_lambertT0;
00207 Lambert_Origin_Latitude = parameters->_lambertOriginLatitude;
00208 Lambert_False_Northing = parameters->_lambertFalseNorthing;
00209 Lambert_Scale_Factor = parameters->_lambertScaleFactor;
00210
00211 Lambert_2_Origin_Lat = parameters->_lambertOriginLatitude;
00212
00213 double sinOlat = sin(Lambert_Origin_Latitude);
00214 double esSinOlat = es * sinOlat;
00215 double w0 = sqrt(1 - es2 * sinOlat * sinOlat);
00216 double f0 = cos(Lambert_Origin_Latitude) / (w0 * pow(Lambert_1_t0, Lambert_1_n));
00217 double c = Lambert_Scale_Factor * f0;
00218
00219 double phi = 0.0;
00220 double tempPhi = 1.0;
00221 Lambert_2_Std_Parallel_1 = calculateLambert2StandardParallel(es2, phi, tempPhi, c);
00222
00223 phi = 89.0 * PI_OVER_180;
00224 tempPhi = 0.0;
00225 Lambert_2_Std_Parallel_2 = calculateLambert2StandardParallel(es2, phi, tempPhi, c);
00226
00227 delete parameters;
00228 parameters = 0;
00229 }
00230
00231
00232 LambertConformalConic::LambertConformalConic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double standardParallel1, double standardParallel2, double falseEasting, double falseNorthing ) :
00233 CoordinateSystem(),
00234 coordinateType( CoordinateType::lambertConformalConic2Parallels ),
00235 es( 0.081819190842621 ),
00236 es_OVER_2( 0.040909595421311 ),
00237 Lambert_1_n( 0.70710678118655 ),
00238 Lambert_1_rho0( 6388838.2901212 ),
00239 Lambert_1_rho_olat( 6388838.2901211 ),
00240 Lambert_1_t0( 0.41618115138974 ),
00241 Lambert_Origin_Long( 0.0 ),
00242 Lambert_Origin_Latitude( (45 * PI / 180) ),
00243 Lambert_False_Easting( 0.0 ),
00244 Lambert_False_Northing( 0.0 ),
00245 Lambert_Scale_Factor( 1.0 ),
00246 Lambert_2_Origin_Lat( (45 * PI / 180) ),
00247 Lambert_2_Std_Parallel_1( (40 * PI / 180) ),
00248 Lambert_2_Std_Parallel_2( (50 * PI / 180) ),
00249 Lambert_Delta_Easting( 400000000.0 ),
00250 Lambert_Delta_Northing( 400000000.0 )
00251 {
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272 double es2;
00273 double es_sin;
00274 double t0;
00275 double t1;
00276 double t2;
00277 double t_olat;
00278 double m0;
00279 double m1;
00280 double m2;
00281 double m_olat;
00282 double n;
00283 double const_value;
00284 double Lambert_lat0;
00285 double Lambert_k0;
00286 double Lambert_false_northing;
00287 double inv_f = 1 / ellipsoidFlattening;
00288 char errorStatus[500] = "";
00289
00290 if (ellipsoidSemiMajorAxis <= 0.0)
00291 {
00292 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00293 }
00294 if ((inv_f < 250) || (inv_f > 350))
00295 {
00296 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00297 }
00298 if ((originLatitude < -MAX_LAT) || (originLatitude > MAX_LAT))
00299 {
00300 strcat( errorStatus, ErrorMessages::originLatitude );
00301 }
00302 if ((standardParallel1 < -MAX_LAT) || (standardParallel1 > MAX_LAT))
00303 {
00304 strcat( errorStatus, ErrorMessages::standardParallel1 );
00305 }
00306 if ((standardParallel2 < -MAX_LAT) || (standardParallel2 > MAX_LAT))
00307 {
00308 strcat( errorStatus, ErrorMessages::standardParallel2 );
00309 }
00310 if ((standardParallel1 == 0) && (standardParallel2 == 0))
00311 {
00312 strcat( errorStatus, ErrorMessages::standardParallel1_2 );
00313 }
00314 if (standardParallel1 == -standardParallel2)
00315 {
00316 strcat( errorStatus, ErrorMessages::standardParallelHemisphere );
00317 }
00318 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00319 {
00320 strcat( errorStatus, ErrorMessages::centralMeridian );
00321 }
00322
00323 if( strlen( errorStatus ) > 0)
00324 throw CoordinateConversionException( errorStatus );
00325
00326 semiMajorAxis = ellipsoidSemiMajorAxis;
00327 flattening = ellipsoidFlattening;
00328
00329 Lambert_2_Origin_Lat = originLatitude;
00330 Lambert_2_Std_Parallel_1 = standardParallel1;
00331 Lambert_2_Std_Parallel_2 = standardParallel2;
00332
00333 if (centralMeridian > PI)
00334 centralMeridian -= TWO_PI;
00335 Lambert_Origin_Long = centralMeridian;
00336 Lambert_False_Easting = falseEasting;
00337
00338 es2 = 2 * flattening - flattening * flattening;
00339 es = sqrt(es2);
00340 es_OVER_2 = es / 2.0;
00341
00342 if (fabs(Lambert_2_Std_Parallel_1 - Lambert_2_Std_Parallel_2) > 1.0e-10)
00343 {
00344 es_sin = esSin(sin(originLatitude));
00345 m_olat = lambertM(cos(originLatitude), es_sin);
00346 t_olat = lambertT(originLatitude, es_sin);
00347
00348 es_sin = esSin(sin(Lambert_2_Std_Parallel_1));
00349 m1 = lambertM(cos(Lambert_2_Std_Parallel_1), es_sin);
00350 t1 = lambertT(Lambert_2_Std_Parallel_1, es_sin);
00351
00352 es_sin = esSin(sin(Lambert_2_Std_Parallel_2));
00353 m2 = lambertM(cos(Lambert_2_Std_Parallel_2), es_sin);
00354 t2 = lambertT(Lambert_2_Std_Parallel_2, es_sin);
00355
00356 n = log(m1 / m2) / log(t1 / t2);
00357
00358 Lambert_lat0 = asin(n);
00359
00360 es_sin = esSin(sin(Lambert_lat0));
00361 m0 = lambertM(cos(Lambert_lat0), es_sin);
00362 t0 = lambertT(Lambert_lat0, es_sin);
00363
00364 Lambert_k0 = (m1 / m0) * (pow(t0 / t1, n));
00365
00366 const_value = ((semiMajorAxis * m2) / (n * pow(t2, n)));
00367
00368 Lambert_false_northing = (const_value * pow(t_olat, n)) - (const_value * pow(t0, n)) + falseNorthing;
00369 }
00370 else
00371 {
00372 Lambert_lat0 = Lambert_2_Std_Parallel_1;
00373 Lambert_k0 = 1.0;
00374 Lambert_false_northing = falseNorthing;
00375 }
00376
00377 CommonParameters* parameters = setCommonLambert1StandardParallelParameters(Lambert_lat0, Lambert_false_northing, Lambert_k0);
00378
00379 Lambert_1_n = parameters->_lambertN;
00380 Lambert_1_rho0 = parameters->_lambertRho0;
00381 Lambert_1_rho_olat = parameters->_lambertRhoOlat;
00382 Lambert_1_t0 = parameters->_lambertT0;
00383 Lambert_Origin_Latitude = parameters->_lambertOriginLatitude;
00384 Lambert_False_Northing = parameters->_lambertFalseNorthing;
00385 Lambert_Scale_Factor = parameters->_lambertScaleFactor;
00386
00387 delete parameters;
00388 parameters = 0;
00389 }
00390
00391 LambertConformalConic::LambertConformalConic( const LambertConformalConic &lcc )
00392 {
00393 coordinateType = lcc.coordinateType;
00394 semiMajorAxis = lcc.semiMajorAxis;
00395 flattening = lcc.flattening;
00396 es = lcc.es;
00397 es_OVER_2 = lcc.es_OVER_2;
00398 Lambert_1_n = lcc.Lambert_1_n;
00399 Lambert_1_rho0 = lcc.Lambert_1_rho0;
00400 Lambert_1_rho_olat = lcc.Lambert_1_rho_olat;
00401 Lambert_1_t0 = lcc.Lambert_1_t0;
00402 Lambert_Origin_Long = lcc.Lambert_Origin_Long;
00403 Lambert_Origin_Latitude = lcc.Lambert_Origin_Latitude;
00404 Lambert_False_Easting = lcc.Lambert_False_Easting;
00405 Lambert_False_Northing = lcc.Lambert_False_Northing;
00406 Lambert_Scale_Factor = lcc.Lambert_Scale_Factor;
00407 Lambert_2_Origin_Lat = lcc.Lambert_2_Origin_Lat;
00408 Lambert_2_Std_Parallel_1 = lcc.Lambert_2_Std_Parallel_1;
00409 Lambert_2_Std_Parallel_2 = lcc.Lambert_2_Std_Parallel_2;
00410 Lambert_Delta_Easting = lcc.Lambert_Delta_Easting;
00411 Lambert_Delta_Northing = lcc.Lambert_Delta_Northing;
00412 }
00413
00414
00415 LambertConformalConic::~LambertConformalConic()
00416 {
00417 }
00418
00419
00420 LambertConformalConic& LambertConformalConic::operator=( const LambertConformalConic &lcc )
00421 {
00422 if( this != &lcc )
00423 {
00424 coordinateType = lcc.coordinateType;
00425 semiMajorAxis = lcc.semiMajorAxis;
00426 flattening = lcc.flattening;
00427 es = lcc.es;
00428 es_OVER_2 = lcc.es_OVER_2;
00429 Lambert_1_n = lcc.Lambert_1_n;
00430 Lambert_1_rho0 = lcc.Lambert_1_rho0;
00431 Lambert_1_rho_olat = lcc.Lambert_1_rho_olat;
00432 Lambert_1_t0 = lcc.Lambert_1_t0;
00433 Lambert_Origin_Long = lcc.Lambert_Origin_Long;
00434 Lambert_Origin_Latitude = lcc.Lambert_Origin_Latitude;
00435 Lambert_False_Easting = lcc.Lambert_False_Easting;
00436 Lambert_False_Northing = lcc.Lambert_False_Northing;
00437 Lambert_Scale_Factor = lcc.Lambert_Scale_Factor;
00438 Lambert_2_Origin_Lat = lcc.Lambert_2_Origin_Lat;
00439 Lambert_2_Std_Parallel_1 = lcc.Lambert_2_Std_Parallel_1;
00440 Lambert_2_Std_Parallel_2 = lcc.Lambert_2_Std_Parallel_2;
00441 Lambert_Delta_Easting = lcc.Lambert_Delta_Easting;
00442 Lambert_Delta_Northing = lcc.Lambert_Delta_Northing;
00443 }
00444
00445 return *this;
00446 }
00447
00448
00449 MapProjection5Parameters* LambertConformalConic::get1StandardParallelParameters() const
00450 {
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464 return new MapProjection5Parameters( CoordinateType::lambertConformalConic1Parallel, Lambert_Origin_Long, Lambert_Origin_Latitude, Lambert_Scale_Factor, Lambert_False_Easting, Lambert_False_Northing );
00465 }
00466
00467
00468 MapProjection6Parameters* LambertConformalConic::get2StandardParallelParameters() const
00469 {
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 return new MapProjection6Parameters( CoordinateType::lambertConformalConic2Parallels, Lambert_Origin_Long, Lambert_2_Origin_Lat, Lambert_2_Std_Parallel_1, Lambert_2_Std_Parallel_2, Lambert_False_Easting, Lambert_False_Northing );
00485 }
00486
00487
00488 MSP::CCS::MapProjectionCoordinates* LambertConformalConic::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00489 {
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503 double t;
00504 double rho;
00505 double dlam;
00506 double theta;
00507 char errorStatus[50] = "";
00508
00509 double longitude = geodeticCoordinates->longitude();
00510 double latitude = geodeticCoordinates->latitude();
00511
00512 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00513 {
00514 strcat( errorStatus, ErrorMessages::latitude );
00515 }
00516 if ((longitude < -PI) || (longitude > TWO_PI))
00517 {
00518 strcat( errorStatus, ErrorMessages::longitude );
00519 }
00520
00521 if( strlen( errorStatus ) > 0)
00522 throw CoordinateConversionException( errorStatus );
00523
00524 if (fabs(fabs(latitude) - PI_OVER_2) > 1.0e-10)
00525 {
00526 t = lambertT(latitude, esSin(sin(latitude)));
00527 rho = Lambert_1_rho0 * pow(t / Lambert_1_t0, Lambert_1_n);
00528 }
00529 else
00530 {
00531 if ((latitude * Lambert_1_n) <= 0)
00532 {
00533 throw CoordinateConversionException( ErrorMessages::latitude );
00534 }
00535 rho = 0.0;
00536 }
00537
00538 dlam = longitude - Lambert_Origin_Long;
00539
00540 if (dlam > PI)
00541 {
00542 dlam -= TWO_PI;
00543 }
00544 if (dlam < -PI)
00545 {
00546 dlam += TWO_PI;
00547 }
00548
00549 theta = Lambert_1_n * dlam;
00550
00551 double easting = rho * sin(theta) + Lambert_False_Easting;
00552 double northing = Lambert_1_rho_olat - rho * cos(theta) + Lambert_False_Northing;
00553
00554 return new MapProjectionCoordinates( coordinateType, easting, northing );
00555 }
00556
00557
00558 MSP::CCS::GeodeticCoordinates* LambertConformalConic::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00559 {
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573 double dx;
00574 double dy;
00575 double rho;
00576 double rho_olat_MINUS_dy;
00577 double t;
00578 double PHI;
00579 double es_sin;
00580 double tempPHI = 0.0;
00581 double theta = 0.0;
00582 double tolerance = 4.85e-10;
00583 int count = 30;
00584 double longitude, latitude;
00585 char errorStatus[50] = "";
00586
00587 double easting = mapProjectionCoordinates->easting();
00588 double northing = mapProjectionCoordinates->northing();
00589
00590 if ((easting < (Lambert_False_Easting - Lambert_Delta_Easting))
00591 ||(easting > (Lambert_False_Easting + Lambert_Delta_Easting)))
00592 {
00593 strcat( errorStatus, ErrorMessages::easting );
00594 }
00595 if ((northing < (Lambert_False_Northing - Lambert_Delta_Northing))
00596 || (northing > (Lambert_False_Northing + Lambert_Delta_Northing)))
00597 {
00598 strcat( errorStatus, ErrorMessages::northing );
00599 }
00600
00601 if( strlen( errorStatus ) > 0)
00602 throw CoordinateConversionException( errorStatus );
00603
00604 dy = northing - Lambert_False_Northing;
00605 dx = easting - Lambert_False_Easting;
00606 rho_olat_MINUS_dy = Lambert_1_rho_olat - dy;
00607 rho = sqrt(dx * dx + (rho_olat_MINUS_dy) * (rho_olat_MINUS_dy));
00608
00609 if (Lambert_1_n < 0.0)
00610 {
00611 rho *= -1.0;
00612 dx *= -1.0;
00613 rho_olat_MINUS_dy *= -1.0;
00614 }
00615
00616 if (rho != 0.0)
00617 {
00618 theta = atan2(dx, rho_olat_MINUS_dy) / Lambert_1_n;
00619 t = Lambert_1_t0 * pow(rho / Lambert_1_rho0, 1 / Lambert_1_n);
00620 PHI = PI_OVER_2 - 2.0 * atan(t);
00621 while (fabs(PHI - tempPHI) > tolerance && count)
00622 {
00623 tempPHI = PHI;
00624 es_sin = esSin(sin(PHI));
00625 PHI = PI_OVER_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin), es_OVER_2));
00626 count --;
00627 }
00628
00629 if(!count)
00630 throw CoordinateConversionException( ErrorMessages::northing );
00631
00632 latitude = PHI;
00633 longitude = theta + Lambert_Origin_Long;
00634
00635 if (fabs(latitude) < 2.0e-7)
00636 latitude = 0.0;
00637 if (latitude > PI_OVER_2)
00638 latitude = PI_OVER_2;
00639 else if (latitude < -PI_OVER_2)
00640 latitude = -PI_OVER_2;
00641
00642 if (longitude > PI)
00643 {
00644 if (longitude - PI < 3.5e-6)
00645 longitude = PI;
00646 else
00647 longitude -= TWO_PI;
00648 }
00649 if (longitude < -PI)
00650 {
00651 if (fabs(longitude + PI) < 3.5e-6)
00652 longitude = -PI;
00653 else
00654 longitude += TWO_PI;
00655 }
00656
00657 if (fabs(longitude) < 2.0e-7)
00658 longitude = 0.0;
00659 if (longitude > PI)
00660 longitude = PI;
00661 else if (longitude < -PI)
00662 longitude = -PI;
00663 }
00664 else
00665 {
00666 if (Lambert_1_n > 0.0)
00667 latitude = PI_OVER_2;
00668 else
00669 latitude = -PI_OVER_2;
00670 longitude = Lambert_Origin_Long;
00671 }
00672
00673 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00674 }
00675
00676
00677 LambertConformalConic::CommonParameters* LambertConformalConic::setCommonLambert1StandardParallelParameters( double originLatitude, double falseNorthing, double scaleFactor )
00678 {
00691 double _esSin;
00692 double m0;
00693 char errorStatus[500] = "";
00694
00695 if (((originLatitude < -MAX_LAT) || (originLatitude > MAX_LAT)) ||
00696 (originLatitude > -ONE_SECOND) && (originLatitude < ONE_SECOND))
00697 {
00698 strcat( errorStatus, ErrorMessages::originLatitude );
00699 }
00700 if (scaleFactor < MIN_SCALE_FACTOR)
00701 {
00702 strcat( errorStatus, ErrorMessages::scaleFactor );
00703 }
00704
00705 if(strlen( errorStatus ) > 0)
00706 throw CoordinateConversionException(errorStatus);
00707
00708 CommonParameters* parameters = new CommonParameters();
00709
00710 parameters->_lambertOriginLatitude = originLatitude;
00711 parameters->_lambertFalseNorthing = falseNorthing;
00712 parameters->_lambertScaleFactor = scaleFactor;
00713
00714 parameters->_lambertN = sin(originLatitude);
00715
00716 _esSin = esSin(sin(originLatitude));
00717
00718 m0 = cos(originLatitude) / sqrt(1.0 - _esSin * _esSin);
00719 parameters->_lambertT0 = lambertT(originLatitude, _esSin);
00720
00721 parameters->_lambertRho0 = semiMajorAxis * scaleFactor * m0 / parameters->_lambertN;
00722
00723 parameters->_lambertRhoOlat = parameters->_lambertRho0;
00724
00725 return parameters;
00726 }
00727
00728
00729 double LambertConformalConic::calculateLambert2StandardParallel(double es2, double phi, double tempPhi, double c)
00730 {
00731
00732
00733
00734
00735
00736 double tolerance = 1.0e-11;
00737 int count = 30;
00738 while (fabs(phi - tempPhi) > tolerance && count > 0)
00739 {
00740 tempPhi = phi;
00741 double sinPhi = sin(phi);
00742 double esSinPhi = es * sinPhi;
00743 double tPhi = lambertT(phi, esSinPhi);
00744 double wPhi = sqrt(1 - es2 * sinPhi * sinPhi);
00745 double fPhi = cos(phi) / (wPhi * pow(tPhi, Lambert_1_n));
00746 double fprPhi = ((Lambert_1_n - sinPhi) * (1 - es2)) / (pow(wPhi, 3) * pow(tPhi, Lambert_1_n));
00747
00748 phi = phi + (c - fPhi) / fprPhi;
00749
00750 count--;
00751 }
00752
00753 return phi;
00754 }
00755
00756
00757 double LambertConformalConic::lambertM( double clat, double essin )
00758 {
00759 return clat / sqrt(1.0 - essin * essin);
00760 }
00761
00762
00763 double LambertConformalConic::lambertT( double lat, double essin )
00764 {
00765 return tan(PI_OVER_4 - lat / 2) / pow((1.0 - essin) / (1.0 + essin), es_OVER_2);
00766 }
00767
00768
00769 double LambertConformalConic::esSin( double sinlat )
00770 {
00771 return es * sinlat;
00772 }
00773
00774
00775
00776