28char map_et_radius_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Map/map_et_radius.C,v 1.6 2014/10/13 08:53:05 j_novak Exp $" ;
79#include "utilitaires.h"
83double fonc_invr_map_et_noyau(
double,
const Param&) ;
84double fonc_invr_map_et_coq(
double,
const Param&) ;
85double fonc_invr_map_et_zec(
double,
const Param&) ;
92double Map_et::val_r(
int l,
double xi,
double theta,
double pphi)
const {
95 assert( l<mg->get_nzone() ) ;
104 double xi_2 = xi * xi ;
105 double xi_3 = xi * xi_2 ;
106 double a = xi_2 * xi_2 * (3. - 2.*xi_2) ;
107 double b = ( 2.5 - 1.5 * xi_2 ) * xi_3 ;
108 resu =
alpha[l] * ( xi + a * ftp + b * gtp ) +
beta[l] ;
114 double xm1 = xi - 1. ;
115 double xp1 = xi + 1. ;
116 double a = 0.25* xm1 * xm1 * (xi + 2.) ;
117 double b = 0.25* xp1 * xp1 * (2. - xi) ;
118 resu =
alpha[l] * ( xi + a * ftp + b * gtp ) +
beta[l] ;
123 double xm1 = xi - 1. ;
124 double a = 0.25* xm1 * xm1 * (xi + 2.) ;
125 resu = double(1) / (
alpha[l] * ( xi + a * ftp ) +
beta[l] ) ;
130 cout <<
"Map_et::val_r: unknown type_r ! " << endl ;
147 int& lz,
double& xi)
const {
151 double precis = 1e-14 ;
160 val_lx(rr, theta, pphi, par, lz, xi) ;
170 const Param& par,
int& lz,
double& xi)
const {
184 if (rr == __infinity) {
192 cout.setf(ios::showpoint);
193 cout <<
"Map_et::val_lx: the domain containing r = " << rr <<
194 " has not been found ! "
206 double ftp = double(0) ;
207 double gtp = double(0) ;
208 for (
int l=0; l<nz; l++) {
215 rmax =
alpha[l] * ( double(1) + ftp + gtp ) +
beta[l] ;
222 rmax =
alpha[l] * ( double(1) + gtp ) +
beta[l] ;
229 rmax = double(1) / (
alpha[l] +
beta[l] ) ;
234 cout <<
"Map_et::val_lx: unknown type_r ! " << endl ;
240 if ( rr <= rmax + 1.e-14 ) {
242 if (ftp ==
double(0)) ftp =
ff.
val_point(l, 0, theta, pphi) ;
243 if (gtp ==
double(0)) gtp =
gg.
val_point(l, 0, theta, pphi) ;
250 cout.setf(ios::showpoint);
251 cout <<
"Map_et::val_lx: the domain containing r = " << rr <<
252 " has not been found ! "
254 for (
int l=0; l<nz; l++) {
259 rmax =
alpha[l] * ( double(1) + ftp + gtp ) +
beta[l] ;
263 rmax =
alpha[l] * ( double(1) + gtp ) +
beta[l] ;
267 rmax = double(1) / (
alpha[l] +
beta[l] ) ;
271 cout <<
"Map_et::val_lx: unknown type_r ! " << endl ;
276 cout <<
"domain: " << l <<
" theta = " << theta <<
" phi = "
277 << pphi <<
" : rmax = " << rmax << endl ;
285 if ( (rr >= rmax) && ( rr <= rmax + 1.e-14) ) {
310 xi =
zerosec(fonc_invr_map_et_noyau, parzerosec, xmin, xmax,
311 precis, nitermax, niter) ;
323 xi =
zerosec(fonc_invr_map_et_coq, parzerosec, xmin, xmax,
324 precis, nitermax, niter) ;
331 xi = ( double(1)/rr -
beta[lz] ) /
alpha[lz] ;
336 xi = ( double(1) / rr -
beta[lz] ) /
alpha[lz] ;
341 xi =
zerosec(fonc_invr_map_et_zec, parzerosec, xmin, xmax,
342 precis, nitermax, niter) ;
349 cout <<
"Map_et::val_lx: unknown type_r ! " << endl ;
370 assert( l<mg->get_nzone() ) ;
373 double ftp =
ff(l, k, j, 0) ;
378 double gtp =
gg(l, k, j, 0) ;
379 double xi_2 = xi * xi ;
380 double xi_3 = xi * xi_2 ;
381 double a = xi_2 * xi_2 * (3. - 2.*xi_2) ;
382 double b = ( 2.5 - 1.5 * xi_2 ) * xi_3 ;
383 resu =
alpha[l] * ( xi + a * ftp + b * gtp ) +
beta[l] ;
388 double gtp =
gg(l, k, j, 0) ;
389 double xm1 = xi - 1. ;
390 double xp1 = xi + 1. ;
391 double a = 0.25* xm1 * xm1 * (xi + 2.) ;
392 double b = 0.25* xp1 * xp1 * (2. - xi) ;
393 resu =
alpha[l] * ( xi + a * ftp + b * gtp ) +
beta[l] ;
398 double xm1 = xi - 1. ;
399 double a = 0.25* xm1 * xm1 * (xi + 2.) ;
400 resu = double(1) / (
alpha[l] * ( xi + a * ftp ) +
beta[l] ) ;
405 cout <<
"Map_et::val_r_jk: unknown type_r ! " << endl ;
419 int& lz,
double& xi)
const {
433 if (rr == __infinity) {
441 cout.setf(ios::showpoint);
442 cout <<
"Map_et::val_lx_jk: the domain containing r = " << rr <<
443 " has not been found ! "
455 double ftp = double(0) ;
456 double gtp = double(0) ;
457 for (
int l=0; l<nz; l++) {
462 ftp =
ff(l, k, j, 0) ;
463 gtp =
gg(l, k, j, 0) ;
464 rmax =
alpha[l] * ( double(1) + ftp + gtp ) +
beta[l] ;
470 gtp =
gg(l, k, j, 0) ;
471 rmax =
alpha[l] * ( double(1) + gtp ) +
beta[l] ;
478 rmax = double(1) / (
alpha[l] +
beta[l] ) ;
483 cout <<
"Map_et::val_lx_jk: unknown type_r ! " << endl ;
489 if ( rr <= rmax + 1.e-14 ) {
491 if (ftp ==
double(0)) ftp =
ff(l, k, j, 0) ;
492 if (gtp ==
double(0)) gtp =
gg(l, k, j, 0) ;
499 cout.setf(ios::showpoint);
500 cout <<
"Map_et::val_lx_jk: the domain containing r = " << rr <<
501 " has not been found ! "
503 for (
int l=0; l<nz; l++) {
504 ftp =
ff(l, k, j, 0) ;
505 gtp =
gg(l, k, j, 0) ;
508 rmax =
alpha[l] * ( double(1) + ftp + gtp ) +
beta[l] ;
512 rmax =
alpha[l] * ( double(1) + gtp ) +
beta[l] ;
516 rmax = double(1) / (
alpha[l] +
beta[l] ) ;
520 cout <<
"Map_et::val_lx_jk: unknown type_r ! " << endl ;
525 cout <<
"domain: " << l <<
" j = " << j <<
" k = " << k
526 <<
" : rmax = " << rmax << endl ;
534 if ( (rr >= rmax) && ( rr <= rmax + 1.e-14) ) {
559 xi =
zerosec(fonc_invr_map_et_noyau, parzerosec, xmin, xmax,
560 precis, nitermax, niter) ;
572 xi =
zerosec(fonc_invr_map_et_coq, parzerosec, xmin, xmax,
573 precis, nitermax, niter) ;
580 xi = ( double(1)/rr -
beta[lz] ) /
alpha[lz] ;
585 xi = ( double(1) / rr -
beta[lz] ) /
alpha[lz] ;
590 xi =
zerosec(fonc_invr_map_et_zec, parzerosec, xmin, xmax,
591 precis, nitermax, niter) ;
598 cout <<
"Map_et::val_lx_jk: unknown type_r ! " << endl ;
612double fonc_invr_map_et_noyau(
double x,
const Param& par) {
620 double x_3 = x_2 * x ;
621 double a = x_2 * x_2 * (3. - 2.*x_2) ;
622 double b = ( 2.5 - 1.5 * x_2 ) * x_3 ;
624 return alp * ( x + a * f + b * g ) + bet - r ;
630double fonc_invr_map_et_coq(
double x,
const Param& par) {
632 double r = par.get_double(0) ;
633 double f = par.get_double(1) ;
634 double g = par.get_double(2) ;
635 double alp = par.get_double(3) ;
636 double bet = par.get_double(4) ;
637 double xm1 = x - 1. ;
638 double xp1 = x + 1. ;
639 double a = 0.25* xm1 * xm1 * (x + 2.) ;
640 double b = 0.25* xp1 * xp1 * (2. - x) ;
642 return alp * ( x + a * f + b * g ) + bet - r ;
648double fonc_invr_map_et_zec(
double x,
const Param& par) {
650 double r = par.get_double(0) ;
651 double f = par.get_double(1) ;
652 double alp = par.get_double(3) ;
653 double bet = par.get_double(4) ;
654 double xm1 = x - 1. ;
655 double a = 0.25* xm1 * xm1 * (x + 2.) ;
657 return alp * ( x + a * f ) + bet -
double(1) / r ;
double * beta
Array (size: mg->nzone ) of the values of in each domain.
virtual void val_lx_jk(double rr, int j, int k, const Param &par, int &l, double &xi) const
Computes the domain index l and the value of corresponding to a point of arbitrary r but collocation...
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
Valeur ff
Values of the function at the nt*np angular collocation points in each domain.
virtual void val_lx(double rr, double theta, double pphi, int &l, double &xi) const
Computes the domain index l and the value of corresponding to a point given by its physical coordina...
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Valeur gg
Values of the function at the nt*np angular collocation points in each domain.
virtual double val_r_jk(int l, double xi, int j, int k) const
Returns the value of the radial coordinate r for a given and a given collocation point in in a give...
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
int get_nzone() const
Returns the number of domains.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
int get_etat() const
Gives the logical state.
double val_point(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
int get_etat() const
Returns the logical state.
Mtbl * c
Values of the function at the points of the multi-grid
double zerosec(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter, bool abort=true)
Finding the zero a function.