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