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