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
74namespace Lorene {
75Scalar 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
107Scalar 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
139Scalar 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
170Scalar 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
202Scalar 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
234Scalar 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
266Scalar 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
296Scalar racine_cubique(const Scalar& ci)
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
326Scalar 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
358Scalar Heaviside(const Scalar& ci)
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
388Scalar 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
421Scalar 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
454Scalar 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
490Scalar 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
526Scalar 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
556double 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
585double 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
614Tbl 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
642Tbl 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
670Tbl 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
698Tbl 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
733Tbl 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}
Lorene::Scalar::set_etat_one
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition scalar.C:340
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:359
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:330
Lorene::Scalar::Scalar
Scalar(const Map &mpi)
Constructor from mapping.
Definition scalar.C:210
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:560
Lorene::Scalar::va
Valeur va
The numerical value of the Scalar.
Definition scalar.h:411
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition tbl.h:394
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:375
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::atan
Cmp atan(const Cmp &)
Arctangent.
Definition cmp_math.C:198
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:223
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:325
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:273
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:72
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:507
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:484
Lorene::acos
Cmp acos(const Cmp &)
Arccosine.
Definition cmp_math.C:172
Lorene::asin
Cmp asin(const Cmp &)
Arcsine.
Definition cmp_math.C:147
Lorene::racine_cubique
Cmp racine_cubique(const Cmp &)
Cube root.
Definition cmp_math.C:248
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:461
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:438
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:97
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:413
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:542
Lorene::tan
Cmp tan(const Cmp &)
Tangent.
Definition cmp_math.C:123
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:299
Lorene::Heaviside
Mtbl Heaviside(const Mtbl &)
Heaviside function.
Definition mtbl_math.C:320
Lorene::totalmin
double totalmin(const Mtbl &)
Minimum value of the Mtbl elements in all domain.
Definition mtbl_math.C:525
Lorene::totalmax
double totalmax(const Mtbl &)
Maximum value of the Mtbl elements in all domains.
Definition mtbl_math.C:497
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:874
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:738