]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_Lambert93.cxx
Salome HOME
Merge remote-tracking branch 'origin/BR_IMPROVEMENTS' into BR_v14_rc
[modules/hydro.git] / src / HYDROData / HYDROData_Lambert93.cxx
1 // Copyright (C) 2014-2015  EDF-R&D
2 // This library is free software; you can redistribute it and/or
3 // modify it under the terms of the GNU Lesser General Public
4 // License as published by the Free Software Foundation; either
5 // version 2.1 of the License, or (at your option) any later version.
6 //
7 // This library is distributed in the hope that it will be useful,
8 // but WITHOUT ANY WARRANTY; without even the implied warranty of
9 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
10 // Lesser General Public License for more details.
11 //
12 // You should have received a copy of the GNU Lesser General Public
13 // License along with this library; if not, write to the Free Software
14 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
15 //
16 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
17 //
18
19 #include <HYDROData_Lambert93.h>
20
21 #include <math.h>
22
23
24
25 // Standard values
26
27 const double PI = 3.14159265;
28
29
30
31 // Base values of the Lambert-93
32
33 const double a = 6378137;               //  m  -- le demi-grand axe
34
35 const double f = 1.0 / 298.257222101;   // l'aplatissement
36
37
38
39 const double phi_1_deg = 44;            // deg -- le premier parallèle d'échelle
40
41 const double phi_2_deg = 49;            // deg -- le deuxième parallèle d'échell
42
43
44
45 const double lambda_0_deg = 3;          // deg -- la longitude d'origine donnée par le méridien central de Greenwich
46
47 const double phi_0_deg = 46.5;          // deg -- la latitude d'origine 
48
49
50
51 const double X_0 = 700000;              //  m  -- la coordonnée à l'origine
52
53 const double Y_0 = 6600000;             //  m  -- la coordonnée à l'origine 
54
55
56
57 // Derived values of the Lambert-93
58
59 const double b = a * ( 1 - f );         //  m  -- le demi-petit axe
60
61 const double e = sqrt( a*a - b*b ) / a; // l'excentricité
62
63
64
65 const double phi_0 = HYDROData_Lambert93::toRad( phi_0_deg );
66
67 const double phi_1 = HYDROData_Lambert93::toRad( phi_1_deg );
68
69 const double phi_2 = HYDROData_Lambert93::toRad( phi_2_deg );
70
71 const double lambda_0 = HYDROData_Lambert93::toRad( lambda_0_deg );
72
73
74
75
76
77 double cot( double x )
78
79 {
80
81   return cos( x ) / sin( x );
82
83 }
84
85
86
87 double ln( double x )
88
89 {
90
91   return log( x );
92
93 }
94
95
96
97 const double s1 = sin( phi_1 );
98
99 const double s2 = sin( phi_2 );
100
101 const double c1 = cos( phi_1 );
102
103 const double c2 = cos( phi_2 );
104
105
106
107 const double n1 = ln( c2/c1 ) + 1.0/2.0 * ln( (1-e*e*s1*s1)/(1-e*e*s2*s2) );
108
109 const double n2 = tan( phi_1 / 2 + PI/4 ) * pow( 1-e*s1, e/2 ) * pow( 1+e*s2, e/2 );
110
111 const double n3 = tan( phi_2 / 2 + PI/4 ) * pow( 1+e*s1, e/2 ) * pow( 1-e*s2, e/2 );
112
113 const double n = n1 / ( ln( n2/n3 ) );
114
115
116
117 const double p1 = a * c1 / ( n * sqrt( 1-e*e*s1*s1 ) );
118
119 const double p2 = tan( phi_1 / 2 + PI / 4 );
120
121 const double p3 = pow( (1-e*s1)/(1+e*s1), e/2 );
122
123 const double p_0 = p1 * pow( p2*p3, n );
124
125
126
127 double HYDROData_Lambert93::toRad( double theDeg )
128
129 {
130
131   return theDeg * PI / 180.0;
132
133 }
134
135
136
137 double HYDROData_Lambert93::toDeg( double theRad )
138
139 {
140
141   return theRad / PI * 180.0;
142
143 }
144
145
146
147 double HYDROData_Lambert93::calc_rho( double phi )
148
149 {
150
151   double c1 = cot( phi/2 + PI/4 );
152
153   double c2 = ( 1 + e * sin( phi ) ) / ( 1 - e * sin( phi ) );
154
155   double rho = p_0 * pow( c1 * pow( c2, e/2 ), n );
156
157   return rho;
158
159 }
160
161
162
163 void HYDROData_Lambert93::toXY( double theLatitudeDeg, double theLongitudeDeg,
164
165                                 double& theX, double& theY )
166
167 {
168
169   double phi = toRad( theLatitudeDeg );
170
171   double lambda = toRad( theLongitudeDeg );
172
173
174
175   double rho = calc_rho( phi );
176
177   double rho_0 = calc_rho( phi_0 );
178
179   double theta = n * ( lambda - lambda_0 );
180
181
182
183   theX = X_0 + rho * sin( theta );
184
185   theY = Y_0 + rho_0 - rho * cos( theta );
186
187 }
188
189
190
191 double arctan( double x )
192
193 {
194
195   return atan( x );
196
197 }
198
199
200
201 typedef double (*FUNC)( double );
202
203 double solve( FUNC f, double c, double x1, double x2, double eps )
204
205 {
206
207   double f1 = f( x1 ) - c;
208
209   double f2 = f( x2 ) - c;
210
211   while( fabs( x1 - x2 ) > eps )
212
213   {
214
215     double x = ( x1 + x2 ) / 2;
216
217     double fx = f( x ) - c;
218
219     bool b1 = f1>=0;
220
221     bool b2 = f2>=0;
222
223     bool b = fx>=0;
224
225     if( b==b1 )
226
227     {
228
229       x1 = x;
230
231       f1 = fx;
232
233     }
234
235     else
236
237     {
238
239       x2 = x;
240
241       f2 = fx;
242
243     }
244
245   }
246
247   return ( x1 + x2 ) / 2;
248
249 }
250
251
252
253 double F( double phi )
254
255 {
256
257   double f1 = tan( phi/2 + PI/4 );
258
259   double f2 = ( 1 - e*sin(phi) ) / ( 1 + e*sin(phi) );
260
261   return f1 * pow( f2, e/2 );
262
263 }
264
265
266
267 double Finv( double x, double eps )
268
269 {
270
271   return solve( F, x, 0, PI/2-eps, eps );
272
273 }
274
275
276
277 double HYDROData_Lambert93::calc_phi_inv( double rho, double eps )
278
279 {
280
281   double x = pow( p_0 / rho, 1/n );
282
283   double phi = Finv( x, eps );
284
285   return phi;
286
287 }
288
289
290
291 double HYDROData_Lambert93::calc_phi_ign( double rho, double eps )
292
293 {
294
295   double x = p_0 / rho;
296
297   double y = pow( x, 1/n );
298
299   double phi_i_1, phi_i = 2*arctan( y ) - PI/2;
300
301   while( true )
302
303   {
304
305     phi_i_1 = phi_i;
306
307     double z = y * pow( ( 1 + e*sin(phi_i_1) ) / ( 1 - e*sin(phi_i_1) ), e/2 );
308
309     phi_i = 2*arctan( pow( x, 1/n ) * z ) - PI/2;
310
311     if( fabs( phi_i - phi_i_1 ) < eps )
312
313       return phi_i;
314
315   }
316
317   return -1;
318
319 }
320
321
322
323 void HYDROData_Lambert93::toGeo( double theX, double theY,
324
325                        double& theLatitudeDeg, double& theLongitudeDeg,
326
327                        double theEps )
328
329 {
330
331   double rho_0 = calc_rho( phi_0 );
332
333   double rho = sqrt( pow( theX - X_0, 2 ) + pow( Y_0 - theY + rho_0, 2 ) );
334
335   double theta = 2 * arctan( ( theX - X_0 ) / ( Y_0 - theY + rho_0 + rho ) );
336
337
338
339   double lambda = theta / n + lambda_0;
340
341   double phi = calc_phi_inv( rho, theEps );
342
343
344
345   theLatitudeDeg = toDeg( phi );
346
347   theLongitudeDeg = toDeg( lambda );
348
349 }
350
351 void HYDROData_Lambert93::DMSToDeg( int theDeg,
352                                     int theMin,
353                                     double theSec,
354                                     double& theDegOut )
355 {
356   double aCoef = theDeg > 0 ? 1.0 : -1.0;
357
358   theDegOut = (double)theDeg;
359   theDegOut += aCoef * ( (double)theMin ) / 60.0;
360   theDegOut += aCoef * theSec / 3600.0;
361 }
362
363 void HYDROData_Lambert93::DMSToSec( int theDeg,
364                                     int theMin,
365                                     double theSec,
366                                     double& theSecOut )
367 {
368   double aCoef = theDeg > 0 ? 1.0 : -1.0;
369
370   theSecOut = theDeg * 3600.0;
371   theSecOut += aCoef * theMin * 60.0;
372   theSecOut += aCoef * theSec;
373 }
374
375 void HYDROData_Lambert93::degToDMS( double theDegIn,
376                                     int& theDeg,
377                                     int& theMin,
378                                     double& theSec )
379 {
380   theDeg = int( theDegIn );
381
382   double aRemainder = fabs( theDegIn - theDeg ) * 60.0;
383   theMin = int( aRemainder );
384
385   aRemainder = ( aRemainder - theMin ) * 60.0;
386   theSec = aRemainder;
387 }
388
389 void HYDROData_Lambert93::secToDMS( double theSecIn,
390                                     int& theDeg,
391                                     int& theMin,
392                                     double& theSec )
393 {
394   theDeg = int( theSecIn / 3600.0 );
395
396   double aRemainder = fabs( theSecIn - theDeg * 3600.0 );
397   theMin = int( aRemainder / 60.0 );
398
399   theSec = fabs( aRemainder - theMin * 60.0 );
400 }