LORENE
base_val_theta_funct.C
1 /*
2  * Method of the class Base_val to get the values of the theta basis functions
3  * at the theta collocation points.
4  *
5  * (see file base_val.h for the documentation)
6  */
7 
8 /*
9  * Copyright (c) 1999-2001 Eric Gourgoulhon
10  * Copyright (c) 2001 Jerome Novak
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 char base_val_theta_funct_C[] = "$Header: /cvsroot/Lorene/C++/Source/Base_val/base_val_theta_funct.C,v 1.9 2014/10/13 08:52:39 j_novak Exp $" ;
32 
33 /*
34  * $Id: base_val_theta_funct.C,v 1.9 2014/10/13 08:52:39 j_novak Exp $
35  * $Log: base_val_theta_funct.C,v $
36  * Revision 1.9 2014/10/13 08:52:39 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.8 2014/10/06 15:12:57 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.7 2009/10/08 16:20:13 j_novak
43  * Addition of new bases T_COS and T_SIN.
44  *
45  * Revision 1.6 2007/10/23 17:15:12 j_novak
46  * Added the bases T_COSSIN_C and T_COSSIN_S
47  *
48  * Revision 1.5 2006/05/30 08:30:12 n_vasset
49  * Implementation of sine-like bases (T_SIN_P, T_SIN_I, T_COSSIN_SI, etc...).
50  *
51  * Revision 1.4 2005/05/27 14:54:59 j_novak
52  * Added new bases T_COSSIN_CI and T_COS_I
53  *
54  * Revision 1.3 2004/11/23 15:08:01 m_forot
55  * Added the bases for the cases without any equatorial symmetry
56  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
57  *
58  * Revision 1.2 2002/10/16 14:36:30 j_novak
59  * Reorganization of #include instructions of standard C++, in order to
60  * use experimental version 3 of gcc.
61  *
62  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
63  * LORENE
64  *
65  * Revision 1.3 2001/11/16 09:28:35 novak
66  * The case nt=1 is treated separately
67  *
68  * Revision 1.2 1999/12/29 10:49:47 eric
69  * Mehtode const.
70  *
71  * Revision 1.1 1999/12/28 12:58:35 eric
72  * Initial revision
73  *
74  *
75  * $Header: /cvsroot/Lorene/C++/Source/Base_val/base_val_theta_funct.C,v 1.9 2014/10/13 08:52:39 j_novak Exp $
76  *
77  */
78 
79 // Headers C
80 #include <cstdlib>
81 #include <cmath>
82 
83 
84 // Headers Lorene
85 #include "base_val.h"
86 #include "type_parite.h"
87 #include "tbl.h"
88 
89 // Local prototypes
90 namespace Lorene {
91 void theta_funct_pas_prevu(int, double*) ;
92 void theta_funct_cos(int, double*) ;
93 void theta_funct_sin(int, double*) ;
94 void theta_funct_cos_p(int, double*) ;
95 void theta_funct_cos_i(int, double*) ;
96 void theta_funct_sin_p(int, double*) ;
97 void theta_funct_sin_i(int, double*) ;
98 void theta_funct_cossin_cp(int, double*) ;
99 void theta_funct_cossin_ci(int, double*) ;
100 void theta_funct_cossin_sp(int, double*) ;
101 void theta_funct_cossin_si(int, double*) ;
102 void theta_funct_cossin_c(int, double*) ;
103 void theta_funct_cossin_s(int, double*) ;
104 
105 //************************************************************************
106 // user interface : method Base_val::theta_functions
107 //************************************************************************
108 
109 const Tbl& Base_val::theta_functions(int l, int nt) const {
110 
111  const int nmax = 20 ; // maximum number of couples (base_t, nt)
112  static int nb_done = 0 ; // number of Tbl already computed
113  static int base_t_done[nmax] ; // theta bases already treated
114  static int nt_done[nmax] ; // number of points already treated
115  static Tbl* tab[nmax] ; // result for couples (base_t, nt)
116 
117  static void(*vbasecol[MAX_BASE])(int, double*) ; // computation routines
118  static int dim2[MAX_BASE] ; // dim(2) of the Tbl's
119 
120  static int premier_appel = 1 ;
121 
122  // Initializations at first call
123  // -----------------------------
124  if (premier_appel == 1) {
125 
126  premier_appel = 0 ;
127 
128  for (int i=0 ; i<MAX_BASE ; i++) {
129  vbasecol[i] = theta_funct_pas_prevu ;
130  dim2[i] = 0 ;
131  }
132 
133  vbasecol[T_COS >> TRA_T] = theta_funct_cos ;
134  vbasecol[T_SIN >> TRA_T] = theta_funct_sin ;
135  vbasecol[T_COS_P >> TRA_T] = theta_funct_cos_p ;
136  vbasecol[T_COS_I >> TRA_T] = theta_funct_cos_i ;
137  vbasecol[T_SIN_I >> TRA_T] = theta_funct_sin_i ;
138  vbasecol[T_SIN_P >> TRA_T] = theta_funct_sin_p ;
139  vbasecol[T_COSSIN_CP >> TRA_T] = theta_funct_cossin_cp ;
140  vbasecol[T_COSSIN_CI >> TRA_T] = theta_funct_cossin_ci ;
141  vbasecol[T_COSSIN_SP >> TRA_T] = theta_funct_cossin_sp ;
142  vbasecol[T_COSSIN_SI >> TRA_T] = theta_funct_cossin_si ;
143  vbasecol[T_COSSIN_C >> TRA_T] = theta_funct_cossin_c ;
144  vbasecol[T_COSSIN_S >> TRA_T] = theta_funct_cossin_s ;
145 
146  dim2[T_COS >> TRA_T] = 1 ;
147  dim2[T_SIN >> TRA_T] = 1 ;
148  dim2[T_COS_P >> TRA_T] = 1 ;
149  dim2[T_COS_I >> TRA_T] = 1 ;
150  dim2[T_SIN_P >> TRA_T] = 1 ;
151  dim2[T_SIN_I >> TRA_T] = 1 ;
152  dim2[T_COSSIN_CP >> TRA_T] = 2 ;
153  dim2[T_COSSIN_CI >> TRA_T] = 2 ;
154  dim2[T_COSSIN_SP >> TRA_T] = 2 ;
155  dim2[T_COSSIN_SI >> TRA_T] = 2 ;
156  dim2[T_COSSIN_C >> TRA_T] = 2 ;
157  dim2[T_COSSIN_S >> TRA_T] = 2 ;
158 
159  }
160 
161  // Computation
162  // -----------
163 
164  int base_t = ( b[l] & MSQ_T ) >> TRA_T ;
165 
166  // Has this couple (base_t, nt) been previously considered ?
167  // ---------------------------------------------------------
168  int index = -1 ;
169  for (int i=0; i<nb_done; i++) {
170  if ( (base_t_done[i] == base_t) && (nt_done[i] == nt) ) {
171  index = i ;
172  }
173  }
174 
175  // If not, a new computation must be performed
176  // -------------------------------------------
177  if (index == -1) {
178  if ( nb_done >= nmax ) {
179  cout << "Base_val::theta_functions : nb_done >= nmax ! " << endl ;
180  abort() ;
181  }
182 
183  index = nb_done ;
184 
185  tab[index] = new Tbl( dim2[base_t], nt, nt ) ;
186  (tab[index])->set_etat_qcq() ;
187 
188  vbasecol[base_t](nt, (tab[index])->t ) ;
189 
190  base_t_done[index] = base_t ;
191  nt_done[index] = nt ;
192  nb_done++ ;
193 
194  } // end of the case where the computation had to be done
195 
196 
197  return *(tab[index]) ;
198 
199 }
200 
201 
202 //************************************************************************
203 // computational subroutines
204 //************************************************************************
205 
206 //====================================
207 // Unknown case
208 //====================================
209 
210 void theta_funct_pas_prevu(int, double*) {
211 
212  cout << "Base_val::theta_functions : theta basis not implemented !"
213  << endl ;
214  abort() ;
215 
216 }
217 
218 //==============================================
219 // Basis cos(j*theta) T_COS
220 //==============================================
221 
222 void theta_funct_cos(int nt, double* ff) {
223 
224  double xx = ( nt > 1 ? M_PI / double(nt-1) : 0.) ;
225 
226  for (int i = 0; i < nt ; i++ ) {
227  for (int j = 0; j < nt ; j++ ) {
228  double theta = xx*j ;
229  ff[nt*i+ j] = cos(i * theta);
230  }
231  }
232 
233 }
234 
235 //==============================================
236 // Basis sin(j* theta) T_SIN
237 //==============================================
238 
239 void theta_funct_sin(int nt, double* ff) {
240 
241  double xx = ( nt > 1 ? M_PI / double(nt-1) : 0.) ;
242 
243  for (int i = 0; i < nt ; i++ ) {
244  for (int j = 0; j < nt ; j++ ) {
245  double theta = xx*j ;
246  ff[nt*i+ j] = sin(i * theta);
247  }
248  }
249 
250 }
251 
252 //==============================================
253 // Basis cos(2*j* theta) T_COS_P
254 //==============================================
255 
256 void theta_funct_cos_p(int nt, double* ff) {
257 
258  double xx = ( nt > 1 ? M_PI / double(2*(nt-1)) : 0.) ;
259 
260  for (int i = 0; i < nt ; i++ ) {
261  for (int j = 0; j < nt ; j++ ) {
262  double theta = xx*j ;
263  ff[nt*i+ j] = cos(2*i * theta);
264  }
265  }
266 
267 }
268 
269 //==============================================
270 // Basis cos((2*j+1)* theta) T_COS_I
271 //==============================================
272 
273 void theta_funct_cos_i(int nt, double* ff) {
274 
275  double xx = ( nt > 1 ? M_PI / double(2*(nt-1)) : 0.) ;
276 
277  for (int i = 0; i < nt ; i++ ) {
278  for (int j = 0; j < nt ; j++ ) {
279  double theta = xx*j ;
280  ff[nt*i+ j] = cos((2*i+1) * theta);
281  }
282  }
283 
284 }
285 
286 //==============================================
287 // Basis sin(2*j* theta) T_SIN_P
288 //==============================================
289 
290 void theta_funct_sin_p(int nt, double* ff) {
291 
292  double xx = ( nt > 1 ? M_PI / double(2*(nt-1)) : 0.) ;
293 
294  for (int i = 0; i < nt ; i++ ) {
295  for (int j = 0; j < nt ; j++ ) {
296  double theta = xx*j ;
297  ff[nt*i+ j] = sin(2*i * theta);
298  }
299  }
300 
301 }
302 
303 //==============================================
304 // Basis sin((2*j+1)* theta) T_SIN_I
305 //==============================================
306 
307 void theta_funct_sin_i(int nt, double* ff) {
308 
309  double xx = ( nt > 1 ? M_PI / double(2*(nt-1)) : 0.) ;
310 
311  for (int i = 0; i < nt ; i++ ) {
312  for (int j = 0; j < nt ; j++ ) {
313  double theta = xx*j ;
314  ff[nt*i+ j] = sin((2*i+1) * theta);
315  }
316  }
317 
318 }
319 
320 //===========================================================
321 // Basis cos(2*j* theta)/sin((2*j+1) theta) T_COSSIN_CP
322 //===========================================================
323 
324 void theta_funct_cossin_cp(int nt, double* ff) {
325 
326  int nt2 = nt*nt ;
327 
328  double xx = ( nt > 1 ? M_PI / double(2*(nt-1)) : 0.) ;
329 
330  // cos(2i theta_j)
331  // ---------------
332  for (int i = 0; i < nt ; i++ ) {
333  for (int j = 0; j < nt ; j++ ) {
334  double theta = xx*j ;
335  ff[nt*i+ j] = cos(2*i * theta);
336  }
337  }
338 
339  // sin((2i+1) theta_j)
340  // -------------------
341 
342  for (int i = 0; i < nt-1 ; i++ ) {
343  for (int j = 0; j < nt ; j++ ) {
344  double theta = xx*j ;
345  ff[nt2+nt*i+ j] = sin((2*i+1) * theta);
346  }
347  }
348 
349  for (int j = 0; j < nt ; j++ ) {
350  ff[nt2+nt*(nt-1) + j] = 0 ;
351  }
352 
353 
354 }
355 
356 //===========================================================
357 // Basis cos((2*j+1)* theta)/sin(2*j*theta) T_COSSIN_CI
358 //===========================================================
359 
360 void theta_funct_cossin_ci(int nt, double* ff) {
361 
362  int nt2 = nt*nt ;
363 
364  double xx = ( nt > 1 ? M_PI / double(2*(nt-1)) : 0.) ;
365 
366  // cos((2i+1) theta_j)
367  // ---------------
368  for (int i = 0; i < nt ; i++ ) {
369  for (int j = 0; j < nt ; j++ ) {
370  double theta = xx*j ;
371  ff[nt*i+ j] = cos((2*i+1) * theta);
372  }
373  }
374 
375  // sin(2i theta_j)
376  // -------------------
377 
378  for (int i = 0; i < nt-1 ; i++ ) {
379  for (int j = 0; j < nt ; j++ ) {
380  double theta = xx*j ;
381  ff[nt2+nt*i+ j] = sin(2*i * theta);
382  }
383  }
384 
385  for (int j = 0; j < nt ; j++ ) {
386  ff[nt2+nt*(nt-1) + j] = 0 ;
387  }
388 
389 
390 }
391 
392 
393 //===========================================================
394 // Basis sin(2*j* theta)/cos((2*j+1) theta) T_COSSIN_SP
395 //===========================================================
396 
397 void theta_funct_cossin_sp(int nt, double* ff) {
398 
399  int nt2 = nt*nt ;
400 
401  double xx = ( nt > 1 ? M_PI / double(2*(nt-1)) : 0.) ;
402 
403  // sin(2i theta_j)
404  // ---------------
405  for (int i = 0; i < nt ; i++ ) {
406  for (int j = 0; j < nt ; j++ ) {
407  double theta = xx*j ;
408  ff[nt*i+ j] = sin(2*i * theta);
409  }
410  }
411 
412  // cos((2i+1) theta_j)
413  // -------------------
414 
415  for (int i = 0; i < nt-1 ; i++ ) {
416  for (int j = 0; j < nt ; j++ ) {
417  double theta = xx*j ;
418  ff[nt2+nt*i+ j] = cos((2*i+1) * theta);
419  }
420  }
421 
422  for (int j = 0; j < nt ; j++ ) {
423  ff[nt2+nt*(nt-1) + j] = 0 ;
424  }
425 
426 
427 }
428 
429 //===========================================================
430 // Basis sin((2*j+1)* theta)/cos(2*j*theta) T_COSSIN_SI
431 //===========================================================
432 
433 void theta_funct_cossin_si(int nt, double* ff) {
434 
435  int nt2 = nt*nt ;
436 
437  double xx = ( nt > 1 ? M_PI / double(2*(nt-1)) : 0.) ;
438 
439  // sin((2i+1) theta_j)
440  // ---------------
441  for (int i = 0; i < nt ; i++ ) {
442  for (int j = 0; j < nt ; j++ ) {
443  double theta = xx*j ;
444  ff[nt*i+ j] = sin((2*i+1) * theta);
445  }
446  }
447 
448  // cos(2i theta_j)
449  // -------------------
450 
451  for (int i = 0; i < nt-1 ; i++ ) {
452  for (int j = 0; j < nt ; j++ ) {
453  double theta = xx*j ;
454  ff[nt2+nt*i+ j] = cos(2*i * theta);
455  }
456  }
457 
458  for (int j = 0; j < nt ; j++ ) {
459  ff[nt2+nt*(nt-1) + j] = 0 ;
460  }
461 
462 
463 }
464 
465 //===========================================================
466 // Basis cos(j* theta)/sin(j*theta) T_COSSIN_C
467 //===========================================================
468 
469 void theta_funct_cossin_c(int nt, double* ff) {
470 
471  int nt2 = nt*nt ;
472 
473  double xx = ( nt > 1 ? M_PI / double(nt-1) : 0.) ;
474 
475  // cos(i theta_j)
476  // ---------------
477  for (int i = 0; i < nt ; i++ ) {
478  for (int j = 0; j < nt ; j++ ) {
479  double theta = xx*j ;
480  ff[nt*i+ j] = cos(i * theta);
481  }
482  }
483 
484  // sin(i theta_j)
485  // -------------------
486 
487  for (int i = 0; i < nt-1 ; i++ ) {
488  for (int j = 0; j < nt ; j++ ) {
489  double theta = xx*j ;
490  ff[nt2+nt*i+ j] = sin(i * theta);
491  }
492  }
493 
494  for (int j = 0; j < nt ; j++ ) {
495  ff[nt2+nt*(nt-1) + j] = 0 ;
496  }
497 
498 
499 }
500 
501 //===========================================================
502 // Basis sin(j* theta)/cos(j*theta) T_COSSIN_S
503 //===========================================================
504 
505 void theta_funct_cossin_s(int nt, double* ff) {
506 
507  int nt2 = nt*nt ;
508 
509  double xx = ( nt > 1 ? M_PI / double(nt-1) : 0.) ;
510 
511  // sin(i theta_j)
512  // ---------------
513  for (int i = 0; i < nt-1 ; i++ ) {
514  for (int j = 0; j < nt ; j++ ) {
515  double theta = xx*j ;
516  ff[nt*i+ j] = sin(i * theta);
517  }
518  }
519 
520  for (int j = 0; j < nt ; j++ ) {
521  ff[nt*(nt-1) + j] = 0 ;
522  }
523 
524  // cos(i theta_j)
525  // -------------------
526 
527  for (int i = 0; i < nt ; i++ ) {
528  for (int j = 0; j < nt ; j++ ) {
529  double theta = xx*j ;
530  ff[nt2+nt*i+ j] = cos(i * theta);
531  }
532  }
533 
534 }
535 
536 }
Lorene::Base_val::b
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:331
T_COS
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
T_COSSIN_SP
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
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
T_COS_I
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
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
Lorene::Base_val::theta_functions
const Tbl & theta_functions(int l, int nt) const
Values of the theta basis functions at the theta collocation points.
Definition: base_val_theta_funct.C:109
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
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
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
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
T_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206