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
24char 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//--------------------------------
81namespace Lorene {
82void _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//----------
90void _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//----------
170void _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//------------
236void _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//------------
316void _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//------------
382void _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//------------
448void _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//----------------
514void _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//----------------
603void _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//----------------
692void _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//----------------
781void _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//----------------
870void _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