LORENE
mat_sin_legmi.C
1 /*
2  * Copyright (c) 2003-2009 Jerome Novak
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 char mat_sinp_legmi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/mat_sin_legmi.C,v 1.3 2014/10/13 08:53:14 j_novak Exp $" ;
24 
25 /*
26  * Fournit la matrice de passage pour la transformation des coefficients du
27  * developpement en sin(j*theta) dans les coefficients du developpement en
28  * fonctions associees de Legendre P_l^m(cos(theta)) avec m impair.
29  *
30  * Cette routine n'effectue le calcul de la matrice que si celui-ci n'a pas
31  * deja ete fait, sinon elle renvoie le pointeur sur une valeur precedemment
32  * calculee.
33  *
34  * Entree:
35  * -------
36  * int np : Nombre de degres de liberte en phi
37  * int nt : Nombre de degres de liberte en theta
38  *
39  * Sortie (valeur de retour) :
40  * ---------------------------
41  * double* mat_sin_legmi : pointeur sur le tableau contenant l'ensemble
42  * (pour les np/2 valeurs de m: m=1,3,...,np-1) des
43  * matrices de passage.
44  * La dimension du tableau est (np/2+1)*nt^2
45  * Le stokage est le suivant:
46  *
47  * mat_sin_legmi[ nt*nt* m + nt*l + j] = A_{mlj}
48  *
49  * ou A_{mlj} est defini par
50  *
51  * sin(j*theta) = som_{l=m}^{nt-2} A_{mlj} P_l^m( cos(theta) )
52  * pour 1 <= j <= nt-2
53  *
54  * ou P_n^m(x) represente la fonction de Legendre associee de degre n et
55  * d'ordre m normalisee de facon a ce que
56  *
57  * int_0^pi [ P_n^m(cos(theta)) ]^2 sin(theta) dtheta = 1
58  *
59  *
60  */
61 
62 /*
63  * $Id: mat_sin_legmi.C,v 1.3 2014/10/13 08:53:14 j_novak Exp $
64  * $Log: mat_sin_legmi.C,v $
65  * Revision 1.3 2014/10/13 08:53:14 j_novak
66  * Lorene classes and functions now belong to the namespace Lorene.
67  *
68  * Revision 1.2 2014/10/06 15:16:03 j_novak
69  * Modified #include directives to use c++ syntax.
70  *
71  * Revision 1.1 2009/10/23 12:54:47 j_novak
72  * New base T_LEG_MI
73  *
74  *
75  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/mat_sin_legmi.C,v 1.3 2014/10/13 08:53:14 j_novak Exp $
76  *
77  */
78 
79 // headers du C
80 #include <cstdlib>
81 #include <cmath>
82 #include <cassert>
83 
84 // Prototypage
85 #include "headcpp.h"
86 #include "proto.h"
87 
88 namespace Lorene {
89 //******************************************************************************
90 
91 double* mat_sin_legmi(int np, int nt) {
92 
93 #define NMAX 30 // Nombre maximun de couples(np,nt) differents
94  static double* tab[NMAX] ; // Tableau des pointeurs sur les tableaux
95  static int nb_dejafait = 0 ; // Nombre de tableaux deja initialises
96  static int np_dejafait[NMAX] ; // Valeurs de np pour lesquelles le
97  // calcul a deja ete fait
98  static int nt_dejafait[NMAX] ; // Valeurs de np pour lesquelles le
99  // calcul a deja ete fait
100 
101  int i, indice, j, j2, m, l ;
102 
103  // Les matrices A_{mlj} pour ce couple (np,nt) ont-elles deja ete calculees ?
104  indice = -1 ;
105  for ( i=0 ; i < nb_dejafait ; i++ ) {
106  if ( (np_dejafait[i] == np) && (nt_dejafait[i] == nt) ) indice = i ;
107  }
108 
109 
110 // Si le calcul n'a pas deja ete fait, il faut le faire :
111  if (indice == -1) {
112  if ( nb_dejafait >= NMAX ) {
113  cout << "mat_sinp_legii: nb_dejafait >= NMAX : "
114  << nb_dejafait << " <-> " << NMAX << endl ;
115  abort () ;
116  exit(-1) ;
117  }
118  indice = nb_dejafait ;
119  nb_dejafait++ ;
120  np_dejafait[indice] = np ;
121  nt_dejafait[indice] = nt ;
122 
123  tab[indice] = new double[(np/2+1)*nt*nt] ;
124  for (int qq=0; qq<nt*nt*(np/2+1); qq++)
125  tab[indice][qq] = -1.34 ;
126 
127 //-----------------------
128 // Preparation du calcul
129 //-----------------------
130 
131 // Sur-echantillonnage pour calculer les produits scalaires sans aliasing:
132  int nt2 = 2*nt - 1 ;
133  int nt2m1 = nt2 - 1 ;
134 
135  int deg[3] ;
136  deg[0] = 1 ;
137  deg[1] = 1 ;
138  deg[2] = nt2 ;
139 
140 // Tableaux de travail
141  double* yy = new double[nt2] ;
142  double* sint = new double[nt*nt2];
143 
144 // Calcul des sin( j theta) aux points de collocation
145 // de l'echantillonnage double :
146 
147  double dt = M_PI / double(nt2-1) ;
148  for (j=0; j<nt-1; j++) {
149  for (j2=0; j2<nt2; j2++) {
150  double theta = j2*dt ;
151  sint[nt2*j + j2] = sin( j * theta ) ;
152  }
153  }
154 
155 
156 //-------------------
157 // Boucle sur m
158 //-------------------
159 
160  int m_max = (np == 1) ? 1 : np-1 ;
161 
162  for (m=1; m <= m_max ; m+=2) {
163 
164  // Recherche des fonctions de Legendre associees d'ordre m :
165 
166  double* leg = legendre_norm(m, nt) ;
167 
168  for (l=m; l<nt-1; l++) { // boucle sur les P_l^m
169 
170  int parite = 1 - 2*((l-m)%2) ; // parite du P_l^m par rapport au plan theta=pi/2
171 
172  for (j=0; j<nt-1; j++) { // boucle sur les sin(j theta)
173 
174  //... produit scalaire de sin(j theta) par
175  // P_l^m(cos(theta))
176 
177  for (j2=0; j2<nt; j2++) {
178  yy[nt2m1-j2] = sint[nt2*j + j2] *
179  leg[nt2* (l-m) + 2*j2] ;
180 
181  }
182 
183  for (j2=nt; j2<nt2; j2++) {
184  yy[nt2m1-j2] = sint[nt2*j + j2] *
185  parite * leg[nt2* (l-m) + 2*nt2 -2 - 2*j2] ;
186 
187  }
188 
189 //....... on passe en Tchebyshev vis-a-vis de x=cos(theta) pour calculer
190 // l'integrale (routine int1d_cheb) :
191  cfrcheb(deg, deg, yy, deg, yy) ;
192 // for (int o=0; o<nt2; o++)
193 // cout << yy[o] << endl ;
194 
195  tab[indice][ nt*nt* ((m-1)/2) + nt*l + j] =
196  int1d_cheb(nt2, yy) ;
197 
198  } // fin de la boucle sur j (indice de sin(j theta) )
199 
200  } // fin de la boucle sur l (indice de P_l^m)
201 
202  delete [] leg ;
203 
204  } // fin de la boucle sur m
205 
206 // Liberation espace memoire
207 // -------------------------
208 
209  delete [] yy ;
210  delete [] sint ;
211 
212  } // fin du cas ou le calcul etait necessaire
213 
214  return tab[indice] ;
215 
216 }
217 
218 
219 }
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69
Lorene prototypes.
Definition: app_hor.h:64