/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ParseGrid parse the grid command lines */
#include "cddefines.h"
#include "optimize.h"
#include "grid.h"
#include "parser.h"

/* ParseGrid - called from ParseCommands if GRID command found */
void ParseGrid(
	/* command line, which was changed to all caps in main parsing routine */
	Parser &p)
{
	DEBUG_ENTRY( "ParseGrid()" );

	/* RP fake optimizer to run a grid of calculations, also accepts
	 * keyword XSPEC */
	strcpy( optimize.chOptRtn, "XSPE" );
	grid.lgGrid = true;

	if( p.nMatch("REPE") )
	{
		/* just keep repeating, don't actually change the values in the grid.
		 * useful for debugging unintentional crosstalk */
		grid.lgStrictRepeat = true;
	}

	/* 06 aug 22, change to accept three parameters: lower and upper limit and number of points. */
	/* scan off range for the previously selected variable */
	if( optimize.nparm > 0 )
	{
		ASSERT( optimize.nparm <= LIMPAR );

		grid.paramLimits[optimize.nparm-1][0] = (realnum)p.FFmtRead();
		grid.paramLimits[optimize.nparm-1][1] = (realnum)p.FFmtRead();
		grid.paramIncrements[optimize.nparm-1] = (realnum)p.FFmtRead();
		grid.lgLinearIncrements[optimize.nparm-1] = p.nMatch("LINE") ? true : false ;
		if( grid.paramIncrements[optimize.nparm-1] < realnum(0.) )
			grid.lgNegativeIncrements = true;

		/* the increase step should not be 0 */
		if( grid.paramIncrements[optimize.nparm-1] == 0. )
		{
			fprintf( ioQQQ," The increment (third parameter) should not be zero.\n" );
			fprintf( ioQQQ," Sorry.\n" );
			cdEXIT( EXIT_FAILURE );
		}

		if( abs( grid.paramIncrements[optimize.nparm-1] ) > abs( grid.paramLimits[optimize.nparm-1][1] - grid.paramLimits[optimize.nparm-1][0] ) )
		{
			fprintf( ioQQQ," The increment (third parameter) must not be greater than the difference between the limits (first and second parameters).\n" );
			fprintf( ioQQQ," Sorry.\n" );
			cdEXIT( EXIT_FAILURE );
		}

		if( p.lgEOL() )
		{
			fprintf( ioQQQ," This command has changed since the definition given in Porter et al. 2006, PASP, 118, 920.\n" );
			fprintf( ioQQQ," The grid command now requires three parameters: lower limit, upper limit, and increment.\n" );
			fprintf( ioQQQ," The keywords RANGE and STEPS are no longer necessary.\n" );
			fprintf( ioQQQ," Sorry.\n" );
			cdEXIT( EXIT_FAILURE );
		}
		else
		{
			++optimize.nRangeSet;
		}

		realnum ratio =	(grid.paramLimits[optimize.nparm-1][1] - grid.paramLimits[optimize.nparm-1][0])/
		      grid.paramIncrements[optimize.nparm-1];

		/* Alert if the uplimit and lowlimit are wrong */
		if( ratio < realnum(0.) )
		{
			fprintf( ioQQQ, "The increment (third parameter) has the wrong sign. It doesn't take you from the initial to the final grid value (first and second parameter, resp.).\n" );
			fprintf( ioQQQ," Sorry.\n" );
			cdEXIT( EXIT_FAILURE );
		}

		// this takes care of the blowup in the error due to cancellation in limits[1]-limits[0]
		// it assumes that limits[1]-limits[0] is accurate within 3*eps*(limits[0]+limits[1])/2
		// which should be a very conservative estimate...
		realnum feps = realnum(1.5)*
			(grid.paramLimits[optimize.nparm-1][1] + grid.paramLimits[optimize.nparm-1][0])/
			(grid.paramLimits[optimize.nparm-1][1] - grid.paramLimits[optimize.nparm-1][0]);
		long eps = max(nint(abs(feps)),3);

		// this will blow for pathologically narrow grid ranges
		ASSERT( eps <= INT_MAX );

		// take special care if step is integer fraction of max-min (which is nearly always the case)
		if( fp_equal( ratio, realnum(nint(ratio)), int(eps) ) )
			grid.numParamValues[optimize.nparm-1] = nint(ratio) + 1;
		else
			grid.numParamValues[optimize.nparm-1] = long(ratio) + 1;

		if( grid.numParamValues[optimize.nparm-1] < 2 )
			fprintf( ioQQQ, " NOTE must have at least two grid points\n" );

		grid.numParamValues[optimize.nparm-1] = MAX2( 2, grid.numParamValues[optimize.nparm-1] );

		// Create some buffer area in the allowed range of parameter values to prevent
		// accidentally going over the limit due to roundoff error. The buffer is 1/10th
		// of a step, so should still guard against doing too many steps due to bugs.
		realnum safety = 0.001f*grid.paramIncrements[optimize.nparm-1];

		if( grid.lgLinearIncrements[optimize.nparm-1] )
		{
			if( grid.paramLimits[optimize.nparm-1][0]-safety<=0. )
			{
				fprintf(ioQQQ,"The current implementation of the grid command works with log parameter values even when you specify LINEAR.\n");
				fprintf(ioQQQ,"A non-positive value was entered.  The grid command cannot deal with this.\n");
				cdEXIT( EXIT_FAILURE );
			}
			optimize.varang[optimize.nparm-1][0] = log10(grid.paramLimits[optimize.nparm-1][0]-safety);
			optimize.varang[optimize.nparm-1][1] = log10(grid.paramLimits[optimize.nparm-1][1]+safety);
		}
		else
		{
			optimize.varang[optimize.nparm-1][0] = grid.paramLimits[optimize.nparm-1][0]-safety;
			optimize.varang[optimize.nparm-1][1] = grid.paramLimits[optimize.nparm-1][1]+safety;
		}
	}

	return;
}