LORENE
tbl_math.C
1/*
2 * Mathematical functions for class Tbl
3 *
4 */
5
6/*
7 * Copyright (c) 1999-2000 Jean-Alain Marck
8 * Copyright (c) 1999-2001 Eric Gourgoulhon
9 * Copyright (c) 1999-2001 Philippe Grandclement
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30
31
32/*
33 * $Id: tbl_math.C,v 1.5 2016/12/05 16:18:16 j_novak Exp $
34 * $Log: tbl_math.C,v $
35 * Revision 1.5 2016/12/05 16:18:16 j_novak
36 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37 *
38 * Revision 1.4 2014/10/13 08:53:41 j_novak
39 * Lorene classes and functions now belong to the namespace Lorene.
40 *
41 * Revision 1.3 2014/10/06 15:13:18 j_novak
42 * Modified #include directives to use c++ syntax.
43 *
44 * Revision 1.2 2012/01/17 10:38:48 j_penner
45 * added a Heaviside function
46 *
47 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
48 * LORENE
49 *
50 * Revision 2.4 1999/12/02 17:47:48 phil
51 * ajout de racine_cubique
52 *
53 * Revision 2.3 1999/10/29 15:05:32 eric
54 * Ajout de fonctions mathematiques (abs, norme, etc...).
55 *
56 * Revision 2.2 1999/03/01 15:10:19 eric
57 * *** empty log message ***
58 *
59 * Revision 2.1 1999/02/24 15:26:13 hyc
60 * *** empty log message ***
61 *
62 *
63 * $Header: /cvsroot/Lorene/C++/Source/Tbl/tbl_math.C,v 1.5 2016/12/05 16:18:16 j_novak Exp $
64 *
65 */
66
67// Headers C
68// ---------
69#include <cmath>
70#include <cstdlib>
71
72// Headers Lorene
73// --------------
74#include "tbl.h"
75
76 //-------//
77 // Sinus //
78 //-------//
79
80namespace Lorene {
81Tbl sin(const Tbl& ti)
82{
83 // Protection
84 assert(ti.get_etat() != ETATNONDEF) ;
85
86 // Cas ETATZERO
87 if (ti.get_etat() == ETATZERO) {
88 return ti ;
89 }
90
91 // Cas general
92 assert(ti.get_etat() == ETATQCQ) ; // sinon...
93 Tbl to(ti.dim) ; // Tbl resultat
94 to.set_etat_qcq() ;
95 int taille = ti.get_taille() ;
96 for (int i=0 ; i<taille ; i++) {
97 to.t[i] = sin(ti.t[i]) ;
98 }
99 return to ;
100}
101
102 //---------//
103 // Cosinus //
104 //---------//
105
106Tbl cos(const Tbl& ti)
107{
108 // Protection
109 assert(ti.get_etat() != ETATNONDEF) ;
110
111 Tbl to(ti.dim) ; // Tbl resultat
112 to.set_etat_qcq() ;
113
114 // Cas ETATZERO
115 if (ti.get_etat() == ETATZERO) {
116 int taille = ti.get_taille() ;
117 for (int i=0 ; i<taille ; i++) {
118 to.t[i] = 1 ;
119 }
120 return to ;
121 }
122
123 // Cas general
124 assert(ti.get_etat() == ETATQCQ) ; // sinon...
125 int taille = ti.get_taille() ;
126 for (int i=0 ; i<taille ; i++) {
127 to.t[i] = cos(ti.t[i]) ;
128 }
129 return to ;
130}
131
132 //----------//
133 // Tangente //
134 //----------//
135
136Tbl tan(const Tbl& ti)
137{
138 // Protection
139 assert(ti.get_etat() != ETATNONDEF) ;
140
141 // Cas ETATZERO
142 if (ti.get_etat() == ETATZERO) {
143 return ti ;
144 }
145
146 // Cas general
147 assert(ti.get_etat() == ETATQCQ) ; // sinon...
148 Tbl to(ti.dim) ; // Tbl resultat
149 to.set_etat_qcq() ;
150 int taille = ti.get_taille() ;
151 for (int i=0 ; i<taille ; i++) {
152 to.t[i] = tan(ti.t[i]) ;
153 }
154 return to ;
155}
156
157 //----------//
158 // ArcSinus //
159 //----------//
160
161Tbl asin(const Tbl& ti)
162{
163 // Protection
164 assert(ti.get_etat() != ETATNONDEF) ;
165
166 // Cas ETATZERO
167 if (ti.get_etat() == ETATZERO) {
168 return ti ;
169 }
170
171 // Cas general
172 assert(ti.get_etat() == ETATQCQ) ; // sinon...
173 Tbl to(ti.dim) ; // Tbl resultat
174 to.set_etat_qcq() ;
175 int taille = ti.get_taille() ;
176 for (int i=0 ; i<taille ; i++) {
177 to.t[i] = asin(ti.t[i]) ;
178 }
179 return to ;
180}
181
182 //------------//
183 // ArcCosinus //
184 //------------//
185
186Tbl acos(const Tbl& ti)
187{
188 // Protection
189 assert(ti.get_etat() != ETATNONDEF) ;
190
191 Tbl to(ti.dim) ; // Tbl resultat
192 to.set_etat_qcq() ;
193
194 // Cas ETATZERO
195 if (ti.get_etat() == ETATZERO) {
196 int taille = ti.get_taille() ;
197 for (int i=0 ; i<taille ; i++) {
198 to.t[i] = M_PI * .5 ;
199 }
200 return to ;
201 }
202
203 // Cas general
204 assert(ti.get_etat() == ETATQCQ) ; // sinon...
205 int taille = ti.get_taille() ;
206 for (int i=0 ; i<taille ; i++) {
207 to.t[i] = acos(ti.t[i]) ;
208 }
209 return to ;
210}
211
212 //-------------//
213 // ArcTangente //
214 //-------------//
215
216Tbl atan(const Tbl& ti)
217{
218 // Protection
219 assert(ti.get_etat() != ETATNONDEF) ;
220
221 // Cas ETATZERO
222 if (ti.get_etat() == ETATZERO) {
223 return ti ;
224 }
225
226 // Cas general
227 assert(ti.get_etat() == ETATQCQ) ; // sinon...
228 Tbl to(ti.dim) ; // Tbl resultat
229 to.set_etat_qcq() ;
230 int taille = ti.get_taille() ;
231 for (int i=0 ; i<taille ; i++) {
232 to.t[i] = atan(ti.t[i]) ;
233 }
234 return to ;
235}
236
237 //------//
238 // Sqrt //
239 //------//
240
241Tbl sqrt(const Tbl& ti)
242{
243 // Protection
244 assert(ti.get_etat() != ETATNONDEF) ;
245
246 // Cas ETATZERO
247 if (ti.get_etat() == ETATZERO) {
248 return ti ;
249 }
250
251 // Cas general
252 assert(ti.get_etat() == ETATQCQ) ; // sinon...
253 Tbl to(ti.dim) ; // Tbl resultat
254 to.set_etat_qcq() ;
255 int taille = ti.get_taille() ;
256 for (int i=0 ; i<taille ; i++) {
257 to.t[i] = sqrt(ti.t[i]) ;
258 }
259 return to ;
260}
261
262 //---------------//
263 // Exponantielle //
264 //---------------//
265
266Tbl exp(const Tbl& ti)
267{
268 // Protection
269 assert(ti.get_etat() != ETATNONDEF) ;
270
271 Tbl to(ti.dim) ; // Tbl resultat
272 to.set_etat_qcq() ;
273
274 // Cas ETATZERO
275 if (ti.get_etat() == ETATZERO) {
276 int taille = ti.get_taille() ;
277 for (int i=0 ; i<taille ; i++) {
278 to.t[i] = 1 ;
279 }
280 return to ;
281 }
282
283 // Cas general
284 assert(ti.get_etat() == ETATQCQ) ; // sinon...
285 int taille = ti.get_taille() ;
286 for (int i=0 ; i<taille ; i++) {
287 to.t[i] = exp(ti.t[i]) ;
288 }
289 return to ;
290}
291
292 //--------------------//
293 // Heaviside Function //
294 //--------------------//
295
296Tbl Heaviside(const Tbl& ti)
297{
298 // Protection
299 assert(ti.get_etat() != ETATNONDEF) ;
300
301 Tbl to(ti.dim) ; // Tbl resultat
302 to.set_etat_qcq() ;
303
304 // Cas ETATZERO
305 if (ti.get_etat() == ETATZERO) {
306 int taille = ti.get_taille() ;
307 for (int i=0 ; i<taille ; i++) {
308 to.t[i] = 0 ;
309 }
310 return to ;
311 }
312
313 // Cas general
314 assert(ti.get_etat() == ETATQCQ) ; // Otherwise
315 int taille = ti.get_taille() ;
316 for (int i=0 ; i<taille ; i++) {
317 if(ti.t[i] >= 0)
318 to.t[i] = 1 ;
319 else
320 to.t[i] = 0 ;
321 }
322 return to ;
323}
324
325 //-------------//
326 // Log naturel //
327 //-------------//
328
329Tbl log(const Tbl& ti)
330{
331 // Protection
332 assert(ti.get_etat() != ETATNONDEF) ;
333
334 // Cas ETATZERO
335 if (ti.get_etat() == ETATZERO) {
336 cout << "Tbl log: log(ETATZERO) !" << endl ;
337 abort () ;
338 }
339
340 // Cas general
341 assert(ti.get_etat() == ETATQCQ) ; // sinon...
342 Tbl to(ti.dim) ; // Tbl resultat
343 to.set_etat_qcq() ;
344 int taille = ti.get_taille() ;
345 for (int i=0 ; i<taille ; i++) {
346 to.t[i] = log(ti.t[i]) ;
347 }
348 return to ;
349}
350
351 //-------------//
352 // Log decimal //
353 //-------------//
354
355Tbl log10(const Tbl& ti)
356{
357 // Protection
358 assert(ti.get_etat() != ETATNONDEF) ;
359
360 // Cas ETATZERO
361 if (ti.get_etat() == ETATZERO) {
362 cout << "Tbl log10: log10(ETATZERO) !" << endl ;
363 abort () ;
364 }
365
366 // Cas general
367 assert(ti.get_etat() == ETATQCQ) ; // sinon...
368 Tbl to(ti.dim) ; // Tbl resultat
369 to.set_etat_qcq() ;
370 int taille = ti.get_taille() ;
371 for (int i=0 ; i<taille ; i++) {
372 to.t[i] = log10(ti.t[i]) ;
373 }
374 return to ;
375}
376
377 //--------------//
378 // Power entier //
379 //--------------//
380
381Tbl pow(const Tbl& ti, int n)
382{
383 // Protection
384 assert(ti.get_etat() != ETATNONDEF) ;
385
386 // Cas ETATZERO
387 if (ti.get_etat() == ETATZERO) {
388 if (n > 0) {
389 return ti ;
390 }
391 else {
392 cout << "Tbl pow: ETATZERO^n avec n<=0 ! "<< endl ;
393 abort () ;
394 }
395 }
396
397 // Cas general
398 assert(ti.get_etat() == ETATQCQ) ; // sinon...
399 Tbl to(ti.dim) ; // Tbl resultat
400 to.set_etat_qcq() ;
401 double x = n ;
402 int taille = ti.get_taille() ;
403 for (int i=0 ; i<taille ; i++) {
404 to.t[i] = pow(ti.t[i], x) ;
405 }
406 return to ;
407}
408
409 //--------------//
410 // Power double //
411 //--------------//
412
413Tbl pow(const Tbl& ti, double x)
414{
415 // Protection
416 assert(ti.get_etat() != ETATNONDEF) ;
417
418 // Cas ETATZERO
419 if (ti.get_etat() == ETATZERO) {
420 if (x > 0) {
421 return ti ;
422 }
423 else {
424 cout << "Tbl pow: ETATZERO^x avec x<=0 !" << endl ;
425 abort () ;
426 }
427 }
428
429 // Cas general
430 assert(ti.get_etat() == ETATQCQ) ; // sinon...
431 Tbl to(ti.dim) ; // Tbl resultat
432 to.set_etat_qcq() ;
433 int taille = ti.get_taille() ;
434 for (int i=0 ; i<taille ; i++) {
435 to.t[i] = pow(ti.t[i], x) ;
436 }
437 return to ;
438}
439
440 //----------------//
441 // Absolute value //
442 //----------------//
443
444Tbl abs(const Tbl& ti)
445{
446 // Protection
447 assert(ti.get_etat() != ETATNONDEF) ;
448
449 // Cas ETATZERO
450 if (ti.get_etat() == ETATZERO) {
451 return ti ;
452 }
453
454 // Cas general
455 assert(ti.get_etat() == ETATQCQ) ; // sinon...
456
457 Tbl to(ti.dim) ; // Tbl resultat
458 to.set_etat_qcq() ;
459
460 const double* xi = ti.t ;
461 double* xo = to.t ;
462 int taille = ti.get_taille() ;
463
464 for (int i=0 ; i<taille ; i++) {
465 xo[i] = fabs( xi[i] ) ;
466 }
467
468 return to ;
469}
470 //----------------//
471 // Cubic //
472 //----------------//
473
474Tbl racine_cubique(const Tbl& ti)
475{
476 // Protection
477 assert(ti.get_etat() != ETATNONDEF) ;
478
479 // Cas ETATZERO
480 if (ti.get_etat() == ETATZERO) {
481 return ti ;
482 }
483
484 // Cas general
485 assert(ti.get_etat() == ETATQCQ) ; // sinon...
486
487 Tbl absolute(abs(ti)) ;
488 Tbl res (pow(absolute, 1./3.)) ;
489
490 for (int i=0 ; i<ti.get_taille() ; i++)
491 if (ti.t[i] < 0)
492 res.t[i] *= -1 ;
493
494 return res ;
495}
496
497 //-------------------------------//
498 // max //
499 //-------------------------------//
500
501double max(const Tbl& ti) {
502
503 // Protection
504 assert(ti.get_etat() != ETATNONDEF) ;
505
506 // Cas particulier
507 if (ti.get_etat() == ETATZERO) {
508 return double(0) ;
509 }
510
511 // Cas general
512 assert(ti.get_etat() == ETATQCQ) ; // sinon....
513
514 const double* x = ti.t ;
515 double resu = x[0] ;
516 for (int i=1; i<ti.get_taille(); i++) {
517 if ( x[i] > resu ) resu = x[i] ;
518 }
519
520 return resu ;
521}
522
523 //-------------------------------//
524 // min //
525 //-------------------------------//
526
527double min(const Tbl& ti) {
528
529 // Protection
530 assert(ti.get_etat() != ETATNONDEF) ;
531
532 // Cas particulier
533 if (ti.get_etat() == ETATZERO) {
534 return double(0) ;
535 }
536
537 // Cas general
538 assert(ti.get_etat() == ETATQCQ) ; // sinon....
539
540 const double* x = ti.t ;
541 double resu = x[0] ;
542 for (int i=1; i<ti.get_taille(); i++) {
543 if ( x[i] < resu ) resu = x[i] ;
544 }
545
546 return resu ;
547}
548
549 //-------------------------------//
550 // norme //
551 //-------------------------------//
552
553double norme(const Tbl& ti) {
554
555 // Protection
556 assert(ti.get_etat() != ETATNONDEF) ;
557
558 double resu = 0 ;
559
560 if (ti.get_etat() != ETATZERO) { // on n'effectue la somme que si necessaire
561
562 assert(ti.get_etat() == ETATQCQ) ; // sinon....
563 const double* x = ti.t ;
564 for (int i=0; i<ti.get_taille(); i++) {
565 resu += fabs( x[i] ) ;
566 }
567
568 }
569
570 return resu ;
571}
572
573 //-------------------------------//
574 // diffrel //
575 //-------------------------------//
576
577double diffrel(const Tbl& t1, const Tbl& t2) {
578
579 // Protections
580 assert(t1.get_etat() != ETATNONDEF) ;
581 assert(t2.get_etat() != ETATNONDEF) ;
582
583 double norm2 = norme(t2) ;
584 double normdiff = norme(t1-t2) ;
585 double resu ;
586 if ( norm2 == double(0) ) {
587 resu = normdiff ;
588 }
589 else {
590 resu = normdiff / norm2 ;
591 }
592
593 return resu ;
594
595}
596
597 //-------------------------------//
598 // diffrelmax //
599 //-------------------------------//
600
601double diffrelmax(const Tbl& t1, const Tbl& t2) {
602
603 // Protections
604 assert(t1.get_etat() != ETATNONDEF) ;
605 assert(t2.get_etat() != ETATNONDEF) ;
606
607 double max2 = max(abs(t2)) ;
608 double maxdiff = max(abs(t1-t2)) ;
609 double resu ;
610 if ( max2 == double(0) ) {
611 resu = maxdiff ;
612 }
613 else {
614 resu = maxdiff / max2 ;
615 }
616
617 return resu ;
618
619}
620}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::dim
Dim_tbl dim
Number of dimensions, size,...
Definition tbl.h:172
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition tbl.h:394
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille).
Definition tbl.h:397
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene::atan
Cmp atan(const Cmp &)
Arctangent.
Definition cmp_math.C:198
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:223
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:325
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:273
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:72
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:507
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:484
Lorene::acos
Cmp acos(const Cmp &)
Arccosine.
Definition cmp_math.C:172
Lorene::asin
Cmp asin(const Cmp &)
Arcsine.
Definition cmp_math.C:147
Lorene::racine_cubique
Cmp racine_cubique(const Cmp &)
Cube root.
Definition cmp_math.C:248
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:461
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:438
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:97
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:413
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:542
Lorene::tan
Cmp tan(const Cmp &)
Tangent.
Definition cmp_math.C:123
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:299
Lorene::Heaviside
Mtbl Heaviside(const Mtbl &)
Heaviside function.
Definition mtbl_math.C:320
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:738