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

Last change on this file since 75 was 66, checked in by Marek Kurdej, 11 years ago

Documentation: file info.
  • Property svn:keywords set to Id
File size: 2.5 KB
Line 
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
15double atanh(double z)
16{
17 return (0.5 * log((z+1)/(1-z)));
18}
19
20double abs_me(double x)
21{
22 if(x<0) x=-x;
23 return x;
24}
25
26//transformation des longitudes tatitudes en Lambert93
27void 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
59void 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}
Note: See TracBrowser for help on using the repository browser.