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 char 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 
79 namespace Lorene {
80 Mtbl 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 
105 Mtbl 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 
134 Mtbl 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 
159 Mtbl 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 
184 Mtbl 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 
213 Mtbl 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 
238 Mtbl 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 
288 Mtbl 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 
317 Mtbl Heaviside(const Mtbl& ti)
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 
346 Mtbl 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 
372 Mtbl 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 
398 Mtbl 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 
430 Mtbl 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 
461 Mtbl 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 
494 double 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 
522 double 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 
551 Tbl 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 
580 Tbl 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 
609 Tbl 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 
638 Tbl 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 
671 Tbl 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 }
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Cmp asin(const Cmp &)
Arcsine.
Definition: cmp_math.C:144
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Multi-domain array.
Definition: mtbl.h:118
Lorene prototypes.
Definition: app_hor.h:64
Cmp racine_cubique(const Cmp &)
Cube root.
Definition: cmp_math.C:245
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:458
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int get_etat() const
Gives the logical state.
Definition: mtbl.h:277
Cmp tan(const Cmp &)
Tangent.
Definition: cmp_math.C:120
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
int get_nzone() const
Gives the number of zones (domains)
Definition: mtbl.h:280
Cmp atan(const Cmp &)
Arctangent.
Definition: cmp_math.C:195
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
double totalmax(const Mtbl &)
Maximum value of the Mtbl elements in all domains.
Definition: mtbl_math.C:494
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
double totalmin(const Mtbl &)
Minimum value of the Mtbl elements in all domain.
Definition: mtbl_math.C:522
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
Cmp acos(const Cmp &)
Arccosine.
Definition: cmp_math.C:169
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69
Mtbl Heaviside(const Mtbl &)
Heaviside function.
Definition: mtbl_math.C:317
Basic array class.
Definition: tbl.h:161
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
const Mg3d * get_mg() const
Gives the Mg3d on which the Mtbl is defined.
Definition: mtbl.h:274
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:539