LORENE
cmp_math.C
1/*
2 * Mathematical functions for the Cmp class.
3 *
4 * These functions are not member functions of the Cmp class.
5 *
6 */
7
8/*
9 * Copyright (c) 1999-2001 Eric Gourgoulhon
10 * Copyright (c) 1999-2001 Philippe Grandclement
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
30char cmp_math_C[] = "$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_math.C,v 1.3 2014/10/13 08:52:47 j_novak Exp $" ;
31
32/*
33 * $Id: cmp_math.C,v 1.3 2014/10/13 08:52:47 j_novak Exp $
34 * $Log: cmp_math.C,v $
35 * Revision 1.3 2014/10/13 08:52:47 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.2 2014/10/06 15:13:03 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
42 * LORENE
43 *
44 * Revision 2.2 1999/12/02 17:59:16 phil
45 * /
46 *
47 * Revision 2.1 1999/11/26 14:23:30 eric
48 * Modif commentaires.
49 *
50 * Revision 2.0 1999/11/15 14:13:05 eric
51 * *** empty log message ***
52 *
53 *
54 * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_math.C,v 1.3 2014/10/13 08:52:47 j_novak Exp $
55 *
56 */
57
58// Headers C
59#include <cmath>
60
61// Headers Lorene
62#include "cmp.h"
63
64 //-------//
65 // Sine //
66 //-------//
67
68namespace Lorene {
69Cmp sin(const Cmp& ci)
70{
71 // Protection
72 assert(ci.get_etat() != ETATNONDEF) ;
73
74 // Cas ETATZERO
75 if (ci.get_etat() == ETATZERO) {
76 return ci ;
77 }
78
79 // Cas general
80 assert(ci.get_etat() == ETATQCQ) ; // sinon...
81
82 Cmp co(ci.get_mp()) ; // result
83
84 co.set_etat_qcq() ;
85 co.va = sin( ci.va ) ;
86
87 return co ;
88}
89
90 //---------//
91 // Cosine //
92 //---------//
93
94Cmp cos(const Cmp& ci)
95{
96 // Protection
97 assert(ci.get_etat() != ETATNONDEF) ;
98
99 Cmp co(ci.get_mp()) ; // result
100
101 // Cas ETATZERO
102 if (ci.get_etat() == ETATZERO) {
103 co.set_etat_qcq() ;
104 co.va = double(1) ;
105 }
106 else {
107 // Cas general
108 assert(ci.get_etat() == ETATQCQ) ; // sinon...
109 co.set_etat_qcq() ;
110 co.va = cos( ci.va ) ;
111 }
112
113 return co ;
114}
115
116 //----------//
117 // Tangent //
118 //----------//
119
120Cmp tan(const Cmp& ci)
121{
122 // Protection
123 assert(ci.get_etat() != ETATNONDEF) ;
124
125 // Cas ETATZERO
126 if (ci.get_etat() == ETATZERO) {
127 return ci ;
128 }
129
130 // Cas general
131 assert(ci.get_etat() == ETATQCQ) ; // sinon...
132
133 Cmp co(ci.get_mp()) ; // result
134 co.set_etat_qcq() ;
135 co.va = tan( ci.va ) ;
136
137 return co ;
138}
139
140 //----------//
141 // Arcsine //
142 //----------//
143
144Cmp asin(const Cmp& ci)
145{
146 // Protection
147 assert(ci.get_etat() != ETATNONDEF) ;
148
149 // Cas ETATZERO
150 if (ci.get_etat() == ETATZERO) {
151 return ci ;
152 }
153
154 // Cas general
155 assert(ci.get_etat() == ETATQCQ) ; // sinon...
156
157 Cmp co(ci.get_mp()) ; // result
158
159 co.set_etat_qcq() ;
160 co.va = asin( ci.va ) ;
161
162 return co ;
163}
164
165 //-------------//
166 // Arccossine //
167 //-------------//
168
169Cmp acos(const Cmp& ci)
170{
171 // Protection
172 assert(ci.get_etat() != ETATNONDEF) ;
173
174 Cmp co(ci.get_mp()) ; // result
175
176 // Cas ETATZERO
177 if (ci.get_etat() == ETATZERO) {
178 co.set_etat_qcq() ;
179 co.va = double(0.5) * M_PI ;
180 }
181 else {
182 // Cas general
183 assert(ci.get_etat() == ETATQCQ) ; // sinon...
184 co.set_etat_qcq() ;
185 co.va = acos( ci.va ) ;
186 }
187
188 return co ;
189}
190
191 //-------------//
192 // Arctangent //
193 //-------------//
194
195Cmp atan(const Cmp& ci)
196{
197 // Protection
198 assert(ci.get_etat() != ETATNONDEF) ;
199
200 // Cas ETATZERO
201 if (ci.get_etat() == ETATZERO) {
202 return ci ;
203 }
204
205 // Cas general
206 assert(ci.get_etat() == ETATQCQ) ; // sinon...
207
208 Cmp co(ci.get_mp()) ; // result
209
210 co.set_etat_qcq() ;
211 co.va = atan( ci.va ) ;
212
213 return co ;
214}
215
216 //-------------//
217 // Square root //
218 //-------------//
219
220Cmp sqrt(const Cmp& ci)
221{
222 // Protection
223 assert(ci.get_etat() != ETATNONDEF) ;
224
225 // Cas ETATZERO
226 if (ci.get_etat() == ETATZERO) {
227 return ci ;
228 }
229
230 // Cas general
231 assert(ci.get_etat() == ETATQCQ) ; // sinon...
232
233 Cmp co(ci.get_mp()) ; // result
234
235 co.set_etat_qcq() ;
236 co.va = sqrt( ci.va ) ;
237
238 return co ;
239}
240
241 //-------------//
242 // Cube root //
243 //-------------//
244
246{
247 // Protection
248 assert(ci.get_etat() != ETATNONDEF) ;
249
250 // Cas ETATZERO
251 if (ci.get_etat() == ETATZERO) {
252 return ci ;
253 }
254
255 // Cas general
256 assert(ci.get_etat() == ETATQCQ) ; // sinon...
257
258 Cmp co(ci.get_mp()) ; // result
259
260 co.set_etat_qcq() ;
261 co.va = racine_cubique( ci.va ) ;
262
263 return co ;
264}
265
266 //--------------//
267 // Exponential //
268 //--------------//
269
270Cmp exp(const Cmp& ci)
271{
272 // Protection
273 assert(ci.get_etat() != ETATNONDEF) ;
274
275 Cmp co(ci.get_mp()) ; // result
276
277 // Cas ETATZERO
278 if (ci.get_etat() == ETATZERO) {
279 co.set_etat_qcq() ;
280 co.va = double(1) ;
281 }
282 else {
283 // Cas general
284 assert(ci.get_etat() == ETATQCQ) ; // sinon...
285 co.set_etat_qcq() ;
286 co.va = exp( ci.va ) ;
287 }
288
289 return co ;
290}
291
292 //---------------------//
293 // Neperian logarithm //
294 //---------------------//
295
296Cmp log(const Cmp& ci)
297{
298 // Protection
299 assert(ci.get_etat() != ETATNONDEF) ;
300
301 // Cas ETATZERO
302 if (ci.get_etat() == ETATZERO) {
303 cout << "Argument of log is ZERO in log(Cmp) !" << endl ;
304 abort() ;
305 }
306
307 // Cas general
308 assert(ci.get_etat() == ETATQCQ) ; // sinon...
309
310 Cmp co(ci.get_mp()) ; // result
311
312 co.set_etat_qcq() ;
313 co.va = log( ci.va ) ;
314
315 return co ;
316}
317
318 //---------------------//
319 // Decimal logarithm //
320 //---------------------//
321
322Cmp log10(const Cmp& ci)
323{
324 // Protection
325 assert(ci.get_etat() != ETATNONDEF) ;
326
327 // Cas ETATZERO
328 if (ci.get_etat() == ETATZERO) {
329 cout << "Argument of log10 is ZERO in log10(Cmp) !" << endl ;
330 abort() ;
331 }
332
333 // Cas general
334 assert(ci.get_etat() == ETATQCQ) ; // sinon...
335
336 Cmp co(ci.get_mp()) ; // result
337
338 co.set_etat_qcq() ;
339 co.va = log10( ci.va ) ;
340
341 return co ;
342}
343
344 //------------------//
345 // Power (integer) //
346 //------------------//
347
348Cmp pow(const Cmp& ci, int n)
349{
350 // Protection
351 assert(ci.get_etat() != ETATNONDEF) ;
352
353 // Cas ETATZERO
354 if (ci.get_etat() == ETATZERO) {
355 if (n > 0) {
356 return ci ;
357 }
358 else {
359 cout << "pow(Cmp, int) : ETATZERO^n with n <= 0 !" << endl ;
360 abort() ;
361 }
362 }
363
364 // Cas general
365 assert(ci.get_etat() == ETATQCQ) ; // sinon...
366
367 Cmp co(ci.get_mp()) ; // result
368
369 co.set_etat_qcq() ;
370 co.va = pow(ci.va, n) ;
371
372 return co ;
373}
374
375 //-----------------//
376 // Power (double) //
377 //-----------------//
378
379Cmp pow(const Cmp& ci, double x)
380{
381 // Protection
382 assert(ci.get_etat() != ETATNONDEF) ;
383
384 // Cas ETATZERO
385 if (ci.get_etat() == ETATZERO) {
386 if (x > double(0)) {
387 return ci ;
388 }
389 else {
390 cout << "pow(Cmp, double) : ETATZERO^x with x <= 0 !" << endl ;
391 abort() ;
392 }
393 }
394
395 // Cas general
396 assert(ci.get_etat() == ETATQCQ) ; // sinon...
397
398 Cmp co(ci.get_mp()) ; // result
399
400 co.set_etat_qcq() ;
401 co.va = pow(ci.va, x) ;
402
403 return co ;
404}
405
406 //-----------------//
407 // Absolute value //
408 //-----------------//
409
410Cmp abs(const Cmp& ci)
411{
412 // Protection
413 assert(ci.get_etat() != ETATNONDEF) ;
414
415 // Cas ETATZERO
416 if (ci.get_etat() == ETATZERO) {
417 return ci ;
418 }
419
420 // Cas general
421 assert(ci.get_etat() == ETATQCQ) ; // sinon...
422
423 Cmp co(ci.get_mp()) ; // result
424
425 co.set_etat_qcq() ;
426 co.va = abs( ci.va ) ;
427
428 return co ;
429}
430
431 //-------------------------------//
432 // max //
433 //-------------------------------//
434
435Tbl max(const Cmp& ci) {
436
437 // Protection
438 assert(ci.get_etat() != ETATNONDEF) ;
439
440 Tbl resu( ci.get_mp()->get_mg()->get_nzone() ) ;
441
442 if (ci.get_etat() == ETATZERO) {
443 resu.annule_hard() ;
444 }
445 else {
446 assert(ci.get_etat() == ETATQCQ) ;
447
448 resu = max( ci.va ) ; // max(Valeur)
449 }
450
451 return resu ;
452}
453
454 //-------------------------------//
455 // min //
456 //-------------------------------//
457
458Tbl min(const Cmp& ci) {
459
460 // Protection
461 assert(ci.get_etat() != ETATNONDEF) ;
462
463 Tbl resu( ci.get_mp()->get_mg()->get_nzone() ) ;
464
465 if (ci.get_etat() == ETATZERO) {
466 resu.annule_hard() ;
467 }
468 else {
469 assert(ci.get_etat() == ETATQCQ) ;
470
471 resu = min( ci.va ) ; // min(Valeur)
472 }
473
474 return resu ;
475}
476
477 //-------------------------------//
478 // norme //
479 //-------------------------------//
480
481Tbl norme(const Cmp& ci) {
482
483 // Protection
484 assert(ci.get_etat() != ETATNONDEF) ;
485
486 Tbl resu( ci.get_mp()->get_mg()->get_nzone() ) ;
487
488 if (ci.get_etat() == ETATZERO) {
489 resu.annule_hard() ;
490 }
491 else {
492 assert(ci.get_etat() == ETATQCQ) ;
493
494 resu = norme( ci.va ) ; // norme(Valeur)
495 }
496
497 return resu ;
498}
499
500 //-------------------------------//
501 // diffrel //
502 //-------------------------------//
503
504Tbl diffrel(const Cmp& c1, const Cmp& c2) {
505
506 // Protections
507 assert(c1.get_etat() != ETATNONDEF) ;
508 assert(c2.get_etat() != ETATNONDEF) ;
509
510 int nz = c1.get_mp()->get_mg()->get_nzone() ;
511 Tbl resu(nz) ;
512
513 Cmp diff = c1 - c2 ; // la compatibilite dzpuis est testee a ce niveau
514
515 Tbl normdiff = norme(diff) ;
516 Tbl norme2 = norme(c2) ;
517
518 assert(normdiff.get_etat() == ETATQCQ) ;
519 assert(norme2.get_etat() == ETATQCQ) ;
520
521 resu.set_etat_qcq() ;
522 for (int l=0; l<nz; l++) {
523 if ( norme2(l) == double(0) ) {
524 resu.set(l) = normdiff(l) ;
525 }
526 else{
527 resu.set(l) = normdiff(l) / norme2(l) ;
528 }
529 }
530
531 return resu ;
532
533}
534
535 //-------------------------------//
536 // diffrelmax //
537 //-------------------------------//
538
539Tbl diffrelmax(const Cmp& c1, const Cmp& c2) {
540
541 // Protections
542 assert(c1.get_etat() != ETATNONDEF) ;
543 assert(c2.get_etat() != ETATNONDEF) ;
544
545 int nz = c1.get_mp()->get_mg()->get_nzone() ;
546 Tbl resu(nz) ;
547
548 Tbl max2 = max(abs(c2)) ;
549
550 Cmp diff = c1 - c2 ; // la compatibilite dzpuis est testee a ce niveau
551
552 Tbl maxdiff = max(abs(diff)) ;
553
554 assert(maxdiff.get_etat() == ETATQCQ) ;
555 assert(max2.get_etat() == ETATQCQ) ;
556
557 resu.set_etat_qcq() ;
558 for (int l=0; l<nz; l++) {
559 if ( max2(l) == double(0) ) {
560 resu.set(l) = maxdiff(l) ;
561 }
562 else{
563 resu.set(l) = maxdiff(l) / max2(l) ;
564 }
565 }
566
567 return resu ;
568
569}
570
571}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition cmp.C:304
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Basic array class.
Definition tbl.h:161
int get_etat() const
Gives the logical state.
Definition tbl.h:394
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
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