LORENE
op_d2sdx2.C
1 /*
2  * Copyright (c) 1999-2000 Jean-Alain Marck
3  * Copyright (c) 1999-2001 Eric Gourgoulhon
4  * Copyright (c) 1999-2001 Philippe Grandclement
5  *
6  * This file is part of LORENE.
7  *
8  * LORENE is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * LORENE is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with LORENE; if not, write to the Free Software
20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21  *
22  */
23 
24 
25 
26 
27 /*
28  * Ensemble des routines de base de derivation seconde par rapport a r
29  * (Utilisation interne)
30  *
31  * void _d2sdx2_XXXX(Tbl * t, int & b)
32  * t pointeur sur le Tbl d'entree-sortie
33  * b base spectrale
34  *
35  */
36 
37 /*
38  * $Id: op_d2sdx2.C,v 1.8 2016/12/05 16:18:07 j_novak Exp $
39  * $Log: op_d2sdx2.C,v $
40  * Revision 1.8 2016/12/05 16:18:07 j_novak
41  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
42  *
43  * Revision 1.7 2015/03/05 08:49:32 j_novak
44  * Implemented operators with Legendre bases.
45  *
46  * Revision 1.6 2014/10/13 08:53:24 j_novak
47  * Lorene classes and functions now belong to the namespace Lorene.
48  *
49  * Revision 1.5 2013/06/14 15:54:06 j_novak
50  * Inclusion of Legendre bases.
51  *
52  * Revision 1.4 2008/08/27 08:50:10 jl_cornou
53  * Added Jacobi(0,2) polynomials
54  *
55  * Revision 1.3 2007/12/11 15:28:18 jl_cornou
56  * Jacobi(0,2) polynomials partially implemented
57  *
58  * Revision 1.2 2004/11/23 15:16:01 m_forot
59  *
60  * Added the bases for the cases without any equatorial symmetry
61  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
62  *
63  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
64  * LORENE
65  *
66  * Revision 2.7 2000/02/24 16:40:50 eric
67  * Initialisation a zero du tableau xo avant le calcul.
68  *
69  * Revision 2.6 1999/11/15 16:37:53 eric
70  * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
71  *
72  * Revision 2.5 1999/10/11 15:33:31 phil
73  * *** empty log message ***
74  *
75  * Revision 2.4 1999/10/11 14:25:47 phil
76  * realloc -> delete + new
77  *
78  * Revision 2.3 1999/09/13 11:30:23 phil
79  * gestion du cas k==1
80  *
81  * Revision 2.2 1999/03/01 15:06:31 eric
82  * *** empty log message ***
83  *
84  * Revision 2.1 1999/02/23 10:39:13 hyc
85  * *** empty log message ***
86  *
87  *
88  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_d2sdx2.C,v 1.8 2016/12/05 16:18:07 j_novak Exp $
89  *
90  */
91 
92 // Fichier includes
93 #include "tbl.h"
94 
95 namespace Lorene {
96 void _d2sdx2_1d_r_jaco02(int, double*, double*) ;
97 
98 // Prototypage
99 void c_est_pas_fait(char * ) ;
100 
101 // Routine pour les cas non prevus
102 //--------------------------------
103 void _d2sdx2_pas_prevu(Tbl* , int & b) {
104  cout << "d2sdx2 pas prevu..." << endl ;
105  cout << " base: " << b << endl ;
106  abort () ;
107 }
108 
109 // cas R_CHEBU - dzpuis = 0
110 //-------------------------
111 void _d2sdx2_r_chebu_0(Tbl *tb, int & )
112 {
113 
114  // Peut-etre rien a faire ?
115  if (tb->get_etat() == ETATZERO) {
116  return ;
117  }
118 
119  // Protection
120  assert(tb->get_etat() == ETATQCQ) ;
121 
122  // Pour le confort
123  int nr = (tb->dim).dim[0] ; // Nombre
124  int nt = (tb->dim).dim[1] ; // de points
125  int np = (tb->dim).dim[2] ; // physiques REELS
126  np = np - 2 ; // Nombre de points physiques
127 
128  // Variables statiques
129  static double* cx1 = 0x0 ;
130  static double* cx2 = 0x0 ;
131  static double* cx3 = 0x0 ;
132  static int nr_pre = 0 ;
133 
134  // Test sur np pour initialisation eventuelle
135  if (nr > nr_pre) {
136  nr_pre = nr ;
137  if (cx1 != 0x0) delete [] cx1 ;
138  if (cx2 != 0x0) delete [] cx2 ;
139  if (cx3 != 0x0) delete [] cx3 ;
140  cx1 = new double [nr] ;
141  cx2 = new double [nr] ;
142  cx3 = new double [nr] ;
143  for (int i=0 ; i<nr ; i++) {
144  cx1[i] = (i+2)*(i+2)*(i+2) ;
145  cx2[i] = (i+2) ;
146  cx3[i] = i*i ;
147  }
148  }
149 
150  // pt. sur le tableau de double resultat
151  double* xo = new double[(tb->dim).taille] ;
152 
153  // Initialisation a zero :
154  for (int i=0; i<(tb->dim).taille; i++) {
155  xo[i] = 0 ;
156  }
157 
158  // On y va...
159  double* xi = tb->t ;
160  double* xci = xi ; // Pointeurs
161  double* xco = xo ; // courants
162 
163  for (int k=0 ; k<np+1 ; k++)
164  if (k == 1) {
165  xci += nr*nt ;
166  xco += nr*nt ;
167  }
168  else {
169  for (int j=0 ; j<nt ; j++) {
170 
171  double som1, som2 ;
172 
173  xco[nr-1] = 0 ;
174  som1 = (nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
175  som2 = (nr-1) * xci[nr-1] ;
176  xco[nr-3] = som1 - (nr-3)*(nr-3)*som2 ;
177  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
178  som1 += cx1[i] * xci[i+2] ;
179  som2 += cx2[i] * xci[i+2] ;
180  xco[i] = som1 - cx3[i] * som2 ;
181  } // Fin de la premiere boucle sur r
182  xco[nr-2] = 0 ;
183  som1 = (nr-2)*(nr-2)*(nr-2) * xci[nr-2] ;
184  som2 = (nr-2) * xci[nr-2] ;
185  xco[nr-4] = som1 - (nr-4)*(nr-4)*som2 ;
186  for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
187  som1 += cx1[i] * xci[i+2] ;
188  som2 += cx2[i] * xci[i+2] ;
189  xco[i] = som1 - cx3[i] * som2 ;
190  } // Fin de la deuxieme boucle sur r
191  xco[0] *= .5 ;
192 
193  xci += nr ;
194  xco += nr ;
195  } // Fin de la boucle sur theta
196  } // Fin de la boucle sur phi
197 
198  // On remet les choses la ou il faut
199  delete [] tb->t ;
200  tb->t = xo ;
201 
202  // base de developpement
203  // inchangee
204 }
205 
206 // cas R_CHEB
207 //-----------
208 void _d2sdx2_r_cheb(Tbl *tb, int & )
209 {
210 
211  // Peut-etre rien a faire ?
212  if (tb->get_etat() == ETATZERO) {
213  return ;
214  }
215 
216  // Protection
217  assert(tb->get_etat() == ETATQCQ) ;
218 
219  // Pour le confort
220  int nr = (tb->dim).dim[0] ; // Nombre
221  int nt = (tb->dim).dim[1] ; // de points
222  int np = (tb->dim).dim[2] ; // physiques REELS
223  np = np - 2 ; // Nombre de points physiques
224 
225  // Variables statiques
226  static double* cx1 = 0x0 ;
227  static double* cx2 = 0x0 ;
228  static double* cx3 = 0x0 ;
229  static int nr_pre = 0 ;
230 
231  // Test sur np pour initialisation eventuelle
232  if (nr > nr_pre) {
233  nr_pre = nr ;
234  if (cx1 != 0x0) delete [] cx1 ;
235  if (cx2 != 0x0) delete [] cx2 ;
236  if (cx3 != 0x0) delete [] cx3 ;
237  cx1 = new double [nr] ;
238  cx2 = new double [nr] ;
239  cx3 = new double [nr] ;
240 
241  for (int i=0 ; i<nr ; i++) {
242  cx1[i] = (i+2)*(i+2)*(i+2) ;
243  cx2[i] = (i+2) ;
244  cx3[i] = i*i ;
245  }
246  }
247 
248  // pt. sur le tableau de double resultat
249  double* xo = new double[(tb->dim).taille] ;
250 
251  // Initialisation a zero :
252  for (int i=0; i<(tb->dim).taille; i++) {
253  xo[i] = 0 ;
254  }
255 
256  // On y va...
257  double* xi = tb->t ;
258  double* xci = xi ; // Pointeurs
259  double* xco = xo ; // courants
260 
261  for (int k=0 ; k<np+1 ; k++)
262  if (k == 1) {
263  xci += nr*nt ;
264  xco += nr*nt ;
265  }
266  else {
267  for (int j=0 ; j<nt ; j++) {
268 
269  double som1, som2 ;
270 
271  xco[nr-1] = 0 ;
272  som1 = (nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
273  som2 = (nr-1) * xci[nr-1] ;
274  xco[nr-3] = som1 - (nr-3)*(nr-3)*som2 ;
275  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
276  som1 += cx1[i] * xci[i+2] ;
277  som2 += cx2[i] * xci[i+2] ;
278  xco[i] = som1 - cx3[i] * som2 ;
279  } // Fin de la premiere boucle sur r
280  xco[nr-2] = 0 ;
281  som1 = (nr-2)*(nr-2)*(nr-2) * xci[nr-2] ;
282  som2 = (nr-2) * xci[nr-2] ;
283  xco[nr-4] = som1 - (nr-4)*(nr-4)*som2 ;
284  for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
285  som1 += cx1[i] * xci[i+2] ;
286  som2 += cx2[i] * xci[i+2] ;
287  xco[i] = som1 - cx3[i] * som2 ;
288  } // Fin de la deuxieme boucle sur r
289  xco[0] *= .5 ;
290 
291  xci += nr ;
292  xco += nr ;
293  } // Fin de la boucle sur theta
294  } // Fin de la boucle sur phi
295 
296  // On remet les choses la ou il faut
297  delete [] tb->t ;
298  tb->t = xo ;
299 
300  // base de developpement
301  // inchangee
302 }
303 
304 // cas R_CHEBP
305 //------------
306 void _d2sdx2_r_chebp(Tbl *tb, int & )
307 {
308 
309  // Peut-etre rien a faire ?
310  if (tb->get_etat() == ETATZERO) {
311  return ;
312  }
313 
314  // Protection
315  assert(tb->get_etat() == ETATQCQ) ;
316 
317  // Pour le confort
318  int nr = (tb->dim).dim[0] ; // Nombre
319  int nt = (tb->dim).dim[1] ; // de points
320  int np = (tb->dim).dim[2] ; // physiques REELS
321  np = np - 2 ; // Nombre de points physiques
322 
323  // Variables statiques
324  static double* cx1 = 0x0 ;
325  static double* cx2 = 0x0 ;
326  static double* cx3 = 0x0 ;
327  static int nr_pre = 0 ;
328 
329  // Test sur np pour initialisation eventuelle
330  if (nr > nr_pre) {
331  nr_pre = nr ;
332  if (cx1 != 0x0) delete [] cx1 ;
333  if (cx2 != 0x0) delete [] cx2 ;
334  if (cx3 != 0x0) delete [] cx3 ;
335  cx1 = new double [nr] ;
336  cx2 = new double [nr] ;
337  cx3 = new double [nr] ;
338  for (int i=0 ; i<nr ; i++) {
339  cx1[i] = 8*(i+1)*(i+1)*(i+1) ;
340  cx2[i] = 2*(i+1) ;
341  cx3[i] = 4*i*i ;
342  }
343  }
344  // pt. sur le tableau de double resultat
345  double* xo = new double[(tb->dim).taille] ;
346 
347  // Initialisation a zero :
348  for (int i=0; i<(tb->dim).taille; i++) {
349  xo[i] = 0 ;
350  }
351 
352  // On y va...
353  double* xi = tb->t ;
354  double* xci = xi ; // Pointeurs
355  double* xco = xo ; // courants
356 
357  for (int k=0 ; k<np+1 ; k++)
358  if (k == 1) {
359  xci += nr*nt ;
360  xco += nr*nt ;
361  }
362  else {
363  for (int j=0 ; j<nt ; j++) {
364 
365  double som1, som2 ;
366 
367  xco[nr-1] = 0 ;
368  som1 = 8*(nr-1)*(nr-1)*(nr-1) * xci[nr-1] ;
369  som2 = 2*(nr-1) * xci[nr-1] ;
370  xco[nr-2] = som1 - 4*(nr-2)*(nr-2)*som2 ;
371  for (int i = nr-3 ; i >= 0 ; i-- ) {
372  som1 += cx1[i] * xci[i+1] ;
373  som2 += cx2[i] * xci[i+1] ;
374  xco[i] = som1 - cx3[i] * som2 ;
375  } // Fin de la boucle sur r
376  xco[0] *= .5 ;
377 
378  xci += nr ;
379  xco += nr ;
380  } // Fin de la boucle sur theta
381  } // Fin de la boucle sur phi
382 
383  // On remet les choses la ou il faut
384  delete [] tb->t ;
385  tb->t = xo ;
386 
387  // base de developpement
388  // inchangee
389 }
390 
391 // cas R_CHEBI
392 //------------
393 void _d2sdx2_r_chebi(Tbl *tb, int & )
394 {
395 
396  // Peut-etre rien a faire ?
397  if (tb->get_etat() == ETATZERO) {
398  return ;
399  }
400 
401  // Protection
402  assert(tb->get_etat() == ETATQCQ) ;
403 
404  // Pour le confort
405  int nr = (tb->dim).dim[0] ; // Nombre
406  int nt = (tb->dim).dim[1] ; // de points
407  int np = (tb->dim).dim[2] ; // physiques REELS
408  np = np - 2 ; // Nombre de points physiques
409 
410  // Variables statiques
411  static double* cx1 = 0x0 ;
412  static double* cx2 = 0x0 ;
413  static double* cx3 = 0x0 ;
414  static int nr_pre = 0 ;
415 
416  // Test sur np pour initialisation eventuelle
417  if (nr > nr_pre) {
418  nr_pre = nr ;
419  if (cx1 != 0x0) delete [] cx1 ;
420  if (cx2 != 0x0) delete [] cx2 ;
421  if (cx3 != 0x0) delete [] cx3 ;
422  cx1 = new double [nr] ;
423  cx2 = new double [nr] ;
424  cx3 = new double [nr] ;
425 
426  for (int i=0 ; i<nr ; i++) {
427  cx1[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
428  cx2[i] = (2*i+3) ;
429  cx3[i] = (2*i+1)*(2*i+1) ;
430  }
431  }
432 
433  // pt. sur le tableau de double resultat
434  double* xo = new double[(tb->dim).taille] ;
435 
436  // Initialisation a zero :
437  for (int i=0; i<(tb->dim).taille; i++) {
438  xo[i] = 0 ;
439  }
440 
441  // On y va...
442  double* xi = tb->t ;
443  double* xci = xi ; // Pointeurs
444  double* xco = xo ; // courants
445 
446  for (int k=0 ; k<np+1 ; k++)
447  if (k == 1) {
448  xci += nr*nt ;
449  xco += nr*nt ;
450  }
451  else {
452  for (int j=0 ; j<nt ; j++) {
453 
454  double som1, som2 ;
455 
456  xco[nr-1] = 0 ;
457  som1 = (2*nr-1)*(2*nr-1)*(2*nr-1) * xci[nr-1] ;
458  som2 = (2*nr-1) * xci[nr-1] ;
459  xco[nr-2] = som1 - (2*nr-3)*(2*nr-3)*som2 ;
460  for (int i = nr-3 ; i >= 0 ; i-- ) {
461  som1 += cx1[i] * xci[i+1] ;
462  som2 += cx2[i] * xci[i+1] ;
463  xco[i] = som1 - cx3[i] * som2 ;
464  } // Fin de la boucle sur r
465 
466  xci += nr ;
467  xco += nr ;
468  } // Fin de la boucle sur theta
469  } // Fin de la boucle sur phi
470 
471  // On remet les choses la ou il faut
472  delete [] tb->t ;
473  tb->t = xo ;
474 
475  // base de developpement
476  // inchangee
477 }
478 
479 // cas R_CHEBPIM_P
480 //----------------
481 void _d2sdx2_r_chebpim_p(Tbl *tb, int & )
482 {
483 
484  // Peut-etre rien a faire ?
485  if (tb->get_etat() == ETATZERO) {
486  return ;
487  }
488 
489  // Protection
490  assert(tb->get_etat() == ETATQCQ) ;
491 
492  // Pour le confort
493  int nr = (tb->dim).dim[0] ; // Nombre
494  int nt = (tb->dim).dim[1] ; // de points
495  int np = (tb->dim).dim[2] ; // physiques REELS
496  np = np - 2 ; // Nombre de points physiques
497 
498  // Variables statiques
499  static double* cx1p = 0x0 ;
500  static double* cx2p = 0x0 ;
501  static double* cx3p = 0x0 ;
502  static double* cx1i = 0x0 ;
503  static double* cx2i = 0x0 ;
504  static double* cx3i = 0x0 ;
505  static int nr_pre = 0 ;
506 
507  // Test sur np pour initialisation eventuelle
508  if (nr > nr_pre) {
509  nr_pre = nr ;
510  if (cx1p != 0x0) delete [] cx1p ;
511  if (cx2p != 0x0) delete [] cx2p ;
512  if (cx3p != 0x0) delete [] cx3p ;
513  if (cx1i != 0x0) delete [] cx1i ;
514  if (cx2i != 0x0) delete [] cx2i ;
515  if (cx3i != 0x0) delete [] cx3i ;
516  cx1p = new double[nr] ;
517  cx2p = new double[nr] ;
518  cx3p = new double[nr] ;
519  cx1i = new double[nr] ;
520  cx2i = new double[nr] ;
521  cx3i = new double[nr] ;
522  for (int i=0 ; i<nr ; i++) {
523  cx1p[i] = 8*(i+1)*(i+1)*(i+1) ;
524  cx2p[i] = 2*(i+1) ;
525  cx3p[i] = 4*i*i ;
526 
527  cx1i[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
528  cx2i[i] = (2*i+3) ;
529  cx3i[i] = (2*i+1)*(2*i+1) ;
530  }
531  }
532 
533  double* cx1t[2] ;
534  double* cx2t[2] ;
535  double* cx3t[2] ;
536  cx1t[0] = cx1p ; cx1t[1] = cx1i ;
537  cx2t[0] = cx2p ; cx2t[1] = cx2i ;
538  cx3t[0] = cx3p ; cx3t[1] = cx3i ;
539 
540  // pt. sur le tableau de double resultat
541  double* xo = new double[(tb->dim).taille] ;
542 
543  // Initialisation a zero :
544  for (int i=0; i<(tb->dim).taille; i++) {
545  xo[i] = 0 ;
546  }
547 
548  // On y va...
549  double* xi = tb->t ;
550  double* xci ; // Pointeurs
551  double* xco ; // courants
552 
553  // Position depart
554  xci = xi ;
555  xco = xo ;
556 
557  double *cx1, *cx2, *cx3 ;
558 
559  // k = 0
560  cx1 = cx1t[0] ;
561  cx2 = cx2t[0] ;
562  cx3 = cx3t[0] ;
563  for (int j=0 ; j<nt ; j++) {
564 
565  double som1 = 0 ;
566  double som2 = 0 ;
567 
568  xco[nr-1] = 0 ;
569  for (int i = nr-2 ; i >= 0 ; i-- ) {
570  som1 += cx1[i] * xci[i+1] ;
571  som2 += cx2[i] * xci[i+1] ;
572  xco[i] = som1 - cx3[i] * som2 ;
573  } // Fin de la boucle sur r
574  xco[0] *= .5 ; // normalisation
575  xci += nr ; // au
576  xco += nr ; // suivant
577  } // Fin de la boucle sur theta
578 
579  // k = 1
580  xci += nr*nt ;
581  xco += nr*nt ;
582 
583  // k >= 2
584  for (int k=2 ; k<np+1 ; k++) {
585  int m = (k/2) % 2 ;
586  cx1 = cx1t[m] ;
587  cx2 = cx2t[m] ;
588  cx3 = cx3t[m] ;
589  for (int j=0 ; j<nt ; j++) {
590 
591  double som1 = 0 ;
592  double som2 = 0 ;
593 
594  xco[nr-1] = 0 ;
595  for (int i = nr-2 ; i >= 0 ; i-- ) {
596  som1 += cx1[i] * xci[i+1] ;
597  som2 += cx2[i] * xci[i+1] ;
598  xco[i] = som1 - cx3[i] * som2 ;
599  } // Fin de la boucle sur r
600  if (m == 0) xco[0] *= .5 ; // normalisation eventuelle
601  xci += nr ; // au
602  xco += nr ; // suivant
603  } // Fin de la boucle sur theta
604  } // Fin de la boucle sur phi
605 
606  // On remet les choses la ou il faut
607  delete [] tb->t ;
608  tb->t = xo ;
609 
610  // base de developpement
611  // inchangee
612 }
613 
614 // cas R_CHEBPIM_I
615 //----------------
616 void _d2sdx2_r_chebpim_i(Tbl *tb, int & )
617 {
618 
619  // Peut-etre rien a faire ?
620  if (tb->get_etat() == ETATZERO) {
621  return ;
622  }
623 
624  // Protection
625  assert(tb->get_etat() == ETATQCQ) ;
626 
627  // Pour le confort
628  int nr = (tb->dim).dim[0] ; // Nombre
629  int nt = (tb->dim).dim[1] ; // de points
630  int np = (tb->dim).dim[2] ; // physiques REELS
631  np = np - 2 ; // Nombre de points physiques
632 
633  // Variables statiques
634  static double* cx1p = 0x0 ;
635  static double* cx2p = 0x0 ;
636  static double* cx3p = 0x0 ;
637  static double* cx1i = 0x0 ;
638  static double* cx2i = 0x0 ;
639  static double* cx3i = 0x0 ;
640  static int nr_pre = 0 ;
641 
642  // Test sur np pour initialisation eventuelle
643  if (nr > nr_pre) {
644  nr_pre = nr ;
645  if (cx1p != 0x0) delete [] cx1p ;
646  if (cx2p != 0x0) delete [] cx2p ;
647  if (cx3p != 0x0) delete [] cx3p ;
648  if (cx1i != 0x0) delete [] cx1i ;
649  if (cx2i != 0x0) delete [] cx2i ;
650  if (cx3i != 0x0) delete [] cx3i ;
651  cx1p = new double[nr] ;
652  cx2p = new double[nr] ;
653  cx3p = new double[nr] ;
654  cx1i = new double[nr] ;
655  cx2i = new double[nr] ;
656  cx3i = new double[nr] ;
657  for (int i=0 ; i<nr ; i++) {
658  cx1p[i] = 8*(i+1)*(i+1)*(i+1) ;
659  cx2p[i] = 2*(i+1) ;
660  cx3p[i] = 4*i*i ;
661 
662  cx1i[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
663  cx2i[i] = (2*i+3) ;
664  cx3i[i] = (2*i+1)*(2*i+1) ;
665  }
666  }
667 
668  double* cx1t[2] ;
669  double* cx2t[2] ;
670  double* cx3t[2] ;
671  cx1t[1] = cx1p ; cx1t[0] = cx1i ;
672  cx2t[1] = cx2p ; cx2t[0] = cx2i ;
673  cx3t[1] = cx3p ; cx3t[0] = cx3i ;
674 
675  // pt. sur le tableau de double resultat
676  double* xo = new double[(tb->dim).taille] ;
677 
678  // Initialisation a zero :
679  for (int i=0; i<(tb->dim).taille; i++) {
680  xo[i] = 0 ;
681  }
682 
683  // On y va...
684  double* xi = tb->t ;
685  double* xci ; // Pointeurs
686  double* xco ; // courants
687 
688  // Position depart
689  xci = xi ;
690  xco = xo ;
691 
692  double *cx1, *cx2, *cx3 ;
693 
694  // k = 0
695  cx1 = cx1t[0] ;
696  cx2 = cx2t[0] ;
697  cx3 = cx3t[0] ;
698  for (int j=0 ; j<nt ; j++) {
699 
700  double som1 = 0 ;
701  double som2 = 0 ;
702 
703  xco[nr-1] = 0 ;
704  for (int i = nr-2 ; i >= 0 ; i-- ) {
705  som1 += cx1[i] * xci[i+1] ;
706  som2 += cx2[i] * xci[i+1] ;
707  xco[i] = som1 - cx3[i] * som2 ;
708  } // Fin de la boucle sur r
709  xci += nr ;
710  xco += nr ;
711  } // Fin de la boucle sur theta
712 
713  // k = 1
714  xci += nr*nt ;
715  xco += nr*nt ;
716 
717  // k >= 2
718  for (int k=2 ; k<np+1 ; k++) {
719  int m = (k/2) % 2 ;
720  cx1 = cx1t[m] ;
721  cx2 = cx2t[m] ;
722  cx3 = cx3t[m] ;
723  for (int j=0 ; j<nt ; j++) {
724 
725  double som1 = 0 ;
726  double som2 = 0 ;
727 
728  xco[nr-1] = 0 ;
729  for (int i = nr-2 ; i >= 0 ; i-- ) {
730  som1 += cx1[i] * xci[i+1] ;
731  som2 += cx2[i] * xci[i+1] ;
732  xco[i] = som1 - cx3[i] * som2 ;
733  } // Fin de la boucle sur r
734  if (m == 1) xco[0] *= .5 ;
735  xci += nr ;
736  xco += nr ;
737  } // Fin de la boucle sur theta
738  } // Fin de la boucle sur phi
739 
740  // On remet les choses la ou il faut
741  delete [] tb->t ;
742  tb->t = xo ;
743 
744  // base de developpement
745  // inchangee
746 }
747 
748 // cas R_CHEBPI_P
749 //----------------
750 void _d2sdx2_r_chebpi_p(Tbl *tb, int & )
751 {
752 
753  // Peut-etre rien a faire ?
754  if (tb->get_etat() == ETATZERO) {
755  return ;
756  }
757 
758  // Protection
759  assert(tb->get_etat() == ETATQCQ) ;
760 
761  // Pour le confort
762  int nr = (tb->dim).dim[0] ; // Nombre
763  int nt = (tb->dim).dim[1] ; // de points
764  int np = (tb->dim).dim[2] ; // physiques REELS
765  np = np - 2 ; // Nombre de points physiques
766 
767  // Variables statiques
768  static double* cx1p = 0x0 ;
769  static double* cx2p = 0x0 ;
770  static double* cx3p = 0x0 ;
771  static double* cx1i = 0x0 ;
772  static double* cx2i = 0x0 ;
773  static double* cx3i = 0x0 ;
774  static int nr_pre = 0 ;
775 
776  // Test sur np pour initialisation eventuelle
777  if (nr > nr_pre) {
778  nr_pre = nr ;
779  if (cx1p != 0x0) delete [] cx1p ;
780  if (cx2p != 0x0) delete [] cx2p ;
781  if (cx3p != 0x0) delete [] cx3p ;
782  if (cx1i != 0x0) delete [] cx1i ;
783  if (cx2i != 0x0) delete [] cx2i ;
784  if (cx3i != 0x0) delete [] cx3i ;
785  cx1p = new double[nr] ;
786  cx2p = new double[nr] ;
787  cx3p = new double[nr] ;
788  cx1i = new double[nr] ;
789  cx2i = new double[nr] ;
790  cx3i = new double[nr] ;
791  for (int i=0 ; i<nr ; i++) {
792  cx1p[i] = 8*(i+1)*(i+1)*(i+1) ;
793  cx2p[i] = 2*(i+1) ;
794  cx3p[i] = 4*i*i ;
795 
796  cx1i[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
797  cx2i[i] = (2*i+3) ;
798  cx3i[i] = (2*i+1)*(2*i+1) ;
799  }
800  }
801 
802  double* cx1t[2] ;
803  double* cx2t[2] ;
804  double* cx3t[2] ;
805  cx1t[0] = cx1p ; cx1t[1] = cx1i ;
806  cx2t[0] = cx2p ; cx2t[1] = cx2i ;
807  cx3t[0] = cx3p ; cx3t[1] = cx3i ;
808 
809  // pt. sur le tableau de double resultat
810  double* xo = new double[(tb->dim).taille] ;
811 
812  // Initialisation a zero :
813  for (int i=0; i<(tb->dim).taille; i++) {
814  xo[i] = 0 ;
815  }
816 
817  // On y va...
818  double* xi = tb->t ;
819  double* xci ; // Pointeurs
820  double* xco ; // courants
821 
822  // Position depart
823  xci = xi ;
824  xco = xo ;
825 
826  double *cx1, *cx2, *cx3 ;
827 
828  // k = 0
829  for (int j=0 ; j<nt ; j++) {
830  int l = j % 2 ;
831  cx1 = cx1t[l] ;
832  cx2 = cx2t[l] ;
833  cx3 = cx3t[l] ;
834  double som1 = 0 ;
835  double som2 = 0 ;
836 
837  xco[nr-1] = 0 ;
838  for (int i = nr-2 ; i >= 0 ; i-- ) {
839  som1 += cx1[i] * xci[i+1] ;
840  som2 += cx2[i] * xci[i+1] ;
841  xco[i] = som1 - cx3[i] * som2 ;
842  } // Fin de la boucle sur r
843  if (l == 0) xco[0] *= .5 ; // normalisation
844  xci += nr ; // au
845  xco += nr ; // suivant
846  } // Fin de la boucle sur theta
847 
848  // k = 1
849  xci += nr*nt ;
850  xco += nr*nt ;
851 
852  // k >= 2
853  for (int k=2 ; k<np+1 ; k++) {
854  for (int j=0 ; j<nt ; j++) {
855  int l = j % 2 ;
856  cx1 = cx1t[l] ;
857  cx2 = cx2t[l] ;
858  cx3 = cx3t[l] ;
859  double som1 = 0 ;
860  double som2 = 0 ;
861 
862  xco[nr-1] = 0 ;
863  for (int i = nr-2 ; i >= 0 ; i-- ) {
864  som1 += cx1[i] * xci[i+1] ;
865  som2 += cx2[i] * xci[i+1] ;
866  xco[i] = som1 - cx3[i] * som2 ;
867  } // Fin de la boucle sur r
868  if (l == 0) xco[0] *= .5 ; // normalisation eventuelle
869  xci += nr ; // au
870  xco += nr ; // suivant
871  } // Fin de la boucle sur theta
872  } // Fin de la boucle sur phi
873 
874  // On remet les choses la ou il faut
875  delete [] tb->t ;
876  tb->t = xo ;
877 
878  // base de developpement
879  // inchangee
880 }
881 
882 // cas R_CHEBPI_I
883 //----------------
884 void _d2sdx2_r_chebpi_i(Tbl *tb, int & )
885 {
886 
887  // Peut-etre rien a faire ?
888  if (tb->get_etat() == ETATZERO) {
889  return ;
890  }
891 
892  // Protection
893  assert(tb->get_etat() == ETATQCQ) ;
894 
895  // Pour le confort
896  int nr = (tb->dim).dim[0] ; // Nombre
897  int nt = (tb->dim).dim[1] ; // de points
898  int np = (tb->dim).dim[2] ; // physiques REELS
899  np = np - 2 ; // Nombre de points physiques
900 
901  // Variables statiques
902  static double* cx1p = 0x0 ;
903  static double* cx2p = 0x0 ;
904  static double* cx3p = 0x0 ;
905  static double* cx1i = 0x0 ;
906  static double* cx2i = 0x0 ;
907  static double* cx3i = 0x0 ;
908  static int nr_pre = 0 ;
909 
910  // Test sur np pour initialisation eventuelle
911  if (nr > nr_pre) {
912  nr_pre = nr ;
913  if (cx1p != 0x0) delete [] cx1p ;
914  if (cx2p != 0x0) delete [] cx2p ;
915  if (cx3p != 0x0) delete [] cx3p ;
916  if (cx1i != 0x0) delete [] cx1i ;
917  if (cx2i != 0x0) delete [] cx2i ;
918  if (cx3i != 0x0) delete [] cx3i ;
919  cx1p = new double[nr] ;
920  cx2p = new double[nr] ;
921  cx3p = new double[nr] ;
922  cx1i = new double[nr] ;
923  cx2i = new double[nr] ;
924  cx3i = new double[nr] ;
925  for (int i=0 ; i<nr ; i++) {
926  cx1p[i] = 8*(i+1)*(i+1)*(i+1) ;
927  cx2p[i] = 2*(i+1) ;
928  cx3p[i] = 4*i*i ;
929 
930  cx1i[i] = (2*i+3)*(2*i+3)*(2*i+3) ;
931  cx2i[i] = (2*i+3) ;
932  cx3i[i] = (2*i+1)*(2*i+1) ;
933  }
934  }
935 
936  double* cx1t[2] ;
937  double* cx2t[2] ;
938  double* cx3t[2] ;
939  cx1t[1] = cx1p ; cx1t[0] = cx1i ;
940  cx2t[1] = cx2p ; cx2t[0] = cx2i ;
941  cx3t[1] = cx3p ; cx3t[0] = cx3i ;
942 
943  // pt. sur le tableau de double resultat
944  double* xo = new double[(tb->dim).taille] ;
945 
946  // Initialisation a zero :
947  for (int i=0; i<(tb->dim).taille; i++) {
948  xo[i] = 0 ;
949  }
950 
951  // On y va...
952  double* xi = tb->t ;
953  double* xci ; // Pointeurs
954  double* xco ; // courants
955 
956  // Position depart
957  xci = xi ;
958  xco = xo ;
959 
960  double *cx1, *cx2, *cx3 ;
961 
962  // k = 0
963  for (int j=0 ; j<nt ; j++) {
964  int l = j % 2 ;
965  cx1 = cx1t[l] ;
966  cx2 = cx2t[l] ;
967  cx3 = cx3t[l] ;
968  double som1 = 0 ;
969  double som2 = 0 ;
970 
971  xco[nr-1] = 0 ;
972  for (int i = nr-2 ; i >= 0 ; i-- ) {
973  som1 += cx1[i] * xci[i+1] ;
974  som2 += cx2[i] * xci[i+1] ;
975  xco[i] = som1 - cx3[i] * som2 ;
976  } // Fin de la boucle sur r
977  if (l == 1) xco[0] *= .5 ;
978  xci += nr ;
979  xco += nr ;
980  } // Fin de la boucle sur theta
981 
982  // k = 1
983  xci += nr*nt ;
984  xco += nr*nt ;
985 
986  // k >= 2
987  for (int k=2 ; k<np+1 ; k++) {
988  for (int j=0 ; j<nt ; j++) {
989  int l = j % 2 ;
990  cx1 = cx1t[l] ;
991  cx2 = cx2t[l] ;
992  cx3 = cx3t[l] ;
993  double som1 = 0 ;
994  double som2 = 0 ;
995 
996  xco[nr-1] = 0 ;
997  for (int i = nr-2 ; i >= 0 ; i-- ) {
998  som1 += cx1[i] * xci[i+1] ;
999  som2 += cx2[i] * xci[i+1] ;
1000  xco[i] = som1 - cx3[i] * som2 ;
1001  } // Fin de la boucle sur r
1002  if (l == 1) xco[0] *= .5 ;
1003  xci += nr ;
1004  xco += nr ;
1005  } // Fin de la boucle sur theta
1006  } // Fin de la boucle sur phi
1007 
1008  // On remet les choses la ou il faut
1009  delete [] tb->t ;
1010  tb->t = xo ;
1011 
1012  // base de developpement
1013  // inchangee
1014 }
1015 
1016 
1017 // cas R_LEG
1018 //-----------
1019 void _d2sdx2_r_leg(Tbl *tb, int & )
1020 {
1021 
1022  // Peut-etre rien a faire ?
1023  if (tb->get_etat() == ETATZERO) {
1024  return ;
1025  }
1026 
1027  // Protection
1028  assert(tb->get_etat() == ETATQCQ) ;
1029 
1030  // Pour le confort
1031  int nr = (tb->dim).dim[0] ; // Nombre
1032  int nt = (tb->dim).dim[1] ; // de points
1033  int np = (tb->dim).dim[2] ; // physiques REELS
1034  np = np - 2 ; // Nombre de points physiques
1035 
1036  // Variables statiques
1037  static double* cx1 = 0x0 ;
1038  static double* cx2 = 0x0 ;
1039  static double* cx3 = 0x0 ;
1040  static int nr_pre = 0 ;
1041 
1042  // Test sur np pour initialisation eventuelle
1043  if (nr > nr_pre) {
1044  nr_pre = nr ;
1045  if (cx1 != 0x0) delete [] cx1 ;
1046  if (cx2 != 0x0) delete [] cx2 ;
1047  if (cx3 != 0x0) delete [] cx3 ;
1048  cx1 = new double [nr] ;
1049  cx2 = new double [nr] ;
1050  cx3 = new double [nr] ;
1051 
1052  for (int i=0 ; i<nr ; i++) {
1053  cx1[i] = (i+2)*(i+3) ;
1054  cx2[i] = i*(i+1) ;
1055  cx3[i] = double(i) + 0.5 ;
1056  }
1057  }
1058 
1059  // pt. sur le tableau de double resultat
1060  double* xo = new double[(tb->dim).taille] ;
1061 
1062  // Initialisation a zero :
1063  for (int i=0; i<(tb->dim).taille; i++) {
1064  xo[i] = 0 ;
1065  }
1066 
1067  // On y va...
1068  double* xi = tb->t ;
1069  double* xci = xi ; // Pointeurs
1070  double* xco = xo ; // courants
1071 
1072  for (int k=0 ; k<np+1 ; k++)
1073  if (k == 1) {
1074  xci += nr*nt ;
1075  xco += nr*nt ;
1076  }
1077  else {
1078  for (int j=0 ; j<nt ; j++) {
1079 
1080  double som1, som2 ;
1081 
1082  xco[nr-1] = 0 ;
1083  som1 = (nr-1)*nr * xci[nr-1] ;
1084  som2 = xci[nr-1] ;
1085  if (nr > 2)
1086  xco[nr-3] = (double(nr) -2.5) * (som1 - (nr-3)*(nr-2)*som2) ;
1087  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
1088  som1 += cx1[i] * xci[i+2] ;
1089  som2 += xci[i+2] ;
1090  xco[i] = cx3[i]*(som1 - cx2[i] * som2) ;
1091  } // Fin de la premiere boucle sur r
1092  if (nr > 1) xco[nr-2] = 0 ;
1093  if (nr > 3) {
1094  som1 = (nr-2)*(nr-1)* xci[nr-2] ;
1095  som2 = xci[nr-2] ;
1096  xco[nr-4] = (double(nr) - 3.5) * (som1 - (nr-4)*(nr-3)*som2) ;
1097  }
1098  for (int i = nr-6 ; i >= 0 ; i -= 2 ) {
1099  som1 += cx1[i] * xci[i+2] ;
1100  som2 += xci[i+2] ;
1101  xco[i] = cx3[i]*(som1 - cx2[i] * som2) ;
1102  } // Fin de la deuxieme boucle sur r
1103 
1104  xci += nr ;
1105  xco += nr ;
1106  } // Fin de la boucle sur theta
1107  } // Fin de la boucle sur phi
1108 
1109  // On remet les choses la ou il faut
1110  delete [] tb->t ;
1111  tb->t = xo ;
1112 
1113  // base de developpement
1114  // inchangee
1115 }
1116 
1117 // cas R_LEGP
1118 //------------
1119 void _d2sdx2_r_legp(Tbl *tb, int & )
1120 {
1121 
1122  // Peut-etre rien a faire ?
1123  if (tb->get_etat() == ETATZERO) {
1124  return ;
1125  }
1126 
1127  // Protection
1128  assert(tb->get_etat() == ETATQCQ) ;
1129 
1130  // Pour le confort
1131  int nr = (tb->dim).dim[0] ; // Nombre
1132  int nt = (tb->dim).dim[1] ; // de points
1133  int np = (tb->dim).dim[2] ; // physiques REELS
1134  np = np - 2 ; // Nombre de points physiques
1135 
1136  // Variables statiques
1137  static double* cx1 = 0x0 ;
1138  static double* cx2 = 0x0 ;
1139  static double* cx3 = 0x0 ;
1140  static int nr_pre = 0 ;
1141 
1142  // Test sur np pour initialisation eventuelle
1143  if (nr > nr_pre) {
1144  nr_pre = nr ;
1145  if (cx1 != 0x0) delete [] cx1 ;
1146  if (cx2 != 0x0) delete [] cx2 ;
1147  if (cx3 != 0x0) delete [] cx3 ;
1148  cx1 = new double [nr] ;
1149  cx2 = new double [nr] ;
1150  cx3 = new double [nr] ;
1151  for (int i=0 ; i<nr ; i++) {
1152  cx1[i] = (2*i+2)*(2*i+3) ;
1153  cx2[i] = 2*i*(2*i+1) ;
1154  cx3[i] = double(2*i) + 0.5 ;
1155  }
1156  }
1157  // pt. sur le tableau de double resultat
1158  double* xo = new double[(tb->dim).taille] ;
1159 
1160  // Initialisation a zero :
1161  for (int i=0; i<(tb->dim).taille; i++) {
1162  xo[i] = 0 ;
1163  }
1164 
1165  // On y va...
1166  double* xi = tb->t ;
1167  double* xci = xi ; // Pointeurs
1168  double* xco = xo ; // courants
1169 
1170  for (int k=0 ; k<np+1 ; k++)
1171  if (k == 1) {
1172  xci += nr*nt ;
1173  xco += nr*nt ;
1174  }
1175  else {
1176  for (int j=0 ; j<nt ; j++) {
1177 
1178  double som1, som2 ;
1179 
1180  xco[nr-1] = 0 ;
1181  som1 = (2*nr-2)*(2*nr-1)* xci[nr-1] ;
1182  som2 = xci[nr-1] ;
1183  if (nr > 1)
1184  xco[nr-2] = (double(2*nr) - 1.5)*(som1 - 2*(nr-2)*(2*nr-1)*som2) ;
1185  for (int i = nr-3 ; i >= 0 ; i-- ) {
1186  som1 += cx1[i] * xci[i+1] ;
1187  som2 += xci[i+1] ;
1188  xco[i] = cx3[i]*(som1 - cx2[i] * som2) ;
1189  } // Fin de la boucle sur r
1190 
1191  xci += nr ;
1192  xco += nr ;
1193  } // Fin de la boucle sur theta
1194  } // Fin de la boucle sur phi
1195 
1196  // On remet les choses la ou il faut
1197  delete [] tb->t ;
1198  tb->t = xo ;
1199 
1200  // base de developpement
1201  // inchangee
1202 }
1203 
1204 // cas R_LEGI
1205 //------------
1206 void _d2sdx2_r_legi(Tbl *tb, int & )
1207 {
1208 
1209  // Peut-etre rien a faire ?
1210  if (tb->get_etat() == ETATZERO) {
1211  return ;
1212  }
1213 
1214  // Protection
1215  assert(tb->get_etat() == ETATQCQ) ;
1216 
1217  // Pour le confort
1218  int nr = (tb->dim).dim[0] ; // Nombre
1219  int nt = (tb->dim).dim[1] ; // de points
1220  int np = (tb->dim).dim[2] ; // physiques REELS
1221  np = np - 2 ; // Nombre de points physiques
1222 
1223  // Variables statiques
1224  static double* cx1 = 0x0 ;
1225  static double* cx2 = 0x0 ;
1226  static double* cx3 = 0x0 ;
1227  static int nr_pre = 0 ;
1228 
1229  // Test sur np pour initialisation eventuelle
1230  if (nr > nr_pre) {
1231  nr_pre = nr ;
1232  if (cx1 != 0x0) delete [] cx1 ;
1233  if (cx2 != 0x0) delete [] cx2 ;
1234  if (cx3 != 0x0) delete [] cx3 ;
1235  cx1 = new double [nr] ;
1236  cx2 = new double [nr] ;
1237  cx3 = new double [nr] ;
1238 
1239  for (int i=0 ; i<nr ; i++) {
1240  cx1[i] = (2*i+3)*(2*i+4) ;
1241  cx2[i] = (2*i+1)*(2*i+2) ;
1242  cx3[i] = double(2*i) + 1.5 ;
1243  }
1244  }
1245 
1246  // pt. sur le tableau de double resultat
1247  double* xo = new double[(tb->dim).taille] ;
1248 
1249  // Initialisation a zero :
1250  for (int i=0; i<(tb->dim).taille; i++) {
1251  xo[i] = 0 ;
1252  }
1253 
1254  // On y va...
1255  double* xi = tb->t ;
1256  double* xci = xi ; // Pointeurs
1257  double* xco = xo ; // courants
1258 
1259  for (int k=0 ; k<np+1 ; k++)
1260  if (k == 1) {
1261  xci += nr*nt ;
1262  xco += nr*nt ;
1263  }
1264  else {
1265  for (int j=0 ; j<nt ; j++) {
1266 
1267  double som1, som2 ;
1268 
1269  xco[nr-1] = 0 ;
1270  som1 = (2*nr-1)*(2*nr) * xci[nr-1] ;
1271  som2 = xci[nr-1] ;
1272  if (nr > 1)
1273  xco[nr-2] = (double(nr) - 1.5)*(som1 - (2*nr-3)*(2*nr-2)*som2) ;
1274  for (int i = nr-3 ; i >= 0 ; i-- ) {
1275  som1 += cx1[i] * xci[i+1] ;
1276  som2 += xci[i+1] ;
1277  xco[i] = cx3[i]*(som1 - cx2[i] * som2) ;
1278  } // Fin de la boucle sur r
1279 
1280  xci += nr ;
1281  xco += nr ;
1282  } // Fin de la boucle sur theta
1283  } // Fin de la boucle sur phi
1284 
1285  // On remet les choses la ou il faut
1286  delete [] tb->t ;
1287  tb->t = xo ;
1288 
1289  // base de developpement
1290  // inchangee
1291 }
1292 
1293 
1294 
1295 // cas R_JACO02
1296 //-----------
1297 void _d2sdx2_r_jaco02(Tbl *tb, int & )
1298 {
1299 
1300  // Peut-etre rien a faire ?
1301  if (tb->get_etat() == ETATZERO) {
1302  return ;
1303  }
1304 
1305  // Protection
1306  assert(tb->get_etat() == ETATQCQ) ;
1307 
1308  // Pour le confort
1309  int nr = (tb->dim).dim[0] ; // Nombre
1310  int nt = (tb->dim).dim[1] ; // de points
1311  int np = (tb->dim).dim[2] ; // physiques REELS
1312  np = np - 2 ; // Nombre de points physiques
1313 
1314  // pt. sur le tableau de double resultat
1315  double* xo = new double[(tb->dim).taille] ;
1316 
1317  // Initialisation a zero :
1318  for (int i=0; i<(tb->dim).taille; i++) {
1319  xo[i] = 0 ;
1320  }
1321 
1322  // On y va...
1323  double* xi = tb->t ;
1324  double* xci = xi ; // Pointeurs
1325  double* xco = xo ; // courants
1326 
1327  for (int k=0 ; k<np+1 ; k++)
1328  if (k == 1) {
1329  xci += nr*nt ;
1330  xco += nr*nt ;
1331  }
1332  else {
1333  for (int j=0 ; j<nt ; j++) {
1334 
1335  double* tb1 = new double[nr] ;
1336  for (int m =0 ; m<nr ; m++) {
1337  tb1[m]=xci[m];
1338  }
1339  double* res = new double[nr] ;
1340  _d2sdx2_1d_r_jaco02(nr,tb1,res) ;
1341  for (int i = 0 ; i<nr ; i++ ) {
1342  xco[i] = res[i] ;
1343  }
1344  delete [] res ;
1345  delete [] tb1 ;
1346 
1347  xci += nr ;
1348  xco += nr ;
1349  } // Fin de la boucle sur theta
1350  } // Fin de la boucle sur phi
1351 
1352  // On remet les choses la ou il faut
1353  delete [] tb->t ;
1354  tb->t = xo ;
1355 
1356  // base de developpement
1357  // inchangee
1358 }
1359 
1360 
1361 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::c_est_pas_fait
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Definition: c_est_pas_fait.C:76