LORENE
val_dern_1d.C
1 /*
2  * Copyright (c) 2004 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 val_dern_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/val_dern_1d.C,v 1.2 2014/10/13 08:53:27 j_novak Exp $" ;
24 
25 /*
26  * $Id: val_dern_1d.C,v 1.2 2014/10/13 08:53:27 j_novak Exp $
27  * $Log: val_dern_1d.C,v $
28  * Revision 1.2 2014/10/13 08:53:27 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.1 2004/02/17 09:21:39 j_novak
32  * New functions for calculating values of the derivatives of a function
33  * using its Chebyshev coefficients.
34  *
35  *
36  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/val_dern_1d.C,v 1.2 2014/10/13 08:53:27 j_novak Exp $
37  *
38  */
39 
40 #include "type_parite.h"
41 #include "tbl.h"
42 
43 /*
44  * Functions computing value of f^(n) at boundaries of the interval [-1, 1],
45  * using the Chebyshev expansion of f. Note: n=0 works too.
46  *
47  * Input : 1-dimensional Tbl containing the Chebyshev coefficients of f.
48  * int base : base of spectral expansion.
49  *
50  * Output : double : the value of the n-th derivative of f at x=+/- 1.
51  *
52  */
53 
54 
55 namespace Lorene {
56 double val1_dern_1d(int n, const Tbl& tb, int base_r)
57 {
58 
59  //This function should be OK for any radial base
60  assert ( (base_r == R_CHEB) || (base_r == R_CHEBI) || (base_r == R_CHEBP) ||
61  (base_r == R_CHEBU) ) ;
62 
63  assert (n>=0) ;
64  assert (tb.get_ndim() == 1) ;
65  int nr = tb.get_dim(0) ;
66 
67  double resu = 0. ;
68 
69  int n_ini = ( (base_r == R_CHEBP) || (base_r == R_CHEBI) ) ? n / 2 : n ;
70 
71  double *tbi = &tb.t[n_ini] ;
72  for (int i=n_ini; i<nr; i++) {
73  double fact = 1. ;
74  int ii = i ;
75  if (base_r == R_CHEBP) ii *= 2 ;
76  if (base_r == R_CHEBI) ii = 2*i + 1 ;
77  for (int j=0; j<n; j++)
78  fact *= double(ii*ii - j*j)/double(2*j + 1) ;
79  resu += fact * (*tbi) ;
80  tbi++ ;
81  }
82 
83  return resu ;
84 }
85 
86 double valm1_dern_1d(int n, const Tbl& tb, int base_r)
87 {
88 
89  //This function should be OK for any radial base
90  assert ( (base_r == R_CHEB) || (base_r == R_CHEBI) || (base_r == R_CHEBP) ||
91  (base_r == R_CHEBU) ) ;
92 
93  assert (n>=0) ;
94  assert (tb.get_ndim() == 1) ;
95  int nr = tb.get_dim(0) ;
96 
97  double resu = 0. ;
98  double parite, fac ;
99  int n_ini ;
100  switch (base_r) {
101  case R_CHEBP:
102  n_ini = n / 2 ;
103  parite = 1 ;
104  fac = (n%2 == 0 ? 1 : -1) ;
105  break ;
106  case R_CHEBI:
107  n_ini = n / 2 ;
108  fac = (n%2 == 0 ? -1 : 1) ;
109  parite = 1 ;
110  break ;
111  default:
112  n_ini = n ;
113  parite = -1 ;
114  fac = 1 ;
115  break ;
116  }
117  double *tbi = &tb.t[n_ini] ;
118 
119  for (int i=n_ini; i<nr; i++) {
120  double fact = fac ;
121  int ii = i ;
122  if (base_r == R_CHEBP) ii *= 2 ;
123  if (base_r == R_CHEBI) ii = 2*i + 1 ;
124  for (int j=0; j<n; j++)
125  fact *= double(ii*ii - j*j)/double(2*j + 1) ;
126  resu += fact * (*tbi) ;
127  fac *= parite ;
128  tbi++ ;
129  }
130 
131  return resu ;
132 }
133 }
Lorene prototypes.
Definition: app_hor.h:64
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166