view general/numerical/scalar/hzeta.c @ 61:eff6bddf82e3 tip

Finally implemented perceptual brightness thing.

/*=================================================================
 * hzeta.c 
 *
 * Wrapper function for gsl_sf_hzeta_e
 *
 * This MEX function computes the Hurwitz-zeta-function by calling 
 * the function gsl_sf_hzeta_e from the GNU Scientific Library. For 
 * x>1 and a>=0 the Hurwitz-zeta-function is defined as
 * 
 *                        oo 
 *                       -----
 *                        \         1
 *       hzeta(x, a) =     )     -------- .
 *                        /             x
 *                       -----   (i + a)
 *                       i = 0
 *
 * The MEX function takes two arguments, x and a. If both arguments 
 * are matrices of the same size, then the MEX function will compute
 * the Hurwitz-zeta-function elementwise and the result of the MEX 
 * function will be a matrix of the same size as the inputs. If one 
 * input argument is a scalar s, then the function behaves as if this
 * argument would be a matrix filled with the value s of the same 
 * size as the other argument.
 *
 * Compile this MEX function with "mex hzeta.c -lgsl". See Matlab
 * documentation for additional information on compiling MEX functions.
 *
 * (C) by Heiko Bauke, 2007, [e-mail]
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License in
 * version 2 as published by the Free Software Foundation.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 * 02111-1307, USA.
 *
 *=================================================================*/

#include "mex.h"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_zeta.h>

void mexFunction(int nlhs, mxArray *plhs[], 
                 int nrhs, const mxArray *prhs[]) {
  double *x, *y, *a;
  int nrows_x, ncols_x, i_x, 
      nrows_a, ncols_a, i_a, 
      nrows_y, ncols_y, i_y, 
      status;
  gsl_sf_result result;

  /* check input arguments */
  if (nrhs!=2)
    mexErrMsgTxt("Two input arguments required.\n");
  if (!mxIsDouble(prhs[0]))
    mexErrMsgTxt("1st argument not a matrix of double."); 
  nrows_x=mxGetM(prhs[0]);
  ncols_x=mxGetN(prhs[0]);
  if (!mxIsDouble(prhs[1]))
    mexErrMsgTxt("2nd argument not a matrix of double."); 
  nrows_a=mxGetM(prhs[1]);
  ncols_a=mxGetN(prhs[1]);
  if (!( (nrows_a==nrows_x && ncols_a==ncols_x) ||
	 (nrows_x==1 && ncols_x==1) ||
	 (nrows_a==1 && ncols_a==1) ) ) 
    mexErrMsgTxt("arguments must be matrices of the same size or scalars."); 
  /* create matrix for the return argument */ 
  nrows_y=nrows_x>nrows_a ? nrows_x : nrows_a;
  ncols_y=ncols_x>ncols_a ? ncols_x : ncols_a;
  plhs[0]=mxCreateDoubleMatrix(nrows_y, ncols_y, mxREAL); 
  /* assign pointers to each input and output */ 
  x=mxGetPr(prhs[0]); 
  a=mxGetPr(prhs[1]);
  y=mxGetPr(plhs[0]); 
  /* compute zeta function */
  gsl_set_error_handler_off();
  for (i_x=i_a=i_y=0; i_y<nrows_y*ncols_y; ++i_y) {
    if (gsl_sf_hzeta_e(x[i_x], a[i_a], &result)==GSL_SUCCESS)
      y[i_y]=result.val;
    else
      y[i_y]=GSL_NAN;
    if (!(nrows_x==1 && ncols_x==1))
      ++i_x;
    if (!(nrows_a==1 && ncols_a==1))
      ++i_a;
  }
}