LORENE
map_af_integ_surf.C
1/*
2 * Copyright (c) 2000-2001 Philippe Grandclement
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * $Id: map_af_integ_surf.C,v 1.10 2018/11/16 14:34:35 j_novak Exp $
27 * $Log: map_af_integ_surf.C,v $
28 * Revision 1.10 2018/11/16 14:34:35 j_novak
29 * Changed minor points to avoid some compilation warnings.
30 *
31 * Revision 1.9 2017/02/24 16:50:27 j_novak
32 * *** empty log message ***
33 *
34 * Revision 1.8 2016/12/05 16:17:57 j_novak
35 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36 *
37 * Revision 1.7 2014/10/13 08:53:02 j_novak
38 * Lorene classes and functions now belong to the namespace Lorene.
39 *
40 * Revision 1.6 2014/10/06 15:13:12 j_novak
41 * Modified #include directives to use c++ syntax.
42 *
43 * Revision 1.5 2009/10/08 16:20:47 j_novak
44 * Addition of new bases T_COS and T_SIN.
45 *
46 * Revision 1.4 2007/10/05 15:56:19 j_novak
47 * Addition of new bases for the non-symmetric case in theta.
48 *
49 * Revision 1.3 2004/03/10 12:43:06 jl_jaramillo
50 * Treatment of case ETATUN in surface integrals for Scalar's.
51 *
52 * Revision 1.2 2004/01/29 08:50:03 p_grandclement
53 * Modification of Map::operator==(const Map&) and addition of the surface
54 * integrales using Scalar.
55 *
56 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
57 * LORENE
58 *
59 * Revision 1.7 2001/02/19 11:40:27 phil
60 * correction indices
61 *
62 * Revision 1.6 2001/02/12 14:14:29 phil
63 * *** empty log message ***
64 *
65 * Revision 1.5 2001/02/12 14:00:52 phil
66 * on prends tous les coefficients now
67 *
68 * Revision 1.4 2001/02/12 12:35:34 phil
69 * gestion des bases angulaires plus proprement
70 *
71 * Revision 1.3 2001/01/02 10:52:27 phil
72 * ajout calcul a l'infini
73 *
74 * Revision 1.2 2000/09/19 13:54:14 phil
75 * *** empty log message ***
76 *
77 * Revision 1.1 2000/09/19 13:09:32 phil
78 * Initial revision
79 *
80 *
81 * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_integ_surf.C,v 1.10 2018/11/16 14:34:35 j_novak Exp $
82 *
83 */
84
85
86#include <cstdlib>
87#include <cmath>
88
89#include "map.h"
90#include "cmp.h"
91#include "proto.h"
92#include "scalar.h"
93
94 //=============
95 // Cmp version
96 //===============
97
98namespace Lorene {
99double Map_af::integrale_surface (const Cmp& ci, double rayon) const {
100
101 assert (ci.get_etat() != ETATNONDEF) ;
102 assert (rayon > 0) ;
103 if (ci.get_etat() == ETATZERO)
104 return 0 ;
105
106 assert (ci.get_etat() == ETATQCQ) ;
107
108 int l ;
109 double xi ;
110 val_lx (rayon, 0, 0, l, xi) ;
111
112 if (l == get_mg()->get_nzone()-1) {
113 ci.check_dzpuis(0) ;
114 }
115
116 ci.va.coef() ;
117 int nr = get_mg()->get_nr(l) ;
118 int nt = get_mg()->get_nt(l) ;
119
120 int base_r = ci.va.base.get_base_r(l) ;
121 int base_t = ci.va.base.get_base_t(l) ;
122 int base_p = ci.va.base.get_base_p(l) ;
123
124 double result = 0 ;
125 double* coef = new double [nr] ;
126 double* auxi = new double[1] ;
127
128 bool odd_theta = false ;
129 double c_cos = 2 ;
130 switch (base_t) {
131 case T_COS_P : case T_COSSIN_CP :
132 break ;
133 case T_COS_I : case T_COSSIN_CI :
134 odd_theta = true ;
135 break ;
136 case T_COS : case T_COSSIN_C :
137 c_cos = 1. ;
138 break ;
139 default :
140 cout << "base_t cas non prevu dans Map_af::integrale_surface" << endl ;
141 abort() ;
142 break ;
143 }
144
145 if (!odd_theta) {
146 for (int j=0 ; j<nt ; j++) {
147 for (int i=0 ; i<nr ; i++)
148 coef[i] = (*ci.va.c_cf)(l, 0, j, i) ;
149
150 switch (base_r) {
151
152 case R_CHEB :
153 som_r_cheb (coef, nr, 1, 1, xi, auxi) ;
154 break ;
155 case R_CHEBP :
156 som_r_chebp (coef, nr, 1, 1, xi, auxi) ;
157 break ;
158 case R_CHEBI :
159 som_r_chebi (coef, nr, 1, 1, xi, auxi) ;
160 break ;
161 case R_CHEBU :
162 som_r_chebu (coef, nr, 1, 1, xi, auxi) ;
163 break ;
164 case R_CHEBPI_P :
165 som_r_chebpi_p (coef, nr, 1, 1, xi, auxi) ;
166 break ;
167 case R_CHEBPI_I :
168 som_r_chebpi_i (coef, nr, 1, 1, xi, auxi) ;
169 break ;
170 case R_CHEBPIM_P :
171 som_r_chebpim_p (coef, nr, 1, 1, xi, auxi) ;
172 break ;
173 case R_CHEBPIM_I :
174 som_r_chebpim_i (coef, nr, 1, 1, xi, auxi) ;
175 break ;
176 case R_LEG :
177 som_r_leg (coef, nr, 1, 1, xi, auxi) ;
178 break ;
179 case R_LEGP :
180 som_r_legp (coef, nr, 1, 1, xi, auxi) ;
181 break ;
182 case R_LEGI :
183 som_r_legi (coef, nr, 1, 1, xi, auxi) ;
184 break ;
185 default :
186 som_r_pas_prevu (coef, nr, 1, 1, xi, auxi) ;
187 break ;
188 }
189 result += 2 * (*auxi)/(1-c_cos*c_cos*j*j) ;
190 if (c_cos == 1.) j++ ;
191 }
192 }
193 delete [] auxi ;
194 delete [] coef ;
195
196 switch (base_p) {
197 case P_COSSIN :
198 result *= 2*rayon*rayon*M_PI ;
199 break ;
200 case P_COSSIN_P :
201 result *= 2*rayon*rayon*M_PI ;
202 break ;
203 default :
204 cout << "base_p cas non prevu dans Map_af::integrale_surface" << endl ;
205 abort() ;
206 break ;
207 }
208
209 return result ;
210}
211
212
213// Integrale a l'infini
214double Map_af::integrale_surface_infini (const Cmp& ci) const {
215
216 assert (ci.get_etat() != ETATNONDEF) ;
217 assert (ci.check_dzpuis(2));
218
219 if (ci.get_etat() == ETATZERO)
220 return 0 ;
221
222 assert (ci.get_etat() == ETATQCQ) ;
223
224 int nz = ci.get_mp()->get_mg()->get_nzone() ;
225
226 ci.va.coef() ;
227 int nr = get_mg()->get_nr(nz-1) ;
228 int nt = get_mg()->get_nt(nz-1) ;
229
230 int base_r = ci.va.base.get_base_r(nz-1) ;
231 int base_t = ci.va.base.get_base_t(nz-1) ;
232 int base_p = ci.va.base.get_base_p(nz-1) ;
233
234 double result = 0 ;
235 double* coef = new double [nr] ;
236 double* auxi = new double[1] ;
237
238 bool odd_theta = false ;
239 double c_cos = 2. ;
240 switch (base_t) {
241 case T_COS_P : case T_COSSIN_CP :
242 break ;
243 case T_COS_I : case T_COSSIN_CI :
244 odd_theta = true ;
245 break ;
246 case T_COS : case T_COSSIN_C :
247 c_cos = 1. ;
248 break ;
249 default :
250 cout << "base_t cas non prevu dans Map_af::integrale_surface" << endl ;
251 abort() ;
252 break ;
253 }
254
255 if (!odd_theta) {
256 for (int j=0 ; j<nt ; j++) {
257 for (int i=0 ; i<nr ; i++)
258 coef[i] = (*ci.va.c_cf)(nz-1, 0, j, i) ;
259
260 switch (base_r) {
261 case R_CHEBU :
262 som_r_chebu (coef, nr, 1, 1, 1, auxi) ;
263 break ;
264 default :
265 som_r_pas_prevu (coef, nr, 1, 1, 1, auxi) ;
266 break ;
267 }
268 result += 2 * (*auxi)/(1-c_cos*c_cos*j*j) ;
269 if (c_cos == 1.) j++ ;
270 }
271 }
272 delete [] auxi ;
273 delete [] coef ;
274
275 switch (base_p) {
276 case P_COSSIN :
277 result *= 2*M_PI ;
278 break ;
279 case P_COSSIN_P :
280 result *= 2*M_PI ;
281 break ;
282 default :
283 cout << "base_p cas non prevu dans Map_af::integrale_surface_infini" << endl ;
284 abort() ;
285 break ;
286 }
287
288 return result ;
289}
290
291 //=============
292 // Scalar version
293 //===============
294
295double Map_af::integrale_surface (const Scalar& ci, double rayon) const {
296
297 assert (ci.get_etat() != ETATNONDEF) ;
298 assert (rayon > 0) ;
299 if (ci.get_etat() == ETATZERO)
300 return 0 ;
301
302 assert ( (ci.get_etat() == ETATQCQ) || (ci.get_etat() == ETATUN) ) ;
303
304 int l ;
305 double xi ;
306 val_lx (rayon, 0, 0, l, xi) ;
307
308 if (l == get_mg()->get_nzone()-1) {
309 ci.check_dzpuis(0) ;
310 }
311
312 ci.get_spectral_va().coef() ;
313 int nr = get_mg()->get_nr(l) ;
314 int nt = get_mg()->get_nt(l) ;
315
316 int base_r = ci.get_spectral_va().base.get_base_r(l) ;
317 int base_t = ci.get_spectral_va().base.get_base_t(l) ;
318 int base_p = ci.get_spectral_va().base.get_base_p(l) ;
319
320 double result = 0 ;
321 double* coef = new double [nr] ;
322 double* auxi = new double[1] ;
323
324 bool odd_theta = false ;
325 double c_cos = 2. ;
326
327 switch (base_t) {
328 case T_COS_P : case T_COSSIN_CP :
329 break ;
330 case T_COS_I: case T_COSSIN_CI :
331 odd_theta = true ;
332 break ;
333 case T_COS : case T_COSSIN_C :
334 c_cos = 1. ;
335 break ;
336 default :
337 cout << "base_t cas non prevu dans Map_af::integrale_surface" << endl ;
338 abort() ;
339 break ;
340 }
341
342 if (!odd_theta) {
343 for (int j=0 ; j<nt ; j++) {
344 for (int i=0 ; i<nr ; i++)
345 coef[i] = (*ci.get_spectral_va().c_cf)(l, 0, j, i) ;
346
347 switch (base_r) {
348
349 case R_CHEB :
350 som_r_cheb (coef, nr, 1, 1, xi, auxi) ;
351 break ;
352 case R_CHEBP :
353 som_r_chebp (coef, nr, 1, 1, xi, auxi) ;
354 break ;
355 case R_CHEBI :
356 som_r_chebi (coef, nr, 1, 1, xi, auxi) ;
357 break ;
358 case R_CHEBU :
359 som_r_chebu (coef, nr, 1, 1, xi, auxi) ;
360 break ;
361 case R_CHEBPI_P :
362 som_r_chebpi_p (coef, nr, 1, 1, xi, auxi) ;
363 break ;
364 case R_CHEBPI_I :
365 som_r_chebpi_i (coef, nr, 1, 1, xi, auxi) ;
366 break ;
367 case R_CHEBPIM_P :
368 som_r_chebpim_p (coef, nr, 1, 1, xi, auxi) ;
369 break ;
370 case R_CHEBPIM_I :
371 som_r_chebpim_i (coef, nr, 1, 1, xi, auxi) ;
372 break ;
373 case R_LEG :
374 som_r_leg (coef, nr, 1, 1, xi, auxi) ;
375 break ;
376 case R_LEGP :
377 som_r_legp (coef, nr, 1, 1, xi, auxi) ;
378 break ;
379 case R_LEGI :
380 som_r_legi (coef, nr, 1, 1, xi, auxi) ;
381 break ;
382 default :
383 som_r_pas_prevu (coef, nr, 1, 1, xi, auxi) ;
384 break ;
385 }
386 result += 2 * (*auxi)/(1-c_cos*c_cos*j*j) ;
387 if (c_cos == 1.) j++ ;
388 }
389 }
390
391 delete [] auxi ;
392 delete [] coef ;
393
394 switch (base_p) {
395 case P_COSSIN :
396 result *= 2*rayon*rayon*M_PI ;
397 break ;
398 case P_COSSIN_P :
399 result *= 2*rayon*rayon*M_PI ;
400 break ;
401 default :
402 cout << "base_p cas non prevu dans Map_af::integrale_surface" << endl ;
403 abort() ;
404 break ;
405 }
406
407 return result ;
408}
409
410
411// Integrale a l'infini
412double Map_af::integrale_surface_infini (const Scalar& ci) const {
413
414 assert (ci.get_etat() != ETATNONDEF) ;
415 assert (ci.check_dzpuis(2));
416
417 if (ci.get_etat() == ETATZERO)
418 return 0 ;
419
420 assert ( (ci.get_etat() == ETATQCQ) || (ci.get_etat() == ETATUN) ) ;
421
422 int nz = ci.get_mp().get_mg()->get_nzone() ;
423
424 ci.get_spectral_va().coef() ;
425 int nr = get_mg()->get_nr(nz-1) ;
426 int nt = get_mg()->get_nt(nz-1) ;
427
428 int base_r = ci.get_spectral_va().base.get_base_r(nz-1) ;
429 int base_t = ci.get_spectral_va().base.get_base_t(nz-1) ;
430 int base_p = ci.get_spectral_va().base.get_base_p(nz-1) ;
431
432 double result = 0 ;
433 double* coef = new double [nr] ;
434 double* auxi = new double[1] ;
435
436 bool odd_theta = false ;
437 double c_cos = 2. ;
438
439 switch (base_t) {
440 case T_COS_P : case T_COSSIN_CP :
441 break ;
442 case T_COS_I: case T_COSSIN_CI :
443 odd_theta = true ;
444 break ;
445 case T_COS : case T_COSSIN_C :
446 c_cos = 1. ;
447 break ;
448 default :
449 cout << "base_t cas non prevu dans Map_af::integrale_surface" << endl ;
450 abort() ;
451 break ;
452 }
453
454 if (!odd_theta) {
455 for (int j=0 ; j<nt ; j++) {
456 for (int i=0 ; i<nr ; i++)
457 coef[i] = (*ci.get_spectral_va().c_cf)(nz-1, 0, j, i) ;
458
459 switch (base_r) {
460 case R_CHEBU :
461 som_r_chebu (coef, nr, 1, 1, 1, auxi) ;
462 break ;
463 default :
464 som_r_pas_prevu (coef, nr, 1, 1, 1, auxi) ;
465 break ;
466 }
467 result += 2 * (*auxi)/(1-c_cos*c_cos*j*j) ;
468 if (c_cos == 1.) j++ ;
469 }
470 }
471
472 delete [] auxi ;
473 delete [] coef ;
474
475 switch (base_p) {
476 case P_COSSIN :
477 result *= 2*M_PI ;
478 break ;
479 case P_COSSIN_P :
480 result *= 2*M_PI ;
481 break ;
482 default :
483 cout << "base_p cas non prevu dans Map_af::integrale_surface_infini" << endl ;
484 abort() ;
485 break ;
486 }
487
488 return result ;
489}
490}
Lorene::Base_val::get_base_p
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition base_val.h:425
Lorene::Base_val::get_base_r
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition base_val.h:403
Lorene::Base_val::get_base_t
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition base_val.h:414
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::check_dzpuis
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition cmp.C:718
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Lorene::Map_af::val_lx
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...
Definition map_af_radius.C:131
Lorene::Map_af::integrale_surface_infini
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
Definition map_af_integ_surf.C:214
Lorene::Map_af::integrale_surface
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Definition map_af_integ_surf.C:99
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:474
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:469
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
Lorene::Scalar::check_dzpuis
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition scalar.C:879
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:607
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:560
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:312
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:151
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:315
P_COSSIN_P
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition type_parite.h:247
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition type_parite.h:184
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition type_parite.h:180
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition type_parite.h:186
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition type_parite.h:178
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition type_parite.h:174
R_LEG
#define R_LEG
base de Legendre ordinaire (fin)
Definition type_parite.h:182
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition type_parite.h:200
T_COSSIN_CI
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition type_parite.h:212
P_COSSIN
#define P_COSSIN
dev. standart
Definition type_parite.h:245
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition type_parite.h:176
T_COSSIN_CP
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition type_parite.h:208
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
T_COS
#define T_COS
dev. cos seulement
Definition type_parite.h:196
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
T_COS_I
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition type_parite.h:204
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition type_parite.h:192
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition type_parite.h:172
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:874
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:777