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 
82 namespace Lorene {
83 Mtbl 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 
108 Mtbl 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 
137 Mtbl 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 
162 Mtbl 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 
187 Mtbl 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 
216 Mtbl 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 
241 Mtbl 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 
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 
291 Mtbl 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 
320 Mtbl 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 
349 Mtbl 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 
375 Mtbl 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 
401 Mtbl 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 
433 Mtbl 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 
464 Mtbl 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 
497 double 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 
525 double 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 
554 Tbl 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 
583 Tbl 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 
612 Tbl 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 
641 Tbl 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 
674 Tbl 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 }
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Cmp asin(const Cmp &)
Arcsine.
Definition: cmp_math.C:147
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Multi-domain array.
Definition: mtbl.h:118
Lorene prototypes.
Definition: app_hor.h:67
Cmp racine_cubique(const Cmp &)
Cube root.
Definition: cmp_math.C:248
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
int get_etat() const
Gives the logical state.
Definition: tbl.h:414
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:461
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
int get_etat() const
Gives the logical state.
Definition: mtbl.h:277
Cmp tan(const Cmp &)
Tangent.
Definition: cmp_math.C:123
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
int get_nzone() const
Gives the number of zones (domains)
Definition: mtbl.h:280
Cmp atan(const Cmp &)
Arctangent.
Definition: cmp_math.C:198
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
double totalmax(const Mtbl &)
Maximum value of the Mtbl elements in all domains.
Definition: mtbl_math.C:497
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:302
double totalmin(const Mtbl &)
Minimum value of the Mtbl elements in all domain.
Definition: mtbl_math.C:525
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
Cmp acos(const Cmp &)
Arccosine.
Definition: cmp_math.C:172
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:413
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:72
Mtbl Heaviside(const Mtbl &)
Heaviside function.
Definition: mtbl_math.C:320
Basic array class.
Definition: tbl.h:164
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:375
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542