[62] | 1 | // This file is part of the PACPUS framework distributed under the
|
---|
| 2 | // CECILL-C License, Version 1.0.
|
---|
| 3 | //
|
---|
| 4 | /// @author Firstname Surname <firstname.surname@utc.fr>
|
---|
| 5 | /// @date Month, Year
|
---|
| 6 | /// @version $Id: transf_lamb_93.h 64 2013-01-09 16:41:12Z kurdejma $
|
---|
| 7 | /// @copyright Copyright (c) UTC/CNRS Heudiasyc 2006 - 2013. All rights reserved.
|
---|
| 8 | /// @brief Brief description.
|
---|
| 9 | ///
|
---|
| 10 | /// Detailed description.
|
---|
| 11 |
|
---|
[3] | 12 | #include <math.h>
|
---|
| 13 | #define pi 3.141592653589793
|
---|
| 14 | double atanh(double z)
|
---|
| 15 | {
|
---|
| 16 | return (0.5 * log((z+1)/(1-z)));
|
---|
| 17 | }
|
---|
| 18 |
|
---|
| 19 | double abs_me(double x)
|
---|
| 20 | {
|
---|
| 21 | if(x<0) x=-x;
|
---|
| 22 | return x;
|
---|
| 23 | }
|
---|
| 24 |
|
---|
| 25 | //transformation des longitudes tatitudes en Lambert93
|
---|
| 26 | void lonlattolam(double lon, double lat, double & lam93x, double & lam93y)
|
---|
| 27 | {
|
---|
| 28 | double GRS_a = 6378137;
|
---|
| 29 | double GRS_f = 1./298.257222101;
|
---|
| 30 |
|
---|
| 31 | double GRS_b = GRS_a*(1-GRS_f);
|
---|
| 32 | double GRS_bb= GRS_b*GRS_b;
|
---|
| 33 | double GRS_aa= 40680631590769.;
|
---|
| 34 | double GRS_e = sqrt((GRS_aa - GRS_bb) / (GRS_aa));
|
---|
| 35 |
|
---|
| 36 | double n = 0.725607765053267;
|
---|
| 37 | double C = 11754255.4261;
|
---|
| 38 | double XS = 700000;
|
---|
| 39 | double YS = 12655612.0499;
|
---|
| 40 |
|
---|
| 41 | double latiso;
|
---|
| 42 | latiso = atanh(sin(lat)) - GRS_e*atanh(GRS_e*sin(lat));
|
---|
| 43 | double gamma;
|
---|
| 44 | gamma = (lon - 0.0523598775598299)*n;
|
---|
| 45 | double R;
|
---|
| 46 | R = C * exp(-n*latiso);
|
---|
| 47 |
|
---|
| 48 | lam93x = R *sin(gamma)+XS;
|
---|
| 49 | lam93y = -R *cos(gamma)+YS;
|
---|
| 50 | //printf("x: %.20lf\n y: %.20lf\n",lam93x,lam93y);
|
---|
| 51 |
|
---|
| 52 | }
|
---|
| 53 |
|
---|
| 54 |
|
---|
| 55 |
|
---|
| 56 | //Transformation des coordonnees de Lambert93 en Longitude lattitude
|
---|
| 57 |
|
---|
| 58 | void lamtolonlat(double lamx, double lamy, double & lon, double & lat)
|
---|
| 59 | {
|
---|
| 60 |
|
---|
| 61 | double GRS_a = 6378137;
|
---|
| 62 | double GRS_f = 1./298.257222101;
|
---|
| 63 | double GRS_b = GRS_a*(1-GRS_f);
|
---|
| 64 | double GRS_bb= GRS_b*GRS_b;
|
---|
| 65 | double GRS_aa= 40680631590769.;
|
---|
| 66 | double GRS_e = sqrt((GRS_aa - GRS_bb) / (GRS_aa));
|
---|
| 67 |
|
---|
| 68 | //double n = 0.725607765053267;
|
---|
| 69 | //double C = 11754255.4261;
|
---|
| 70 | //double XS = 700000;
|
---|
| 71 | //double YS = 12655612.0499;
|
---|
| 72 |
|
---|
| 73 |
|
---|
| 74 | //lamx = lamx-700000;
|
---|
| 75 | //lamy = lamy-12655612.0499;
|
---|
| 76 |
|
---|
| 77 | double gamma;
|
---|
| 78 | gamma = atan(-(lamx-700000)/(lamy-12655612.0499));
|
---|
| 79 |
|
---|
| 80 |
|
---|
| 81 | lon = gamma/0.725607765053267 + 0.0523598775598299;
|
---|
| 82 | double R;
|
---|
| 83 | R = sqrt((lamx-700000) * (lamx-700000) + (lamy-12655612.0499) * (lamy-12655612.0499));
|
---|
| 84 |
|
---|
| 85 | double latiso;
|
---|
| 86 | latiso = log((11754255.4261)/R)/(0.725607765053267);
|
---|
| 87 |
|
---|
| 88 | double phiNew, phiOld;
|
---|
| 89 | phiOld =1;
|
---|
| 90 | phiNew= asin (tanh ( latiso + GRS_e * atanh(GRS_e * sin(phiOld))));
|
---|
| 91 | //printf("\nphiNew: %.20lf",phiNew);
|
---|
| 92 | while (abs_me(phiOld-phiNew) > 1e-10)
|
---|
| 93 | {
|
---|
| 94 |
|
---|
| 95 | if(abs_me(phiOld-phiNew) > 1e-10)
|
---|
| 96 | {
|
---|
| 97 |
|
---|
| 98 | phiOld = phiNew;
|
---|
| 99 | phiNew = asin(tanh(latiso+GRS_e*atanh(GRS_e*sin(phiOld))));
|
---|
| 100 | }
|
---|
| 101 | else
|
---|
| 102 | phiOld = phiNew;
|
---|
| 103 | }
|
---|
| 104 |
|
---|
| 105 | lat = phiNew;
|
---|
| 106 |
|
---|
| 107 | }
|
---|