htswanalysis/MACS/lib/gsl/gsl-1.11/dht/dht.c

   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 <config.h>
  23 #include <stdlib.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_sf_bessel.h>
  27 #include <gsl/gsl_dht.h>
  28
  29
  30 gsl_dht *
  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
  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 *
  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
 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
 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
 174 double gsl_dht_k_sample(const gsl_dht * t, int n)
 175 {
 176   return t->j[n+1] / t->xmax;
 177 }
 178
 179
 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
 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