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