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
32char 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
71namespace Lorene {
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
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
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
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
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
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
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
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
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
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
451Scalar 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
487Scalar 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
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
553double 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
582double 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
611Tbl 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
639Tbl 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
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
695Tbl 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
730Tbl 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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Basic array class.
Definition tbl.h:161
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
Lorene prototypes.
Definition app_hor.h:64