LORENE
op_primr.C
1 /*
2  * Computation of primitive in a single domain
3  *
4  *
5  */
6 
7 /*
8  * Copyright (c) 2004 Eric Gourgoulhon.
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License version 2
14  * as published by the Free Software Foundation.
15  *
16  * LORENE is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with LORENE; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24  *
25  */
26 
27 
28 
29 /*
30  * $Id: op_primr.C,v 1.13 2017/02/24 14:55:55 j_novak Exp $
31  * $Log: op_primr.C,v $
32  * Revision 1.13 2017/02/24 14:55:55 j_novak
33  * Corrected error in the primitive formula
34  *
35  * Revision 1.12 2017/02/22 17:11:33 j_novak
36  * Addition of new Legendre basis.
37  *
38  * Revision 1.11 2016/12/05 16:18:08 j_novak
39  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
40  *
41  * Revision 1.10 2014/10/13 08:53:26 j_novak
42  * Lorene classes and functions now belong to the namespace Lorene.
43  *
44  * Revision 1.9 2014/10/06 15:16:06 j_novak
45  * Modified #include directives to use c++ syntax.
46  *
47  * Revision 1.8 2013/04/25 15:46:06 j_novak
48  * Added special treatment in the case np = 1, for type_p = NONSYM.
49  *
50  * Revision 1.7 2007/12/21 13:59:02 j_novak
51  * Suppression of call to pow(-1, something).
52  *
53  * Revision 1.6 2007/12/11 15:28:18 jl_cornou
54  * Jacobi(0,2) polynomials partially implemented
55  *
56  * Revision 1.5 2006/05/17 15:01:16 j_novak
57  * Treatment of the case nr = 1 and R_CHEB
58  *
59  * Revision 1.4 2004/11/23 15:16:01 m_forot
60  *
61  * Added the bases for the cases without any equatorial symmetry
62  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
63  *
64  * Revision 1.3 2004/10/12 09:58:24 j_novak
65  * Better memory management.
66  *
67  * Revision 1.2 2004/06/14 15:24:57 e_gourgoulhon
68  * First operationnal version (tested).
69  *
70  * Revision 1.1 2004/06/13 21:33:13 e_gourgoulhon
71  * First version.
72  *
73  *
74  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_primr.C,v 1.13 2017/02/24 14:55:55 j_novak Exp $
75  *
76  */
77 
78 
79 // C headers
80 #include <cstdlib>
81 #include <cmath>
82 
83 // Lorene headers
84 #include "tbl.h"
85 
86 // Unexpected case
87 //----------------
88 namespace Lorene {
89 void _primr_pas_prevu(const Tbl&, int bin, const Tbl&, Tbl&, int&, Tbl& ) {
90 
91  cout << "Unexpected basis in primr : basis = " << hex << bin << endl ;
92  abort() ;
93 
94 }
95 
96 // case R_CHEB
97 //------------
98 void _primr_r_cheb(const Tbl& tin, int bin, const Tbl& valm1, Tbl& tout,
99  int& bout, Tbl& valp1) {
100 
101  assert(tin.dim == tout.dim) ;
102 
103  // Output spectral basis
104  bout = bin ;
105 
106  // Number of coefficients
107  int nr = tin.get_dim(0) ;
108  int nt = tin.get_dim(1) ;
109  int np = tin.get_dim(2) - 2 ;
110  int borne_phi = np + 1 ;
111  if (np == 1) borne_phi = 1 ;
112 
113  // Case of a zero input or pure angular grid
114  // -----------------------------------------
115  if ((tin.get_etat() == ETATZERO)||(nr == 1)) {
116  if (valm1.get_etat() == ETATZERO) {
117  tout.set_etat_zero() ;
118  valp1.set_etat_zero() ;
119  return ;
120  }
121  else {
122  assert(valm1.get_etat() == ETATQCQ) ;
123  tout.set_etat_qcq() ;
124  valp1.set_etat_qcq() ;
125  double* xco = tout.t ;
126  for (int k=0 ; k< borne_phi ; k++) {
127  if (k==1) { // jump over the coefficient of sin(0*phi)
128  xco += nr*nt ;
129  }
130  else {
131  for (int j=0 ; j<nt ; j++) {
132  xco[0] = valm1(k,j) ; // constant value = boundary value
133  for (int i=1; i<nr; i++) xco[i] = 0 ;
134  valp1.set(k,j) = xco[0] ;
135  xco += nr ;
136  }
137  }
138  }
139  return ;
140  }
141  }
142 
143  // Case of a non-zero input
144  // ------------------------
145 
146  assert(tin.get_etat() == ETATQCQ ) ;
147  tout.annule_hard() ;
148  valp1.annule_hard() ;
149 
150  const double* xci = tin.t ;
151  double* xco = tout.t ;
152 
153  for (int k=0 ; k< borne_phi ; k++) {
154  if (k==1) { // jump over the coefficient of sin(0*phi)
155  xci += nr*nt ;
156  xco += nr*nt ;
157  }
158  else {
159  for (int j=0 ; j<nt ; j++) {
160 
161  xco[1] = xci[0] - 0.5 * xci[2] ; // special case i = 1
162 
163  for (int i=2; i<nr-2; i++) {
164  xco[i] = (xci[i-1] - xci[i+1]) / double(2*i) ;
165  }
166 
167  xco[nr-2] = xci[nr-3] / double(2*nr - 4) ;
168  xco[nr-1] = xci[nr-2] / double(2*nr - 2) ;
169 
170  // Determination of the T_0 coefficient by matching with
171  // provided value at xi = - 1 :
172  double som = - xco[1] ;
173  for (int i=2; i<nr; i+=2) som += xco[i] ;
174  for (int i=3; i<nr; i+=2) som -= xco[i] ;
175  xco[0] = valm1(k,j) - som ;
176 
177  // Value of primitive at xi = + 1 :
178  som = xco[0] ;
179  for (int i=1; i<nr; i++) som += xco[i] ;
180  valp1.set(k,j) = som ;
181 
182  xci += nr ;
183  xco += nr ;
184  } // end of theta loop
185  }
186  } // end of phi loop
187 
188 }
189 
190 
191 
192 // case R_CHEBP
193 //-------------
194 void _primr_r_chebp(const Tbl& tin, int bin, const Tbl&, Tbl& tout,
195  int& bout, Tbl& valp1) {
196 
197  assert(tin.dim == tout.dim) ;
198 
199  // Output spectral basis
200  int base_t = bin & MSQ_T ;
201  int base_p = bin & MSQ_P ;
202  bout = base_p | base_t | R_CHEBI ;
203 
204  // Number of coefficients
205  int nr = tin.get_dim(0) ;
206  int nt = tin.get_dim(1) ;
207  int np = tin.get_dim(2) - 2 ;
208  int borne_phi = np + 1 ;
209  if (np == 1) borne_phi = 1 ;
210 
211  // Case of a zero input
212  // --------------------
213  if (tin.get_etat() == ETATZERO) {
214  tout.set_etat_zero() ;
215  valp1.set_etat_zero() ;
216  return ;
217  }
218 
219  // Case of a non-zero input
220  // ------------------------
221 
222  assert(tin.get_etat() == ETATQCQ ) ;
223  tout.annule_hard() ;
224  valp1.annule_hard() ;
225 
226  const double* xci = tin.t ;
227  double* xco = tout.t ;
228 
229  for (int k=0 ; k< borne_phi ; k++) {
230  if (k==1) { // jump over the coefficient of sin(0*phi)
231  xci += nr*nt ;
232  xco += nr*nt ;
233  }
234  else {
235  for (int j=0 ; j<nt ; j++) {
236 
237  xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
238 
239  for (int i=1; i<nr-2; i++) {
240  xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
241  }
242 
243  xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
244  xco[nr-1] = 0 ;
245 
246  // Value of primitive at xi = + 1 :
247  double som = xco[0] ;
248  for (int i=1; i<nr; i++) som += xco[i] ;
249  valp1.set(k,j) = som ;
250 
251  xci += nr ;
252  xco += nr ;
253  } // end of theta loop
254  }
255  } // end of phi loop
256 
257 }
258 
259 
260 // case R_CHEBI
261 //-------------
262 void _primr_r_chebi(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
263  int& bout, Tbl& valp1) {
264 
265  assert(tin.dim == tout.dim) ;
266 
267  // Output spectral basis
268  int base_t = bin & MSQ_T ;
269  int base_p = bin & MSQ_P ;
270  bout = base_p | base_t | R_CHEBP ;
271 
272  // Number of coefficients
273  int nr = tin.get_dim(0) ;
274  int nt = tin.get_dim(1) ;
275  int np = tin.get_dim(2) - 2 ;
276  int borne_phi = np + 1 ;
277  if (np == 1) borne_phi = 1 ;
278 
279 
280  // Case of a zero input
281  // --------------------
282  if (tin.get_etat() == ETATZERO) {
283  if (val0.get_etat() == ETATZERO) {
284  tout.set_etat_zero() ;
285  valp1.set_etat_zero() ;
286  return ;
287  }
288  else {
289  assert(val0.get_etat() == ETATQCQ) ;
290  tout.annule_hard() ;
291  valp1.annule_hard() ;
292  double* xco = tout.t ;
293  for (int k=0 ; k< borne_phi ; k++) {
294  if (k==1) { // jump over the coefficient of sin(0*phi)
295  xco += nr*nt ;
296  }
297  else {
298  for (int j=0 ; j<nt ; j++) {
299  xco[0] = val0(k,j) ; // constant value = boundary value
300  for (int i=1; i<nr; i++) xco[i] = 0 ;
301  valp1.set(k,j) = xco[0] ;
302  xco += nr ;
303  }
304  }
305  }
306  return ;
307  }
308  }
309 
310  // Case of a non-zero input
311  // ------------------------
312 
313  assert(tin.get_etat() == ETATQCQ ) ;
314  tout.annule_hard() ;
315  valp1.annule_hard() ;
316 
317  const double* xci = tin.t ;
318  double* xco = tout.t ;
319 
320  for (int k=0 ; k< borne_phi ; k++) {
321  if (k==1) { // jump over the coefficient of sin(0*phi)
322  xci += nr*nt ;
323  xco += nr*nt ;
324  }
325  else {
326  for (int j=0 ; j<nt ; j++) {
327 
328  for (int i=1; i<nr-1; i++) {
329  xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
330  }
331 
332  xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
333 
334  // Determination of the T_0 coefficient by maching with
335  // provided value at xi = 0 :
336  double som = - xco[1] ;
337  for (int i=2; i<nr; i+=2) som += xco[i] ;
338  for (int i=3; i<nr; i+=2) som -= xco[i] ;
339  xco[0] = val0(k,j) - som ;
340 
341  // Value of primitive at xi = + 1 :
342  som = xco[0] ;
343  for (int i=1; i<nr; i++) som += xco[i] ;
344  valp1.set(k,j) = som ;
345 
346  xci += nr ;
347  xco += nr ;
348  } // end of theta loop
349  }
350  } // end of phi loop
351 
352 }
353 
354 
355 
356 // case R_CHEBPIM_P
357 //-----------------
358 void _primr_r_chebpim_p(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
359  int& bout, Tbl& valp1) {
360 
361  assert(tin.dim == tout.dim) ;
362 
363  // Output spectral basis
364  int base_t = bin & MSQ_T ;
365  int base_p = bin & MSQ_P ;
366  bout = base_p | base_t | R_CHEBPIM_I ;
367 
368  // Number of coefficients
369  int nr = tin.get_dim(0) ;
370  int nt = tin.get_dim(1) ;
371  int np = tin.get_dim(2) - 2 ;
372  int borne_phi = np + 1 ;
373  if (np == 1) borne_phi = 1 ;
374 
375 
376  // Case of a zero input
377  // --------------------
378  if (tin.get_etat() == ETATZERO) {
379  if (val0.get_etat() == ETATZERO) {
380  tout.set_etat_zero() ;
381  valp1.set_etat_zero() ;
382  return ;
383  }
384  else {
385  assert(val0.get_etat() == ETATQCQ) ;
386  tout.annule_hard() ;
387  valp1.annule_hard() ;
388  double* xco = tout.t ;
389 
390  // m even part
391  for (int k=0 ; k<borne_phi ; k += 4) {
392  int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
393  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
394  if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
395  xco += nr*nt ;
396  }
397  else {
398  for (int j=0 ; j<nt ; j++) {
399  assert( val0(k+kmod,j) == double(0) ) ;
400  for (int i=0; i<nr; i++) xco[i] = 0 ;
401  valp1.set(k+kmod,j) = 0. ;
402  xco += nr ;
403  }
404  }
405  }
406  xco += 2*nr*nt ; // next even m
407  }
408 
409  // m odd part
410  xco = tout.t + 2*nr*nt ;
411  for (int k=2 ; k<borne_phi ; k += 4) {
412  int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
413  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
414  for (int j=0 ; j<nt ; j++) {
415  xco[0] = val0(k+kmod,j) ; // constant value = boundary value
416  for (int i=1; i<nr; i++) xco[i] = 0 ;
417  valp1.set(k+kmod,j) = xco[0] ;
418  xco += nr ;
419  }
420  }
421  xco += 2*nr*nt ; // next odd m
422  }
423  return ;
424  }
425  }
426 
427  // Case of a non-zero input
428  // ------------------------
429 
430  assert(tin.get_etat() == ETATQCQ ) ;
431  tout.annule_hard() ;
432  valp1.annule_hard() ;
433 
434  const double* xci = tin.t ;
435  double* xco = tout.t ;
436 
437  // m even part
438  // -----------
439  for (int k=0 ; k<borne_phi ; k += 4) {
440  int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
441  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
442  if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
443  xci += nr*nt ;
444  xco += nr*nt ;
445  }
446  else {
447  for (int j=0 ; j<nt ; j++) {
448  xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
449 
450  for (int i=1; i<nr-2; i++) {
451  xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
452  }
453 
454  xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
455  xco[nr-1] = 0 ;
456 
457  // Value of primitive at xi = + 1 :
458  double som = xco[0] ;
459  for (int i=1; i<nr; i++) som += xco[i] ;
460  valp1.set(k+kmod,j) = som ;
461 
462  xci += nr ;
463  xco += nr ;
464  } // end of theta loop
465 
466  }
467  }
468  xci += 2*nr*nt ; // next even m
469  xco += 2*nr*nt ; //
470  }
471 
472  // m odd part
473  // ----------
474  xci = tin.t + 2*nr*nt ;
475  xco = tout.t + 2*nr*nt ;
476  for (int k=2 ; k<borne_phi ; k += 4) {
477  int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
478  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
479  for (int j=0 ; j<nt ; j++) {
480 
481  for (int i=1; i<nr-1; i++) {
482  xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
483  }
484 
485  xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
486 
487  // Determination of the T_0 coefficient by maching with
488  // provided value at xi = 0 :
489  double som = - xco[1] ;
490  for (int i=2; i<nr; i+=2) som += xco[i] ;
491  for (int i=3; i<nr; i+=2) som -= xco[i] ;
492  xco[0] = val0(k+kmod,j) - som ;
493 
494  // Value of primitive at xi = + 1 :
495  som = xco[0] ;
496  for (int i=1; i<nr; i++) som += xco[i] ;
497  valp1.set(k+kmod,j) = som ;
498 
499  xci += nr ;
500  xco += nr ;
501  } // end of theta loop
502  }
503  xci += 2*nr*nt ; // next odd m
504  xco += 2*nr*nt ; //
505  }
506 
507 
508 }
509 
510 
511 
512 // case R_CHEBPIM_I
513 //-----------------
514 void _primr_r_chebpim_i(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
515  int& bout, Tbl& valp1) {
516 
517  assert(tin.dim == tout.dim) ;
518 
519  // Output spectral basis
520  int base_t = bin & MSQ_T ;
521  int base_p = bin & MSQ_P ;
522  bout = base_p | base_t | R_CHEBPIM_P ;
523 
524  // Number of coefficients
525  int nr = tin.get_dim(0) ;
526  int nt = tin.get_dim(1) ;
527  int np = tin.get_dim(2) - 2 ;
528  int borne_phi = np + 1 ;
529  if (np == 1) borne_phi = 1 ;
530 
531  // Case of a zero input
532  // --------------------
533  if (tin.get_etat() == ETATZERO) {
534  if (val0.get_etat() == ETATZERO) {
535  tout.set_etat_zero() ;
536  valp1.set_etat_zero() ;
537  return ;
538  }
539  else {
540  assert(val0.get_etat() == ETATQCQ) ;
541  tout.annule_hard() ;
542  valp1.annule_hard() ;
543  double* xco = tout.t ;
544 
545  // m odd part
546  for (int k=0 ; k<borne_phi ; k += 4) {
547  int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
548  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
549  if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
550  xco += nr*nt ;
551  }
552  else {
553  for (int j=0 ; j<nt ; j++) {
554  xco[0] = val0(k+kmod,j) ; // constant value = boundary value
555  for (int i=1; i<nr; i++) xco[i] = 0 ;
556  valp1.set(k+kmod,j) = xco[0] ;
557  xco += nr ;
558  }
559  }
560  }
561  xco += 2*nr*nt ; // next odd m
562  }
563 
564  // m even part
565  xco = tout.t + 2*nr*nt ;
566  for (int k=2 ; k<borne_phi ; k += 4) {
567  int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
568  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
569  for (int j=0 ; j<nt ; j++) {
570  assert( val0(k+kmod,j) == double(0) ) ;
571  for (int i=0; i<nr; i++) xco[i] = 0 ;
572  valp1.set(k+kmod,j) = 0. ;
573  xco += nr ;
574  }
575  }
576  xco += 2*nr*nt ; // next even m
577  }
578  return ;
579  }
580  }
581 
582  // Case of a non-zero input
583  // ------------------------
584 
585  assert(tin.get_etat() == ETATQCQ ) ;
586  tout.annule_hard() ;
587  valp1.annule_hard() ;
588 
589  const double* xci = tin.t ;
590  double* xco = tout.t ;
591 
592  // m odd part
593  // ----------
594  for (int k=0 ; k<borne_phi ; k += 4) {
595  int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
596  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
597  if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
598  xci += nr*nt ;
599  xco += nr*nt ;
600  }
601  else {
602  for (int j=0 ; j<nt ; j++) {
603 
604  for (int i=1; i<nr-1; i++) {
605  xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
606  }
607 
608  xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
609 
610  // Determination of the T_0 coefficient by maching with
611  // provided value at xi = 0 :
612  double som = - xco[1] ;
613  for (int i=2; i<nr; i+=2) som += xco[i] ;
614  for (int i=3; i<nr; i+=2) som -= xco[i] ;
615  xco[0] = val0(k+kmod,j) - som ;
616 
617  // Value of primitive at xi = + 1 :
618  som = xco[0] ;
619  for (int i=1; i<nr; i++) som += xco[i] ;
620  valp1.set(k+kmod,j) = som ;
621 
622  xci += nr ;
623  xco += nr ;
624  } // end of theta loop
625  }
626  }
627  xci += 2*nr*nt ; // next odd m
628  xco += 2*nr*nt ; //
629  }
630 
631  // m even part
632  // -----------
633  xci = tin.t + 2*nr*nt ;
634  xco = tout.t + 2*nr*nt ;
635  for (int k=2 ; k<borne_phi ; k += 4) {
636  int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
637  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
638 
639  for (int j=0 ; j<nt ; j++) {
640  xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
641 
642  for (int i=1; i<nr-2; i++) {
643  xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
644  }
645 
646  xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
647  xco[nr-1] = 0 ;
648 
649  // Value of primitive at xi = + 1 :
650  double som = xco[0] ;
651  for (int i=1; i<nr; i++) som += xco[i] ;
652  valp1.set(k+kmod,j) = som ;
653 
654  xci += nr ;
655  xco += nr ;
656  } // end of theta loop
657 
658  }
659  xci += 2*nr*nt ; // next even m
660  xco += 2*nr*nt ; //
661  }
662 
663 
664 }
665 
666 // case R_CHEBPI_P
667 //-------------
668 void _primr_r_chebpi_p(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
669  int& bout, Tbl& valp1) {
670 
671  assert(tin.dim == tout.dim) ;
672 
673  // Output spectral basis
674  int base_t = bin & MSQ_T ;
675  int base_p = bin & MSQ_P ;
676  bout = base_p | base_t | R_CHEBPI_I ;
677 
678  // Number of coefficients
679  int nr = tin.get_dim(0) ;
680  int nt = tin.get_dim(1) ;
681  int np = tin.get_dim(2) - 2 ;
682  int borne_phi = np + 1 ;
683  if (np == 1) borne_phi = 1 ;
684 
685 
686  // Case of a zero input
687  // --------------------
688  if (tin.get_etat() == ETATZERO) {
689  if (val0.get_etat() == ETATZERO) {
690  tout.set_etat_zero() ;
691  valp1.set_etat_zero() ;
692  return ;
693  }
694  else {
695  assert(val0.get_etat() == ETATQCQ) ;
696  tout.annule_hard() ;
697  valp1.annule_hard() ;
698  double* xco = tout.t ;
699  for (int k=0 ; k< borne_phi ; k++) {
700  if (k==1) { // jump over the coefficient of sin(0*phi)
701  xco += nr*nt ;
702  }
703  else {
704  for (int j=0 ; j<nt ; j++) {
705  int l = j%2;
706  if(l==0){
707  for (int i=0; i<nr; i++) xco[i] = 0 ;
708  valp1.set(k,j) = 0. ;
709  } else {
710  xco[0] = val0(k,j) ; // constant value = boundary value
711  for (int i=1; i<nr; i++) xco[i] = 0 ;
712  valp1.set(k,j) = xco[0] ;
713  }
714  xco += nr ;
715  }
716  }
717  }
718  return ;
719  }
720  }
721 
722  // Case of a non-zero input
723  // ------------------------
724 
725  assert(tin.get_etat() == ETATQCQ ) ;
726  tout.annule_hard() ;
727  valp1.annule_hard() ;
728 
729  const double* xci = tin.t ;
730  double* xco = tout.t ;
731 
732  for (int k=0 ; k< borne_phi ; k++) {
733  if (k==1) { // jump over the coefficient of sin(0*phi)
734  xci += nr*nt ;
735  xco += nr*nt ;
736  }
737  else {
738  for (int j=0 ; j<nt ; j++) {
739 
740  int l = j%2;
741  if(l==0){
742  xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
743 
744  for (int i=1; i<nr-2; i++) {
745  xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
746  }
747 
748  xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
749  xco[nr-1] = 0 ;
750 
751  // Value of primitive at xi = + 1 :
752  double som = xco[0] ;
753  for (int i=1; i<nr; i++) som += xco[i] ;
754  valp1.set(k,j) = som ;
755  } else {
756  for (int i=1; i<nr-1; i++) {
757  xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
758  }
759 
760  xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
761 
762  // Determination of the T_0 coefficient by maching with
763  // provided value at xi = 0 :
764  double som = - xco[1] ;
765  for (int i=2; i<nr; i+=2) som += xco[i] ;
766  for (int i=3; i<nr; i+=2) som -= xco[i] ;
767  xco[0] = val0(k,j) - som ;
768 
769  // Value of primitive at xi = + 1 :
770  som = xco[0] ;
771  for (int i=1; i<nr; i++) som += xco[i] ;
772  valp1.set(k,j) = som ;
773  }
774  xci += nr ;
775  xco += nr ;
776  } // end of theta loop
777  }
778  } // end of phi loop
779 
780 }
781 // case R_CHEBPI_I
782 //-------------
783 void _primr_r_chebpi_i(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
784  int& bout, Tbl& valp1) {
785 
786  assert(tin.dim == tout.dim) ;
787 
788  // Output spectral basis
789  int base_t = bin & MSQ_T ;
790  int base_p = bin & MSQ_P ;
791  bout = base_p | base_t | R_CHEBPI_P ;
792 
793  // Number of coefficients
794  int nr = tin.get_dim(0) ;
795  int nt = tin.get_dim(1) ;
796  int np = tin.get_dim(2) - 2 ;
797  int borne_phi = np + 1 ;
798  if (np == 1) borne_phi = 1 ;
799 
800 
801  // Case of a zero input
802  // --------------------
803  if (tin.get_etat() == ETATZERO) {
804  if (val0.get_etat() == ETATZERO) {
805  tout.set_etat_zero() ;
806  valp1.set_etat_zero() ;
807  return ;
808  }
809  else {
810  assert(val0.get_etat() == ETATQCQ) ;
811  tout.annule_hard() ;
812  valp1.annule_hard() ;
813  double* xco = tout.t ;
814  for (int k=0 ; k< borne_phi ; k++) {
815  if (k==1) { // jump over the coefficient of sin(0*phi)
816  xco += nr*nt ;
817  }
818  else {
819  for (int j=0 ; j<nt ; j++) {
820  int l = j%2;
821  if(l==1){
822  for (int i=0; i<nr; i++) xco[i] = 0 ;
823  valp1.set(k,j) = 0. ;
824  } else {
825  xco[0] = val0(k,j) ; // constant value = boundary value
826  for (int i=1; i<nr; i++) xco[i] = 0 ;
827  valp1.set(k,j) = xco[0] ;
828  }
829  xco += nr ;
830  }
831  }
832  }
833  return ;
834  }
835  }
836 
837  // Case of a non-zero input
838  // ------------------------
839 
840  assert(tin.get_etat() == ETATQCQ ) ;
841  tout.annule_hard() ;
842  valp1.annule_hard() ;
843 
844  const double* xci = tin.t ;
845  double* xco = tout.t ;
846 
847  for (int k=0 ; k< borne_phi ; k++) {
848  if (k==1) { // jump over the coefficient of sin(0*phi)
849  xci += nr*nt ;
850  xco += nr*nt ;
851  }
852  else {
853  for (int j=0 ; j<nt ; j++) {
854 
855  int l = j%2;
856  if(l==1){
857  xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
858 
859  for (int i=1; i<nr-2; i++) {
860  xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
861  }
862 
863  xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
864  xco[nr-1] = 0 ;
865 
866  // Value of primitive at xi = + 1 :
867  double som = xco[0] ;
868  for (int i=1; i<nr; i++) som += xco[i] ;
869  valp1.set(k,j) = som ;
870  } else {
871  for (int i=1; i<nr-1; i++) {
872  xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
873  }
874 
875  xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
876 
877  // Determination of the T_0 coefficient by maching with
878  // provided value at xi = 0 :
879  double som = - xco[1] ;
880  for (int i=2; i<nr; i+=2) som += xco[i] ;
881  for (int i=3; i<nr; i+=2) som -= xco[i] ;
882  xco[0] = val0(k,j) - som ;
883 
884  // Value of primitive at xi = + 1 :
885  som = xco[0] ;
886  for (int i=1; i<nr; i++) som += xco[i] ;
887  valp1.set(k,j) = som ;
888  }
889  xci += nr ;
890  xco += nr ;
891  } // end of theta loop
892  }
893  } // end of phi loop
894 
895 }
896 
897 
898 // case R_JACO02
899 //------------
900 void _primr_r_jaco02(const Tbl& tin, int bin, const Tbl& valm1, Tbl& tout,
901  int& bout, Tbl& valp1) {
902 
903  assert(tin.dim == tout.dim) ;
904 
905  // Output spectral basis
906  bout = bin ;
907 
908  // Number of coefficients
909  int nr = tin.get_dim(0) ;
910  int nt = tin.get_dim(1) ;
911  int np = tin.get_dim(2) - 2 ;
912  int borne_phi = np + 1 ;
913  if (np == 1) borne_phi = 1 ;
914 
915  // Case of a zero input or pure angular grid
916  // -----------------------------------------
917  if ((tin.get_etat() == ETATZERO)||(nr == 1)) {
918  if (valm1.get_etat() == ETATZERO) {
919  tout.set_etat_zero() ;
920  valp1.set_etat_zero() ;
921  return ;
922  }
923  else {
924  assert(valm1.get_etat() == ETATQCQ) ;
925  tout.set_etat_qcq() ;
926  valp1.set_etat_qcq() ;
927  double* xco = tout.t ;
928  for (int k=0 ; k< borne_phi ; k++) {
929  if (k==1) { // jump over the coefficient of sin(0*phi)
930  xco += nr*nt ;
931  }
932  else {
933  for (int j=0 ; j<nt ; j++) {
934  xco[0] = valm1(k,j) ; // constant value = boundary value
935  for (int i=1; i<nr; i++) xco[i] = 0 ;
936  valp1.set(k,j) = xco[0] ;
937  xco += nr ;
938  }
939  }
940  }
941  return ;
942  }
943  }
944 
945  // Case of a non-zero input
946  // ------------------------
947 
948  assert(tin.get_etat() == ETATQCQ ) ;
949  tout.annule_hard() ;
950  valp1.annule_hard() ;
951 
952  const double* xci = tin.t ;
953  double* xco = tout.t ;
954 
955  for (int k=0 ; k< borne_phi ; k++) {
956  if (k==1) { // jump over the coefficient of sin(0*phi)
957  xci += nr*nt ;
958  xco += nr*nt ;
959  }
960  else {
961  for (int j=0 ; j<nt ; j++) {
962 
963  for (int i=1; i<nr-1; i++) {
964  xco[i] = (i+2)/double((i+1)*(2*i+1))*xci[i-1] - xci[i]/double((i+1)*(i+2)) - (i+1)/double((i+2)*(2*i+5))*xci[i+1] ;
965  }
966  xco[nr-1] = (nr+1)/double((nr)*(2*nr-1))*xci[nr-2] - xci[nr-1]/double((nr)*(nr+1));
967 
968  // Determination of the J_0 coefficient by matching with
969  // provided value at xi = - 1 :
970 
971  double som = -3*xco[1] ;
972  for (int i=2; i<nr; i++) {
973  int signe = (i%2 == 0 ? 1 : -1) ;
974  som += xco[i]*signe*(i+1)*(i+2)/double(2) ;
975  }
976  xco[0] = valm1(k,j) - som ;
977 
978  // Value of primitive at xi = + 1 :
979  som = xco[0] ;
980  for (int i=1; i<nr; i++) som += xco[i] ;
981  valp1.set(k,j) = som ;
982 
983  xci += nr ;
984  xco += nr ;
985  } // end of theta loop
986  }
987  } // end of phi loop
988 
989 }
990 
991 
992 // case R_LEG
993 //-----------
994 void _primr_r_leg(const Tbl& tin, int bin, const Tbl& valm1, Tbl& tout,
995  int& bout, Tbl& valp1) {
996 
997  assert(tin.dim == tout.dim) ;
998 
999  // Output spectral basis
1000  bout = bin ;
1001 
1002  // Number of coefficients
1003  int nr = tin.get_dim(0) ;
1004  int nt = tin.get_dim(1) ;
1005  int np = tin.get_dim(2) - 2 ;
1006  int borne_phi = np + 1 ;
1007  if (np == 1) borne_phi = 1 ;
1008 
1009  // Case of a zero input or pure angular grid
1010  // -----------------------------------------
1011  if ((tin.get_etat() == ETATZERO)||(nr == 1)) {
1012  if (valm1.get_etat() == ETATZERO) {
1013  tout.set_etat_zero() ;
1014  valp1.set_etat_zero() ;
1015  return ;
1016  }
1017  else {
1018  assert(valm1.get_etat() == ETATQCQ) ;
1019  tout.set_etat_qcq() ;
1020  valp1.set_etat_qcq() ;
1021  double* xco = tout.t ;
1022  for (int k=0 ; k< borne_phi ; k++) {
1023  if (k==1) { // jump over the coefficient of sin(0*phi)
1024  xco += nr*nt ;
1025  }
1026  else {
1027  for (int j=0 ; j<nt ; j++) {
1028  xco[0] = valm1(k,j) ; // constant value = boundary value
1029  for (int i=1; i<nr; i++) xco[i] = 0 ;
1030  valp1.set(k,j) = xco[0] ;
1031  xco += nr ;
1032  }
1033  }
1034  }
1035  return ;
1036  }
1037  }
1038 
1039  // Case of a non-zero input
1040  // ------------------------
1041 
1042  assert(tin.get_etat() == ETATQCQ ) ;
1043  tout.annule_hard() ;
1044  valp1.annule_hard() ;
1045 
1046  const double* xci = tin.t ;
1047  double* xco = tout.t ;
1048 
1049  for (int k=0 ; k< borne_phi ; k++) {
1050  if (k==1) { // jump over the coefficient of sin(0*phi)
1051  xci += nr*nt ;
1052  xco += nr*nt ;
1053  }
1054  else {
1055  for (int j=0 ; j<nt ; j++) {
1056 
1057  for (int i=1; i<nr-2; i++) {
1058  xco[i] = xci[i-1] / double(2*i-1) - xci[i+1] / double(2*i+3) ;
1059  }
1060 
1061  xco[nr-2] = xci[nr-3] / double(2*nr - 5) ;
1062  xco[nr-1] = xci[nr-2] / double(2*nr - 3) ;
1063 
1064  // Determination of the T_0 coefficient by matching with
1065  // provided value at xi = - 1 :
1066  double som = - xco[1] ;
1067  for (int i=2; i<nr; i+=2) som += xco[i] ;
1068  for (int i=3; i<nr; i+=2) som -= xco[i] ;
1069  xco[0] = valm1(k,j) - som ;
1070 
1071  // Value of primitive at xi = + 1 :
1072  som = xco[0] ;
1073  for (int i=1; i<nr; i++) som += xco[i] ;
1074  valp1.set(k,j) = som ;
1075 
1076  xci += nr ;
1077  xco += nr ;
1078  } // end of theta loop
1079  }
1080  } // end of phi loop
1081 
1082 }
1083 
1084 
1085 
1086 // case R_LEGP
1087 //-------------
1088 void _primr_r_legp(const Tbl& tin, int bin, const Tbl&, Tbl& tout,
1089  int& bout, Tbl& valp1) {
1090 
1091  assert(tin.dim == tout.dim) ;
1092 
1093  // Output spectral basis
1094  int base_t = bin & MSQ_T ;
1095  int base_p = bin & MSQ_P ;
1096  bout = base_p | base_t | R_LEGI ;
1097 
1098  // Number of coefficients
1099  int nr = tin.get_dim(0) ;
1100  int nt = tin.get_dim(1) ;
1101  int np = tin.get_dim(2) - 2 ;
1102  int borne_phi = np + 1 ;
1103  if (np == 1) borne_phi = 1 ;
1104 
1105  // Case of a zero input
1106  // --------------------
1107  if ( (tin.get_etat() == ETATZERO) || (nr == 1) ){
1108  tout.set_etat_zero() ;
1109  valp1.set_etat_zero() ;
1110  return ;
1111  }
1112 
1113  // Case of a non-zero input
1114  // ------------------------
1115 
1116  assert(tin.get_etat() == ETATQCQ ) ;
1117  tout.annule_hard() ;
1118  valp1.annule_hard() ;
1119 
1120  const double* xci = tin.t ;
1121  double* xco = tout.t ;
1122 
1123  for (int k=0 ; k< borne_phi ; k++) {
1124  if (k==1) { // jump over the coefficient of sin(0*phi)
1125  xci += nr*nt ;
1126  xco += nr*nt ;
1127  }
1128  else {
1129  for (int j=0 ; j<nt ; j++) {
1130 
1131  for (int i=0; i<nr-2; i++) {
1132  xco[i] = xci[i]/ double(4*i+1) - xci[i+1]/double(4*i+5) ;
1133  }
1134 
1135  xco[nr-2] = xci[nr-2] / double(4*nr - 7) ;
1136  xco[nr-1] = 0 ;
1137 
1138  // Value of primitive at xi = + 1 :
1139  double som = xco[0] ;
1140  for (int i=1; i<nr; i++) som += xco[i] ;
1141  valp1.set(k,j) = som ;
1142 
1143  xci += nr ;
1144  xco += nr ;
1145  } // end of theta loop
1146  }
1147  } // end of phi loop
1148 
1149 }
1150 
1151 
1152 // case R_LEGI
1153 //------------
1154 void _primr_r_legi(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
1155  int& bout, Tbl& valp1) {
1156 
1157  assert(tin.dim == tout.dim) ;
1158 
1159  // Output spectral basis
1160  int base_t = bin & MSQ_T ;
1161  int base_p = bin & MSQ_P ;
1162  bout = base_p | base_t | R_LEGP ;
1163 
1164  // Number of coefficients
1165  int nr = tin.get_dim(0) ;
1166  int nt = tin.get_dim(1) ;
1167  int np = tin.get_dim(2) - 2 ;
1168  int borne_phi = np + 1 ;
1169  if (np == 1) borne_phi = 1 ;
1170 
1171 
1172  // Case of a zero input
1173  // --------------------
1174  if ( (tin.get_etat() == ETATZERO) || (nr == 1) ){
1175  if (val0.get_etat() == ETATZERO) {
1176  tout.set_etat_zero() ;
1177  valp1.set_etat_zero() ;
1178  return ;
1179  }
1180  else {
1181  assert(val0.get_etat() == ETATQCQ) ;
1182  tout.annule_hard() ;
1183  valp1.annule_hard() ;
1184  double* xco = tout.t ;
1185  for (int k=0 ; k< borne_phi ; k++) {
1186  if (k==1) { // jump over the coefficient of sin(0*phi)
1187  xco += nr*nt ;
1188  }
1189  else {
1190  for (int j=0 ; j<nt ; j++) {
1191  xco[0] = val0(k,j) ; // constant value = boundary value
1192  for (int i=1; i<nr; i++) xco[i] = 0 ;
1193  valp1.set(k,j) = xco[0] ;
1194  xco += nr ;
1195  }
1196  }
1197  }
1198  return ;
1199  }
1200  }
1201 
1202  // Case of a non-zero input
1203  // ------------------------
1204 
1205  assert(tin.get_etat() == ETATQCQ ) ;
1206  tout.annule_hard() ;
1207  valp1.annule_hard() ;
1208 
1209  const double* xci = tin.t ;
1210  double* xco = tout.t ;
1211 
1212  for (int k=0 ; k< borne_phi ; k++) {
1213  if (k==1) { // jump over the coefficient of sin(0*phi)
1214  xci += nr*nt ;
1215  xco += nr*nt ;
1216  }
1217  else {
1218  for (int j=0 ; j<nt ; j++) {
1219 
1220  for (int i=1; i<nr-1; i++) {
1221  xco[i] = xci[i-1]/ double(4*i-1) - xci[i]/double(4*i+3) ;
1222  }
1223 
1224  xco[nr-1] = xci[nr-2] / double(4*nr - 5) ;
1225 
1226  // Determination of the T_0 coefficient by matching with
1227  // provided value at xi = 0 :
1228  double val = -0.5 ;
1229  double som = val*xco[1] ;
1230  for (int i=2; i<nr; i++) {
1231  val *= -double(2*i-1) / double(2*i) ;
1232  som += val*xco[i] ;
1233  }
1234  xco[0] = val0(k,j) - som ;
1235 
1236  // Value of primitive at xi = + 1 :
1237  som = xco[0] ;
1238  for (int i=1; i<nr; i++) som += xco[i] ;
1239  valp1.set(k,j) = som ;
1240 
1241  xci += nr ;
1242  xco += nr ;
1243  } // end of theta loop
1244  }
1245  } // end of phi loop
1246 
1247 }
1248 
1249 
1250 
1251 
1252 }
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene
Lorene prototypes.
Definition: app_hor.h:67
MSQ_P
#define MSQ_P
Extraction de l&#39;info sur Phi.
Definition: type_parite.h:156
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
MSQ_T
#define MSQ_T
Extraction de l&#39;info sur Theta.
Definition: type_parite.h:154
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172