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