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