LORENE
valeur_ylm_i.C
1 /*
2  * Copyright (c) 1999-2003 Eric Gourgoulhon
3  * Copyright (c) 1999-2001 Philippe Grandclement
4  *
5  * This file is part of LORENE.
6  *
7  * LORENE is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * LORENE is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with LORENE; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20  *
21  */
22 
23 
24 char valeur_ylm_i_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_ylm_i.C,v 1.12 2014/10/13 08:53:51 j_novak Exp $" ;
25 
26 /*
27  * Fonction membre de la classe Valeur qui calcule les coefficients
28  * de la decomposition en cos(l*theta) / sin(l*theta)
29  * a partir des coefficients de la decomp. en harmoniques spheriques
30  *
31  */
32 
33 
34 /*
35  * $Id: valeur_ylm_i.C,v 1.12 2014/10/13 08:53:51 j_novak Exp $
36  * $Log: valeur_ylm_i.C,v $
37  * Revision 1.12 2014/10/13 08:53:51 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.11 2014/10/06 15:13:25 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.10 2009/10/23 12:56:29 j_novak
44  * New base T_LEG_MI
45  *
46  * Revision 1.9 2009/10/13 13:49:58 j_novak
47  * New base T_LEG_MP.
48  *
49  * Revision 1.8 2005/11/22 10:13:44 p_grandclement
50  * Correction of the normalization in the case nt=1
51  *
52  * Revision 1.7 2005/02/16 15:33:14 m_forot
53  * Correct the case T_LEG
54  *
55  * Revision 1.6 2004/11/23 15:17:20 m_forot
56  * Added the bases for the cases without any equatorial symmetry
57  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
58  *
59  * Revision 1.5 2003/10/13 20:52:58 e_gourgoulhon
60  * Loop variables i and l have now local scope.
61  *
62  * Revision 1.4 2003/09/17 12:30:22 j_novak
63  * New checks for changing to T_LEG* bases.
64  *
65  * Revision 1.3 2003/09/16 08:54:09 j_novak
66  * Addition of the T_LEG_II base (odd in theta, only for odd m) and the
67  * transformation functions to and from the T_SIN_P base.
68  *
69  * Revision 1.2 2002/10/16 14:37:16 j_novak
70  * Reorganization of #include instructions of standard C++, in order to
71  * use experimental version 3 of gcc.
72  *
73  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
74  * LORENE
75  *
76  * Revision 2.11 2000/09/29 16:10:37 eric
77  * Ajout des bases T_LEG_IP et T_LEG_PI.
78  *
79  * Revision 2.10 2000/03/31 15:58:26 phil
80  * changement des bases meme si etat est zero
81  *
82  * Revision 2.9 1999/12/22 16:25:52 eric
83  * Traitement du cas ETATZERO
84  * Test sur c_cf avant d'appeler coef().
85  *
86  * Revision 2.8 1999/12/16 16:41:43 phil
87  * *** empty log message ***
88  *
89  * Revision 2.7 1999/12/16 16:10:07 phil
90  * correction cas nt = 1.
91  *
92  * Revision 2.6 1999/11/30 12:46:40 eric
93  * Valeur::base est desormais du type Base_val et non plus Base_val*.
94  *
95  * Revision 2.5 1999/11/22 11:36:14 eric
96  * Suppression de #include "tenseur.h"
97  *
98  * Revision 2.4 1999/10/18 14:12:39 eric
99  * Les bases sont desormais membres des Mtbl_cf.
100  *
101  * Revision 2.3 1999/04/14 10:19:52 phil
102  * *** empty log message ***
103  *
104  * Revision 2.2 1999/04/14 09:53:21 phil
105  * remplacement de malloc en new
106  *
107  * Revision 2.1 1999/04/14 09:38:47 phil
108  * Changement liberation memoire : free -> delete
109  *
110  * Revision 2.0 1999/04/13 16:47:28 phil
111  * *** empty log message ***
112  *
113  *
114  * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_ylm_i.C,v 1.12 2014/10/13 08:53:51 j_novak Exp $
115  *
116  */
117 
118 // headers C
119 #include <cassert>
120 #include <cstdlib>
121 #include <cmath>
122 
123 // headers Lorene
124 #include "type_parite.h"
125 #include "valeur.h"
126 #include "proto.h"
127 
128 namespace Lorene {
129 void ylm_i_pasprevu(const int*, const double*, double*) ;
130 
132 
133  static void (*chbase_t[MAX_BASE])(const int*, const double*,
134  double*) ;
135  static int nouv_base_t[MAX_BASE] ;
136  static int premier_appel = 1 ;
137 
138  int deg[3] ;
139 
140  if (premier_appel==1) {
141  premier_appel = 0 ;
142 
143  for (int i=0; i<MAX_BASE; i++) {
144  chbase_t[i] = ylm_i_pasprevu ;
145  nouv_base_t[i] = NONDEF ;
146  }
147 
148  chbase_t[T_LEG_P >> TRA_T] = chb_legp_cossincp ;
149  nouv_base_t[T_LEG_P >> TRA_T] = T_COSSIN_CP ;
150 
151  chbase_t[T_LEG_I >> TRA_T] = chb_legi_cossinci ;
152  nouv_base_t[T_LEG_I >> TRA_T] = T_COSSIN_CI ;
153 
154  chbase_t[T_LEG >> TRA_T] = chb_leg_cossinc ;
155  nouv_base_t[T_LEG >> TRA_T] = T_COSSIN_C ;
156 
157  chbase_t[T_LEG_MP >> TRA_T] = chb_legmp_cos ;
158  nouv_base_t[T_LEG_MP >> TRA_T] = T_COS ;
159 
160  chbase_t[T_LEG_MI >> TRA_T] = chb_legmi_sin ;
161  nouv_base_t[T_LEG_MI >> TRA_T] = T_SIN ;
162 
163  chbase_t[T_LEG_PP >> TRA_T] = chb_legpp_cosp ;
164  nouv_base_t[T_LEG_PP >> TRA_T] = T_COS_P ;
165 
166  chbase_t[T_LEG_IP >> TRA_T] = chb_legip_cosi ;
167  nouv_base_t[T_LEG_IP >> TRA_T] = T_COS_I ;
168 
169  chbase_t[T_LEG_PI >> TRA_T] = chb_legpi_sini ;
170  nouv_base_t[T_LEG_PI >> TRA_T] = T_SIN_I ;
171 
172  chbase_t[T_LEG_II >> TRA_T] = chb_legii_sinp ;
173  nouv_base_t[T_LEG_II >> TRA_T] = T_SIN_P ;
174  }
175 
176 //---------------------------------------------------------------------------
177 // fin des operation de premier appel
178 //---------------------------------------------------------------------------
179 
180  // Tout null ?
181  int nzone = get_mg()->get_nzone() ;
182  if (etat == ETATZERO) {
183  for (int l=0 ; l<nzone ; l++) {
184  int vbase_r = base.b[l] & MSQ_R ;
185  int vbase_t = base.b[l] & MSQ_T ;
186  int vbase_p = base.b[l] & MSQ_P ;
187 
188  int vbase_t_tra = vbase_t >> TRA_T ;
189  base.b[l] = ( vbase_p | nouv_base_t[vbase_t_tra] ) | vbase_r ;
190  }
191  return ;
192  }
193 
194  // Protection
195  assert(etat != ETATNONDEF) ;
196 
197  if (c_cf == 0x0) {
198  coef() ; // The coefficients are required
199  }
200 
201 // Boucle sur les differentes zones
202 
203 
204  for (int l=0; l<nzone; l++) {
205 
206 // On recupere les anciennes bases en r, phi et theta :
207  int vbase_r = base.b[l] & MSQ_R ;
208  int vbase_t = base.b[l] & MSQ_T ;
209  int vbase_p = base.b[l] & MSQ_P ;
210 
211  if ((vbase_t != T_COSSIN_CP) && (vbase_t != T_COS_I) &&
212  (vbase_t != T_COS_P) && (vbase_t != T_COSSIN_CI) &&
213  (vbase_t != T_SIN_P) && (vbase_t != T_SIN_I) &&
214  (vbase_t != T_COSSIN_C) && (vbase_t != T_COSSIN_S) &&
215  (vbase_t != T_COS) && (vbase_t != T_SIN) )
216  { // cas ou le calcul est necessaire
217 
218  int vbase_t_tra = vbase_t >> TRA_T ;
219  assert(vbase_t_tra < MAX_BASE) ;
220  bool pair = ( (vbase_t == T_LEG_PP) || (vbase_t == T_LEG_IP) ||
221  (vbase_t == T_LEG_MP) ) ;
222  bool impair = ( (vbase_t == T_LEG_PI) || (vbase_t == T_LEG_II) ||
223  (vbase_t == T_LEG_MI) ) ;
224 
225 // Nouvelle base :
226  base.b[l] = ( vbase_p | nouv_base_t[vbase_t_tra] ) | vbase_r ;
227  if (get_mg()->get_nt(l)==1) {
228  *c_cf->t[l] /= sqrt(2.) ;
229  continue ;
230  }
231 //... tbl contenant les coefficients dans la zone l :
232  Tbl* cf = c_cf->t[l] ;
233 
234  if (cf->get_etat() == ETATZERO) continue ; // On ne fait rien si le tbl = 0
235 
236 //... resultat du calcul :
237  double* resu = new double [cf->get_taille()] ;
238 
239  deg[0] = get_mg()->get_np(l) ; // nb. de degres de liberte en phi
240  deg[1] = get_mg()->get_nt(l) ; // nb. de degres de liberte en theta
241  deg[2] = get_mg()->get_nr(l) ; // nb. de degres de liberte en r
242 
243 
244 // Transformation en theta:
245 //-------------------------
246 
247  if ((pair && (vbase_p == P_COSSIN_I)) ||
248  (impair && (vbase_p == P_COSSIN_P)) )
249  ylm_i_pasprevu(deg, (cf->t), resu ) ;
250  else
251  chbase_t[vbase_t_tra](deg, (cf->t), resu ) ;
252 
253  // On branche le tbl contenant les coef sur resu :
254  delete [] cf->t ; // les anciens coef. sont oublies
255  cf->t = resu ; // nouveaux coef.
256 
257  } // fin du cas ou la transformation devait etre effectuee
258 
259  } // fin de la boucle sur les differentes zones
260 
261  // On met les bonnes bases dans c_cf :
262  c_cf->base = base ;
263 
264 }
265 
266 //******************************************************************************
267 
268 void ylm_i_pasprevu(const int*, const double*, double*) {
269 
270  cout <<
271  "Valeur::ylm_i: change of basis not implemented yet !"
272  << endl ;
273  abort() ;
274 }
275 
276 
277 }
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:331
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
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition: mtbl_cf.h:205
Basic array class.
Definition: tbl.h:161
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
double * t
The array of double.
Definition: tbl.h:173
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:729
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
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: valeur.h:295
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
#define TRA_T
Translation en Theta, used for a bitwise shift (in hex)
Definition: type_parite.h:160
#define T_LEG_MP
fct. de Legendre associees avec m pair
Definition: type_parite.h:238
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition: type_parite.h:247
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define T_LEG_PI
fct. de Legendre associees paires avec m impair
Definition: type_parite.h:224
#define T_LEG
fct. de Legendre associees
Definition: type_parite.h:236
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition: type_parite.h:202
#define MSQ_R
Extraction de l'info sur R.
Definition: type_parite.h:152
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
#define T_LEG_P
fct. de Legendre associees paires
Definition: type_parite.h:216
#define T_LEG_IP
fct. de Legendre associees impaires avec m pair
Definition: type_parite.h:222
#define NONDEF
base inconnue
Definition: type_parite.h:150
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
#define T_LEG_MI
fct. de Legendre associees avec m impair
Definition: type_parite.h:240
#define MSQ_T
Extraction de l'info sur Theta.
Definition: type_parite.h:154
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
#define T_LEG_II
fct. de Legendre associees impaires avec m impair
Definition: type_parite.h:226
#define T_LEG_I
fct. de Legendre associees impaires
Definition: type_parite.h:220
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
#define MSQ_P
Extraction de l'info sur Phi.
Definition: type_parite.h:156
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
#define T_LEG_PP
fct. de Legendre associees paires avec m pair
Definition: type_parite.h:218
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
Lorene prototypes.
Definition: app_hor.h:64