/*****************************************************************************
 * Copyright 2005 Kevin Hobbs <hobbsk@ohiou.edu>                             *
 *                                                                           *
 * This file is part of Balloon.                                             *
 *                                                                           *
 * Balloon 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 2 of the License, or         *
 * (at your option) any later version.                                       *
 *                                                                           *
 * Balloon 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 Balloon; if not, write to the Free Software                    *
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA *
 *****************************************************************************/


#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "mesh.h"
#include "utility.h"
#include "balloon.h"


int inflate( 
		struct point * p, 
		double * s,
	       	const double * s0,
		const size_t ndims,
		gsl_rng * r)
{
	size_t i;
	double * x;
	double y;
	int big;
	int stepped_out;
	int status;
	
	/*********************************/
	/* The initial point must be in! */
	/*********************************/ 

	x = xmalloc( sizeof( double ) * ndims );
	
	big = 1; /* so we go the first time */
	stepped_out = 0; /* so we go the first time */
	/*
	for ( i = 0; i < ndims; i++)
		if ( fabs(s[i]) > s0[i] / 10)
			big++;
	*/
	
	/*****************************************/
	/* We go till the step size is small and */
	/* we have crossed to the outside at     */
	/* least once.                           */
	/*****************************************/
	while ((big != 0) || (stepped_out == 0))
	{
		for ( i = 0; i < ndims; i++)
			x[i] = p -> x_in[i] + s[i];
		
		status = in_out( x, &y, ndims );
		if ( status == FUNCTION_INSIDE )
		{
			for ( i = 0; i < ndims; i++)
			{
				p -> x_in[i] = x[i];
				s[i] = s[i] * gsl_ran_flat(r, 1.0, 2.0);
			}
			p -> y_in = y;
		}
		
		if ( status == FUNCTION_OUTSIDE )
		{
			stepped_out = 1;
			for ( i = 0; i < ndims; i++)
				s[i] = s[i] / gsl_ran_flat(r, 1.0, 2.0);
		}
		
		big = 0;
		for ( i = 0; i < ndims; i++)
			if ( fabs(s[i]) > s0[i] / 10)
				big++;
	}

	for ( i = 0; i < ndims; i++)
		p -> s[i] = s[i];

	free(x);

	return 0;
}