samer@4: /*=================================================================
samer@4:  * hzeta.c 
samer@4:  *
samer@4:  * Wrapper function for gsl_sf_hzeta_e
samer@4:  *
samer@4:  * This MEX function computes the Hurwitz-zeta-function by calling 
samer@4:  * the function gsl_sf_hzeta_e from the GNU Scientific Library. For 
samer@4:  * x>1 and a>=0 the Hurwitz-zeta-function is defined as
samer@4:  * 
samer@4:  *                        oo 
samer@4:  *                       -----
samer@4:  *                        \         1
samer@4:  *       hzeta(x, a) =     )     -------- .
samer@4:  *                        /             x
samer@4:  *                       -----   (i + a)
samer@4:  *                       i = 0
samer@4:  *
samer@4:  * The MEX function takes two arguments, x and a. If both arguments 
samer@4:  * are matrices of the same size, then the MEX function will compute
samer@4:  * the Hurwitz-zeta-function elementwise and the result of the MEX 
samer@4:  * function will be a matrix of the same size as the inputs. If one 
samer@4:  * input argument is a scalar s, then the function behaves as if this
samer@4:  * argument would be a matrix filled with the value s of the same 
samer@4:  * size as the other argument.
samer@4:  *
samer@4:  * Compile this MEX function with "mex hzeta.c -lgsl". See Matlab
samer@4:  * documentation for additional information on compiling MEX functions.
samer@4:  *
samer@4:  * (C) by Heiko Bauke, 2007, heiko.bauke@physics.ox.ac.uk
samer@4:  *
samer@4:  * This program is free software; you can redistribute it and/or
samer@4:  * modify it under the terms of the GNU General Public License in
samer@4:  * version 2 as published by the Free Software Foundation.
samer@4:  *
samer@4:  * This program is distributed in the hope that it will be useful, but
samer@4:  * WITHOUT ANY WARRANTY; without even the implied warranty of
samer@4:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
samer@4:  * General Public License for more details.
samer@4:  *
samer@4:  * You should have received a copy of the GNU General Public License
samer@4:  * along with this program; if not, write to the Free Software
samer@4:  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
samer@4:  * 02111-1307, USA.
samer@4:  *
samer@4:  *=================================================================*/
samer@4: 
samer@4: #include "mex.h"
samer@4: #include <gsl/gsl_errno.h>
samer@4: #include <gsl/gsl_math.h>
samer@4: #include <gsl/gsl_sf_zeta.h>
samer@4: 
samer@4: void mexFunction(int nlhs, mxArray *plhs[], 
samer@4:                  int nrhs, const mxArray *prhs[]) {
samer@4:   double *x, *y, *a;
samer@4:   int nrows_x, ncols_x, i_x, 
samer@4:       nrows_a, ncols_a, i_a, 
samer@4:       nrows_y, ncols_y, i_y, 
samer@4:       status;
samer@4:   gsl_sf_result result;
samer@4: 
samer@4:   /* check input arguments */
samer@4:   if (nrhs!=2)
samer@4:     mexErrMsgTxt("Two input arguments required.\n");
samer@4:   if (!mxIsDouble(prhs[0]))
samer@4:     mexErrMsgTxt("1st argument not a matrix of double."); 
samer@4:   nrows_x=mxGetM(prhs[0]);
samer@4:   ncols_x=mxGetN(prhs[0]);
samer@4:   if (!mxIsDouble(prhs[1]))
samer@4:     mexErrMsgTxt("2nd argument not a matrix of double."); 
samer@4:   nrows_a=mxGetM(prhs[1]);
samer@4:   ncols_a=mxGetN(prhs[1]);
samer@4:   if (!( (nrows_a==nrows_x && ncols_a==ncols_x) ||
samer@4: 	 (nrows_x==1 && ncols_x==1) ||
samer@4: 	 (nrows_a==1 && ncols_a==1) ) ) 
samer@4:     mexErrMsgTxt("arguments must be matrices of the same size or scalars."); 
samer@4:   /* create matrix for the return argument */ 
samer@4:   nrows_y=nrows_x>nrows_a ? nrows_x : nrows_a;
samer@4:   ncols_y=ncols_x>ncols_a ? ncols_x : ncols_a;
samer@4:   plhs[0]=mxCreateDoubleMatrix(nrows_y, ncols_y, mxREAL); 
samer@4:   /* assign pointers to each input and output */ 
samer@4:   x=mxGetPr(prhs[0]); 
samer@4:   a=mxGetPr(prhs[1]);
samer@4:   y=mxGetPr(plhs[0]); 
samer@4:   /* compute zeta function */
samer@4:   gsl_set_error_handler_off();
samer@4:   for (i_x=i_a=i_y=0; i_y<nrows_y*ncols_y; ++i_y) {
samer@4:     if (gsl_sf_hzeta_e(x[i_x], a[i_a], &result)==GSL_SUCCESS)
samer@4:       y[i_y]=result.val;
samer@4:     else
samer@4:       y[i_y]=GSL_NAN;
samer@4:     if (!(nrows_x==1 && ncols_x==1))
samer@4:       ++i_x;
samer@4:     if (!(nrows_a==1 && ncols_a==1))
samer@4:       ++i_a;
samer@4:   }
samer@4: }