LORENE
tbl_math.C
1/*
2 * Mathematical functions for class Tbl
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 tbl_math_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tbl/tbl_math.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $" ;
31
32/*
33 * $Id: tbl_math.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $
34 * $Log: tbl_math.C,v $
35 * Revision 1.4 2014/10/13 08:53:41 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.3 2014/10/06 15:13:18 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.2 2012/01/17 10:38:48 j_penner
42 * added a Heaviside function
43 *
44 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
45 * LORENE
46 *
47 * Revision 2.4 1999/12/02 17:47:48 phil
48 * ajout de racine_cubique
49 *
50 * Revision 2.3 1999/10/29 15:05:32 eric
51 * Ajout de fonctions mathematiques (abs, norme, etc...).
52 *
53 * Revision 2.2 1999/03/01 15:10:19 eric
54 * *** empty log message ***
55 *
56 * Revision 2.1 1999/02/24 15:26:13 hyc
57 * *** empty log message ***
58 *
59 *
60 * $Header: /cvsroot/Lorene/C++/Source/Tbl/tbl_math.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $
61 *
62 */
63
64// Headers C
65// ---------
66#include <cmath>
67#include <cstdlib>
68
69// Headers Lorene
70// --------------
71#include "tbl.h"
72
73 //-------//
74 // Sinus //
75 //-------//
76
77namespace Lorene {
78Tbl sin(const Tbl& ti)
79{
80 // Protection
81 assert(ti.get_etat() != ETATNONDEF) ;
82
83 // Cas ETATZERO
84 if (ti.get_etat() == ETATZERO) {
85 return ti ;
86 }
87
88 // Cas general
89 assert(ti.get_etat() == ETATQCQ) ; // sinon...
90 Tbl to(ti.dim) ; // Tbl resultat
91 to.set_etat_qcq() ;
92 int taille = ti.get_taille() ;
93 for (int i=0 ; i<taille ; i++) {
94 to.t[i] = sin(ti.t[i]) ;
95 }
96 return to ;
97}
98
99 //---------//
100 // Cosinus //
101 //---------//
102
103Tbl cos(const Tbl& ti)
104{
105 // Protection
106 assert(ti.get_etat() != ETATNONDEF) ;
107
108 Tbl to(ti.dim) ; // Tbl resultat
109 to.set_etat_qcq() ;
110
111 // Cas ETATZERO
112 if (ti.get_etat() == ETATZERO) {
113 int taille = ti.get_taille() ;
114 for (int i=0 ; i<taille ; i++) {
115 to.t[i] = 1 ;
116 }
117 return to ;
118 }
119
120 // Cas general
121 assert(ti.get_etat() == ETATQCQ) ; // sinon...
122 int taille = ti.get_taille() ;
123 for (int i=0 ; i<taille ; i++) {
124 to.t[i] = cos(ti.t[i]) ;
125 }
126 return to ;
127}
128
129 //----------//
130 // Tangente //
131 //----------//
132
133Tbl tan(const Tbl& ti)
134{
135 // Protection
136 assert(ti.get_etat() != ETATNONDEF) ;
137
138 // Cas ETATZERO
139 if (ti.get_etat() == ETATZERO) {
140 return ti ;
141 }
142
143 // Cas general
144 assert(ti.get_etat() == ETATQCQ) ; // sinon...
145 Tbl to(ti.dim) ; // Tbl resultat
146 to.set_etat_qcq() ;
147 int taille = ti.get_taille() ;
148 for (int i=0 ; i<taille ; i++) {
149 to.t[i] = tan(ti.t[i]) ;
150 }
151 return to ;
152}
153
154 //----------//
155 // ArcSinus //
156 //----------//
157
158Tbl asin(const Tbl& ti)
159{
160 // Protection
161 assert(ti.get_etat() != ETATNONDEF) ;
162
163 // Cas ETATZERO
164 if (ti.get_etat() == ETATZERO) {
165 return ti ;
166 }
167
168 // Cas general
169 assert(ti.get_etat() == ETATQCQ) ; // sinon...
170 Tbl to(ti.dim) ; // Tbl resultat
171 to.set_etat_qcq() ;
172 int taille = ti.get_taille() ;
173 for (int i=0 ; i<taille ; i++) {
174 to.t[i] = asin(ti.t[i]) ;
175 }
176 return to ;
177}
178
179 //------------//
180 // ArcCosinus //
181 //------------//
182
183Tbl acos(const Tbl& ti)
184{
185 // Protection
186 assert(ti.get_etat() != ETATNONDEF) ;
187
188 Tbl to(ti.dim) ; // Tbl resultat
189 to.set_etat_qcq() ;
190
191 // Cas ETATZERO
192 if (ti.get_etat() == ETATZERO) {
193 int taille = ti.get_taille() ;
194 for (int i=0 ; i<taille ; i++) {
195 to.t[i] = M_PI * .5 ;
196 }
197 return to ;
198 }
199
200 // Cas general
201 assert(ti.get_etat() == ETATQCQ) ; // sinon...
202 int taille = ti.get_taille() ;
203 for (int i=0 ; i<taille ; i++) {
204 to.t[i] = acos(ti.t[i]) ;
205 }
206 return to ;
207}
208
209 //-------------//
210 // ArcTangente //
211 //-------------//
212
213Tbl atan(const Tbl& ti)
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 Tbl to(ti.dim) ; // Tbl resultat
226 to.set_etat_qcq() ;
227 int taille = ti.get_taille() ;
228 for (int i=0 ; i<taille ; i++) {
229 to.t[i] = atan(ti.t[i]) ;
230 }
231 return to ;
232}
233
234 //------//
235 // Sqrt //
236 //------//
237
238Tbl sqrt(const Tbl& ti)
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 Tbl to(ti.dim) ; // Tbl resultat
251 to.set_etat_qcq() ;
252 int taille = ti.get_taille() ;
253 for (int i=0 ; i<taille ; i++) {
254 to.t[i] = sqrt(ti.t[i]) ;
255 }
256 return to ;
257}
258
259 //---------------//
260 // Exponantielle //
261 //---------------//
262
263Tbl exp(const Tbl& ti)
264{
265 // Protection
266 assert(ti.get_etat() != ETATNONDEF) ;
267
268 Tbl to(ti.dim) ; // Tbl resultat
269 to.set_etat_qcq() ;
270
271 // Cas ETATZERO
272 if (ti.get_etat() == ETATZERO) {
273 int taille = ti.get_taille() ;
274 for (int i=0 ; i<taille ; i++) {
275 to.t[i] = 1 ;
276 }
277 return to ;
278 }
279
280 // Cas general
281 assert(ti.get_etat() == ETATQCQ) ; // sinon...
282 int taille = ti.get_taille() ;
283 for (int i=0 ; i<taille ; i++) {
284 to.t[i] = exp(ti.t[i]) ;
285 }
286 return to ;
287}
288
289 //--------------------//
290 // Heaviside Function //
291 //--------------------//
292
294{
295 // Protection
296 assert(ti.get_etat() != ETATNONDEF) ;
297
298 Tbl to(ti.dim) ; // Tbl resultat
299 to.set_etat_qcq() ;
300
301 // Cas ETATZERO
302 if (ti.get_etat() == ETATZERO) {
303 int taille = ti.get_taille() ;
304 for (int i=0 ; i<taille ; i++) {
305 to.t[i] = 0 ;
306 }
307 return to ;
308 }
309
310 // Cas general
311 assert(ti.get_etat() == ETATQCQ) ; // Otherwise
312 int taille = ti.get_taille() ;
313 for (int i=0 ; i<taille ; i++) {
314 if(ti.t[i] >= 0)
315 to.t[i] = 1 ;
316 else
317 to.t[i] = 0 ;
318 }
319 return to ;
320}
321
322 //-------------//
323 // Log naturel //
324 //-------------//
325
326Tbl log(const Tbl& ti)
327{
328 // Protection
329 assert(ti.get_etat() != ETATNONDEF) ;
330
331 // Cas ETATZERO
332 if (ti.get_etat() == ETATZERO) {
333 cout << "Tbl log: log(ETATZERO) !" << endl ;
334 abort () ;
335 }
336
337 // Cas general
338 assert(ti.get_etat() == ETATQCQ) ; // sinon...
339 Tbl to(ti.dim) ; // Tbl resultat
340 to.set_etat_qcq() ;
341 int taille = ti.get_taille() ;
342 for (int i=0 ; i<taille ; i++) {
343 to.t[i] = log(ti.t[i]) ;
344 }
345 return to ;
346}
347
348 //-------------//
349 // Log decimal //
350 //-------------//
351
353{
354 // Protection
355 assert(ti.get_etat() != ETATNONDEF) ;
356
357 // Cas ETATZERO
358 if (ti.get_etat() == ETATZERO) {
359 cout << "Tbl log10: log10(ETATZERO) !" << endl ;
360 abort () ;
361 }
362
363 // Cas general
364 assert(ti.get_etat() == ETATQCQ) ; // sinon...
365 Tbl to(ti.dim) ; // Tbl resultat
366 to.set_etat_qcq() ;
367 int taille = ti.get_taille() ;
368 for (int i=0 ; i<taille ; i++) {
369 to.t[i] = log10(ti.t[i]) ;
370 }
371 return to ;
372}
373
374 //--------------//
375 // Power entier //
376 //--------------//
377
378Tbl pow(const Tbl& ti, int n)
379{
380 // Protection
381 assert(ti.get_etat() != ETATNONDEF) ;
382
383 // Cas ETATZERO
384 if (ti.get_etat() == ETATZERO) {
385 if (n > 0) {
386 return ti ;
387 }
388 else {
389 cout << "Tbl pow: ETATZERO^n avec n<=0 ! "<< endl ;
390 abort () ;
391 }
392 }
393
394 // Cas general
395 assert(ti.get_etat() == ETATQCQ) ; // sinon...
396 Tbl to(ti.dim) ; // Tbl resultat
397 to.set_etat_qcq() ;
398 double x = n ;
399 int taille = ti.get_taille() ;
400 for (int i=0 ; i<taille ; i++) {
401 to.t[i] = pow(ti.t[i], x) ;
402 }
403 return to ;
404}
405
406 //--------------//
407 // Power double //
408 //--------------//
409
410Tbl pow(const Tbl& ti, double x)
411{
412 // Protection
413 assert(ti.get_etat() != ETATNONDEF) ;
414
415 // Cas ETATZERO
416 if (ti.get_etat() == ETATZERO) {
417 if (x > 0) {
418 return ti ;
419 }
420 else {
421 cout << "Tbl pow: ETATZERO^x avec x<=0 !" << endl ;
422 abort () ;
423 }
424 }
425
426 // Cas general
427 assert(ti.get_etat() == ETATQCQ) ; // sinon...
428 Tbl to(ti.dim) ; // Tbl resultat
429 to.set_etat_qcq() ;
430 int taille = ti.get_taille() ;
431 for (int i=0 ; i<taille ; i++) {
432 to.t[i] = pow(ti.t[i], x) ;
433 }
434 return to ;
435}
436
437 //----------------//
438 // Absolute value //
439 //----------------//
440
441Tbl abs(const Tbl& ti)
442{
443 // Protection
444 assert(ti.get_etat() != ETATNONDEF) ;
445
446 // Cas ETATZERO
447 if (ti.get_etat() == ETATZERO) {
448 return ti ;
449 }
450
451 // Cas general
452 assert(ti.get_etat() == ETATQCQ) ; // sinon...
453
454 Tbl to(ti.dim) ; // Tbl resultat
455 to.set_etat_qcq() ;
456
457 const double* xi = ti.t ;
458 double* xo = to.t ;
459 int taille = ti.get_taille() ;
460
461 for (int i=0 ; i<taille ; i++) {
462 xo[i] = fabs( xi[i] ) ;
463 }
464
465 return to ;
466}
467 //----------------//
468 // Cubic //
469 //----------------//
470
472{
473 // Protection
474 assert(ti.get_etat() != ETATNONDEF) ;
475
476 // Cas ETATZERO
477 if (ti.get_etat() == ETATZERO) {
478 return ti ;
479 }
480
481 // Cas general
482 assert(ti.get_etat() == ETATQCQ) ; // sinon...
483
484 Tbl absolute(abs(ti)) ;
485 Tbl res (pow(absolute, 1./3.)) ;
486
487 for (int i=0 ; i<ti.get_taille() ; i++)
488 if (ti.t[i] < 0)
489 res.t[i] *= -1 ;
490
491 return res ;
492}
493
494 //-------------------------------//
495 // max //
496 //-------------------------------//
497
498double max(const Tbl& ti) {
499
500 // Protection
501 assert(ti.get_etat() != ETATNONDEF) ;
502
503 // Cas particulier
504 if (ti.get_etat() == ETATZERO) {
505 return double(0) ;
506 }
507
508 // Cas general
509 assert(ti.get_etat() == ETATQCQ) ; // sinon....
510
511 const double* x = ti.t ;
512 double resu = x[0] ;
513 for (int i=1; i<ti.get_taille(); i++) {
514 if ( x[i] > resu ) resu = x[i] ;
515 }
516
517 return resu ;
518}
519
520 //-------------------------------//
521 // min //
522 //-------------------------------//
523
524double min(const Tbl& ti) {
525
526 // Protection
527 assert(ti.get_etat() != ETATNONDEF) ;
528
529 // Cas particulier
530 if (ti.get_etat() == ETATZERO) {
531 return double(0) ;
532 }
533
534 // Cas general
535 assert(ti.get_etat() == ETATQCQ) ; // sinon....
536
537 const double* x = ti.t ;
538 double resu = x[0] ;
539 for (int i=1; i<ti.get_taille(); i++) {
540 if ( x[i] < resu ) resu = x[i] ;
541 }
542
543 return resu ;
544}
545
546 //-------------------------------//
547 // norme //
548 //-------------------------------//
549
550double norme(const Tbl& ti) {
551
552 // Protection
553 assert(ti.get_etat() != ETATNONDEF) ;
554
555 double resu = 0 ;
556
557 if (ti.get_etat() != ETATZERO) { // on n'effectue la somme que si necessaire
558
559 assert(ti.get_etat() == ETATQCQ) ; // sinon....
560 const double* x = ti.t ;
561 for (int i=0; i<ti.get_taille(); i++) {
562 resu += fabs( x[i] ) ;
563 }
564
565 }
566
567 return resu ;
568}
569
570 //-------------------------------//
571 // diffrel //
572 //-------------------------------//
573
574double diffrel(const Tbl& t1, const Tbl& t2) {
575
576 // Protections
577 assert(t1.get_etat() != ETATNONDEF) ;
578 assert(t2.get_etat() != ETATNONDEF) ;
579
580 double norm2 = norme(t2) ;
581 double normdiff = norme(t1-t2) ;
582 double resu ;
583 if ( norm2 == double(0) ) {
584 resu = normdiff ;
585 }
586 else {
587 resu = normdiff / norm2 ;
588 }
589
590 return resu ;
591
592}
593
594 //-------------------------------//
595 // diffrelmax //
596 //-------------------------------//
597
598double diffrelmax(const Tbl& t1, const Tbl& t2) {
599
600 // Protections
601 assert(t1.get_etat() != ETATNONDEF) ;
602 assert(t2.get_etat() != ETATNONDEF) ;
603
604 double max2 = max(abs(t2)) ;
605 double maxdiff = max(abs(t1-t2)) ;
606 double resu ;
607 if ( max2 == double(0) ) {
608 resu = maxdiff ;
609 }
610 else {
611 resu = maxdiff / max2 ;
612 }
613
614 return resu ;
615
616}
617}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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
Lorene prototypes.
Definition app_hor.h:64