scalar_math.C
/*
* Mathematical functions for the Scalar class.
*
* These functions are not member functions of the Scalar class.
*
*/
/*
* Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
*
* Copyright (c) 1999-2003 Eric Gourgoulhon
* Copyright (c) 1999-2001 Philippe Grandclement
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char scalar_math_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_math.C,v 1.6 2014/10/13 08:53:46 j_novak Exp $" ;
/*
* $Id: scalar_math.C,v 1.6 2014/10/13 08:53:46 j_novak Exp $
* $Log: scalar_math.C,v $
* Revision 1.6 2014/10/13 08:53:46 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.5 2014/10/06 15:16:16 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.4 2012/01/17 10:27:46 j_penner
* added a Heaviside function
*
* Revision 1.3 2003/10/10 15:57:29 j_novak
* Added the state one (ETATUN) to the class Scalar
*
* Revision 1.2 2003/10/01 13:04:44 e_gourgoulhon
* The method Tensor::get_mp() returns now a reference (and not
* a pointer) onto a mapping.
*
* Revision 1.1 2003/09/25 08:06:56 e_gourgoulhon
* First versions (use Cmp as intermediate quantities).
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_math.C,v 1.6 2014/10/13 08:53:46 j_novak Exp $
*
*/
// Headers C
#include <cmath>
// Headers Lorene
#include "tensor.h"
//-------//
// Sine //
//-------//
namespace Lorene {
Scalar sin(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
return ci ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
Scalar resu(ci.get_mp()) ;
resu = sin(double(1)) ;
return resu ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = sin( ci.va ) ;
return co ;
}
//---------//
// Cosine //
//---------//
Scalar cos(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
Scalar co(ci.get_mp()) ; // result
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
co.set_etat_qcq() ;
co.va = double(1) ;
}
else {
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
co = cos(double(1)) ;
}
else {
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
co.set_etat_qcq() ;
co.va = cos( ci.va ) ;
}
}
return co ;
}
//----------//
// Tangent //
//----------//
Scalar tan(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
return ci ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
Scalar resu(ci.get_mp()) ;
resu = tan(double(1)) ;
return resu ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = tan( ci.va ) ;
return co ;
}
//----------//
// Arcsine //
//----------//
Scalar asin(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
return ci ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
Scalar resu(ci.get_mp()) ;
resu = M_PI_2 ;
return resu ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = asin( ci.va ) ;
return co ;
}
//-------------//
// Arccosine //
//-------------//
Scalar acos(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
Scalar co(ci.get_mp()) ; // result
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
co.set_etat_qcq() ;
co.va = double(0.5) * M_PI ;
}
else {
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
co.set_etat_zero() ;
}
else {
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
co.set_etat_qcq() ;
co.va = acos( ci.va ) ;
}
}
return co ;
}
//-------------//
// Arctangent //
//-------------//
Scalar atan(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
return ci ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
Scalar resu(ci.get_mp()) ;
resu = 0.25*M_PI ;
return resu ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = atan( ci.va ) ;
return co ;
}
//-------------//
// Square root //
//-------------//
Scalar sqrt(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
return ci ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
return ci ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = sqrt( ci.va ) ;
return co ;
}
//-------------//
// Cube root //
//-------------//
Scalar racine_cubique(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
return ci ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
return ci ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = racine_cubique( ci.va ) ;
return co ;
}
//--------------//
// Exponential //
//--------------//
Scalar exp(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
Scalar co(ci.get_mp()) ; // result
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
co.set_etat_one() ;
}
else {
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
co.set_etat_qcq() ;
co = M_E ;
}
else {
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
co.set_etat_qcq() ;
co.va = exp( ci.va ) ;
}
}
return co ;
}
//---------------------//
// Heaviside Function //
//---------------------//
Scalar Heaviside(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
Scalar co(ci.get_mp()) ; // make output a copy, to ensure the same structure
// if input state is zero, return zero
if (ci.get_etat() == ETATZERO) {
co.set_etat_zero() ;
}
else {
// if input state is one, return one
if (ci.get_etat() == ETATUN) {
co.set_etat_one() ;
}
else {
// In general return the Heaviside function
assert(ci.get_etat() == ETATQCQ) ; // otherwise
co.set_etat_qcq() ;
co.va = Heaviside( ci.va ) ;
}
}
return co ;
}
//---------------------//
// Neperian logarithm //
//---------------------//
Scalar log(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
cout << "Argument of log is ZERO in log(Scalar) !" << endl ;
abort() ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
Scalar resu(ci.get_mp()) ;
resu.set_etat_zero() ;
return resu ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = log( ci.va ) ;
return co ;
}
//---------------------//
// Decimal logarithm //
//---------------------//
Scalar log10(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
cout << "Argument of log10 is ZERO in log10(Scalar) !" << endl ;
abort() ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
Scalar resu(ci.get_mp()) ;
resu.set_etat_zero() ;
return resu ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = log10( ci.va ) ;
return co ;
}
//------------------//
// Power (integer) //
//------------------//
Scalar pow(const Scalar& ci, int n)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
if (n > 0) {
return ci ;
}
else {
cout << "pow(Scalar, int) : ETATZERO^n with n <= 0 !" << endl ;
abort() ;
}
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
return ci ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = pow(ci.va, n) ;
return co ;
}
//-----------------//
// Power (double) //
//-----------------//
Scalar pow(const Scalar& ci, double x)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
if (x > double(0)) {
return ci ;
}
else {
cout << "pow(Scalar, double) : ETATZERO^x with x <= 0 !" << endl ;
abort() ;
}
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
return ci ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = pow(ci.va, x) ;
return co ;
}
//-----------------//
// Absolute value //
//-----------------//
Scalar abs(const Scalar& ci)
{
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Cas ETATZERO
if (ci.get_etat() == ETATZERO) {
return ci ;
}
// Cas ETATUN
if (ci.get_etat() == ETATUN) {
return ci ;
}
// Cas general
assert(ci.get_etat() == ETATQCQ) ; // sinon...
Scalar co(ci.get_mp()) ; // result
co.set_etat_qcq() ;
co.va = abs( ci.va ) ;
return co ;
}
//-------------------------------//
// totalmax //
//-------------------------------//
double totalmax(const Scalar& ci) {
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
double resu ;
if (ci.get_etat() == ETATZERO) {
resu = 0 ;
}
else {
if (ci.get_etat() == ETATUN) {
resu = 1 ;
}
else {
assert(ci.get_etat() == ETATQCQ) ;
resu = totalmax( ci.va ) ; // max(Valeur)
}
}
return resu ;
}
//-------------------------------//
// totalmin //
//-------------------------------//
double totalmin(const Scalar& ci) {
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
// Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
double resu ;
if (ci.get_etat() == ETATZERO) {
resu = 0 ;
}
else {
if (ci.get_etat() == ETATUN) {
resu = 1 ;
}
else {
assert(ci.get_etat() == ETATQCQ) ;
resu = totalmin( ci.va ) ; // min(Valeur)
}
}
return resu ;
}
//-------------------------------//
// max //
//-------------------------------//
Tbl max(const Scalar& ci) {
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
if (ci.get_etat() == ETATZERO) {
resu.annule_hard() ;
}
else {
if (ci.get_etat() == ETATUN) {
resu = 1 ;
}
else {
assert(ci.get_etat() == ETATQCQ) ;
resu = max( ci.va ) ; // max(Valeur)
}
}
return resu ;
}
//-------------------------------//
// min //
//-------------------------------//
Tbl min(const Scalar& ci) {
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
if (ci.get_etat() == ETATZERO) {
resu.annule_hard() ;
}
else {
if (ci.get_etat() == ETATUN) {
resu = 1 ;
}
else {
assert(ci.get_etat() == ETATQCQ) ;
resu = min( ci.va ) ; // min(Valeur)
}
}
return resu ;
}
//-------------------------------//
// norme //
//-------------------------------//
Tbl norme(const Scalar& ci) {
// Protection
assert(ci.get_etat() != ETATNONDEF) ;
Tbl resu( ci.get_mp().get_mg()->get_nzone() ) ;
if (ci.get_etat() == ETATZERO) {
resu.annule_hard() ;
}
else {
if (ci.get_etat() == ETATUN) {
resu = 1 ;
}
else {
assert(ci.get_etat() == ETATQCQ) ;
resu = norme( ci.va ) ; // norme(Valeur)
}
}
return resu ;
}
//-------------------------------//
// diffrel //
//-------------------------------//
Tbl diffrel(const Scalar& c1, const Scalar& c2) {
// Protections
assert(c1.get_etat() != ETATNONDEF) ;
assert(c2.get_etat() != ETATNONDEF) ;
int nz = c1.get_mp().get_mg()->get_nzone() ;
Tbl resu(nz) ;
Scalar diff = c1 - c2 ; // la compatibilite dzpuis est testee a ce niveau
Tbl normdiff = norme(diff) ;
Tbl norme2 = norme(c2) ;
assert(normdiff.get_etat() == ETATQCQ) ;
assert(norme2.get_etat() == ETATQCQ) ;
resu.set_etat_qcq() ;
for (int l=0; l<nz; l++) {
if ( norme2(l) == double(0) ) {
resu.set(l) = normdiff(l) ;
}
else{
resu.set(l) = normdiff(l) / norme2(l) ;
}
}
return resu ;
}
//-------------------------------//
// diffrelmax //
//-------------------------------//
Tbl diffrelmax(const Scalar& c1, const Scalar& c2) {
// Protections
assert(c1.get_etat() != ETATNONDEF) ;
assert(c2.get_etat() != ETATNONDEF) ;
int nz = c1.get_mp().get_mg()->get_nzone() ;
Tbl resu(nz) ;
Tbl max2 = max(abs(c2)) ;
Scalar diff = c1 - c2 ; // la compatibilite dzpuis est testee a ce niveau
Tbl maxdiff = max(abs(diff)) ;
assert(maxdiff.get_etat() == ETATQCQ) ;
assert(max2.get_etat() == ETATQCQ) ;
resu.set_etat_qcq() ;
for (int l=0; l<nz; l++) {
if ( max2(l) == double(0) ) {
resu.set(l) = maxdiff(l) ;
}
else{
resu.set(l) = maxdiff(l) / max2(l) ;
}
}
return resu ;
}
}
Lorene::Scalar::set_etat_one
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition scalar.C:334
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::Scalar
Scalar(const Map &mpi)
Constructor from mapping.
Definition scalar.C:204
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Scalar::va
Valeur va
The numerical value of the Scalar.
Definition scalar.h:405
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:372
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
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:195
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:322
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
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:481
Lorene::acos
Cmp acos(const Cmp &)
Arccosine.
Definition cmp_math.C:169
Lorene::asin
Cmp asin(const Cmp &)
Arcsine.
Definition cmp_math.C:144
Lorene::racine_cubique
Cmp racine_cubique(const Cmp &)
Cube root.
Definition cmp_math.C:245
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:458
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:410
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
Lorene::tan
Cmp tan(const Cmp &)
Tangent.
Definition cmp_math.C:120
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene::Heaviside
Mtbl Heaviside(const Mtbl &)
Heaviside function.
Definition mtbl_math.C:317
Lorene::totalmin
double totalmin(const Mtbl &)
Minimum value of the Mtbl elements in all domain.
Definition mtbl_math.C:522
Lorene::totalmax
double totalmax(const Mtbl &)
Maximum value of the Mtbl elements in all domains.
Definition mtbl_math.C:494
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:726