1 /* dht/dht.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* Author: G. Jungman
21 */
22 #include "gsl__config.h"
23 #include <stdlib.h>
24 #include "gsl_errno.h"
25 #include "gsl_math.h"
26 #include "gsl_sf_bessel.h"
27 #include "gsl_dht.h"
28
29
30 gsl_dht *
gsl_dht_alloc(size_t size)31 gsl_dht_alloc (size_t size)
32 {
33 gsl_dht * t;
34
35 if(size == 0) {
36 GSL_ERROR_VAL("size == 0", GSL_EDOM, 0);
37 }
38
39 t = (gsl_dht *)malloc(sizeof(gsl_dht));
40
41 if(t == 0) {
42 GSL_ERROR_VAL("out of memory", GSL_ENOMEM, 0);
43 }
44
45 t->size = size;
46
47 t->xmax = -1.0; /* Make it clear that this needs to be calculated. */
48 t->nu = -1.0;
49
50 t->j = (double *)malloc((size+2)*sizeof(double));
51
52 if(t->j == 0) {
53 free(t);
54 GSL_ERROR_VAL("could not allocate memory for j", GSL_ENOMEM, 0);
55 }
56
57 t->Jjj = (double *)malloc(size*(size+1)/2 * sizeof(double));
58
59 if(t->Jjj == 0) {
60 free(t->j);
61 free(t);
62 GSL_ERROR_VAL("could not allocate memory for Jjj", GSL_ENOMEM, 0);
63 }
64
65 t->J2 = (double *)malloc((size+1)*sizeof(double));
66
67 if(t->J2 == 0) {
68 free(t->Jjj);
69 free(t->j);
70 free(t);
71 GSL_ERROR_VAL("could not allocate memory for J2", GSL_ENOMEM, 0);
72 }
73
74 return t;
75 }
76
77 /* Handle internal calculation of Bessel zeros. */
78 static int
dht_bessel_zeros(gsl_dht * t)79 dht_bessel_zeros(gsl_dht * t)
80 {
81 unsigned int s;
82 gsl_sf_result z;
83 int stat_z = 0;
84 t->j[0] = 0.0;
85 for(s=1; s < t->size + 2; s++) {
86 stat_z += gsl_sf_bessel_zero_Jnu_e(t->nu, s, &z);
87 t->j[s] = z.val;
88 }
89 if(stat_z != 0) {
90 GSL_ERROR("could not compute bessel zeroes", GSL_EFAILED);
91 }
92 else {
93 return GSL_SUCCESS;
94 }
95 }
96
97 gsl_dht *
gsl_dht_new(size_t size,double nu,double xmax)98 gsl_dht_new (size_t size, double nu, double xmax)
99 {
100 int status;
101
102 gsl_dht * dht = gsl_dht_alloc (size);
103
104 if (dht == 0)
105 return 0;
106
107 status = gsl_dht_init(dht, nu, xmax);
108
109 if (status)
110 return 0;
111
112 return dht;
113 }
114
115 int
gsl_dht_init(gsl_dht * t,double nu,double xmax)116 gsl_dht_init(gsl_dht * t, double nu, double xmax)
117 {
118 if(xmax <= 0.0) {
119 GSL_ERROR ("xmax is not positive", GSL_EDOM);
120 } else if(nu < 0.0) {
121 GSL_ERROR ("nu is negative", GSL_EDOM);
122 }
123 else {
124 size_t n, m;
125 int stat_bz = GSL_SUCCESS;
126 int stat_J = 0;
127 double jN;
128
129 if(nu != t->nu) {
130 /* Recalculate Bessel zeros if necessary. */
131 t->nu = nu;
132 stat_bz = dht_bessel_zeros(t);
133 }
134
135 jN = t->j[t->size+1];
136
137 t->xmax = xmax;
138 t->kmax = jN / xmax;
139
140 t->J2[0] = 0.0;
141 for(m=1; m<t->size+1; m++) {
142 gsl_sf_result J;
143 stat_J += gsl_sf_bessel_Jnu_e(nu + 1.0, t->j[m], &J);
144 t->J2[m] = J.val * J.val;
145 }
146
147 /* J_nu(j[n] j[m] / j[N]) = Jjj[n(n-1)/2 + m - 1], 1 <= n,m <= size
148 */
149 for(n=1; n<t->size+1; n++) {
150 for(m=1; m<=n; m++) {
151 double arg = t->j[n] * t->j[m] / jN;
152 gsl_sf_result J;
153 stat_J += gsl_sf_bessel_Jnu_e(nu, arg, &J);
154 t->Jjj[n*(n-1)/2 + m - 1] = J.val;
155 }
156 }
157
158 if(stat_J != 0) {
159 GSL_ERROR("error computing bessel function", GSL_EFAILED);
160 }
161 else {
162 return stat_bz;
163 }
164 }
165 }
166
167
gsl_dht_x_sample(const gsl_dht * t,int n)168 double gsl_dht_x_sample(const gsl_dht * t, int n)
169 {
170 return t->j[n+1]/t->j[t->size+1] * t->xmax;
171 }
172
173
gsl_dht_k_sample(const gsl_dht * t,int n)174 double gsl_dht_k_sample(const gsl_dht * t, int n)
175 {
176 return t->j[n+1] / t->xmax;
177 }
178
179
gsl_dht_free(gsl_dht * t)180 void gsl_dht_free(gsl_dht * t)
181 {
182 free(t->J2);
183 free(t->Jjj);
184 free(t->j);
185 free(t);
186 }
187
188
189 int
gsl_dht_apply(const gsl_dht * t,double * f_in,double * f_out)190 gsl_dht_apply(const gsl_dht * t, double * f_in, double * f_out)
191 {
192 const double jN = t->j[t->size + 1];
193 const double r = t->xmax / jN;
194 size_t m;
195 size_t i;
196
197 for(m=0; m<t->size; m++) {
198 double sum = 0.0;
199 double Y;
200 for(i=0; i<t->size; i++) {
201 /* Need to find max and min so that we
202 * address the symmetric Jjj matrix properly.
203 * FIXME: we can presumably optimize this
204 * by just running over the elements of Jjj
205 * in a deterministic manner.
206 */
207 size_t m_local;
208 size_t n_local;
209 if(i < m) {
210 m_local = i;
211 n_local = m;
212 }
213 else {
214 m_local = m;
215 n_local = i;
216 }
217 Y = t->Jjj[n_local*(n_local+1)/2 + m_local] / t->J2[i+1];
218 sum += Y * f_in[i];
219 }
220 f_out[m] = sum * 2.0 * r*r;
221 }
222
223 return GSL_SUCCESS;
224 }
225
226