LORENE
valeur_math.C
1/*
2 * Mathematical functions for class Valeur
3 *
4 */
5
6/*
7 * Copyright (c) 1999-2000 Jean-Alain Marck
8 * Copyright (c) 1999-2001 Eric Gourgoulhon
9 * Copyright (c) 1999-2001 Philippe Grandclement
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30char valeur_math_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_math.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $" ;
31
32/*
33 * $Id: valeur_math.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $
34 * $Log: valeur_math.C,v $
35 * Revision 1.5 2014/10/13 08:53:50 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.4 2014/10/06 15:13:23 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.3 2012/01/17 10:39:27 j_penner
42 * added a Heaviside function
43 *
44 * Revision 1.2 2002/10/16 14:37:15 j_novak
45 * Reorganization of #include instructions of standard C++, in order to
46 * use experimental version 3 of gcc.
47 *
48 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
49 * LORENE
50 *
51 * Revision 2.4 1999/12/02 17:57:42 phil
52 * *** empty log message ***
53 *
54 * Revision 2.3 1999/11/09 16:18:02 phil
55 * ajout de racine_cubique
56 *
57 * Revision 2.2 1999/10/29 15:14:50 eric
58 * Ajout de fonctions mathematiques (abs, norme, etc...).
59 *
60 * Revision 2.1 1999/03/01 15:11:03 eric
61 * *** empty log message ***
62 *
63 * Revision 2.0 1999/02/24 15:40:32 hyc
64 * *** empty log message ***
65 *
66 *
67 * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_math.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $
68 *
69 */
70
71// Fichiers include
72// ----------------
73#include <cmath>
74#include <cstdlib>
75
76#include "valeur.h"
77
78
79 //-------//
80 // Sinus //
81 //-------//
82
83namespace Lorene {
85{
86 // Protection
87 assert(ti.get_etat() != ETATNONDEF) ;
88
89 // Cas ETATZERO
90 if (ti.get_etat() == ETATZERO) {
91 return ti ;
92 }
93
94 // Cas general
95 assert(ti.get_etat() == ETATQCQ) ; // sinon...
96 if (ti.c == 0x0) { // Il faut la valeur physique
97 ti.coef_i() ;
98 }
99 Valeur to(ti.get_mg()) ; // Valeur resultat
100 to.set_etat_c_qcq() ;
101 *(to.c) = sin( *(ti.c) ) ;
102 return to ;
103}
104
105 //---------//
106 // Cosinus //
107 //---------//
108
110{
111 // Protection
112 assert(ti.get_etat() != ETATNONDEF) ;
113
114 Valeur to(ti.get_mg()) ; // Valeur resultat
115 to.set_etat_c_qcq() ;
116
117 // Cas ETATZERO
118 if (ti.get_etat() == ETATZERO) {
119 *(to.c) = 1. ;
120 return to ;
121 }
122
123 // Cas general
124 assert(ti.get_etat() == ETATQCQ) ; // sinon...
125 if (ti.c == 0x0) { // Il faut la valeur physique
126 ti.coef_i() ;
127 }
128 *(to.c) = cos( *(ti.c) ) ;
129 return to ;
130}
131
132 //----------//
133 // Tangente //
134 //----------//
135
137{
138 // Protection
139 assert(ti.get_etat() != ETATNONDEF) ;
140
141 // Cas ETATZERO
142 if (ti.get_etat() == ETATZERO) {
143 return ti ;
144 }
145
146 // Cas general
147 assert(ti.get_etat() == ETATQCQ) ; // sinon...
148 if (ti.c == 0x0) { // Il faut la valeur physique
149 ti.coef_i() ;
150 }
151 Valeur to(ti.get_mg()) ; // Valeur resultat
152 to.set_etat_c_qcq() ;
153 *(to.c) = tan( *(ti.c) ) ;
154 return to ;
155}
156
157 //----------//
158 // ArcSinus //
159 //----------//
160
162{
163 // Protection
164 assert(ti.get_etat() != ETATNONDEF) ;
165
166 // Cas ETATZERO
167 if (ti.get_etat() == ETATZERO) {
168 return ti ;
169 }
170
171 // Cas general
172 assert(ti.get_etat() == ETATQCQ) ; // sinon...
173 if (ti.c == 0x0) { // Il faut la valeur physique
174 ti.coef_i() ;
175 }
176 Valeur to(ti.get_mg()) ; // Valeur resultat
177 to.set_etat_c_qcq() ;
178 *(to.c) = asin( *(ti.c) ) ;
179 return to ;
180}
181
182 //------------//
183 // ArcCosinus //
184 //------------//
185
187{
188 // Protection
189 assert(ti.get_etat() != ETATNONDEF) ;
190
191 Valeur to(ti.get_mg()) ; // Valeur resultat
192 to.set_etat_c_qcq() ;
193
194 // Cas ETATZERO
195 if (ti.get_etat() == ETATZERO) {
196 *(to.c) = M_PI * .5 ;
197 return to ;
198 }
199
200 // Cas general
201 assert(ti.get_etat() == ETATQCQ) ; // sinon...
202 if (ti.c == 0x0) { // Il faut la valeur physique
203 ti.coef_i() ;
204 }
205 *(to.c) = acos( *(ti.c) ) ;
206 return to ;
207}
208
209 //-------------//
210 // ArcTangente //
211 //-------------//
212
214{
215 // Protection
216 assert(ti.get_etat() != ETATNONDEF) ;
217
218 // Cas ETATZERO
219 if (ti.get_etat() == ETATZERO) {
220 return ti ;
221 }
222
223 // Cas general
224 assert(ti.get_etat() == ETATQCQ) ; // sinon...
225 if (ti.c == 0x0) { // Il faut la valeur physique
226 ti.coef_i() ;
227 }
228 Valeur to(ti.get_mg()) ; // Valeur resultat
229 to.set_etat_c_qcq() ;
230 *(to.c) = atan( *(ti.c) ) ;
231 return to ;
232}
233
234 //------//
235 // Sqrt //
236 //------//
237
239{
240 // Protection
241 assert(ti.get_etat() != ETATNONDEF) ;
242
243 // Cas ETATZERO
244 if (ti.get_etat() == ETATZERO) {
245 return ti ;
246 }
247
248 // Cas general
249 assert(ti.get_etat() == ETATQCQ) ; // sinon...
250 if (ti.c == 0x0) { // Il faut la valeur physique
251 ti.coef_i() ;
252 }
253 Valeur to(ti.get_mg()) ; // Valeur resultat
254 to.set_etat_c_qcq() ;
255 *(to.c) = sqrt( *(ti.c) ) ;
256 return to ;
257}
258
259 //---------------//
260 // Exponantielle //
261 //---------------//
262
264{
265 // Protection
266 assert(ti.get_etat() != ETATNONDEF) ;
267
268 Valeur to(ti.get_mg()) ; // Valeur resultat
269 to.set_etat_c_qcq() ;
270
271 // Cas ETATZERO
272 if (ti.get_etat() == ETATZERO) {
273 *(to.c) = 1. ;
274 return to ;
275 }
276
277 // Cas general
278 assert(ti.get_etat() == ETATQCQ) ; // sinon...
279 if (ti.c == 0x0) { // Il faut la valeur physique
280 ti.coef_i() ;
281 }
282 *(to.c) = exp( *(ti.c) ) ;
283 return to ;
284}
285
286 //--------------------//
287 // Heaviside Function //
288 //--------------------//
289
291{
292 // Protection
293 assert(ti.get_etat() != ETATNONDEF) ;
294
295 Valeur to(ti.get_mg()) ; // Valeur resultat
296 to.set_etat_c_qcq() ;
297
298 // Cas ETATZERO
299 if (ti.get_etat() == ETATZERO) {
300 *(to.c) = 0. ;
301 return to ;
302 }
303
304 // Cas general
305 assert(ti.get_etat() == ETATQCQ) ; // otherwise
306 if (ti.c == 0x0) { // Use the physical value << What?? (used Google translate)
307 ti.coef_i() ;
308 }
309
310 *(to.c) = Heaviside( *(ti.c) ) ;
311
312 return to ;
313}
314
315 //-------------//
316 // Log naturel //
317 //-------------//
318
320{
321 // Protection
322 assert(ti.get_etat() != ETATNONDEF) ;
323
324 // Cas ETATZERO
325 if (ti.get_etat() == ETATZERO) {
326 cout << "Valeur log: log(ETATZERO) !" << endl ;
327 abort () ;
328 }
329
330 // Cas general
331 assert(ti.get_etat() == ETATQCQ) ; // sinon...
332 if (ti.c == 0x0) { // Il faut la valeur physique
333 ti.coef_i() ;
334 }
335 Valeur to(ti.get_mg()) ; // Valeur resultat
336 to.set_etat_c_qcq() ;
337 *(to.c) = log( *(ti.c) ) ;
338 return to ;
339}
340
341 //-------------//
342 // Log decimal //
343 //-------------//
344
346{
347 // Protection
348 assert(ti.get_etat() != ETATNONDEF) ;
349
350 // Cas ETATZERO
351 if (ti.get_etat() == ETATZERO) {
352 cout << "Valeur log10: log10(ETATZERO) !" << endl ;
353 abort () ;
354 }
355
356 // Cas general
357 assert(ti.get_etat() == ETATQCQ) ; // sinon...
358 if (ti.c == 0x0) { // Il faut la valeur physique
359 ti.coef_i() ;
360 }
361 Valeur to(ti.get_mg()) ; // Valeur resultat
362 to.set_etat_c_qcq() ;
363 *(to.c) = log10( *(ti.c) ) ;
364 return to ;
365}
366
367 //--------------//
368 // Power entier //
369 //--------------//
370
371Valeur pow(const Valeur& ti, int n)
372{
373 // Protection
374 assert(ti.get_etat() != ETATNONDEF) ;
375
376 // Cas ETATZERO
377 if (ti.get_etat() == ETATZERO) {
378 if (n > 0) {
379 return ti ;
380 }
381 else {
382 cout << "Valeur pow: ETATZERO^n with n<=0 ! "<< endl ;
383 abort () ;
384 }
385 }
386
387 // Cas general
388 assert(ti.get_etat() == ETATQCQ) ; // sinon...
389 if (ti.c == 0x0) { // Il faut la valeur physique
390 ti.coef_i() ;
391 }
392 Valeur to(ti.get_mg()) ; // Valeur resultat
393 to.set_etat_c_qcq() ;
394 double x = n ;
395 *(to.c) = pow( *(ti.c), x ) ;
396 return to ;
397}
398
399 //--------------//
400 // Power double //
401 //--------------//
402
403Valeur pow(const Valeur& ti, double x)
404{
405 // Protection
406 assert(ti.get_etat() != ETATNONDEF) ;
407
408 // Cas ETATZERO
409 if (ti.get_etat() == ETATZERO) {
410 if (x > 0) {
411 return ti ;
412 }
413 else {
414 cout << "Valeur pow: ETATZERO^x with x<=0 !" << endl ;
415 abort () ;
416 }
417 }
418
419 // Cas general
420 assert(ti.get_etat() == ETATQCQ) ; // sinon...
421 if (ti.c == 0x0) { // Il faut la valeur physique
422 ti.coef_i() ;
423 }
424 Valeur to(ti.get_mg()) ; // Valeur resultat
425 to.set_etat_c_qcq() ;
426 *(to.c) = pow( *(ti.c), x ) ;
427 return to ;
428}
429
430 //----------------//
431 // Absolute value //
432 //----------------//
433
435{
436 // Protection
437 assert(vi.get_etat() != ETATNONDEF) ;
438
439 // Cas ETATZERO
440 if (vi.get_etat() == ETATZERO) {
441 return vi ;
442 }
443
444 // Cas general
445
446 assert(vi.get_etat() == ETATQCQ) ; // sinon...
447 if (vi.c == 0x0) { // Il faut la valeur physique
448 vi.coef_i() ;
449 }
450
451 Valeur vo(vi.get_mg()) ; // Valeur resultat
452
453 vo.set_etat_c_qcq() ;
454
455 *(vo.c) = abs( *(vi.c) ) ; // abs(Mtbl)
456
457 return vo ;
458}
459 //----------------//
460 // Cube root //
461 //----------------//
462
464{
465// Protection
466 assert(vi.get_etat() != ETATNONDEF) ;
467
468 // Cas ETATZERO
469 if (vi.get_etat() == ETATZERO) {
470 return vi ;
471 }
472
473 // Cas general
474 assert(vi.get_etat() == ETATQCQ) ; // sinon...
475 if (vi.c == 0x0) { // Il faut la valeur physique
476 vi.coef_i() ;
477 }
478 Valeur vo(vi.get_mg()) ; // Valeur resultat
479 vo.set_etat_c_qcq() ;
480 *(vo.c) = racine_cubique( *(vi.c) ) ;
481 return vo ;
482}
483 //-------------------------------//
484 // totalmax //
485 //-------------------------------//
486
487double totalmax(const Valeur& vi) {
488
489 // Protection
490 assert(vi.get_etat() != ETATNONDEF) ;
491
492// Tbl resu(vi.get_mg()->get_nzone()) ;
493 double resu ;
494
495 if (vi.get_etat() == ETATZERO) {
496 resu = 0 ;
497 }
498 else {
499
500 assert(vi.get_etat() == ETATQCQ) ;
501 if (vi.c == 0x0) { // Il faut la valeur physique
502 vi.coef_i() ;
503 }
504
505 resu = totalmax( *(vi.c) ) ; // max(Mtbl)
506
507 }
508
509 return resu ;
510}
511
512 //-------------------------------//
513 // totalmin //
514 //-------------------------------//
515
516double totalmin(const Valeur& vi) {
517
518 // Protection
519 assert(vi.get_etat() != ETATNONDEF) ;
520
521 double resu ;
522
523 if (vi.get_etat() == ETATZERO) {
524 resu = 0 ;
525 }
526 else {
527
528 assert(vi.get_etat() == ETATQCQ) ;
529 if (vi.c == 0x0) { // Il faut la valeur physique
530 vi.coef_i() ;
531 }
532
533 resu = totalmin( *(vi.c) ) ; // min(Mtbl)
534
535 }
536
537 return resu ;
538}
539
540 //-------------------------------//
541 // max //
542 //-------------------------------//
543
544Tbl max(const Valeur& vi) {
545
546 // Protection
547 assert(vi.get_etat() != ETATNONDEF) ;
548
549 Tbl resu(vi.get_mg()->get_nzone()) ;
550
551 if (vi.get_etat() == ETATZERO) {
552 resu.annule_hard() ;
553 }
554 else {
555
556 assert(vi.get_etat() == ETATQCQ) ;
557 if (vi.c == 0x0) { // Il faut la valeur physique
558 vi.coef_i() ;
559 }
560
561 resu = max( *(vi.c) ) ; // max(Mtbl)
562
563 }
564
565 return resu ;
566}
567
568 //-------------------------------//
569 // min //
570 //-------------------------------//
571
572Tbl min(const Valeur& vi) {
573
574 // Protection
575 assert(vi.get_etat() != ETATNONDEF) ;
576
577 Tbl resu(vi.get_mg()->get_nzone()) ;
578
579 if (vi.get_etat() == ETATZERO) {
580 resu.annule_hard() ;
581 }
582 else {
583
584 assert(vi.get_etat() == ETATQCQ) ;
585 if (vi.c == 0x0) { // Il faut la valeur physique
586 vi.coef_i() ;
587 }
588
589 resu = min( *(vi.c) ) ; // min(Mtbl)
590
591 }
592
593 return resu ;
594}
595
596 //-------------------------------//
597 // norme //
598 //-------------------------------//
599
601
602 // Protection
603 assert(vi.get_etat() != ETATNONDEF) ;
604
605 Tbl resu(vi.get_mg()->get_nzone()) ;
606
607 if (vi.get_etat() == ETATZERO) {
608 resu.annule_hard() ;
609 }
610 else {
611
612 assert(vi.get_etat() == ETATQCQ) ;
613 if (vi.c == 0x0) { // Il faut la valeur physique
614 vi.coef_i() ;
615 }
616
617 resu = norme( *(vi.c) ) ; // norme(Mtbl)
618
619 }
620
621 return resu ;
622}
623
624 //-------------------------------//
625 // diffrel //
626 //-------------------------------//
627
628Tbl diffrel(const Valeur& v1, const Valeur& v2) {
629
630 // Protections
631 assert(v1.get_etat() != ETATNONDEF) ;
632 assert(v2.get_etat() != ETATNONDEF) ;
633
634 int nz = v1.get_mg()->get_nzone() ;
635 Tbl resu(nz) ;
636
637 Valeur diff = v1 - v2 ;
638 Tbl normdiff = norme(diff) ;
639 Tbl norme2 = norme(v2) ;
640
641 assert(normdiff.get_etat() == ETATQCQ) ;
642 assert(norme2.get_etat() == ETATQCQ) ;
643
644 resu.set_etat_qcq() ;
645 for (int l=0; l<nz; l++) {
646 if ( norme2(l) == double(0) ) {
647 resu.set(l) = normdiff(l) ;
648 }
649 else{
650 resu.set(l) = normdiff(l) / norme2(l) ;
651 }
652 }
653
654 return resu ;
655
656}
657
658 //-------------------------------//
659 // diffrelmax //
660 //-------------------------------//
661
662Tbl diffrelmax(const Valeur& v1, const Valeur& v2) {
663
664 // Protections
665 assert(v1.get_etat() != ETATNONDEF) ;
666 assert(v2.get_etat() != ETATNONDEF) ;
667
668 int nz = v1.get_mg()->get_nzone() ;
669 Tbl resu(nz) ;
670
671 Tbl max2 = max(abs(v2)) ;
672 Valeur diff = v1 - v2 ;
673 Tbl maxdiff = max(abs(diff)) ;
674
675 assert(maxdiff.get_etat() == ETATQCQ) ;
676 assert(max2.get_etat() == ETATQCQ) ;
677
678 resu.set_etat_qcq() ;
679 for (int l=0; l<nz; l++) {
680 if ( max2(l) == double(0) ) {
681 resu.set(l) = maxdiff(l) ;
682 }
683 else{
684 resu.set(l) = maxdiff(l) / max2(l) ;
685 }
686 }
687
688 return resu ;
689
690}
691
692}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic array class.
Definition tbl.h:161
Values and coefficients of a (real-value) function.
Definition valeur.h:287
int get_etat() const
Returns the logical state.
Definition valeur.h:726
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