LORENE
op_d2sdtet2.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 seconde par rapport a theta
28  * (Utilisation interne)
29  *
30  * void _d2sdtet2_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_d2sdtet2.C,v 1.6 2016/12/05 16:18:07 j_novak Exp $
38  * $Log: op_d2sdtet2.C,v $
39  * Revision 1.6 2016/12/05 16:18:07 j_novak
40  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
41  *
42  * Revision 1.5 2014/10/13 08:53:24 j_novak
43  * Lorene classes and functions now belong to the namespace Lorene.
44  *
45  * Revision 1.4 2009/10/09 14:00:54 j_novak
46  * New bases T_cos and T_SIN.
47  *
48  * Revision 1.3 2006/03/10 12:45:38 j_novak
49  * Use of C++-style cast.
50  *
51  * Revision 1.2 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.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
57  * LORENE
58  *
59  * Revision 2.4 2000/10/04 11:50:20 eric
60  * Ajout d' abort() dans le cas non prevu.
61  *
62  * Revision 2.3 2000/02/24 16:40:42 eric
63  * Initialisation a zero du tableau xo avant le calcul.
64  *
65  * Revision 2.2 1999/11/15 16:37:44 eric
66  * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
67  *
68  * Revision 2.1 1999/03/01 15:06:00 eric
69  * *** empty log message ***
70  *
71  * Revision 2.0 1999/02/23 11:23:09 hyc
72  * *** empty log message ***
73  *
74  *
75  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_d2sdtet2.C,v 1.6 2016/12/05 16:18:07 j_novak Exp $
76  *
77  */
78 
79 // Headers Lorene
80 #include "tbl.h"
81 
82 // Routine pour les cas non prevus
83 //--------------------------------
84 namespace Lorene {
85 void _d2sdtet2_pas_prevu(Tbl* , int & b) {
86  cout << "Unknown theta basis in Mtbl_cf::d2sdt2() !" << endl ;
87  cout << " basis: " << hex << b << endl ;
88  abort() ;
89 }
90 
91 // cas T_COS
92 //----------
93 void _d2sdtet2_t_cos(Tbl* tb, int &)
94 {
95 
96  // Peut-etre rien a faire ?
97  if (tb->get_etat() == ETATZERO) {
98  return ;
99  }
100 
101  // Protection
102  assert(tb->get_etat() == ETATQCQ) ;
103 
104  // Pour le confort
105  int nr = (tb->dim).dim[0] ; // Nombre
106  int nt = (tb->dim).dim[1] ; // de points
107  int np = (tb->dim).dim[2] ; // physiques REELS
108  np = np - 2 ; // Nombre de points physiques
109 
110  // Variables statiques
111  static double* cx = 0 ;
112  static int nt_pre =0 ;
113 
114  // Test sur nt pour initialisation eventuelle
115  if (nt > nt_pre) {
116  nt_pre = nt ;
117  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
118  for (int i=0 ; i<nt ; i++) {
119  cx[i] = - double(i * i) ;
120  }
121  }
122 
123  // pt. sur le tableau de double resultat
124  double* xo = new double[(tb->dim).taille] ;
125 
126  // Initialisation a zero :
127  for (int i=0; i<(tb->dim).taille; i++) {
128  xo[i] = 0 ;
129  }
130 
131  // On y va...
132  double* xi = tb->t ;
133  double* xci = xi ; // Pointeurs
134  double* xco = xo ; // courants
135 
136  // k = 0
137  for (int j=0 ; j<nt ; j++) {
138  for (int i=0 ; i<nr ; i++ ) {
139  *xco = cx[j] * (*xci) ;
140  xci++ ;
141  xco++ ;
142  } // Fin de la boucle sur r
143  } // Fin de la boucle sur theta
144 
145  // k = 1
146  xci += nr*nt ;
147  xco += nr*nt ;
148 
149  // k >= 2
150  int borne_phi = np + 1 ;
151  if (np == 1) borne_phi = 1 ;
152 
153  for (int k=2 ; k<borne_phi ; k++) {
154  for (int j=0 ; j<nt ; j++) {
155  for (int i=0 ; i<nr ; i++ ) {
156  *xco = cx[j] * (*xci) ;
157  xci++ ;
158  xco++ ;
159  } // Fin de la boucle sur r
160  } // Fin de la boucle sur theta
161  } // Fin de la boucle sur phi
162 
163  // On remet les choses la ou il faut
164  delete [] tb->t ;
165  tb->t = xo ;
166 
167  // base de developpement
168  // inchangee
169 }
170 
171 // cas T_SIN
172 //----------
173 void _d2sdtet2_t_sin(Tbl* tb, int &)
174 {
175 
176  // Peut-etre rien a faire ?
177  if (tb->get_etat() == ETATZERO) {
178  return ;
179  }
180 
181  // Protection
182  assert(tb->get_etat() == ETATQCQ) ;
183 
184  // Pour le confort
185  int nr = (tb->dim).dim[0] ; // Nombre
186  int nt = (tb->dim).dim[1] ; // de points
187  int np = (tb->dim).dim[2] ; // physiques REELS
188  np = np - 2 ; // Nombre de points physiques
189 
190  // Variables statiques
191  static double* cx = 0 ;
192  static int nt_pre =0 ;
193 
194  // Test sur nt pour initialisation eventuelle
195  if (nt > nt_pre) {
196  nt_pre = nt ;
197  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
198  for (int i=0 ; i<nt ; i++) {
199  cx[i] = - double(i * i) ;
200  }
201  }
202 
203  // pt. sur le tableau de double resultat
204  double* xo = new double[(tb->dim).taille] ;
205 
206  // Initialisation a zero :
207  for (int i=0; i<(tb->dim).taille; i++) {
208  xo[i] = 0 ;
209  }
210 
211  // On y va...
212  double* xi = tb->t ;
213  double* xci = xi ; // Pointeurs
214  double* xco = xo ; // courants
215 
216  int borne_phi = np + 1 ;
217  if (np == 1) borne_phi = 1 ;
218 
219  for (int k=0 ; k< borne_phi ; k++) {
220  for (int j=0 ; j<nt ; j++) {
221  for (int i=0 ; i<nr ; i++ ) {
222  *xco = cx[j] * (*xci) ;
223  xci++ ;
224  xco++ ;
225  } // Fin de la boucle sur r
226  } // Fin de la boucle sur theta
227  } // Fin de la boucle sur phi
228 
229  // On remet les choses la ou il faut
230  delete [] tb->t ;
231  tb->t = xo ;
232 
233  // base de developpement
234  // inchangee
235 }
236 
237 // cas T_COS_P
238 //------------
239 void _d2sdtet2_t_cos_p(Tbl* tb, int &)
240 {
241 
242  // Peut-etre rien a faire ?
243  if (tb->get_etat() == ETATZERO) {
244  return ;
245  }
246 
247  // Protection
248  assert(tb->get_etat() == ETATQCQ) ;
249 
250  // Pour le confort
251  int nr = (tb->dim).dim[0] ; // Nombre
252  int nt = (tb->dim).dim[1] ; // de points
253  int np = (tb->dim).dim[2] ; // physiques REELS
254  np = np - 2 ; // Nombre de points physiques
255 
256  // Variables statiques
257  static double* cx = 0 ;
258  static int nt_pre =0 ;
259 
260  // Test sur nt pour initialisation eventuelle
261  if (nt > nt_pre) {
262  nt_pre = nt ;
263  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
264  for (int i=0 ; i<nt ; i++) {
265  cx[i] = - (2*i) * (2*i) ;
266  }
267  }
268 
269  // pt. sur le tableau de double resultat
270  double* xo = new double[(tb->dim).taille] ;
271 
272  // Initialisation a zero :
273  for (int i=0; i<(tb->dim).taille; i++) {
274  xo[i] = 0 ;
275  }
276 
277  // On y va...
278  double* xi = tb->t ;
279  double* xci = xi ; // Pointeurs
280  double* xco = xo ; // courants
281 
282  // k = 0
283  for (int j=0 ; j<nt ; j++) {
284  for (int i=0 ; i<nr ; i++ ) {
285  *xco = cx[j] * (*xci) ;
286  xci++ ;
287  xco++ ;
288  } // Fin de la boucle sur r
289  } // Fin de la boucle sur theta
290 
291  // k = 1
292  xci += nr*nt ;
293  xco += nr*nt ;
294 
295  // k >= 2
296  int borne_phi = np + 1 ;
297  if (np == 1) borne_phi = 1 ;
298 
299  for (int k=2 ; k<borne_phi ; k++) {
300  for (int j=0 ; j<nt ; j++) {
301  for (int i=0 ; i<nr ; i++ ) {
302  *xco = cx[j] * (*xci) ;
303  xci++ ;
304  xco++ ;
305  } // Fin de la boucle sur r
306  } // Fin de la boucle sur theta
307  } // Fin de la boucle sur phi
308 
309  // On remet les choses la ou il faut
310  delete [] tb->t ;
311  tb->t = xo ;
312 
313  // base de developpement
314  // inchangee
315 }
316 
317 // cas T_SIN_P
318 //------------
319 void _d2sdtet2_t_sin_p(Tbl* tb, int &)
320 {
321 
322  // Peut-etre rien a faire ?
323  if (tb->get_etat() == ETATZERO) {
324  return ;
325  }
326 
327  // Protection
328  assert(tb->get_etat() == ETATQCQ) ;
329 
330  // Pour le confort
331  int nr = (tb->dim).dim[0] ; // Nombre
332  int nt = (tb->dim).dim[1] ; // de points
333  int np = (tb->dim).dim[2] ; // physiques REELS
334  np = np - 2 ; // Nombre de points physiques
335 
336  // Variables statiques
337  static double* cx = 0 ;
338  static int nt_pre =0 ;
339 
340  // Test sur nt pour initialisation eventuelle
341  if (nt > nt_pre) {
342  nt_pre = nt ;
343  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
344  for (int i=0 ; i<nt ; i++) {
345  cx[i] = - (2*i) * (2*i) ;
346  }
347  }
348 
349  // pt. sur le tableau de double resultat
350  double* xo = new double[(tb->dim).taille] ;
351 
352  // Initialisation a zero :
353  for (int i=0; i<(tb->dim).taille; i++) {
354  xo[i] = 0 ;
355  }
356 
357  // On y va...
358  double* xi = tb->t ;
359  double* xci = xi ; // Pointeurs
360  double* xco = xo ; // courants
361 
362  int borne_phi = np + 1 ;
363  if (np == 1) borne_phi = 1 ;
364 
365  for (int k=0 ; k< borne_phi ; k++) {
366  for (int j=0 ; j<nt ; j++) {
367  for (int i=0 ; i<nr ; i++ ) {
368  *xco = cx[j] * (*xci) ;
369  xci++ ;
370  xco++ ;
371  } // Fin de la boucle sur r
372  } // Fin de la boucle sur theta
373  } // Fin de la boucle sur phi
374 
375  // On remet les choses la ou il faut
376  delete [] tb->t ;
377  tb->t = xo ;
378 
379  // base de developpement
380  // inchangee
381 }
382 
383 // cas T_SIN_I
384 //------------
385 void _d2sdtet2_t_sin_i(Tbl* tb, int &)
386 {
387 
388  // Peut-etre rien a faire ?
389  if (tb->get_etat() == ETATZERO) {
390  return ;
391  }
392 
393  // Protection
394  assert(tb->get_etat() == ETATQCQ) ;
395 
396  // Pour le confort
397  int nr = (tb->dim).dim[0] ; // Nombre
398  int nt = (tb->dim).dim[1] ; // de points
399  int np = (tb->dim).dim[2] ; // physiques REELS
400  np = np - 2 ; // Nombre de points physiques
401 
402  // Variables statiques
403  static double* cx = 0 ;
404  static int nt_pre =0 ;
405 
406  // Test sur nt pour initialisation eventuelle
407  if (nt > nt_pre) {
408  nt_pre = nt ;
409  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
410  for (int i=0 ; i<nt ; i++) {
411  cx[i] = - (2*i+1) * (2*i+1) ;
412  }
413  }
414 
415  // pt. sur le tableau de double resultat
416  double* xo = new double[(tb->dim).taille] ;
417 
418  // Initialisation a zero :
419  for (int i=0; i<(tb->dim).taille; i++) {
420  xo[i] = 0 ;
421  }
422 
423  // On y va...
424  double* xi = tb->t ;
425  double* xci = xi ; // Pointeurs
426  double* xco = xo ; // courants
427 
428  int borne_phi = np + 1 ;
429  if (np == 1) borne_phi = 1 ;
430 
431  for (int k=0 ; k< borne_phi ; k++) {
432  for (int j=0 ; j<nt ; j++) {
433  for (int i=0 ; i<nr ; i++ ) {
434  *xco = cx[j] * (*xci) ;
435  xci++ ;
436  xco++ ;
437  } // Fin de la boucle sur r
438  } // Fin de la boucle sur theta
439  } // Fin de la boucle sur phi
440 
441  // On remet les choses la ou il faut
442  delete [] tb->t ;
443  tb->t = xo ;
444 
445  // base de developpement
446  // inchangee
447 }
448 
449 // cas T_COS_I
450 //------------
451 void _d2sdtet2_t_cos_i(Tbl* tb, int &)
452 {
453 
454  // Peut-etre rien a faire ?
455  if (tb->get_etat() == ETATZERO) {
456  return ;
457  }
458 
459  // Protection
460  assert(tb->get_etat() == ETATQCQ) ;
461 
462  // Pour le confort
463  int nr = (tb->dim).dim[0] ; // Nombre
464  int nt = (tb->dim).dim[1] ; // de points
465  int np = (tb->dim).dim[2] ; // physiques REELS
466  np = np - 2 ; // Nombre de points physiques
467 
468  // Variables statiques
469  static double* cx = 0 ;
470  static int nt_pre =0 ;
471 
472  // Test sur nt pour initialisation eventuelle
473  if (nt > nt_pre) {
474  nt_pre = nt ;
475  cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
476  for (int i=0 ; i<nt ; i++) {
477  cx[i] = - (2*i+1) * (2*i+1) ;
478  }
479  }
480 
481  // pt. sur le tableau de double resultat
482  double* xo = new double[(tb->dim).taille] ;
483 
484  // Initialisation a zero :
485  for (int i=0; i<(tb->dim).taille; i++) {
486  xo[i] = 0 ;
487  }
488 
489  // On y va...
490  double* xi = tb->t ;
491  double* xci = xi ; // Pointeurs
492  double* xco = xo ; // courants
493 
494  int borne_phi = np + 1 ;
495  if (np == 1) borne_phi = 1 ;
496 
497  for (int k=0 ; k< borne_phi ; k++) {
498  for (int j=0 ; j<nt ; j++) {
499  for (int i=0 ; i<nr ; i++ ) {
500  *xco = cx[j] * (*xci) ;
501  xci++ ;
502  xco++ ;
503  } // Fin de la boucle sur r
504  } // Fin de la boucle sur theta
505  } // Fin de la boucle sur phi
506 
507  // On remet les choses la ou il faut
508  delete [] tb->t ;
509  tb->t = xo ;
510 
511  // base de developpement
512  // inchangee
513 }
514 
515 // cas T_COSSIN_CP
516 //----------------
517 void _d2sdtet2_t_cossin_cp(Tbl* tb, int &)
518 {
519 
520  // Peut-etre rien a faire ?
521  if (tb->get_etat() == ETATZERO) {
522  return ;
523  }
524 
525  // Protection
526  assert(tb->get_etat() == ETATQCQ) ;
527 
528  // Pour le confort
529  int nr = (tb->dim).dim[0] ; // Nombre
530  int nt = (tb->dim).dim[1] ; // de points
531  int np = (tb->dim).dim[2] ; // physiques REELS
532  np = np - 2 ; // Nombre de points physiques
533 
534  // Variables statiques
535  static double* cxp = 0 ;
536  static double* cxi = 0 ;
537  static int nt_pre =0 ;
538 
539  // Test sur nt pour initialisation eventuelle
540  if (nt > nt_pre) {
541  nt_pre = nt ;
542  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
543  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
544  for (int i=0 ; i<nt ; i++) {
545  cxp[i] = - (2*i) * (2*i) ;
546  cxi[i] = - (2*i+1) * (2*i+1) ;
547  }
548  }
549 
550  // pt. sur le tableau de double resultat
551  double* xo = new double[(tb->dim).taille] ;
552 
553  // Initialisation a zero :
554  for (int i=0; i<(tb->dim).taille; i++) {
555  xo[i] = 0 ;
556  }
557 
558  // On y va...
559  double* xi = tb->t ;
560  double* xci = xi ; // Pointeurs
561  double* xco = xo ; // courants
562 
563  // Partie cos(pair)
564  int k ;
565  for (k=0 ; k<np+1 ; k += 4) {
566  for (int m=0 ; m<2 ; m++) {
567  for (int j=0 ; j<nt ; j++) {
568  for (int i=0 ; i<nr ; i++ ) {
569  *xco = cxp[j] * (*xci) ;
570  xci++ ;
571  xco++ ;
572  } // Fin de la boucle sur r
573  } // Fin de la boucle sur theta
574  } // Fin de la boucle intermediaire
575  xci += 2*nr*nt ;
576  xco += 2*nr*nt ;
577  } // Fin de la boucle sur phi
578 
579  // Partie sin(impair)
580  xci = xi + 2*nr*nt ;
581  xco = xo + 2*nr*nt ;
582  for (k=2 ; k<np+1 ; k += 4) {
583  for (int m=0 ; m<2 ; m++) {
584  for (int j=0 ; j<nt ; j++) {
585  for (int i=0 ; i<nr ; i++ ) {
586  *xco = cxi[j] * (*xci) ;
587  xci++ ;
588  xco++ ;
589  } // Fin de la boucle sur r
590  } // Fin de la boucle sur theta
591  } // Fin de la boucle intermediaire
592  xci += 2*nr*nt ;
593  xco += 2*nr*nt ;
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  // inchangee
602 }
603 
604 // cas T_COSSIN_SP
605 //----------------
606 void _d2sdtet2_t_cossin_sp(Tbl* tb, int &)
607 {
608 
609  // Peut-etre rien a faire ?
610  if (tb->get_etat() == ETATZERO) {
611  return ;
612  }
613 
614  // Protection
615  assert(tb->get_etat() == ETATQCQ) ;
616 
617  // Pour le confort
618  int nr = (tb->dim).dim[0] ; // Nombre
619  int nt = (tb->dim).dim[1] ; // de points
620  int np = (tb->dim).dim[2] ; // physiques REELS
621  np = np - 2 ; // Nombre de points physiques
622 
623  // Variables statiques
624  static double* cxp = 0 ;
625  static double* cxi = 0 ;
626  static int nt_pre =0 ;
627 
628  // Test sur nt pour initialisation eventuelle
629  if (nt > nt_pre) {
630  nt_pre = nt ;
631  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
632  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
633  for (int i=0 ; i<nt ; i++) {
634  cxp[i] = - (2*i) * (2*i) ;
635  cxi[i] = - (2*i+1) * (2*i+1) ;
636  }
637  }
638 
639  // pt. sur le tableau de double resultat
640  double* xo = new double[(tb->dim).taille] ;
641 
642  // Initialisation a zero :
643  for (int i=0; i<(tb->dim).taille; i++) {
644  xo[i] = 0 ;
645  }
646 
647  // On y va...
648  double* xi = tb->t ;
649  double* xci = xi ; // Pointeurs
650  double* xco = xo ; // courants
651 
652  // Partie sin(pair)
653  int k ;
654  for (k=0 ; k<np+1 ; k += 4) {
655  for (int m=0 ; m<2 ; m++) {
656  for (int j=0 ; j<nt ; j++) {
657  for (int i=0 ; i<nr ; i++ ) {
658  *xco = cxp[j] * (*xci) ;
659  xci++ ;
660  xco++ ;
661  } // Fin de la boucle sur r
662  } // Fin de la boucle sur theta
663  } // Fin de la boucle intermediaire
664  xci += 2*nr*nt ;
665  xco += 2*nr*nt ;
666  } // Fin de la boucle sur phi
667 
668  // Partie cos(impair)
669  xci = xi + 2*nr*nt ;
670  xco = xo + 2*nr*nt ;
671  for (k=2 ; k<np+1 ; k += 4) {
672  for (int m=0 ; m<2 ; m++) {
673  for (int j=0 ; j<nt ; j++) {
674  for (int i=0 ; i<nr ; i++ ) {
675  *xco = cxi[j] * (*xci) ;
676  xci++ ;
677  xco++ ;
678  } // Fin de la boucle sur r
679  } // Fin de la boucle sur theta
680  } // Fin de la boucle intermediaire
681  xci += 2*nr*nt ;
682  xco += 2*nr*nt ;
683  } // Fin de la boucle sur phi
684 
685  // On remet les choses la ou il faut
686  delete [] tb->t ;
687  tb->t = xo ;
688 
689  // base de developpement
690  // inchangee
691 }
692 
693 // cas T_COSSIN_SI
694 //----------------
695 void _d2sdtet2_t_cossin_si(Tbl* tb, int &)
696 {
697 
698  // Peut-etre rien a faire ?
699  if (tb->get_etat() == ETATZERO) {
700  return ;
701  }
702 
703  // Protection
704  assert(tb->get_etat() == ETATQCQ) ;
705 
706  // Pour le confort
707  int nr = (tb->dim).dim[0] ; // Nombre
708  int nt = (tb->dim).dim[1] ; // de points
709  int np = (tb->dim).dim[2] ; // physiques REELS
710  np = np - 2 ; // Nombre de points physiques
711 
712  // Variables statiques
713  static double* cxp = 0 ;
714  static double* cxi = 0 ;
715  static int nt_pre =0 ;
716 
717  // Test sur nt pour initialisation eventuelle
718  if (nt > nt_pre) {
719  nt_pre = nt ;
720  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
721  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
722  for (int i=0 ; i<nt ; i++) {
723  cxp[i] = - (2*i) * (2*i) ;
724  cxi[i] = - (2*i+1) * (2*i+1) ;
725  }
726  }
727 
728  // pt. sur le tableau de double resultat
729  double* xo = new double[(tb->dim).taille] ;
730 
731  // Initialisation a zero :
732  for (int i=0; i<(tb->dim).taille; i++) {
733  xo[i] = 0 ;
734  }
735 
736  // On y va...
737  double* xi = tb->t ;
738  double* xci = xi ; // Pointeurs
739  double* xco = xo ; // courants
740 
741  // Partie sin(impair)
742  int k ;
743  for (k=0 ; k<np+1 ; k += 4) {
744  for (int m=0 ; m<2 ; m++) {
745  for (int j=0 ; j<nt ; j++) {
746  for (int i=0 ; i<nr ; i++ ) {
747  *xco = cxi[j] * (*xci) ;
748  xci++ ;
749  xco++ ;
750  } // Fin de la boucle sur r
751  } // Fin de la boucle sur theta
752  } // Fin de la boucle intermediaire
753  xci += 2*nr*nt ;
754  xco += 2*nr*nt ;
755  } // Fin de la boucle sur phi
756 
757  // Partie cos(pair)
758  xci = xi + 2*nr*nt ;
759  xco = xo + 2*nr*nt ;
760  for (k=2 ; k<np+1 ; k += 4) {
761  for (int m=0 ; m<2 ; m++) {
762  for (int j=0 ; j<nt ; j++) {
763  for (int i=0 ; i<nr ; i++ ) {
764  *xco = cxp[j] * (*xci) ;
765  xci++ ;
766  xco++ ;
767  } // Fin de la boucle sur r
768  } // Fin de la boucle sur theta
769  } // Fin de la boucle intermediaire
770  xci += 2*nr*nt ;
771  xco += 2*nr*nt ;
772  } // Fin de la boucle sur phi
773 
774  // On remet les choses la ou il faut
775  delete [] tb->t ;
776  tb->t = xo ;
777 
778  // base de developpement
779  // inchangee
780 }
781 
782 // cas T_COSSIN_C
783 //----------------
784 void _d2sdtet2_t_cossin_c(Tbl* tb, int &)
785 {
786 
787  // Peut-etre rien a faire ?
788  if (tb->get_etat() == ETATZERO) {
789  return ;
790  }
791 
792  // Protection
793  assert(tb->get_etat() == ETATQCQ) ;
794 
795  // Pour le confort
796  int nr = (tb->dim).dim[0] ; // Nombre
797  int nt = (tb->dim).dim[1] ; // de points
798  int np = (tb->dim).dim[2] ; // physiques REELS
799  np = np - 2 ; // Nombre de points physiques
800 
801  // Variables statiques
802  static double* cxp = 0 ;
803  static double* cxi = 0 ;
804  static int nt_pre =0 ;
805 
806  // Test sur nt pour initialisation eventuelle
807  if (nt > nt_pre) {
808  nt_pre = nt ;
809  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
810  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
811  for (int i=0 ; i<nt ; i++) {
812  cxp[i] = - i*i ;
813  cxi[i] = - i*i ;
814  }
815  }
816 
817  // pt. sur le tableau de double resultat
818  double* xo = new double[(tb->dim).taille] ;
819 
820  // Initialisation a zero :
821  for (int i=0; i<(tb->dim).taille; i++) {
822  xo[i] = 0 ;
823  }
824 
825  // On y va...
826  double* xi = tb->t ;
827  double* xci = xi ; // Pointeurs
828  double* xco = xo ; // courants
829 
830  // Partie cos(pair)
831  int k ;
832  for (k=0 ; k<np+1 ; k += 4) {
833  for (int m=0 ; m<2 ; m++) {
834  for (int j=0 ; j<nt ; j++) {
835  for (int i=0 ; i<nr ; i++ ) {
836  *xco = cxp[j] * (*xci) ;
837  xci++ ;
838  xco++ ;
839  } // Fin de la boucle sur r
840  } // Fin de la boucle sur theta
841  } // Fin de la boucle intermediaire
842  xci += 2*nr*nt ;
843  xco += 2*nr*nt ;
844  } // Fin de la boucle sur phi
845 
846  // Partie sin(impair)
847  xci = xi + 2*nr*nt ;
848  xco = xo + 2*nr*nt ;
849  for (k=2 ; k<np+1 ; k += 4) {
850  for (int m=0 ; m<2 ; m++) {
851  for (int j=0 ; j<nt ; j++) {
852  for (int i=0 ; i<nr ; i++ ) {
853  *xco = cxi[j] * (*xci) ;
854  xci++ ;
855  xco++ ;
856  } // Fin de la boucle sur r
857  } // Fin de la boucle sur theta
858  } // Fin de la boucle intermediaire
859  xci += 2*nr*nt ;
860  xco += 2*nr*nt ;
861  } // Fin de la boucle sur phi
862 
863  // On remet les choses la ou il faut
864  delete [] tb->t ;
865  tb->t = xo ;
866 
867  // base de developpement
868  // inchangee
869 }
870 
871 // cas T_COSSIN_S
872 //----------------
873 void _d2sdtet2_t_cossin_s(Tbl* tb, int &)
874 {
875 
876  // Peut-etre rien a faire ?
877  if (tb->get_etat() == ETATZERO) {
878  return ;
879  }
880 
881  // Protection
882  assert(tb->get_etat() == ETATQCQ) ;
883 
884  // Pour le confort
885  int nr = (tb->dim).dim[0] ; // Nombre
886  int nt = (tb->dim).dim[1] ; // de points
887  int np = (tb->dim).dim[2] ; // physiques REELS
888  np = np - 2 ; // Nombre de points physiques
889 
890  // Variables statiques
891  static double* cxp = 0 ;
892  static double* cxi = 0 ;
893  static int nt_pre =0 ;
894 
895  // Test sur nt pour initialisation eventuelle
896  if (nt > nt_pre) {
897  nt_pre = nt ;
898  cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
899  cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
900  for (int i=0 ; i<nt ; i++) {
901  cxp[i] = - i*i ;
902  cxi[i] = - i*i ;
903  }
904  }
905 
906  // pt. sur le tableau de double resultat
907  double* xo = new double[(tb->dim).taille] ;
908 
909  // Initialisation a zero :
910  for (int i=0; i<(tb->dim).taille; i++) {
911  xo[i] = 0 ;
912  }
913 
914  // On y va...
915  double* xi = tb->t ;
916  double* xci = xi ; // Pointeurs
917  double* xco = xo ; // courants
918 
919  // Partie sin(pair)
920  int k ;
921  for (k=0 ; k<np+1 ; k += 4) {
922  for (int m=0 ; m<2 ; m++) {
923  for (int j=0 ; j<nt ; j++) {
924  for (int i=0 ; i<nr ; i++ ) {
925  *xco = cxp[j] * (*xci) ;
926  xci++ ;
927  xco++ ;
928  } // Fin de la boucle sur r
929  } // Fin de la boucle sur theta
930  } // Fin de la boucle intermediaire
931  xci += 2*nr*nt ;
932  xco += 2*nr*nt ;
933  } // Fin de la boucle sur phi
934 
935  // Partie cos(impair)
936  xci = xi + 2*nr*nt ;
937  xco = xo + 2*nr*nt ;
938  for (k=2 ; k<np+1 ; k += 4) {
939  for (int m=0 ; m<2 ; m++) {
940  for (int j=0 ; j<nt ; j++) {
941  for (int i=0 ; i<nr ; i++ ) {
942  *xco = cxi[j] * (*xci) ;
943  xci++ ;
944  xco++ ;
945  } // Fin de la boucle sur r
946  } // Fin de la boucle sur theta
947  } // Fin de la boucle intermediaire
948  xci += 2*nr*nt ;
949  xco += 2*nr*nt ;
950  } // Fin de la boucle sur phi
951 
952  // On remet les choses la ou il faut
953  delete [] tb->t ;
954  tb->t = xo ;
955 
956  // base de developpement
957  // inchangee
958 }
959 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67