LORENE
sx_1d.C
1 /*
2  * Copyright (c) 2006 Jerome Novak
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char sx_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/sx_1d.C,v 1.4 2015/03/05 08:49:32 j_novak Exp $" ;
24 
25 /*
26  * $Id: sx_1d.C,v 1.4 2015/03/05 08:49:32 j_novak Exp $
27  * $Log: sx_1d.C,v $
28  * Revision 1.4 2015/03/05 08:49:32 j_novak
29  * Implemented operators with Legendre bases.
30  *
31  * Revision 1.3 2014/10/13 08:53:27 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.2 2014/10/06 15:16:07 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.1 2006/04/10 15:19:20 j_novak
38  * New definition of 1D operators dsdx and sx in the nucleus (bases R_CHEBP and
39  * R_CHEBI).
40  *
41  *
42  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/sx_1d.C,v 1.4 2015/03/05 08:49:32 j_novak Exp $
43  *
44  */
45 
46 
47  // Includes
48 #include <cstdlib>
49 
50 #include "headcpp.h"
51 #include "type_parite.h"
52 #include "proto.h"
53 
54 
55 /*
56  * 1/x operator (division by x)
57  *
58  * Only for the bases R_CHEBP and R_CHEBI
59  *
60  * Input :
61  *
62  * int nr : number of coefficients
63  * tb : the array of coefficients of the input function
64  *
65  * Output :
66  * tb : the array of coefficients of the result
67  */
68 
69 
70  //----------------------------
71  // Bases not implemented -----
72  //----------------------------
73 
74 namespace Lorene {
75 void _sx_1d_pas_prevu(int , double* , double *) {
76  cout << "sx_1d : base not implemented..." << endl ;
77  abort() ;
78  exit(-1) ;
79 }
80 
81 
82  //-------------------
83  // case R_CHEBP ---
84  //-------------------
85 
86 void _sx_1d_r_chebp (int nr, double* tb, double* res) {
87 
88  double som ;
89  int sign = 1 ;
90  res[nr-1] = 0 ;
91  som = 2 * sign * tb[nr-1] ;
92  res[nr-2] = som ;
93  for (int i=nr-3 ; i>=0 ; i--) {
94  sign = - sign ;
95  som += 2 * sign * tb[i+1] ;
96  res[i] = (i%2 == 0 ? -1 : 1) * som ;
97  }
98 
99 }
100 
101 
102  //-------------------
103  // case R_CHEBI ---
104  //-------------------
105 
106 void _sx_1d_r_chebi (int nr, double* tb, double* res) {
107 
108  double som ;
109  int sign = 1 ;
110  res[nr-1] = 0 ;
111  som = 2 * sign * tb[nr-2] ;
112  res[nr-2] = som ;
113  for (int i=nr-3 ; i>=0 ; i--) {
114  sign = - sign ;
115  som += 2 * sign * tb[i] ;
116  res[i] = (i%2 == 0 ? -1 : 1) * som ;
117  }
118  res[0] *= 0.5 ;
119 }
120 
121  //-------------------
122  // case R_CHEBU ---
123  //-------------------
124 
125 void _sx_1d_r_chebu (int nr, double* tb, double* res) {
126 
127  for (int i=0; i<nr; i++)
128  res[i] = tb[i] ;
129 
130  sxm1_1d_cheb(nr, res) ;
131 
132 }
133 
134  //------------------
135  // case R_LEGP ---
136  //------------------
137 
138 void _sx_1d_r_legp (int nr, double* tb, double* res) {
139 
140  double term = 0 ;
141  res[nr-1] = 0 ;
142  for (int i=nr-2 ; i>=0 ; i--) {
143  term += tb[i+1] ;
144  res[i] = double(4*i+3)/double(2*i+2)*term ;
145  term *= -double(2*i+1)/double(2*i+2) ;
146  }
147 
148 }
149 
150 
151  //------------------
152  // case R_LEGI ---
153  //------------------
154 
155 void _sx_1d_r_legi (int nr, double* tb, double* res) {
156 
157  double term = 0 ;
158  res[nr-1] = 0 ;
159  for (int i=nr-2; i>=0; i--) {
160  term += tb[i] ;
161  res[i] = double(4*i+1)/double(2*i+1)*term ;
162  term *= -double(2*i)/double(2*i+1) ;
163  }
164 }
165 
166 
167 
168  // ----------------------
169  // The calling function
170  //-----------------------
171 
172 void sx_1d(int nr, double **tb, int base_r)
173 {
174 
175  static void (*sx_1d[MAX_BASE])(int, double *, double *) ;
176  static int nap = 0 ;
177 
178  // Premier appel
179  if (nap==0) {
180  nap = 1 ;
181  for (int i=0 ; i<MAX_BASE ; i++) {
182  sx_1d[i] = _sx_1d_pas_prevu ;
183  }
184  // Les routines existantes
185  sx_1d[R_CHEBP >> TRA_R] = _sx_1d_r_chebp ;
186  sx_1d[R_CHEBI >> TRA_R] = _sx_1d_r_chebi ;
187  sx_1d[R_LEGP >> TRA_R] = _sx_1d_r_legp ;
188  sx_1d[R_LEGI >> TRA_R] = _sx_1d_r_legi ;
189  sx_1d[R_CHEBU >> TRA_R] = _sx_1d_r_chebu ;
190  }
191 
192 
193  double *result = new double[nr] ;
194  sx_1d[base_r](nr, *tb, result) ;
195  delete [] (*tb) ;
196  (*tb) = result ;
197 }
198 
199 }
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene prototypes.
Definition: app_hor.h:64