LORENE
scalar_exp_filter.C
1 /*
2  * Function applying an exponential filter to a Scalar:
3  * sigma(n/N) = exp(alpha*(n/N)^(2p)). See scalar.h for documentation.
4  */
5 
6 /*
7  * Copyright (c) 2007 Jerome Novak
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License version 2
13  * as published by the Free Software Foundation.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 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 $" ;
27 
28 /*
29  * $Id: scalar_exp_filter.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $
30  * $Log: scalar_exp_filter.C,v $
31  * Revision 1.5 2014/10/13 08:53:46 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.4 2014/10/06 15:16:15 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.3 2012/01/17 10:29:27 j_penner
38  * added two routines to handle generalized exponential filtering
39  *
40  * Revision 1.2 2007/10/31 10:50:16 j_novak
41  * Testing whether the coefficients are zero in a given domain.
42  *
43  * Revision 1.1 2007/10/31 10:33:13 j_novak
44  * Added exponential filters to smooth Gibbs-type phenomena.
45  *
46  *
47  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_exp_filter.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $
48  *
49  */
50 
51 // C headers
52 #include <cassert>
53 #include <cmath>
54 
55 // Lorene headers
56 #include "tensor.h"
57 #include "proto.h"
58 
59 namespace Lorene {
60 void Scalar::exponential_filter_r(int lzmin, int lzmax, int p,
61  double alpha) {
62  assert(lzmin >= 0) ;
63  const Mg3d& mgrid = *mp->get_mg() ;
64 #ifndef NDEBUG
65  int nz = mgrid.get_nzone() ;
66 #endif
67  assert(lzmax < nz) ;
68  assert(etat != ETATNONDEF) ;
69  if (etat == ETATZERO) return ;
70  va.coef() ;
71  assert(va.c_cf != 0x0) ;
72  assert(alpha < 0.) ;
73  double alp = log(pow(10., alpha)) ;
74 
75  for (int lz=lzmin; lz<=lzmax; lz++) {
76  if ((*va.c_cf)(lz).get_etat() == ETATQCQ)
77  for (int k=0; k<mgrid.get_np(lz); k++)
78  for (int j=0; j<mgrid.get_nt(lz); j++) {
79  int nr = mgrid.get_nr(lz) ;
80  for (int i=0; i<nr; i++) {
81  double eta = double(i)/double(nr) ;
82  va.c_cf->set(lz, k, j, i) *= exp(alp*pow(eta, 2*p)) ;
83  }
84  }
85  }
86  if (va.c != 0x0) {
87  delete va.c ;
88  va.c = 0x0 ;
89  }
90  va.del_deriv() ;
91  del_deriv() ;
92 
93  return ;
94 }
95 
96 void Scalar::sarra_filter_r(int lzmin, int lzmax, double p,
97  double alpha) {
98  assert(lzmin >= 0) ;
99  const Mg3d& mgrid = *mp->get_mg() ;
100 #ifndef NDEBUG
101  int nz = mgrid.get_nzone() ;
102 #endif
103  assert(lzmax < nz) ;
104  assert(etat != ETATNONDEF) ;
105  if (etat == ETATZERO) return ;
106  va.coef() ;
107  assert(va.c_cf != 0x0) ;
108  assert(alpha < 0.) ;
109 
110  for (int lz=lzmin; lz<=lzmax; lz++) {
111  if ((*va.c_cf)(lz).get_etat() == ETATQCQ)
112  for (int k=0; k<mgrid.get_np(lz); k++)
113  for (int j=0; j<mgrid.get_nt(lz); j++) {
114  int nr = mgrid.get_nr(lz) ;
115  for (int i=0; i<nr; i++) {
116  double eta = double(i)/double(nr) ;
117  va.c_cf->set(lz, k, j, i) *= exp(alpha*pow(eta, p)) ;
118  }
119  }
120  }
121  if (va.c != 0x0) {
122  delete va.c ;
123  va.c = 0x0 ;
124  }
125  va.del_deriv() ;
126  del_deriv() ;
127 
128  return ;
129 }
130 void exp_filter_r_all_domains( Scalar& ss, int p, double alpha ) {
131  int nz = ss.get_mp().get_mg()->get_nzone() ;
132  ss.exponential_filter_r(0, nz-1, p, alpha) ;
133  return ;
134 }
135 
136 void Scalar::sarra_filter_r_all_domains( double p, double alpha ) {
137  int nz = get_mp().get_mg()->get_nzone() ;
138  sarra_filter_r(0, nz-1, p, alpha) ;
139  return ;
140 }
141 
142 void Scalar::exponential_filter_ylm(int lzmin, int lzmax, int p,
143  double alpha) {
144  assert(lzmin >= 0) ;
145  const Mg3d& mgrid = *mp->get_mg() ;
146 #ifndef NDEBUG
147  int nz = mgrid.get_nzone() ;
148 #endif
149  assert(lzmax < nz) ;
150  assert(etat != ETATNONDEF) ;
151  if (etat == ETATZERO) return ;
152  double alp = log(pow(10., alpha)) ;
153  va.ylm() ;
154  assert(va.c_cf != 0x0) ;
155  const Base_val& base = va.base ;
156  int l_q, m_q, base_r ;
157 
158  for (int lz=lzmin; lz<=lzmax; lz++)
159  if ((*va.c_cf)(lz).get_etat() == ETATQCQ) {
160  int np = mgrid.get_np(lz) ;
161  int nt = mgrid.get_nt(lz) ;
162  int nr = mgrid.get_nr(lz) ;
163  int lmax = base.give_lmax(mgrid, lz) ;
164  for (int k=0; k<np; k++)
165  for (int j=0; j<nt; j++) {
166  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
167  if (nullite_plm(j, nt, k, np, base) == 1 ) {
168  double eta = double(l_q) / double(lmax) ;
169  double sigma = exp(alp*pow(eta, 2*p)) ;
170  for (int i=0; i<nr; i++)
171  va.c_cf->set(lz, k, j, i) *= sigma ;
172  }
173  }
174  }
175 
176  va.ylm_i() ;
177  if (va.c != 0x0) {
178  delete va.c ;
179  va.c = 0x0 ;
180  }
181  va.del_deriv() ;
182  del_deriv() ;
183  return ;
184 }
185 
186 void exp_filter_ylm_all_domains(Scalar& ss, int p, double alpha ) {
187  int nz = ss.get_mp().get_mg()->get_nzone() ;
188  ss.exponential_filter_ylm(0, nz-1, p, alpha) ;
189  return ;
190 }
191 
192 }
Bases of the spectral expansions.
Definition: base_val.h:322
int give_lmax(const Mg3d &mgrid, int lz) const
Returns the highest multipole for a given grid.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Multi-domain grid.
Definition: grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
friend Scalar exp(const Scalar &)
Exponential.
Definition: scalar_math.C:323
virtual void del_deriv() const
Logical destructor of the derivatives.
Definition: scalar.C:287
friend Scalar log(const Scalar &)
Neperian logarithm.
Definition: scalar_math.C:385
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.
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.
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition: scalar.h:396
friend Scalar pow(const Scalar &, int)
Power .
Definition: scalar_math.C:451
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...
Valeur va
The numerical value of the Scalar
Definition: scalar.h:405
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.
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Mtbl * c
Values of the function at the points of the multi-grid
Definition: valeur.h:299
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
void del_deriv()
Logical destructor of the derivatives.
Definition: valeur.C:638
void exp_filter_ylm_all_domains(Scalar &ss, int p, double alpha=-16.)
Applies an exponential filter in angular directions in all domains.
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
Lorene prototypes.
Definition: app_hor.h:64