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