26char tbl_val_smooth_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Valencia/tbl_val_smooth.C,v 1.4 2014/10/13 08:53:49 j_novak Exp $" ;
54void radial_smoothing(
double* ,
const double* ,
int ,
double) ;
57void Tbl_val::smooth_atmosphere(
double atmosphere_thr) {
59 const Gval_spher* gspher =
dynamic_cast<const Gval_spher*
>(
gval) ;
60 assert(gspher != 0x0) ;
61 int ndim = gspher->get_ndim() ;
62 int fant = gspher->get_fantome() ;
67 radial_smoothing(
t, gspher->zr->t, nr, atmosphere_thr) ;
72 for (
int j=0; j<nt; j++)
73 radial_smoothing(
t+j*nr, gspher->zr->t, nr, atmosphere_thr) ;
79 for (
int j=0; j<nt; j++)
80 for (
int k=0; k<np; k++)
81 radial_smoothing(
t+k*nt*nr+j*nr, gspher->zr->t, nr, atmosphere_thr) ;
86 cerr <<
"Tbl_val::smooth_atmosphere : strange number of dimensions!"
95void radial_smoothing(
double* tab,
const double* rr,
int n,
double rho) {
97 assert((tab!= 0x0)&&(rr!=0x0)) ;
100 if (fabs(tab[n-1]) > rho)
103 double* t = tab + (n-1) ;
107 for (
int i=0; ((i<n)&&(atmos)); i++) {
108 if (atmos) atmos = ( fabs(*t) < rho) ;
111 jump = ( fabs(*t) > rho ) ;
116 if (indice == -1) return ;
117 int np = 2*(n-indice-2)/3 ;
118 int nm = indice / 100 + 3 ;
120 if (indice < n - np+1) {
125 int ileft = indice - nm + 2 ;
126 int iright = indice + np - 1 ;
127 double alpha = ( rr[ileft - 2] - rr[ileft - 1]) /
128 ( rr[ileft -1] - rr[ileft]) ;
129 double der_l = ( alpha*(alpha+2.)*tab[ileft]
130 - (1.+alpha)*(1.+alpha)*tab[ileft-1]
132 ( (1.+alpha)*(rr[ileft - 1] - rr[ileft - 2]) ) ;
133 double f_l = tab[ileft] ;
134 double f_r = tab[iright] ;
135 double tau = rr[ileft] - rr[iright] ;
136 double alp = der_l / (tau*tau) + 2.*(f_r - f_l)/(tau*tau*tau) ;
137 double bet = 0.5*(der_l + alp*tau*(rr[iright] - 3*rr[ileft])) / tau ;
138 for (
int i=ileft; i<iright; i++) {
139 tab[i] = f_r + (alp*rr[i]+bet)*(rr[i] - rr[iright])*(rr[i] - rr[iright]);
144 double alpha = ( rr[ileft - 2] - rr[ileft - 1]) /
145 ( rr[ileft -1] - rr[ileft]) ;
146 double der_l = ( alpha*(alpha+2.)*tab[ileft]
147 - (1.+alpha)*(1.+alpha)*tab[ileft-1]
149 ( (1.+alpha)*(rr[ileft - 1] - rr[ileft - 2]) ) ;
150 for (
int i=ileft; i<n; i++) {
151 tab[i] = tab[ileft] + (rr[i] - rr[ileft])*der_l ;
const Grille_val * gval
The Grille_val (cartesian or spherical) on which the array is defined.
int get_dim(int i) const
Gives the i th dimension (ie dim->dim[i] , without hidden cells)
double * t
The array of double at the nodes.