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