LORENE
scalar_math.C
1 /*
2  * Mathematical functions for the Scalar class.
3  *
4  * These functions are not member functions of the Scalar class.
5  *
6  */
7 
8 /*
9  * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
10  *
11  * Copyright (c) 1999-2003 Eric Gourgoulhon
12  * Copyright (c) 1999-2001 Philippe Grandclement
13  *
14  * This file is part of LORENE.
15  *
16  * LORENE is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation; either version 2 of the License, or
19  * (at your option) any later version.
20  *
21  * LORENE is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with LORENE; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  *
30  */
31 
32 
33 
34 /*
35  * $Id: scalar_math.C,v 1.7 2016/12/05 16:18:18 j_novak Exp $
36  * $Log: scalar_math.C,v $
37  * Revision 1.7 2016/12/05 16:18:18 j_novak
38  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
39  *
40  * Revision 1.6 2014/10/13 08:53:46 j_novak
41  * Lorene classes and functions now belong to the namespace Lorene.
42  *
43  * Revision 1.5 2014/10/06 15:16:16 j_novak
44  * Modified #include directives to use c++ syntax.
45  *
46  * Revision 1.4 2012/01/17 10:27:46 j_penner
47  * added a Heaviside function
48  *
49  * Revision 1.3 2003/10/10 15:57:29 j_novak
50  * Added the state one (ETATUN) to the class Scalar
51  *
52  * Revision 1.2 2003/10/01 13:04:44 e_gourgoulhon
53  * The method Tensor::get_mp() returns now a reference (and not
54  * a pointer) onto a mapping.
55  *
56  * Revision 1.1 2003/09/25 08:06:56 e_gourgoulhon
57  * First versions (use Cmp as intermediate quantities).
58  *
59  *
60  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_math.C,v 1.7 2016/12/05 16:18:18 j_novak Exp $
61  *
62  */
63 
64 // Headers C
65 #include <cmath>
66 
67 // Headers Lorene
68 #include "tensor.h"
69 
70  //-------//
71  // Sine //
72  //-------//
73 
74 namespace Lorene {
75 Scalar sin(const Scalar& ci)
76 {
77  // Protection
78  assert(ci.get_etat() != ETATNONDEF) ;
79 
80  // Cas ETATZERO
81  if (ci.get_etat() == ETATZERO) {
82  return ci ;
83  }
84 
85  // Cas ETATUN
86  if (ci.get_etat() == ETATUN) {
87  Scalar resu(ci.get_mp()) ;
88  resu = sin(double(1)) ;
89  return resu ;
90  }
91 
92  // Cas general
93  assert(ci.get_etat() == ETATQCQ) ; // sinon...
94 
95  Scalar co(ci.get_mp()) ; // result
96 
97  co.set_etat_qcq() ;
98  co.va = sin( ci.va ) ;
99 
100  return co ;
101 }
102 
103  //---------//
104  // Cosine //
105  //---------//
106 
107 Scalar cos(const Scalar& ci)
108 {
109  // Protection
110  assert(ci.get_etat() != ETATNONDEF) ;
111 
112  Scalar co(ci.get_mp()) ; // result
113 
114  // Cas ETATZERO
115  if (ci.get_etat() == ETATZERO) {
116  co.set_etat_qcq() ;
117  co.va = double(1) ;
118  }
119  else {
120  // Cas ETATUN
121  if (ci.get_etat() == ETATUN) {
122  co = cos(double(1)) ;
123  }
124  else {
125  // Cas general
126  assert(ci.get_etat() == ETATQCQ) ; // sinon...
127  co.set_etat_qcq() ;
128  co.va = cos( ci.va ) ;
129  }
130  }
131 
132  return co ;
133 }
134 
135  //----------//
136  // Tangent //
137  //----------//
138 
139 Scalar tan(const Scalar& ci)
140 {
141  // Protection
142  assert(ci.get_etat() != ETATNONDEF) ;
143 
144  // Cas ETATZERO
145  if (ci.get_etat() == ETATZERO) {
146  return ci ;
147  }
148 
149  // Cas ETATUN
150  if (ci.get_etat() == ETATUN) {
151  Scalar resu(ci.get_mp()) ;
152  resu = tan(double(1)) ;
153  return resu ;
154  }
155 
156  // Cas general
157  assert(ci.get_etat() == ETATQCQ) ; // sinon...
158 
159  Scalar co(ci.get_mp()) ; // result
160  co.set_etat_qcq() ;
161  co.va = tan( ci.va ) ;
162 
163  return co ;
164 }
165 
166  //----------//
167  // Arcsine //
168  //----------//
169 
170 Scalar asin(const Scalar& ci)
171 {
172  // Protection
173  assert(ci.get_etat() != ETATNONDEF) ;
174 
175  // Cas ETATZERO
176  if (ci.get_etat() == ETATZERO) {
177  return ci ;
178  }
179 
180  // Cas ETATUN
181  if (ci.get_etat() == ETATUN) {
182  Scalar resu(ci.get_mp()) ;
183  resu = M_PI_2 ;
184  return resu ;
185  }
186 
187  // Cas general
188  assert(ci.get_etat() == ETATQCQ) ; // sinon...
189 
190  Scalar co(ci.get_mp()) ; // result
191 
192  co.set_etat_qcq() ;
193  co.va = asin( ci.va ) ;
194 
195  return co ;
196 }
197 
198  //-------------//
199  // Arccosine //
200  //-------------//
201 
202 Scalar acos(const Scalar& ci)
203 {
204  // Protection
205  assert(ci.get_etat() != ETATNONDEF) ;
206 
207  Scalar co(ci.get_mp()) ; // result
208 
209  // Cas ETATZERO
210  if (ci.get_etat() == ETATZERO) {
211  co.set_etat_qcq() ;
212  co.va = double(0.5) * M_PI ;
213  }
214  else {
215  // Cas ETATUN
216  if (ci.get_etat() == ETATUN) {
217  co.set_etat_zero() ;
218  }
219  else {
220  // Cas general
221  assert(ci.get_etat() == ETATQCQ) ; // sinon...
222  co.set_etat_qcq() ;
223  co.va = acos( ci.va ) ;
224  }
225  }
226 
227  return co ;
228 }
229 
230  //-------------//
231  // Arctangent //
232  //-------------//
233 
234 Scalar atan(const Scalar& ci)
235 {
236  // Protection
237  assert(ci.get_etat() != ETATNONDEF) ;
238 
239  // Cas ETATZERO
240  if (ci.get_etat() == ETATZERO) {
241  return ci ;
242  }
243 
244  // Cas ETATUN
245  if (ci.get_etat() == ETATUN) {
246  Scalar resu(ci.get_mp()) ;
247  resu = 0.25*M_PI ;
248  return resu ;
249  }
250 
251  // Cas general
252  assert(ci.get_etat() == ETATQCQ) ; // sinon...
253 
254  Scalar co(ci.get_mp()) ; // result
255 
256  co.set_etat_qcq() ;
257  co.va = atan( ci.va ) ;
258 
259  return co ;
260 }
261 
262  //-------------//
263  // Square root //
264  //-------------//
265 
266 Scalar sqrt(const Scalar& ci)
267 {
268  // Protection
269  assert(ci.get_etat() != ETATNONDEF) ;
270 
271  // Cas ETATZERO
272  if (ci.get_etat() == ETATZERO) {
273  return ci ;
274  }
275 
276  // Cas ETATUN
277  if (ci.get_etat() == ETATUN) {
278  return ci ;
279  }
280 
281  // Cas general
282  assert(ci.get_etat() == ETATQCQ) ; // sinon...
283 
284  Scalar co(ci.get_mp()) ; // result
285 
286  co.set_etat_qcq() ;
287  co.va = sqrt( ci.va ) ;
288 
289  return co ;
290 }
291 
292  //-------------//
293  // Cube root //
294  //-------------//
295 
297 {
298  // Protection
299  assert(ci.get_etat() != ETATNONDEF) ;
300 
301  // Cas ETATZERO
302  if (ci.get_etat() == ETATZERO) {
303  return ci ;
304  }
305 
306  // Cas ETATUN
307  if (ci.get_etat() == ETATUN) {
308  return ci ;
309  }
310 
311  // Cas general
312  assert(ci.get_etat() == ETATQCQ) ; // sinon...
313 
314  Scalar co(ci.get_mp()) ; // result
315 
316  co.set_etat_qcq() ;
317  co.va = racine_cubique( ci.va ) ;
318 
319  return co ;
320 }
321 
322  //--------------//
323  // Exponential //
324  //--------------//
325 
326 Scalar exp(const Scalar& ci)
327 {
328  // Protection
329  assert(ci.get_etat() != ETATNONDEF) ;
330 
331  Scalar co(ci.get_mp()) ; // result
332 
333  // Cas ETATZERO
334  if (ci.get_etat() == ETATZERO) {
335  co.set_etat_one() ;
336  }
337  else {
338  // Cas ETATUN
339  if (ci.get_etat() == ETATUN) {
340  co.set_etat_qcq() ;
341  co = M_E ;
342  }
343  else {
344  // Cas general
345  assert(ci.get_etat() == ETATQCQ) ; // sinon...
346  co.set_etat_qcq() ;
347  co.va = exp( ci.va ) ;
348  }
349  }
350 
351  return co ;
352 }
353 
354  //---------------------//
355  // Heaviside Function //
356  //---------------------//
357 
359 {
360  // Protection
361  assert(ci.get_etat() != ETATNONDEF) ;
362 
363  Scalar co(ci.get_mp()) ; // make output a copy, to ensure the same structure
364 
365  // if input state is zero, return zero
366  if (ci.get_etat() == ETATZERO) {
367  co.set_etat_zero() ;
368  }
369  else {
370  // if input state is one, return one
371  if (ci.get_etat() == ETATUN) {
372  co.set_etat_one() ;
373  }
374  else {
375  // In general return the Heaviside function
376  assert(ci.get_etat() == ETATQCQ) ; // otherwise
377  co.set_etat_qcq() ;
378  co.va = Heaviside( ci.va ) ;
379  }
380  }
381 
382  return co ;
383 }
384  //---------------------//
385  // Neperian logarithm //
386  //---------------------//
387 
388 Scalar log(const Scalar& ci)
389 {
390  // Protection
391  assert(ci.get_etat() != ETATNONDEF) ;
392 
393  // Cas ETATZERO
394  if (ci.get_etat() == ETATZERO) {
395  cout << "Argument of log is ZERO in log(Scalar) !" << endl ;
396  abort() ;
397  }
398 
399  // Cas ETATUN
400  if (ci.get_etat() == ETATUN) {
401  Scalar resu(ci.get_mp()) ;
402  resu.set_etat_zero() ;
403  return resu ;
404  }
405 
406  // Cas general
407  assert(ci.get_etat() == ETATQCQ) ; // sinon...
408 
409  Scalar co(ci.get_mp()) ; // result
410 
411  co.set_etat_qcq() ;
412  co.va = log( ci.va ) ;
413 
414  return co ;
415 }
416 
417  //---------------------//
418  // Decimal logarithm //
419  //---------------------//
420 
421 Scalar log10(const Scalar& ci)
422 {
423  // Protection
424  assert(ci.get_etat() != ETATNONDEF) ;
425 
426  // Cas ETATZERO
427  if (ci.get_etat() == ETATZERO) {
428  cout << "Argument of log10 is ZERO in log10(Scalar) !" << endl ;
429  abort() ;
430  }
431 
432  // Cas ETATUN
433  if (ci.get_etat() == ETATUN) {
434  Scalar resu(ci.get_mp()) ;
435  resu.set_etat_zero() ;
436  return resu ;
437  }
438 
439  // Cas general
440  assert(ci.get_etat() == ETATQCQ) ; // sinon...
441 
442  Scalar co(ci.get_mp()) ; // result
443 
444  co.set_etat_qcq() ;
445  co.va = log10( ci.va ) ;
446 
447  return co ;
448 }
449 
450  //------------------//
451  // Power (integer) //
452  //------------------//
453 
454 Scalar pow(const Scalar& ci, int n)
455 {
456  // Protection
457  assert(ci.get_etat() != ETATNONDEF) ;
458 
459  // Cas ETATZERO
460  if (ci.get_etat() == ETATZERO) {
461  if (n > 0) {
462  return ci ;
463  }
464  else {
465  cout << "pow(Scalar, int) : ETATZERO^n with n <= 0 !" << endl ;
466  abort() ;
467  }
468  }
469 
470  // Cas ETATUN
471  if (ci.get_etat() == ETATUN) {
472  return ci ;
473  }
474 
475  // Cas general
476  assert(ci.get_etat() == ETATQCQ) ; // sinon...
477 
478  Scalar co(ci.get_mp()) ; // result
479 
480  co.set_etat_qcq() ;
481  co.va = pow(ci.va, n) ;
482 
483  return co ;
484 }
485 
486  //-----------------//
487  // Power (double) //
488  //-----------------//
489 
490 Scalar pow(const Scalar& ci, double x)
491 {
492  // Protection
493  assert(ci.get_etat() != ETATNONDEF) ;
494 
495  // Cas ETATZERO
496  if (ci.get_etat() == ETATZERO) {
497  if (x > double(0)) {
498  return ci ;
499  }
500  else {
501  cout << "pow(Scalar, double) : ETATZERO^x with x <= 0 !" << endl ;
502  abort() ;
503  }
504  }
505 
506  // Cas ETATUN
507  if (ci.get_etat() == ETATUN) {
508  return ci ;
509  }
510 
511  // Cas general
512  assert(ci.get_etat() == ETATQCQ) ; // sinon...
513 
514  Scalar co(ci.get_mp()) ; // result
515 
516  co.set_etat_qcq() ;
517  co.va = pow(ci.va, x) ;
518 
519  return co ;
520 }
521 
522  //-----------------//
523  // Absolute value //
524  //-----------------//
525 
526 Scalar abs(const Scalar& ci)
527 {
528  // Protection
529  assert(ci.get_etat() != ETATNONDEF) ;
530 
531  // Cas ETATZERO
532  if (ci.get_etat() == ETATZERO) {
533  return ci ;
534  }
535 
536  // Cas ETATUN
537  if (ci.get_etat() == ETATUN) {
538  return ci ;
539  }
540 
541  // Cas general
542  assert(ci.get_etat() == ETATQCQ) ; // sinon...
543 
544  Scalar co(ci.get_mp()) ; // result
545 
546  co.set_etat_qcq() ;
547  co.va = abs( ci.va ) ;
548 
549  return co ;
550 }
551 
552  //-------------------------------//
553  // totalmax //
554  //-------------------------------//
555 
556 double totalmax(const Scalar& ci) {
557 
558  // Protection
559  assert(ci.get_etat() != ETATNONDEF) ;
560 
561 // Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
562  double resu ;
563 
564  if (ci.get_etat() == ETATZERO) {
565  resu = 0 ;
566  }
567  else {
568  if (ci.get_etat() == ETATUN) {
569  resu = 1 ;
570  }
571  else {
572  assert(ci.get_etat() == ETATQCQ) ;
573 
574  resu = totalmax( ci.va ) ; // max(Valeur)
575  }
576  }
577 
578  return resu ;
579 }
580 
581  //-------------------------------//
582  // totalmin //
583  //-------------------------------//
584 
585 double totalmin(const Scalar& ci) {
586 
587  // Protection
588  assert(ci.get_etat() != ETATNONDEF) ;
589 
590 // Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
591  double resu ;
592 
593  if (ci.get_etat() == ETATZERO) {
594  resu = 0 ;
595  }
596  else {
597  if (ci.get_etat() == ETATUN) {
598  resu = 1 ;
599  }
600  else {
601  assert(ci.get_etat() == ETATQCQ) ;
602 
603  resu = totalmin( ci.va ) ; // min(Valeur)
604  }
605  }
606 
607  return resu ;
608 }
609 
610  //-------------------------------//
611  // max //
612  //-------------------------------//
613 
614 Tbl max(const Scalar& ci) {
615 
616  // Protection
617  assert(ci.get_etat() != ETATNONDEF) ;
618 
619  Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
620 
621  if (ci.get_etat() == ETATZERO) {
622  resu.annule_hard() ;
623  }
624  else {
625  if (ci.get_etat() == ETATUN) {
626  resu = 1 ;
627  }
628  else {
629  assert(ci.get_etat() == ETATQCQ) ;
630 
631  resu = max( ci.va ) ; // max(Valeur)
632  }
633  }
634 
635  return resu ;
636 }
637 
638  //-------------------------------//
639  // min //
640  //-------------------------------//
641 
642 Tbl min(const Scalar& ci) {
643 
644  // Protection
645  assert(ci.get_etat() != ETATNONDEF) ;
646 
647  Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
648 
649  if (ci.get_etat() == ETATZERO) {
650  resu.annule_hard() ;
651  }
652  else {
653  if (ci.get_etat() == ETATUN) {
654  resu = 1 ;
655  }
656  else {
657  assert(ci.get_etat() == ETATQCQ) ;
658 
659  resu = min( ci.va ) ; // min(Valeur)
660  }
661  }
662 
663  return resu ;
664 }
665 
666  //-------------------------------//
667  // norme //
668  //-------------------------------//
669 
670 Tbl norme(const Scalar& ci) {
671 
672  // Protection
673  assert(ci.get_etat() != ETATNONDEF) ;
674 
675  Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
676 
677  if (ci.get_etat() == ETATZERO) {
678  resu.annule_hard() ;
679  }
680  else {
681  if (ci.get_etat() == ETATUN) {
682  resu = 1 ;
683  }
684  else {
685  assert(ci.get_etat() == ETATQCQ) ;
686 
687  resu = norme( ci.va ) ; // norme(Valeur)
688  }
689  }
690 
691  return resu ;
692 }
693 
694  //-------------------------------//
695  // diffrel //
696  //-------------------------------//
697 
698 Tbl diffrel(const Scalar& c1, const Scalar& c2) {
699 
700  // Protections
701  assert(c1.get_etat() != ETATNONDEF) ;
702  assert(c2.get_etat() != ETATNONDEF) ;
703 
704  int nz = c1.get_mp().get_mg()->get_nzone() ;
705  Tbl resu(nz) ;
706 
707  Scalar diff = c1 - c2 ; // la compatibilite dzpuis est testee a ce niveau
708 
709  Tbl normdiff = norme(diff) ;
710  Tbl norme2 = norme(c2) ;
711 
712  assert(normdiff.get_etat() == ETATQCQ) ;
713  assert(norme2.get_etat() == ETATQCQ) ;
714 
715  resu.set_etat_qcq() ;
716  for (int l=0; l<nz; l++) {
717  if ( norme2(l) == double(0) ) {
718  resu.set(l) = normdiff(l) ;
719  }
720  else{
721  resu.set(l) = normdiff(l) / norme2(l) ;
722  }
723  }
724 
725  return resu ;
726 
727 }
728 
729  //-------------------------------//
730  // diffrelmax //
731  //-------------------------------//
732 
733 Tbl diffrelmax(const Scalar& c1, const Scalar& c2) {
734 
735  // Protections
736  assert(c1.get_etat() != ETATNONDEF) ;
737  assert(c2.get_etat() != ETATNONDEF) ;
738 
739  int nz = c1.get_mp().get_mg()->get_nzone() ;
740  Tbl resu(nz) ;
741 
742  Tbl max2 = max(abs(c2)) ;
743 
744  Scalar diff = c1 - c2 ; // la compatibilite dzpuis est testee a ce niveau
745 
746  Tbl maxdiff = max(abs(diff)) ;
747 
748  assert(maxdiff.get_etat() == ETATQCQ) ;
749  assert(max2.get_etat() == ETATQCQ) ;
750 
751  resu.set_etat_qcq() ;
752  for (int l=0; l<nz; l++) {
753  if ( max2(l) == double(0) ) {
754  resu.set(l) = maxdiff(l) ;
755  }
756  else{
757  resu.set(l) = maxdiff(l) / max2(l) ;
758  }
759  }
760 
761  return resu ;
762 
763 }
764 
765 }
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
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
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
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
int get_etat() const
Gives the logical state.
Definition: tbl.h:414
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:461
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
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
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
Cmp atan(const Cmp &)
Arctangent.
Definition: cmp_math.C:198
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition: scalar.C:340
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Valeur va
The numerical value of the Scalar.
Definition: scalar.h:411
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
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
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
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