/* 
* 
* Copyright (C) 2013 UMPC
* 
* This file is part of CLaRiNet.
*  
* CLaRiNet is free software: you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* CLaRiNet 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 CLaRiNet.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <stddef.h>
#include "Kriging.h"
#include "Erreur.h"
#include <math.h>
#define EPS 3.0e-14
#define PIM4 0.7511255444649425
#define MAXIT 10

/* GaussHermite : Renvoie les abscisses et les poids pour la quadrature de
 * Gauss-Hermite a l'ordre n (source : Numerical Recipes in C) */

int GaussHermite (double ** _x , double ** _w, int n)
{
	int i,its,j,m;
	double p1,p2,p3,pp,z,z1;

        if(! ((*_x) = (double *) malloc ( n * sizeof(double) )) )
        {
            SetErreur("GaussHermite : could not allocate double pointer for abscissa");
            return(-1);
        }

        if(! ((*_w) = (double *) malloc ( n * sizeof(double) )) )
        {
            SetErreur("GaussHermite : could not allocate double pointer for weights");
            return(-1);
        }

	m=(n+1)/2;
        
	for (i=1;i<=m;i++) {
		if (i == 1) {
			z=sqrt((double)(2*n+1))-1.85575*pow((double)(2*n+1),-0.16667);
		} else if (i == 2) {
			z -= 1.14*pow((double)n,0.426)/z;
		} else if (i == 3) {
			z=1.86*z-0.86*(*_x)[0];
		} else if (i == 4) {
			z=1.91*z-0.91*(*_x)[1];
		} else {
			z=2.0*z-(*_x)[i-3];
		}
		for (its=1;its<=MAXIT;its++) {
			p1=PIM4;
			p2=0.0;
			for (j=1;j<=n;j++) {
				p3=p2;
				p2=p1;
				p1=z*sqrt(2.0/j)*p2-sqrt(((double)(j-1))/j)*p3;
			}
			pp=sqrt((double)(2*n))*p2;
			z1=z;
			z=z1-p1/pp;
			if (fabs(z-z1) <= EPS) break;
		}
		if (its > MAXIT)
                {
                    SetErreur("GaussHermite : max number of iterations exceeded");
                    return(-1);
                }
		(*_x)[i-1]=z;
		(*_x)[n-i] = -z;
		(*_w)[i-1]=2.0/(pp*pp);
		(*_w)[n-i]=(*_w)[i-1];
	}

        return(0);

        /* (C) Copr. 1986-92 Numerical Recipes Software ?421.1-9. */
}
