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