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 char 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 
68 namespace Lorene {
69 Cmp 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 
94 Cmp 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 
120 Cmp 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 
144 Cmp 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 
169 Cmp 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 
195 Cmp 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 
220 Cmp 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 
270 Cmp 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 
296 Cmp 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 
322 Cmp 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 
348 Cmp 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 
379 Cmp 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 
410 Cmp 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 
435 Tbl 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 
458 Tbl 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 
481 Tbl 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 
504 Tbl 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 
539 Tbl 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
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
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