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 
71 namespace Lorene {
72 Cmp 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 
97 Cmp 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 
123 Cmp 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 
147 Cmp 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 
172 Cmp 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 
198 Cmp 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 
223 Cmp 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 
248 Cmp 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 
273 Cmp 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 
299 Cmp 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 
325 Cmp 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 
351 Cmp 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 
382 Cmp 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 
413 Cmp 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 
438 Tbl 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 
461 Tbl 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 
484 Tbl 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 
507 Tbl 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 
542 Tbl 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::get_mp
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Lorene::asin
Cmp asin(const Cmp &)
Arcsine.
Definition: cmp_math.C:147
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Lorene::racine_cubique
Cmp racine_cubique(const Cmp &)
Cube root.
Definition: cmp_math.C:248
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition: tbl.h:414
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:461
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Lorene::tan
Cmp tan(const Cmp &)
Tangent.
Definition: cmp_math.C:123
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::atan
Cmp atan(const Cmp &)
Arctangent.
Definition: cmp_math.C:198
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
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::Cmp::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:307
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
Lorene::acos
Cmp acos(const Cmp &)
Arccosine.
Definition: cmp_math.C:172
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:413
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:72
Lorene::Tbl
Basic array class.
Definition: tbl.h:164
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542