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