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