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