LORENE
op_dsdtet.C
1 /*
2  * Copyright (c) 1999-2000 Jean-Alain Marck
3  * Copyright (c) 1999-2001 Eric Gourgoulhon
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 
25 
26 /*
27  * Ensemble des routines de base pour la derivation par rapport a theta
28  * (Utilisation interne)
29  *
30  * void _dsdtet_XXXX(Tbl * t, int & b)
31  * t pointeur sur le Tbl d'entree-sortie
32  * b base spectrale
33  *
34  */
35 
36 /*
37  * $Id: op_dsdtet.C,v 1.7 2016/12/05 16:18:07 j_novak Exp $
38  * $Log: op_dsdtet.C,v $
39  * Revision 1.7 2016/12/05 16:18:07 j_novak
40  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
41  *
42  * Revision 1.6 2014/10/13 08:53:25 j_novak
43  * Lorene classes and functions now belong to the namespace Lorene.
44  *
45  * Revision 1.5 2009/10/08 16:22:19 j_novak
46  * Addition of new bases T_COS and T_SIN.
47  *
48  * Revision 1.4 2006/03/10 12:45:38 j_novak
49  * Use of C++-style cast.
50  *
51  * Revision 1.3 2004/11/23 15:16:01 m_forot
52  *
53  * Added the bases for the cases without any equatorial symmetry
54  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
55  *
56  * Revision 1.2 2003/12/19 16:21:48 j_novak
57  * Shadow hunt
58  *
59  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
60  * LORENE
61  *
62  * Revision 2.6 2000/10/04 11:50:44 eric
63  * Ajout d' abort() dans le cas non prevu.
64  *
65  * Revision 2.5 2000/02/24 16:41:17 eric
66  * Initialisation a zero du tableau xo avant le calcul.
67  *
68  * Revision 2.4 1999/11/15 16:38:13 eric
69  * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
70  *
71  * Revision 2.3 1999/03/01 15:07:01 eric
72  * *** empty log message ***
73  *
74  * Revision 2.2 1999/02/23 11:04:52 hyc
75  * *** empty log message ***
76  *
77  * Revision 2.1 1999/02/22 17:12:06 hyc
78  * *** empty log message ***
79  *
80  * Revision 2.0 1999/02/22 15:51:02 hyc
81  * *** empty log message ***
82  *
83  *
84  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_dsdtet.C,v 1.7 2016/12/05 16:18:07 j_novak Exp $
85  *
86  */
87 
88 // Headers Lorene
89 #include "tbl.h"
90 
91 
92 // Routine pour les cas non prevus
93 //--------------------------------
94 namespace Lorene {
95 void _dsdtet_pas_prevu(Tbl* , int & b) {
96  cout << "Unknown theta basis in Mtbl_cf::dsdt() !" << endl ;
97  cout << " basis: " << hex << b << endl ;
98  abort() ;
99 }
100 
101 // cas T_COS
102 //------------
103 void _dsdtet_t_cos(Tbl* tb, int & b)
104 {
105 
106  // Peut-etre rien a faire ?
107  if (tb->get_etat() == ETATZERO) {
108  int base_r = b & MSQ_R ;
109  int base_p = b & MSQ_P ;
110  b = base_r | base_p | T_SIN ;
111  return ;
112  }
113 
114  // Protection
115  assert(tb->get_etat() == ETATQCQ) ;
116 
117  // Pour le confort
118  int nr = (tb->dim).dim[0] ; // Nombre
119  int nt = (tb->dim).dim[1] ; // de points
120  int np = (tb->dim).dim[2] ; // physiques REELS
121  np = np - 2 ; // Nombre de points physiques
122 
123  // Variables statiques
124  static double* cx = 0 ;
125  static int nt_pre =0 ;
126 
127  // Test sur np pour initialisation eventuelle
128  {
129  if (nt > nt_pre) {
130  nt_pre = nt ;
131  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
132  for (int i=0 ; i<nt ; i++) {
133  cx[i] = - double(i) ;
134  }
135  }
136  } // Fin de region critique
137 
138  // pt. sur le tableau de double resultat
139  double* xo = new double[(tb->dim).taille] ;
140 
141  // Initialisation a zero :
142  for (int i=0; i<(tb->dim).taille; i++) {
143  xo[i] = 0 ;
144  }
145 
146  // On y va...
147  double* xi = tb->t ;
148  double* xci = xi ; // Pointeurs
149  double* xco = xo ; // courants
150 
151  // k = 0
152  for (int j=0 ; j<nt ; j++) {
153  for (int i=0 ; i<nr ; i++ ) {
154  *xco = cx[j] * (*xci) ;
155  xci++ ;
156  xco++ ;
157  } // Fin de la boucle sur r
158  } // Fin de la boucle sur theta
159 
160  // k = 1
161  xci += nr*nt ;
162  xco += nr*nt ;
163 
164  // k >=2
165  for (int k=2 ; k<np+1 ; k++) {
166  for (int j=0 ; j<nt ; j++) {
167  for (int i=0 ; i<nr ; i++ ) {
168  *xco = cx[j] * (*xci) ;
169  xci++ ;
170  xco++ ;
171  } // Fin de la boucle sur r
172  } // Fin de la boucle sur theta
173  } // Fin de la boucle sur phi
174 
175  // On remet les choses la ou il faut
176  delete [] tb->t ;
177  tb->t = xo ;
178 
179  // base de developpement
180  int base_r = b & MSQ_R ;
181  int base_p = b & MSQ_P ;
182  b = base_r | base_p | T_SIN ;
183 }
184 
185 // cas T_SIN
186 //------------
187 void _dsdtet_t_sin(Tbl* tb, int & b)
188 {
189 
190  // Peut-etre rien a faire ?
191  if (tb->get_etat() == ETATZERO) {
192  int base_r = b & MSQ_R ;
193  int base_p = b & MSQ_P ;
194  b = base_r | base_p | T_COS ;
195  return ;
196  }
197 
198  // Protection
199  assert(tb->get_etat() == ETATQCQ) ;
200 
201  // Pour le confort
202  int nr = (tb->dim).dim[0] ; // Nombre
203  int nt = (tb->dim).dim[1] ; // de points
204  int np = (tb->dim).dim[2] ; // physiques REELS
205  np = np - 2 ; // Nombre de points physiques
206 
207  // Variables statiques
208  static double* cx = 0 ;
209  static int nt_pre = 0 ;
210 
211  // Test sur np pour initialisation eventuelle
212  {
213  if (nt > nt_pre) {
214  nt_pre = nt ;
215  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
216  for (int i=0 ; i<nt ; i++) {
217  cx[i] = double(i) ;
218  }
219  }
220  } // Fin de region critique
221 
222  // pt. sur le tableau de double resultat
223  double* xo = new double[(tb->dim).taille] ;
224 
225  // Initialisation a zero :
226  for (int i=0; i<(tb->dim).taille; i++) {
227  xo[i] = 0 ;
228  }
229 
230  // On y va...
231  double* xi = tb->t ;
232  double* xci = xi ; // Pointeurs
233  double* xco = xo ; // courants
234 
235  // k = 0
236  for (int j=0 ; j<nt ; j++) {
237  for (int i=0 ; i<nr ; i++ ) {
238  *xco = cx[j] * (*xci) ;
239  xci++ ;
240  xco++ ;
241  } // Fin de la boucle sur r
242  } // Fin de la boucle sur theta
243 
244  // k = 1
245  xci += nr*nt ;
246  xco += nr*nt ;
247 
248  // k >=2
249  for (int k=2 ; k<np+1 ; k++) {
250  for (int j=0 ; j<nt ; j++) {
251  for (int i=0 ; i<nr ; i++ ) {
252  *xco = cx[j] * (*xci) ;
253  xci++ ;
254  xco++ ;
255  } // Fin de la boucle sur r
256  } // Fin de la boucle sur theta
257  } // Fin de la boucle sur phi
258 
259  // On remet les choses la ou il faut
260  delete [] tb->t ;
261  tb->t = xo ;
262 
263  // base de developpement
264  int base_r = b & MSQ_R ;
265  int base_p = b & MSQ_P ;
266  b = base_r | base_p | T_COS ;
267 }
268 
269 // cas T_COS_P
270 //------------
271 void _dsdtet_t_cos_p(Tbl* tb, int & b)
272 {
273 
274  // Peut-etre rien a faire ?
275  if (tb->get_etat() == ETATZERO) {
276  int base_r = b & MSQ_R ;
277  int base_p = b & MSQ_P ;
278  b = base_r | base_p | T_SIN_P ;
279  return ;
280  }
281 
282  // Protection
283  assert(tb->get_etat() == ETATQCQ) ;
284 
285  // Pour le confort
286  int nr = (tb->dim).dim[0] ; // Nombre
287  int nt = (tb->dim).dim[1] ; // de points
288  int np = (tb->dim).dim[2] ; // physiques REELS
289  np = np - 2 ; // Nombre de points physiques
290 
291  // Variables statiques
292  static double* cx = 0 ;
293  static int nt_pre =0 ;
294 
295  // Test sur np pour initialisation eventuelle
296  {
297  if (nt > nt_pre) {
298  nt_pre = nt ;
299  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
300  for (int i=0 ; i<nt ; i++) {
301  cx[i] = - (2*i) ;
302  }
303  }
304  } // Fin de region critique
305 
306  // pt. sur le tableau de double resultat
307  double* xo = new double[(tb->dim).taille] ;
308 
309  // Initialisation a zero :
310  for (int i=0; i<(tb->dim).taille; i++) {
311  xo[i] = 0 ;
312  }
313 
314  // On y va...
315  double* xi = tb->t ;
316  double* xci = xi ; // Pointeurs
317  double* xco = xo ; // courants
318 
319  // k = 0
320  for (int j=0 ; j<nt ; j++) {
321  for (int i=0 ; i<nr ; i++ ) {
322  *xco = cx[j] * (*xci) ;
323  xci++ ;
324  xco++ ;
325  } // Fin de la boucle sur r
326  } // Fin de la boucle sur theta
327 
328  // k = 1
329  xci += nr*nt ;
330  xco += nr*nt ;
331 
332  // k >=2
333  for (int k=2 ; k<np+1 ; k++) {
334  for (int j=0 ; j<nt ; j++) {
335  for (int i=0 ; i<nr ; i++ ) {
336  *xco = cx[j] * (*xci) ;
337  xci++ ;
338  xco++ ;
339  } // Fin de la boucle sur r
340  } // Fin de la boucle sur theta
341  } // Fin de la boucle sur phi
342 
343  // On remet les choses la ou il faut
344  delete [] tb->t ;
345  tb->t = xo ;
346 
347  // base de developpement
348  int base_r = b & MSQ_R ;
349  int base_p = b & MSQ_P ;
350  b = base_r | base_p | T_SIN_P ;
351 }
352 
353 // cas T_SIN_P
354 //------------
355 void _dsdtet_t_sin_p(Tbl* tb, int & b)
356 {
357 
358  // Peut-etre rien a faire ?
359  if (tb->get_etat() == ETATZERO) {
360  int base_r = b & MSQ_R ;
361  int base_p = b & MSQ_P ;
362  b = base_r | base_p | T_COS_P ;
363  return ;
364  }
365 
366  // Protection
367  assert(tb->get_etat() == ETATQCQ) ;
368 
369  // Pour le confort
370  int nr = (tb->dim).dim[0] ; // Nombre
371  int nt = (tb->dim).dim[1] ; // de points
372  int np = (tb->dim).dim[2] ; // physiques REELS
373  np = np - 2 ; // Nombre de points physiques
374 
375  // Variables statiques
376  static double* cx = 0 ;
377  static int nt_pre = 0 ;
378 
379  // Test sur np pour initialisation eventuelle
380  {
381  if (nt > nt_pre) {
382  nt_pre = nt ;
383  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
384  for (int i=0 ; i<nt ; i++) {
385  cx[i] = (2*i) ;
386  }
387  }
388  } // Fin de region critique
389 
390  // pt. sur le tableau de double resultat
391  double* xo = new double[(tb->dim).taille] ;
392 
393  // Initialisation a zero :
394  for (int i=0; i<(tb->dim).taille; i++) {
395  xo[i] = 0 ;
396  }
397 
398  // On y va...
399  double* xi = tb->t ;
400  double* xci = xi ; // Pointeurs
401  double* xco = xo ; // courants
402 
403  // k = 0
404  for (int j=0 ; j<nt ; j++) {
405  for (int i=0 ; i<nr ; i++ ) {
406  *xco = cx[j] * (*xci) ;
407  xci++ ;
408  xco++ ;
409  } // Fin de la boucle sur r
410  } // Fin de la boucle sur theta
411 
412  // k = 1
413  xci += nr*nt ;
414  xco += nr*nt ;
415 
416  // k >=2
417  for (int k=2 ; k<np+1 ; k++) {
418  for (int j=0 ; j<nt ; j++) {
419  for (int i=0 ; i<nr ; i++ ) {
420  *xco = cx[j] * (*xci) ;
421  xci++ ;
422  xco++ ;
423  } // Fin de la boucle sur r
424  } // Fin de la boucle sur theta
425  } // Fin de la boucle sur phi
426 
427  // On remet les choses la ou il faut
428  delete [] tb->t ;
429  tb->t = xo ;
430 
431  // base de developpement
432  int base_r = b & MSQ_R ;
433  int base_p = b & MSQ_P ;
434  b = base_r | base_p | T_COS_P ;
435 }
436 
437 // cas T_SIN_I
438 //------------
439 void _dsdtet_t_sin_i(Tbl* tb, int & b)
440 {
441 
442  // Peut-etre rien a faire ?
443  if (tb->get_etat() == ETATZERO) {
444  int base_r = b & MSQ_R ;
445  int base_p = b & MSQ_P ;
446  b = base_r | base_p | T_COS_I ;
447  return ;
448  }
449 
450  // Protection
451  assert(tb->get_etat() == ETATQCQ) ;
452 
453  // Pour le confort
454  int nr = (tb->dim).dim[0] ; // Nombre
455  int nt = (tb->dim).dim[1] ; // de points
456  int np = (tb->dim).dim[2] ; // physiques REELS
457  np = np - 2 ; // Nombre de points physiques
458 
459  // Variables statiques
460  static double* cx = 0 ;
461  static int nt_pre =0 ;
462 
463  // Test sur np pour initialisation eventuelle
464  {
465  if (nt > nt_pre) {
466  nt_pre = nt ;
467  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
468  for (int i=0 ; i<nt ; i++) {
469  cx[i] = (2*i+1) ;
470  }
471  }
472  } // Fin de region critique
473 
474  // pt. sur le tableau de double resultat
475  double* xo = new double[(tb->dim).taille] ;
476 
477  // Initialisation a zero :
478  for (int i=0; i<(tb->dim).taille; i++) {
479  xo[i] = 0 ;
480  }
481 
482  // On y va...
483  double* xi = tb->t ;
484  double* xci = xi ; // Pointeurs
485  double* xco = xo ; // courants
486 
487  // k = 0
488  for (int j=0 ; j<nt ; j++) {
489  for (int i=0 ; i<nr ; i++ ) {
490  *xco = cx[j] * (*xci) ;
491  xci++ ;
492  xco++ ;
493  } // Fin de la boucle sur r
494  } // Fin de la boucle sur theta
495 
496  // k = 1
497  xci += nr*nt ;
498  xco += nr*nt ;
499 
500  // k >=2
501  for (int k=2 ; k<np+1 ; k++) {
502  for (int j=0 ; j<nt ; j++) {
503  for (int i=0 ; i<nr ; i++ ) {
504  *xco = cx[j] * (*xci) ;
505  xci++ ;
506  xco++ ;
507  } // Fin de la boucle sur r
508  } // Fin de la boucle sur theta
509  } // Fin de la boucle sur phi
510 
511  // On remet les choses la ou il faut
512  delete [] tb->t ;
513  tb->t = xo ;
514 
515  // base de developpement
516  int base_r = b & MSQ_R ;
517  int base_p = b & MSQ_P ;
518  b = base_r | base_p | T_COS_I ;
519 }
520 
521 // cas T_COS_I
522 //------------
523 void _dsdtet_t_cos_i(Tbl* tb, int & b)
524 {
525 
526  // Peut-etre rien a faire ?
527  if (tb->get_etat() == ETATZERO) {
528  int base_r = b & MSQ_R ;
529  int base_p = b & MSQ_P ;
530  b = base_r | base_p | T_SIN_I ;
531  return ;
532  }
533 
534  // Protection
535  assert(tb->get_etat() == ETATQCQ) ;
536 
537  // Pour le confort
538  int nr = (tb->dim).dim[0] ; // Nombre
539  int nt = (tb->dim).dim[1] ; // de points
540  int np = (tb->dim).dim[2] ; // physiques REELS
541  np = np - 2 ; // Nombre de points physiques
542 
543  // Variables statiques
544  static double* cx = 0 ;
545  static int nt_pre =0 ;
546 
547  // Test sur np pour initialisation eventuelle
548  {
549  if (nt > nt_pre) {
550  nt_pre = nt ;
551  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
552  for (int i=0 ; i<nt ; i++) {
553  cx[i] = - (2*i+1) ;
554  }
555  }
556  } // Fin de region critique
557 
558  // pt. sur le tableau de double resultat
559  double* xo = new double[(tb->dim).taille] ;
560 
561  // Initialisation a zero :
562  for (int i=0; i<(tb->dim).taille; i++) {
563  xo[i] = 0 ;
564  }
565 
566  // On y va...
567  double* xi = tb->t ;
568  double* xci = xi ; // Pointeurs
569  double* xco = xo ; // courants
570 
571  // k = 0
572  for (int j=0 ; j<nt ; j++) {
573  for (int i=0 ; i<nr ; i++ ) {
574  *xco = cx[j] * (*xci) ;
575  xci++ ;
576  xco++ ;
577  } // Fin de la boucle sur r
578  } // Fin de la boucle sur theta
579 
580  // k = 1
581  xci += nr*nt ;
582  xco += nr*nt ;
583 
584  // k >=2
585  for (int k=2 ; k<np+1 ; k++) {
586  for (int j=0 ; j<nt ; j++) {
587  for (int i=0 ; i<nr ; i++ ) {
588  *xco = cx[j] * (*xci) ;
589  xci++ ;
590  xco++ ;
591  } // Fin de la boucle sur r
592  } // Fin de la boucle sur theta
593  } // Fin de la boucle sur phi
594 
595  // On remet les choses la ou il faut
596  delete [] tb->t ;
597  tb->t = xo ;
598 
599  // base de developpement
600  int base_r = b & MSQ_R ;
601  int base_p = b & MSQ_P ;
602  b = base_r | base_p | T_SIN_I ;
603 }
604 
605 // cas T_COSSIN_CP
606 //----------------
607 void _dsdtet_t_cossin_cp(Tbl* tb, int & b)
608 {
609 
610  // Peut-etre rien a faire ?
611  if (tb->get_etat() == ETATZERO) {
612  int base_r = b & MSQ_R ;
613  int base_p = b & MSQ_P ;
614  b = base_r | base_p | T_COSSIN_SP ;
615  return ;
616  }
617 
618  // Protection
619  assert(tb->get_etat() == ETATQCQ) ;
620 
621  // Pour le confort
622  int nr = (tb->dim).dim[0] ; // Nombre
623  int nt = (tb->dim).dim[1] ; // de points
624  int np = (tb->dim).dim[2] ; // physiques REELS
625  np = np - 2 ; // Nombre de points physiques
626 
627  // Variables statiques
628  static double* cxp = 0 ;
629  static double* cxi = 0 ;
630  static int nt_pre = 0 ;
631 
632  // Test sur np pour initialisation eventuelle
633  {
634  if (nt > nt_pre) {
635  nt_pre = nt ;
636  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
637  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
638  for (int i=0 ; i<nt ; i++) {
639  cxp[i] = - (2*i) ;
640  cxi[i] = (2*i+1) ;
641  }
642  }
643  } // Fin de region critique
644 
645  // pt. sur le tableau de double resultat
646  double* xo = new double[(tb->dim).taille] ;
647 
648  // Initialisation a zero :
649  for (int i=0; i<(tb->dim).taille; i++) {
650  xo[i] = 0 ;
651  }
652 
653  // On y va...
654  double* xi = tb->t ;
655  double* xci = xi ; // Pointeurs
656  double* xco = xo ; // courants
657  double* cx[2] ; // Tableau des Pointeur de coefficient
658 
659  // Initialisation des pointeurs de coefficients
660  cx[0] = cxp ; // cos pairs pour m pair
661  cx[1] = cxi ; // sin impair pour m impair
662 
663  // k = 0
664  // Choix de la parite
665  double* cxl = cx[0] ; // Pointeur de coefficients local
666  for (int j=0 ; j<nt ; j++) {
667  for (int i=0 ; i<nr ; i++) {
668  *xco = cxl[j] * (*xci) ;
669  xci++ ;
670  xco++ ;
671  } // Fin de la Boucle sur r
672  } // Fin de la boucle sur theta
673 
674  // k = 1
675  xci += nr*nt ;
676  xco += nr*nt ;
677 
678  // k >= 2
679  for (int k=2 ; k<np+1 ; k++) {
680  // Choix de la parite
681  int m = (k/2) % 2 ;
682  cxl = cx[m] ; // Pointeur de coefficients local
683  for (int j=0 ; j<nt ; j++) {
684  for (int i=0 ; i<nr ; i++) {
685  *xco = cxl[j] * (*xci) ;
686  xci++ ;
687  xco++ ;
688  } // Fin de la Boucle sur r
689  } // Fin de la boucle sur theta
690  } // Fin de la boucle sur phi
691 
692  // On remet les choses la ou il faut
693  delete [] tb->t ;
694  tb->t = xo ;
695 
696  // base de developpement
697  int base_r = b & MSQ_R ;
698  int base_p = b & MSQ_P ;
699  b = base_r | base_p | T_COSSIN_SP ;
700 }
701 
702 // cas T_COSSIN_SP
703 //----------------
704 void _dsdtet_t_cossin_sp(Tbl* tb, int & b)
705 {
706 
707  // Peut-etre rien a faire ?
708  if (tb->get_etat() == ETATZERO) {
709  int base_r = b & MSQ_R ;
710  int base_p = b & MSQ_P ;
711  b = base_r | base_p | T_COSSIN_CP ;
712  return ;
713  }
714 
715  // Protection
716  assert(tb->get_etat() == ETATQCQ) ;
717 
718  // Pour le confort
719  int nr = (tb->dim).dim[0] ; // Nombre
720  int nt = (tb->dim).dim[1] ; // de points
721  int np = (tb->dim).dim[2] ; // physiques REELS
722  np = np - 2 ; // Nombre de points physiques
723 
724  // Variables statiques
725  static double* cxp = 0 ;
726  static double* cxi = 0 ;
727  static int nt_pre =0 ;
728 
729  // Test sur np pour initialisation eventuelle
730  {
731  if (nt > nt_pre) {
732  nt_pre = nt ;
733  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
734  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
735  for (int i=0 ; i<nt ; i++) {
736  cxp[i] = (2*i) ;
737  cxi[i] = - (2*i+1) ;
738  }
739  }
740  } // Fin de region critique
741 
742  // pt. sur le tableau de double resultat
743  double* xo = new double[(tb->dim).taille] ;
744 
745  // Initialisation a zero :
746  for (int i=0; i<(tb->dim).taille; i++) {
747  xo[i] = 0 ;
748  }
749 
750  // On y va...
751  double* xi = tb->t ;
752  double* xci = xi ; // Pointeurs
753  double* xco = xo ; // courants
754  double* cx[2] ; // Tableau des Pointeur de coefficient
755 
756  // Initialisation des pointeurs de coefficients
757  cx[0] = cxp ; // sin pairs pour m pair
758  cx[1] = cxi ; // cos impair pour m impair
759 
760  // k = 0
761  // Choix de la parite
762  double* cxl = cx[0] ; // Pointeur de coefficients local
763  for (int j=0 ; j<nt ; j++) {
764  for (int i=0 ; i<nr ; i++) {
765  *xco = cxl[j] * (*xci) ;
766  xci++ ;
767  xco++ ;
768  } // Fin de la Boucle sur r
769  } // Fin de la boucle sur theta
770 
771  // k = 1
772  xci += nr*nt ;
773  xco += nr*nt ;
774 
775  // k >= 2
776  for (int k=2 ; k<np+1 ; k++) {
777  // Choix de la parite
778  int m = (k/2) % 2 ;
779  cxl = cx[m] ; // Pointeur de coefficients local
780  for (int j=0 ; j<nt ; j++) {
781  for (int i=0 ; i<nr ; i++) {
782  *xco = cxl[j] * (*xci) ;
783  xci++ ;
784  xco++ ;
785  } // Fin de la Boucle sur r
786  } // Fin de la boucle sur theta
787  } // Fin de la boucle sur phi
788 
789  // On remet les choses la ou il faut
790  delete [] tb->t ;
791  tb->t = xo ;
792 
793  // base de developpement
794  int base_r = b & MSQ_R ;
795  int base_p = b & MSQ_P ;
796  b = base_r | base_p | T_COSSIN_CP ;
797 }
798 
799 // cas T_COSSIN_CI
800 //----------------
801 void _dsdtet_t_cossin_ci(Tbl* tb, int & b)
802 {
803 
804  // Peut-etre rien a faire ?
805  if (tb->get_etat() == ETATZERO) {
806  int base_r = b & MSQ_R ;
807  int base_p = b & MSQ_P ;
808  b = base_r | base_p | T_COSSIN_SI ;
809  return ;
810  }
811 
812  // Protection
813  assert(tb->get_etat() == ETATQCQ) ;
814 
815  // Pour le confort
816  int nr = (tb->dim).dim[0] ; // Nombre
817  int nt = (tb->dim).dim[1] ; // de points
818  int np = (tb->dim).dim[2] ; // physiques REELS
819  np = np - 2 ; // Nombre de points physiques
820 
821  // Variables statiques
822  static double* cxp = 0 ;
823  static double* cxi = 0 ;
824  static int nt_pre =0 ;
825 
826  // Test sur np pour initialisation eventuelle
827  {
828  if (nt > nt_pre) {
829  nt_pre = nt ;
830  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
831  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
832  for (int i=0 ; i<nt ; i++) {
833  cxp[i] = (2*i) ;
834  cxi[i] = - (2*i+1) ;
835  }
836  }
837  } // Fin de region critique
838 
839  // pt. sur le tableau de double resultat
840  double* xo = new double[(tb->dim).taille] ;
841 
842  // Initialisation a zero :
843  for (int i=0; i<(tb->dim).taille; i++) {
844  xo[i] = 0 ;
845  }
846 
847  // On y va...
848  double* xi = tb->t ;
849  double* xci = xi ; // Pointeurs
850  double* xco = xo ; // courants
851  double* cx[2] ; // Tableau des Pointeur de coefficient
852 
853  // Initialisation des pointeurs de coefficients
854  cx[0] = cxi ; // cos impairs pour m pair
855  cx[1] = cxp ; // sin pair pour m impair
856 
857  // k = 0
858  // Choix de la parite
859  double* cxl = cx[0] ; // Pointeur de coefficients local
860  for (int j=0 ; j<nt ; j++) {
861  for (int i=0 ; i<nr ; i++) {
862  *xco = cxl[j] * (*xci) ;
863  xci++ ;
864  xco++ ;
865  } // Fin de la Boucle sur r
866  } // Fin de la boucle sur theta
867 
868  // k = 1
869  xci += nr*nt ;
870  xco += nr*nt ;
871 
872  // k >= 2
873  for (int k=2 ; k<np+1 ; k++) {
874  // Choix de la parite
875  int m = (k/2) % 2 ;
876  cxl = cx[m] ; // Pointeur de coefficients local
877  for (int j=0 ; j<nt ; j++) {
878  for (int i=0 ; i<nr ; i++) {
879  *xco = cxl[j] * (*xci) ;
880  xci++ ;
881  xco++ ;
882  } // Fin de la Boucle sur r
883  } // Fin de la boucle sur theta
884  } // Fin de la boucle sur phi
885 
886  // On remet les choses la ou il faut
887  delete [] tb->t ;
888  tb->t = xo ;
889 
890  // base de developpement
891  int base_r = b & MSQ_R ;
892  int base_p = b & MSQ_P ;
893  b = base_r | base_p | T_COSSIN_SI ;
894 }
895 
896 // cas T_COSSIN_SI
897 //----------------
898 void _dsdtet_t_cossin_si(Tbl* tb, int & b)
899 {
900 
901  // Peut-etre rien a faire ?
902  if (tb->get_etat() == ETATZERO) {
903  int base_r = b & MSQ_R ;
904  int base_p = b & MSQ_P ;
905  b = base_r | base_p | T_COSSIN_CI ;
906  return ;
907  }
908 
909  // Protection
910  assert(tb->get_etat() == ETATQCQ) ;
911 
912  // Pour le confort
913  int nr = (tb->dim).dim[0] ; // Nombre
914  int nt = (tb->dim).dim[1] ; // de points
915  int np = (tb->dim).dim[2] ; // physiques REELS
916  np = np - 2 ; // Nombre de points physiques
917 
918  // Variables statiques
919  static double* cxp = 0 ;
920  static double* cxi = 0 ;
921  static int nt_pre =0 ;
922 
923  // Test sur np pour initialisation eventuelle
924  {
925  if (nt > nt_pre) {
926  nt_pre = nt ;
927  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
928  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
929  for (int i=0 ; i<nt ; i++) {
930  cxp[i] = - (2*i) ;
931  cxi[i] = (2*i+1) ;
932  }
933  }
934  } // Fin de region critique
935 
936  // pt. sur le tableau de double resultat
937  double* xo = new double[(tb->dim).taille] ;
938 
939  // Initialisation a zero :
940  for (int i=0; i<(tb->dim).taille; i++) {
941  xo[i] = 0 ;
942  }
943 
944  // On y va...
945  double* xi = tb->t ;
946  double* xci = xi ; // Pointeurs
947  double* xco = xo ; // courants
948  double* cx[2] ; // Tableau des Pointeur de coefficient
949 
950  // Initialisation des pointeurs de coefficients
951  cx[0] = cxi ; // sin impair pour m pair
952  cx[1] = cxp ; // cos pairs pour m impair
953 
954  // k = 0
955  // Choix de la parite
956  double* cxl = cx[0] ; // Pointeur de coefficients local
957  for (int j=0 ; j<nt ; j++) {
958  for (int i=0 ; i<nr ; i++) {
959  *xco = cxl[j] * (*xci) ;
960  xci++ ;
961  xco++ ;
962  } // Fin de la Boucle sur r
963  } // Fin de la boucle sur theta
964 
965  // k = 1
966  xci += nr*nt ;
967  xco += nr*nt ;
968 
969  // k >= 2
970  for (int k=2 ; k<np+1 ; k++) {
971  // Choix de la parite
972  int m = (k/2) % 2 ;
973  cxl = cx[m] ; // Pointeur de coefficients local
974  for (int j=0 ; j<nt ; j++) {
975  for (int i=0 ; i<nr ; i++) {
976  *xco = cxl[j] * (*xci) ;
977  xci++ ;
978  xco++ ;
979  } // Fin de la Boucle sur r
980  } // Fin de la boucle sur theta
981  } // Fin de la boucle sur phi
982 
983  // On remet les choses la ou il faut
984  delete [] tb->t ;
985  tb->t = xo ;
986 
987  // base de developpement
988  int base_r = b & MSQ_R ;
989  int base_p = b & MSQ_P ;
990  b = base_r | base_p | T_COSSIN_CI ;
991 }
992 
993 // cas T_COSSIN_C
994 //----------------
995 void _dsdtet_t_cossin_c(Tbl* tb, int & b)
996 {
997 
998  // Peut-etre rien a faire ?
999  if (tb->get_etat() == ETATZERO) {
1000  int base_r = b & MSQ_R ;
1001  int base_p = b & MSQ_P ;
1002  b = base_r | base_p | T_COSSIN_S ;
1003  return ;
1004  }
1005 
1006  // Protection
1007  assert(tb->get_etat() == ETATQCQ) ;
1008 
1009  // Pour le confort
1010  int nr = (tb->dim).dim[0] ; // Nombre
1011  int nt = (tb->dim).dim[1] ; // de points
1012  int np = (tb->dim).dim[2] ; // physiques REELS
1013  np = np - 2 ; // Nombre de points physiques
1014 
1015  // Variables statiques
1016  static double* cxp = 0 ;
1017  static double* cxi = 0 ;
1018  static int nt_pre = 0 ;
1019 
1020  // Test sur np pour initialisation eventuelle
1021  {
1022  if (nt > nt_pre) {
1023  nt_pre = nt ;
1024  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
1025  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
1026  for (int i=0 ; i<nt ; i++) {
1027  cxp[i] = - i ;
1028  cxi[i] = i ;
1029  }
1030  }
1031  } // Fin de region critique
1032 
1033  // pt. sur le tableau de double resultat
1034  double* xo = new double[(tb->dim).taille] ;
1035 
1036  // Initialisation a zero :
1037  for (int i=0; i<(tb->dim).taille; i++) {
1038  xo[i] = 0 ;
1039  }
1040 
1041  // On y va...
1042  double* xi = tb->t ;
1043  double* xci = xi ; // Pointeurs
1044  double* xco = xo ; // courants
1045  double* cx[2] ; // Tableau des Pointeur de coefficient
1046 
1047  // Initialisation des pointeurs de coefficients
1048  cx[0] = cxp ; // cos pairs pour m pair
1049  cx[1] = cxi ; // sin impair pour m impair
1050 
1051  // k = 0
1052  // Choix de la parite
1053  double* cxl = cx[0] ; // Pointeur de coefficients local
1054  for (int j=0 ; j<nt ; j++) {
1055  for (int i=0 ; i<nr ; i++) {
1056  *xco = cxl[j] * (*xci) ;
1057  xci++ ;
1058  xco++ ;
1059  } // Fin de la Boucle sur r
1060  } // Fin de la boucle sur theta
1061 
1062  // k = 1
1063  xci += nr*nt ;
1064  xco += nr*nt ;
1065 
1066  // k >= 2
1067  for (int k=2 ; k<np+1 ; k++) {
1068  // Choix de la parite
1069  int m = (k/2) % 2 ;
1070  cxl = cx[m] ; // Pointeur de coefficients local
1071  for (int j=0 ; j<nt ; j++) {
1072  for (int i=0 ; i<nr ; i++) {
1073  *xco = cxl[j] * (*xci) ;
1074  xci++ ;
1075  xco++ ;
1076  } // Fin de la Boucle sur r
1077  } // Fin de la boucle sur theta
1078  } // Fin de la boucle sur phi
1079 
1080  // On remet les choses la ou il faut
1081  delete [] tb->t ;
1082  tb->t = xo ;
1083 
1084  // base de developpement
1085  int base_r = b & MSQ_R ;
1086  int base_p = b & MSQ_P ;
1087  b = base_r | base_p | T_COSSIN_S ;
1088 }
1089 
1090 // cas T_COSSIN_S
1091 //----------------
1092 void _dsdtet_t_cossin_s(Tbl* tb, int & b)
1093 {
1094 
1095  // Peut-etre rien a faire ?
1096  if (tb->get_etat() == ETATZERO) {
1097  int base_r = b & MSQ_R ;
1098  int base_p = b & MSQ_P ;
1099  b = base_r | base_p | T_COSSIN_C ;
1100  return ;
1101  }
1102 
1103  // Protection
1104  assert(tb->get_etat() == ETATQCQ) ;
1105 
1106  // Pour le confort
1107  int nr = (tb->dim).dim[0] ; // Nombre
1108  int nt = (tb->dim).dim[1] ; // de points
1109  int np = (tb->dim).dim[2] ; // physiques REELS
1110  np = np - 2 ; // Nombre de points physiques
1111 
1112  // Variables statiques
1113  static double* cxp = 0 ;
1114  static double* cxi = 0 ;
1115  static int nt_pre =0 ;
1116 
1117  // Test sur np pour initialisation eventuelle
1118  {
1119  if (nt > nt_pre) {
1120  nt_pre = nt ;
1121  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
1122  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
1123  for (int i=0 ; i<nt ; i++) {
1124  cxp[i] = i ;
1125  cxi[i] = - i ;
1126  }
1127  }
1128  } // Fin de region critique
1129 
1130  // pt. sur le tableau de double resultat
1131  double* xo = new double[(tb->dim).taille] ;
1132 
1133  // Initialisation a zero :
1134  for (int i=0; i<(tb->dim).taille; i++) {
1135  xo[i] = 0 ;
1136  }
1137 
1138  // On y va...
1139  double* xi = tb->t ;
1140  double* xci = xi ; // Pointeurs
1141  double* xco = xo ; // courants
1142  double* cx[2] ; // Tableau des Pointeur de coefficient
1143 
1144  // Initialisation des pointeurs de coefficients
1145  cx[0] = cxp ; // sin pairs pour m pair
1146  cx[1] = cxi ; // cos impair pour m impair
1147 
1148  // k = 0
1149  // Choix de la parite
1150  double* cxl = cx[0] ; // Pointeur de coefficients local
1151  for (int j=0 ; j<nt ; j++) {
1152  for (int i=0 ; i<nr ; i++) {
1153  *xco = cxl[j] * (*xci) ;
1154  xci++ ;
1155  xco++ ;
1156  } // Fin de la Boucle sur r
1157  } // Fin de la boucle sur theta
1158 
1159  // k = 1
1160  xci += nr*nt ;
1161  xco += nr*nt ;
1162 
1163  // k >= 2
1164  for (int k=2 ; k<np+1 ; k++) {
1165  // Choix de la parite
1166  int m = (k/2) % 2 ;
1167  cxl = cx[m] ; // Pointeur de coefficients local
1168  for (int j=0 ; j<nt ; j++) {
1169  for (int i=0 ; i<nr ; i++) {
1170  *xco = cxl[j] * (*xci) ;
1171  xci++ ;
1172  xco++ ;
1173  } // Fin de la Boucle sur r
1174  } // Fin de la boucle sur theta
1175  } // Fin de la boucle sur phi
1176 
1177  // On remet les choses la ou il faut
1178  delete [] tb->t ;
1179  tb->t = xo ;
1180 
1181  // base de developpement
1182  int base_r = b & MSQ_R ;
1183  int base_p = b & MSQ_P ;
1184  b = base_r | base_p | T_COSSIN_C ;
1185 }
1186 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67
MSQ_P
#define MSQ_P
Extraction de l&#39;info sur Phi.
Definition: type_parite.h:156
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
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
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
MSQ_R
#define MSQ_R
Extraction de l&#39;info sur R.
Definition: type_parite.h:152
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition: type_parite.h:202
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
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
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
T_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194