source: pacpusframework/trunk/include/Pacpus/PacpusTools/transf_lamb_93.h@ 17

Last change on this file since 17 was 3, checked in by sgosseli, 12 years ago
  • Add the existing Pacpus files from pacpusdev and pacpuscore.
  • Provide a clean build system based on multiple CMake files.
File size: 2.1 KB
Line 
1#include <math.h>
2#define pi 3.141592653589793
3double atanh(double z)
4{
5 return (0.5 * log((z+1)/(1-z)));
6}
7
8double abs_me(double x)
9{
10 if(x<0) x=-x;
11 return x;
12}
13
14//transformation des longitudes tatitudes en Lambert93
15void lonlattolam(double lon, double lat, double & lam93x, double & lam93y)
16 {
17 double GRS_a = 6378137;
18 double GRS_f = 1./298.257222101;
19
20 double GRS_b = GRS_a*(1-GRS_f);
21 double GRS_bb= GRS_b*GRS_b;
22 double GRS_aa= 40680631590769.;
23 double GRS_e = sqrt((GRS_aa - GRS_bb) / (GRS_aa));
24
25 double n = 0.725607765053267;
26 double C = 11754255.4261;
27 double XS = 700000;
28 double YS = 12655612.0499;
29
30 double latiso;
31 latiso = atanh(sin(lat)) - GRS_e*atanh(GRS_e*sin(lat));
32 double gamma;
33 gamma = (lon - 0.0523598775598299)*n;
34 double R;
35 R = C * exp(-n*latiso);
36
37 lam93x = R *sin(gamma)+XS;
38 lam93y = -R *cos(gamma)+YS;
39 //printf("x: %.20lf\n y: %.20lf\n",lam93x,lam93y);
40
41 }
42
43
44
45//Transformation des coordonnees de Lambert93 en Longitude lattitude
46
47void lamtolonlat(double lamx, double lamy, double & lon, double & lat)
48{
49
50 double GRS_a = 6378137;
51 double GRS_f = 1./298.257222101;
52 double GRS_b = GRS_a*(1-GRS_f);
53 double GRS_bb= GRS_b*GRS_b;
54 double GRS_aa= 40680631590769.;
55 double GRS_e = sqrt((GRS_aa - GRS_bb) / (GRS_aa));
56
57 //double n = 0.725607765053267;
58 //double C = 11754255.4261;
59 //double XS = 700000;
60 //double YS = 12655612.0499;
61
62
63 //lamx = lamx-700000;
64 //lamy = lamy-12655612.0499;
65
66 double gamma;
67 gamma = atan(-(lamx-700000)/(lamy-12655612.0499));
68
69
70 lon = gamma/0.725607765053267 + 0.0523598775598299;
71 double R;
72 R = sqrt((lamx-700000) * (lamx-700000) + (lamy-12655612.0499) * (lamy-12655612.0499));
73
74 double latiso;
75 latiso = log((11754255.4261)/R)/(0.725607765053267);
76
77 double phiNew, phiOld;
78 phiOld =1;
79 phiNew= asin (tanh ( latiso + GRS_e * atanh(GRS_e * sin(phiOld))));
80 //printf("\nphiNew: %.20lf",phiNew);
81 while (abs_me(phiOld-phiNew) > 1e-10)
82 {
83
84 if(abs_me(phiOld-phiNew) > 1e-10)
85 {
86
87 phiOld = phiNew;
88 phiNew = asin(tanh(latiso+GRS_e*atanh(GRS_e*sin(phiOld))));
89 }
90 else
91 phiOld = phiNew;
92 }
93
94 lat = phiNew;
95
96}
Note: See TracBrowser for help on using the repository browser.