tbl_val_smooth.C
/*
* Smoothes the junction with an eventual atmosphere.
*
*/
/*
* Copyright (c) 2004 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 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 $" ;
/*
* $Id: tbl_val_smooth.C,v 1.4 2014/10/13 08:53:49 j_novak Exp $
* $Log: tbl_val_smooth.C,v $
* Revision 1.4 2014/10/13 08:53:49 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.3 2004/12/30 16:14:01 j_novak
* Changed the name of a shadowed variable.
*
* Revision 1.2 2004/12/03 13:24:01 j_novak
* Minor modif.
*
* Revision 1.1 2004/11/26 17:02:19 j_novak
* Added a function giving a smooth transition to the atmosphere.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Valencia/tbl_val_smooth.C,v 1.4 2014/10/13 08:53:49 j_novak Exp $
*
*/
// Lorene headers
#include "tbl_val.h"
//Local prototypes
namespace Lorene {
void radial_smoothing(double* , const double* , int , double) ;
//****************************************************************************
void Tbl_val::smooth_atmosphere(double atmosphere_thr) {
const Gval_spher* gspher = dynamic_cast<const Gval_spher*>(gval) ;
assert(gspher != 0x0) ;
int ndim = gspher->get_ndim() ;
int fant = gspher->get_fantome() ;
int nr = get_dim(0) + 2*fant;
switch (ndim) {
case 1: {
radial_smoothing(t, gspher->zr->t, nr, atmosphere_thr) ;
break ;
}
case 2: {
int nt = get_dim(1) + 2*fant ;
for (int j=0; j<nt; j++)
radial_smoothing(t+j*nr, gspher->zr->t, nr, atmosphere_thr) ;
break ;
}
case 3: {
int nt = get_dim(1) + 2*fant ;
int np = get_dim(2) + 2*fant ;
for (int j=0; j<nt; j++)
for (int k=0; k<np; k++)
radial_smoothing(t+k*nt*nr+j*nr, gspher->zr->t, nr, atmosphere_thr) ;
break ;
}
default: {
cerr << "Tbl_val::smooth_atmosphere : strange number of dimensions!"
<< endl ;
abort() ;
break ;
}
}
return ;
}
void radial_smoothing(double* tab, const double* rr, int n, double rho) {
assert((tab!= 0x0)&&(rr!=0x0)) ;
assert (rho >= 0.) ;
if (fabs(tab[n-1]) > rho) // no atmosphere here
return ;
double* t = tab + (n-1) ;
int indice = -1 ;
bool atmos = true ;
bool jump = false ;
for (int i=0; ((i<n)&&(atmos)); i++) {
if (atmos) atmos = ( fabs(*t) < rho) ;
t-- ;
if (atmos) {
jump = ( fabs(*t) > rho ) ;
if (jump) // discontinuity found
indice = n - i - 2 ;
}
}
if (indice == -1) return ;
int np = 2*(n-indice-2)/3 ;
int nm = indice / 100 + 3 ;
assert(n > nm+np) ;
if (indice < n - np+1) { // enough points to interpolate
// The inteprolation is done using a cubic polynomial
//---------------------------------------------------
int ileft = indice - nm + 2 ;
int iright = indice + np - 1 ;
double alpha = ( rr[ileft - 2] - rr[ileft - 1]) /
( rr[ileft -1] - rr[ileft]) ;
double der_l = ( alpha*(alpha+2.)*tab[ileft]
- (1.+alpha)*(1.+alpha)*tab[ileft-1]
+ tab[ileft-2] ) /
( (1.+alpha)*(rr[ileft - 1] - rr[ileft - 2]) ) ;
double f_l = tab[ileft] ;
double f_r = tab[iright] ;
double tau = rr[ileft] - rr[iright] ;
double alp = der_l / (tau*tau) + 2.*(f_r - f_l)/(tau*tau*tau) ;
double bet = 0.5*(der_l + alp*tau*(rr[iright] - 3*rr[ileft])) / tau ;
for (int i=ileft; i<iright; i++) {
tab[i] = f_r + (alp*rr[i]+bet)*(rr[i] - rr[iright])*(rr[i] - rr[iright]);
}
}
else { // too few points to interpolate -> linear extrapolation ...
int ileft = indice ;
double alpha = ( rr[ileft - 2] - rr[ileft - 1]) /
( rr[ileft -1] - rr[ileft]) ;
double der_l = ( alpha*(alpha+2.)*tab[ileft]
- (1.+alpha)*(1.+alpha)*tab[ileft-1]
+ tab[ileft-2] ) /
( (1.+alpha)*(rr[ileft - 1] - rr[ileft - 2]) ) ;
for (int i=ileft; i<n; i++) {
tab[i] = tab[ileft] + (rr[i] - rr[ileft])*der_l ;
}
}
return ;
}
}
Lorene::Tbl_val::gval
const Grille_val * gval
The Grille_val (cartesian or spherical) on which the array is defined.
Definition tbl_val.h:110
Lorene::Tbl_val::get_dim
int get_dim(int i) const
Gives the i th dimension (ie dim->dim[i] , without hidden cells).
Definition tbl_val.h:485
Lorene::Tbl_val::t
double * t
The array of double at the nodes.
Definition tbl_val.h:114
Lorene
Lorene prototypes.
Definition app_hor.h:64