LORENE
valeur_ylm.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_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_ylm.C,v 1.12 2014/10/13 08:53:51 j_novak Exp $" ;
25 
26 
27 /*
28  * Fonction membre de la classe Valeur qui calcule les coefficients
29  * de la decomposition spectrale en harmoniques spheriques
30  * a partir des coefficients en cos(l*theta) / sin(l*theta)
31  *
32  */
33 
34 /*
35  * $Id: valeur_ylm.C,v 1.12 2014/10/13 08:53:51 j_novak Exp $
36  * $Log: valeur_ylm.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:38 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.12 2000/09/29 16:10:21 eric
77  * Ajout des bases T_LEG_IP et T_LEG_PI.
78  *
79  * Revision 2.11 2000/03/31 15:58:40 phil
80  * changement des bases meme si etat est zero
81  *
82  * Revision 2.10 1999/12/22 16:25:26 eric
83  * Traitement du cas ETATZERO
84  * Test sur c_cf avant d'appeler coef().
85  *
86  * Revision 2.9 1999/12/16 16:41:38 phil
87  * *** empty log message ***
88  *
89  * Revision 2.8 1999/12/16 16:09:56 phil
90  * correction cas nt=1
91  *
92  * Revision 2.7 1999/11/30 12:46:33 eric
93  * Valeur::base est desormais du type Base_val et non plus Base_val*.
94  *
95  * Revision 2.6 1999/11/19 09:25:57 eric
96  * *** empty log message ***
97  *
98  * Revision 2.5 1999/10/18 14:12:09 eric
99  * Les bases sont desormais membres des Mtbl_cf.
100  *
101  * Revision 2.4 1999/04/14 10:19:46 phil
102  * *** empty log message ***
103  *
104  * Revision 2.3 1999/04/14 09:52:28 phil
105  * remplacement de malloc en new
106  *
107  * Revision 2.2 1999/04/14 09:36:03 phil
108  * Changement liberation memoire : free -> delete
109  *
110  * Revision 2.1 1999/04/13 16:44:07 phil
111  * *** empty log message ***
112  *
113  * Revision 2.0 1999/04/13 16:36:57 phil
114  * *** empty log message ***
115  *
116  * Revision 1.1 1999/04/13 16:36:14 phil
117  * Initial revision
118  *
119  *
120  * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_ylm.C,v 1.12 2014/10/13 08:53:51 j_novak Exp $
121  *
122  */
123 
124 // headers du C
125 #include <cassert>
126 #include <cstdlib>
127 #include <cmath>
128 
129 // headers Lorene
130 #include "type_parite.h"
131 #include "valeur.h"
132 #include "proto.h"
133 
134 namespace Lorene {
135 void ylm_pasprevu(const int*, const double*, double*) ;
136 
137 //******************************************************************************
138 void Valeur::ylm() {
139 
140  static void (*chbase_t[MAX_BASE])(const int*, const double*,
141  double*) ;
142  static int nouv_base_t[MAX_BASE] ;
143  static int premier_appel = 1 ;
144 
145  int deg[3] ;
146 
147  if (premier_appel==1) {
148  premier_appel = 0 ;
149 
150  for (int i=0; i<MAX_BASE; i++) {
151  chbase_t[i] = ylm_pasprevu ;
152  nouv_base_t[i] = NONDEF ;
153  }
154 
155  chbase_t[T_COSSIN_CP >> TRA_T] = chb_cossincp_legp ;
156  nouv_base_t[T_COSSIN_CP >> TRA_T] = T_LEG_P ;
157 
158  chbase_t[T_COSSIN_CI >> TRA_T] = chb_cossinci_legi ;
159  nouv_base_t[T_COSSIN_CI >> TRA_T] = T_LEG_I ;
160 
161  chbase_t[T_COSSIN_C >> TRA_T] = chb_cossinc_leg ;
162  nouv_base_t[T_COSSIN_C >> TRA_T] = T_LEG ;
163 
164  chbase_t[T_COS_P >> TRA_T] = chb_cosp_legpp ;
165  nouv_base_t[T_COS_P >> TRA_T] = T_LEG_PP ;
166 
167  chbase_t[T_COS_I >> TRA_T] = chb_cosi_legip ;
168  nouv_base_t[T_COS_I >> TRA_T] = T_LEG_IP ;
169 
170  chbase_t[T_SIN_I >> TRA_T] = chb_sini_legpi ;
171  nouv_base_t[T_SIN_I >> TRA_T] = T_LEG_PI ;
172 
173  chbase_t[T_SIN_P >> TRA_T] = chb_sinp_legii ;
174  nouv_base_t[T_SIN_P >> TRA_T] = T_LEG_II ;
175 
176  chbase_t[T_COS >> TRA_T] = chb_cos_legmp ;
177  nouv_base_t[T_COS >> TRA_T] = T_LEG_MP ;
178 
179  chbase_t[T_SIN >> TRA_T] = chb_sin_legmi ;
180  nouv_base_t[T_SIN >> TRA_T] = T_LEG_MI ;
181  }
182 
183 //---------------------------------------------------------------------------
184 // fin des operation de premier appel
185 //---------------------------------------------------------------------------
186 
187  // Tout null ?
188  int nzone = get_mg()->get_nzone() ;
189  if (etat == ETATZERO) {
190  for (int l=0 ; l<nzone ; l++) {
191  int vbase_r = base.b[l] & MSQ_R ;
192  int vbase_t = base.b[l] & MSQ_T ;
193  int vbase_p = base.b[l] & MSQ_P ;
194  int vbase_t_tra = vbase_t >> TRA_T ;
195  base.b[l] = ( vbase_p | nouv_base_t[vbase_t_tra] ) | vbase_r ;
196  }
197  return ;
198  }
199 
200  // Protection
201  assert(etat != ETATNONDEF) ;
202 
203  if (c_cf == 0x0) {
204  coef() ; // The coefficients are required
205  }
206 
207 // Boucle sur les differentes zones
208 
209  for (int l=0; l<nzone; l++) {
210 
211 // On recupere les anciennes bases en r, phi et theta :
212  int vbase_r = base.b[l] & MSQ_R ;
213  int vbase_t = base.b[l] & MSQ_T ;
214  int vbase_p = base.b[l] & MSQ_P ;
215 
216  if ((vbase_t != T_LEG_P) && (vbase_t != T_LEG_IP) &&
217  (vbase_t != T_LEG_PP) && (vbase_t != T_LEG_I) &&
218  (vbase_t != T_LEG_II) && (vbase_t != T_LEG_PI) &&
219  (vbase_t != T_LEG) && (vbase_t != T_LEG_MP) &&
220  (vbase_t != T_LEG_MI) )
221  { // cas ou le calcul est necessaire
222 
223  int vbase_t_tra = vbase_t >> TRA_T ;
224  assert(vbase_t_tra < MAX_BASE) ;
225  bool pair = ( (vbase_t == T_COS_P) || (vbase_t == T_COS_I) ||
226  (vbase_t == T_COS) ) ;
227  bool impair = ( (vbase_t == T_SIN_P) || (vbase_t == T_SIN_I) ||
228  (vbase_t == T_SIN) ) ;
229 
230 // Nouvelle base :
231  base.b[l] = ( vbase_p | nouv_base_t[vbase_t_tra] ) | vbase_r ;
232  if (get_mg()->get_nt(l)==1) {
233  *c_cf->t[l] *=sqrt(2.) ;
234  continue ;
235  }
236 
237 //... tbl contenant les coefficients dans la zone l :
238  Tbl* cf = c_cf->t[l] ;
239 
240  if (cf->get_etat() == ETATZERO) continue ; // On ne fait rien si le tbl = 0
241 
242  deg[0] = get_mg()->get_np(l) ; // nb. de degres de liberte en phi
243  deg[1] = get_mg()->get_nt(l) ; // nb. de degres de liberte en theta
244  deg[2] = get_mg()->get_nr(l) ; // nb. de degres de liberte en r
245 
246 //... resultat du calcul :
247  double* resu = new double [cf->get_taille()] ;
248 
249 // Transformation en theta:
250 //-------------------------
251 
252  if ((pair && (vbase_p == P_COSSIN_I)) ||
253  (impair && (vbase_p == P_COSSIN_P)) )
254  ylm_pasprevu(deg, (cf->t), resu ) ;
255  else
256  chbase_t[vbase_t_tra](deg, (cf->t), resu ) ;
257 
258 // On branche le tbl contenant les coef sur resu :
259  delete [] cf->t ; // les anciens coef. sont oublies
260  cf->t = resu ; // nouveaux coef.
261 
262  } // fin du cas ou la transformation devait etre effectuee
263 
264  } // fin de la boucle sur les differentes zones
265 
266  // On met les bonnes bases dans c_cf :
267  c_cf->base = base ;
268 
269 }
270 
271 //******************************************************************************
272 
273 void ylm_pasprevu(const int* , const double*, double* ) {
274 
275  cout <<
276  "Valeur::ylm: change of basis not implemented yet !"
277  << endl ;
278  abort() ;
279 }
280 
281 }
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Base_val::b
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:331
T_LEG_MI
#define T_LEG_MI
fct. de Legendre associees avec m impair
Definition: type_parite.h:240
T_LEG_MP
#define T_LEG_MP
fct. de Legendre associees avec m pair
Definition: type_parite.h:238
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition: type_parite.h:156
Lorene::Mtbl_cf::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition: mtbl_cf.h:205
T_LEG_PI
#define T_LEG_PI
fct. de Legendre associees paires avec m impair
Definition: type_parite.h:224
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
T_COS
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
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:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
TRA_T
#define TRA_T
Translation en Theta, used for a bitwise shift (in hex)
Definition: type_parite.h:160
T_SIN
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
Lorene::Valeur::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:729
T_COS_I
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
T_LEG_I
#define T_LEG_I
fct. de Legendre associees impaires
Definition: type_parite.h:220
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Lorene::Tbl::t
double * t
The array of double.
Definition: tbl.h:173
MSQ_T
#define MSQ_T
Extraction de l'info sur Theta.
Definition: type_parite.h:154
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition: type_parite.h:202
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
Lorene::Valeur::etat
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: valeur.h:295
MSQ_R
#define MSQ_R
Extraction de l'info sur R.
Definition: type_parite.h:152
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
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:452
T_LEG_P
#define T_LEG_P
fct. de Legendre associees paires
Definition: type_parite.h:216
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
T_LEG_IP
#define T_LEG_IP
fct. de Legendre associees impaires avec m pair
Definition: type_parite.h:222
Lorene::Mtbl_cf::base
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
P_COSSIN_I
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
NONDEF
#define NONDEF
base inconnue
Definition: type_parite.h:150
T_COSSIN_CI
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
T_COSSIN_CP
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
P_COSSIN_P
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition: type_parite.h:247
T_LEG_PP
#define T_LEG_PP
fct. de Legendre associees paires avec m pair
Definition: type_parite.h:218
T_LEG_II
#define T_LEG_II
fct. de Legendre associees impaires avec m impair
Definition: type_parite.h:226
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
T_LEG
#define T_LEG
fct. de Legendre associees
Definition: type_parite.h:236
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206