scalar_exp_filter.C
/*
* Function applying an exponential filter to a Scalar:
* sigma(n/N) = exp(alpha*(n/N)^(2p)). See scalar.h for documentation.
*/
/*
* Copyright (c) 2007 Jerome Novak
*
* 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 version 2
* as published by the Free Software Foundation.
*
* 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_exp_filter_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_exp_filter.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $" ;
/*
* $Id: scalar_exp_filter.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $
* $Log: scalar_exp_filter.C,v $
* Revision 1.5 2014/10/13 08:53:46 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.4 2014/10/06 15:16:15 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.3 2012/01/17 10:29:27 j_penner
* added two routines to handle generalized exponential filtering
*
* Revision 1.2 2007/10/31 10:50:16 j_novak
* Testing whether the coefficients are zero in a given domain.
*
* Revision 1.1 2007/10/31 10:33:13 j_novak
* Added exponential filters to smooth Gibbs-type phenomena.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_exp_filter.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $
*
*/
// C headers
#include <cassert>
#include <cmath>
// Lorene headers
#include "tensor.h"
#include "proto.h"
namespace Lorene {
void Scalar::exponential_filter_r(int lzmin, int lzmax, int p,
double alpha) {
assert(lzmin >= 0) ;
const Mg3d& mgrid = *mp->get_mg() ;
#ifndef NDEBUG
int nz = mgrid.get_nzone() ;
#endif
assert(lzmax < nz) ;
assert(etat != ETATNONDEF) ;
if (etat == ETATZERO) return ;
va.coef() ;
assert(va.c_cf != 0x0) ;
assert(alpha < 0.) ;
double alp = log(pow(10., alpha)) ;
for (int lz=lzmin; lz<=lzmax; lz++) {
if ((*va.c_cf)(lz).get_etat() == ETATQCQ)
for (int k=0; k<mgrid.get_np(lz); k++)
for (int j=0; j<mgrid.get_nt(lz); j++) {
int nr = mgrid.get_nr(lz) ;
for (int i=0; i<nr; i++) {
double eta = double(i)/double(nr) ;
va.c_cf->set(lz, k, j, i) *= exp(alp*pow(eta, 2*p)) ;
}
}
}
if (va.c != 0x0) {
delete va.c ;
va.c = 0x0 ;
}
va.del_deriv() ;
del_deriv() ;
return ;
}
void Scalar::sarra_filter_r(int lzmin, int lzmax, double p,
double alpha) {
assert(lzmin >= 0) ;
const Mg3d& mgrid = *mp->get_mg() ;
#ifndef NDEBUG
int nz = mgrid.get_nzone() ;
#endif
assert(lzmax < nz) ;
assert(etat != ETATNONDEF) ;
if (etat == ETATZERO) return ;
va.coef() ;
assert(va.c_cf != 0x0) ;
assert(alpha < 0.) ;
for (int lz=lzmin; lz<=lzmax; lz++) {
if ((*va.c_cf)(lz).get_etat() == ETATQCQ)
for (int k=0; k<mgrid.get_np(lz); k++)
for (int j=0; j<mgrid.get_nt(lz); j++) {
int nr = mgrid.get_nr(lz) ;
for (int i=0; i<nr; i++) {
double eta = double(i)/double(nr) ;
va.c_cf->set(lz, k, j, i) *= exp(alpha*pow(eta, p)) ;
}
}
}
if (va.c != 0x0) {
delete va.c ;
va.c = 0x0 ;
}
va.del_deriv() ;
del_deriv() ;
return ;
}
void exp_filter_r_all_domains( Scalar& ss, int p, double alpha ) {
int nz = ss.get_mp().get_mg()->get_nzone() ;
ss.exponential_filter_r(0, nz-1, p, alpha) ;
return ;
}
void Scalar::sarra_filter_r_all_domains( double p, double alpha ) {
int nz = get_mp().get_mg()->get_nzone() ;
sarra_filter_r(0, nz-1, p, alpha) ;
return ;
}
void Scalar::exponential_filter_ylm(int lzmin, int lzmax, int p,
double alpha) {
assert(lzmin >= 0) ;
const Mg3d& mgrid = *mp->get_mg() ;
#ifndef NDEBUG
int nz = mgrid.get_nzone() ;
#endif
assert(lzmax < nz) ;
assert(etat != ETATNONDEF) ;
if (etat == ETATZERO) return ;
double alp = log(pow(10., alpha)) ;
va.ylm() ;
assert(va.c_cf != 0x0) ;
const Base_val& base = va.base ;
int l_q, m_q, base_r ;
for (int lz=lzmin; lz<=lzmax; lz++)
if ((*va.c_cf)(lz).get_etat() == ETATQCQ) {
int np = mgrid.get_np(lz) ;
int nt = mgrid.get_nt(lz) ;
int nr = mgrid.get_nr(lz) ;
int lmax = base.give_lmax(mgrid, lz) ;
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if (nullite_plm(j, nt, k, np, base) == 1 ) {
double eta = double(l_q) / double(lmax) ;
double sigma = exp(alp*pow(eta, 2*p)) ;
for (int i=0; i<nr; i++)
va.c_cf->set(lz, k, j, i) *= sigma ;
}
}
}
va.ylm_i() ;
if (va.c != 0x0) {
delete va.c ;
va.c = 0x0 ;
}
va.del_deriv() ;
del_deriv() ;
return ;
}
void exp_filter_ylm_all_domains(Scalar& ss, int p, double alpha ) {
int nz = ss.get_mp().get_mg()->get_nzone() ;
ss.exponential_filter_ylm(0, nz-1, p, alpha) ;
return ;
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::give_lmax
int give_lmax(const Mg3d &mgrid, int lz) const
Returns the highest multipole for a given grid.
Definition base_val_quantum.C:294
Lorene::Base_val::give_quant_numbers
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Definition base_val_quantum.C:81
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::del_deriv
virtual void del_deriv() const
Logical destructor of the derivatives.
Definition scalar.C:287
Lorene::Scalar::sarra_filter_r
void sarra_filter_r(int lzmin, int lzmax, double p, double alpha=-1E-16)
Applies an exponential filter to the spectral coefficients in the radial direction.
Definition scalar_exp_filter.C:96
Lorene::Scalar::exp
friend Scalar exp(const Scalar &)
Exponential.
Definition scalar_math.C:323
Lorene::Scalar::exponential_filter_ylm
virtual void exponential_filter_ylm(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the angular directions.
Definition scalar_exp_filter.C:142
Lorene::Scalar::etat
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
Lorene::Scalar::log
friend Scalar log(const Scalar &)
Neperian logarithm.
Definition scalar_math.C:385
Lorene::Scalar::sarra_filter_r_all_domains
void sarra_filter_r_all_domains(double p, double alpha=1E-16)
Applies an exponential filter in radial direction in all domains for the case where p is a double (se...
Definition scalar_exp_filter.C:136
Lorene::Scalar::va
Valeur va
The numerical value of the Scalar.
Definition scalar.h:405
Lorene::Scalar::pow
friend Scalar pow(const Scalar &, int)
Power .
Definition scalar_math.C:451
Lorene::Scalar::exponential_filter_r
virtual void exponential_filter_r(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the radial direction.
Definition scalar_exp_filter.C:60
Lorene::exp_filter_ylm_all_domains
void exp_filter_ylm_all_domains(Scalar &ss, int p, double alpha=-16.)
Applies an exponential filter in angular directions in all domains.
Definition scalar_exp_filter.C:186
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene
Lorene prototypes.
Definition app_hor.h:64