/* This file is part of Cloudy and is copyright (C) 1978-2006 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
#include "cddefines.h"
/*
#include "conv.h"
#include "input.h"
#include "called.h"
#include "version.h"
#include "prt.h"
#include "trace.h"
#include "grains.h"
#include "parse.h"
*/
#include "punch.h"
#include "optimize.h"
#include "cddrive.h"
#include "physconst.h"
#include "zero.h"
#include "rfield.h"
#include "cloudy.h"
#include "lines.h"
#include "grid.h"
/* #include "opacity.h" */

static int option = 6;
static int FirstRun = TRUE;

void gridFunc(float param[] );

void gridXspec(float xc[], float del[], long int nInterpVars, long nParVals)
{
	long int i, j;

#	ifdef DEBUG_FUN
	fputs( "<+>grid_xspec()\n", debug_fp );
#	endif
	
	if( nInterpVars > LIMPAR )
	{
		fprintf( ioQQQ, "grid_do: too many parameters are varied, increase LIMPAR\n" );
		puts( "[Stop]" );
		cdEXIT(EXIT_FAILURE);
	}

	optimize.nOptimiz = 0;
	
	grid.nintparm = nInterpVars;
	grid.naddparm = 0;
	grid.totNumModels = (long)pow((double)nParVals, (double)nInterpVars);
	ASSERT( grid.totNumModels > 1 );

	if( (grid.paramNames = (char**)MALLOC(sizeof(char*)*(unsigned)(nInterpVars+grid.naddparm) ) )==NULL )
		BadMalloc();

	if( (grid.paramMethods = (long*)MALLOC(sizeof(long)*(unsigned)(nInterpVars+grid.naddparm) ) )==NULL )
		BadMalloc();

	if( (grid.paramRange = (float**)MALLOC(sizeof(float*)*(unsigned)(nInterpVars+grid.naddparm) ) )==NULL )
		BadMalloc();

	if( (grid.paramData = (float**)MALLOC(sizeof(float*)*(unsigned)(nInterpVars+grid.naddparm) ) )==NULL )
		BadMalloc();

	if( (grid.interpParameters = (float**)MALLOC(sizeof(float*)*(unsigned)(grid.totNumModels) ) )==NULL )
		BadMalloc();

	for( i=0; i<nInterpVars+grid.naddparm; i++ )
	{
		if( (grid.paramNames[i] = (char*)MALLOC(sizeof(char)*(unsigned)(12) ) )==NULL )
			BadMalloc();

		if( (grid.paramRange[i] = (float*)MALLOC(sizeof(float)*(unsigned)(6) ) )==NULL )
			BadMalloc();

		if( (grid.paramData[i] = (float*)MALLOC(sizeof(float)*(unsigned)(nParVals) ) )==NULL )
			BadMalloc();

		sprintf( grid.paramNames[i],	"%s%ld", "PARAM", i+1 );
		/* Method is 0 for linear, 1 for logarithmic, for now all linear. */
		grid.paramMethods[i] = 0;
		/* Initial */
		grid.paramRange[i][0] = xc[i];
		/* Delta */
		grid.paramRange[i][1] = del[i];
		/* Minimum */
		grid.paramRange[i][2] = xc[i]*0.5f;
		/* Bottom */
		grid.paramRange[i][3] = xc[i]*0.5f;
		/* Top */
		grid.paramRange[i][4] = (xc[i]+del[i]*(nParVals-1))*2.f;
		/* Maximum */
		grid.paramRange[i][5] = (xc[i]+del[i]*(nParVals-1))*2.f;
		
		for( j=0; j<nParVals; j++ )
		{
			grid.paramData[i][j] = xc[i]+del[i]*j;
		}
	}

	for( i=0; i<grid.totNumModels; i++ )
	{
		if( (grid.interpParameters[i] = (float*)MALLOC(sizeof(float)*(unsigned)(nInterpVars) ) )==NULL )
			BadMalloc();
	}

	for( i=0; i< grid.totNumModels; i++ )
	{
		float variableVector[LIMPAR];
		
		for( j=0; j<nInterpVars; j++ )
		{
			int index;
			
			index = (int)( (i/(int)pow((double)nParVals,nInterpVars-1.-j))%(nParVals) );

			variableVector[j] = xc[j] + del[j]*index;

			grid.interpParameters[i][j] = variableVector[j];
		}
		for( j=nInterpVars; j<LIMPAR; j++ )
		{
			variableVector[j] = xc[j];
		}

		if( i == grid.totNumModels - 1 )
			punch.lgOpenUnits = TRUE;

		gridFunc(variableVector);

		ASSERT( FirstRun == FALSE );
	}
	

	/* Now change these to keV */
	for( i=0; i<grid.numEnergies; i++ )
	{
		grid.Energies[i] *= 0.001f*(float)EVRYD;
	}

#	ifdef DEBUG_FUN
	fputs( " <->grid_xspec()\n", debug_fp );
#	endif
	return;
}

void gridFunc(float param[] )
{
	int lgBAD, lgLimOK;

	long int i;

	double snorm;

#	ifdef DEBUG_FUN
	fputs( "<+>grid_func()\n", debug_fp );
#	endif
	/* limpar defined in optimize.h */

	/* This routine is called by AMOEBA (or POWELL) with values of the
	 * variable parameters for CLOUDY in the array p(i). It returns
	 * the value FUNC = SUM (obs-model)**2/sig**2 for the lines
	 * specified in the observational data file, values held in the
	 * common blocks /OBSLIN/ & /OBSINT/
	 * replacement input strings for CLOUDY READR held in /chCardSav/
	 * parameter information for setting chCardSav held in /parmv/
	 * additional variables
	 * Gary's variables
	 */

	/* initialize the code for this run */
	/* cdInit();*/
	/* zero out lots of variables */
	zero();

	if( optimize.lgOptimFlow )
	{
		fprintf( ioQQQ, " trace, grid_func variables" );
		for( i=0; i < optimize.nvary; i++ )
		{
			fprintf( ioQQQ, "%10.2e", param[i] );
		}
		fprintf( ioQQQ, "\n" );
	}

	for( i=0; i < optimize.nvary; i++ )
	{
		optimize.vparm[0][i] = param[i];
	}

	/* allways increment nOptimiz, even if parameters are out of bounds,
	 * this prevents grid to get stuck in infinite loop */
	++optimize.nOptimiz;

	/* call routine to pack /kardsv/ variable with appropriate
	 * CLOUDY input lines given the array of variable parameters p(i) */
	vary_input(&lgLimOK);

	if( !lgLimOK )
	{
		/* these parameters are not within limits of parameter search
		 * >>chng 96 apr 26, as per Peter van Hoof comment */
		fprintf( ioQQQ, " Iteration %ld not within range.\n", 
		  optimize.nOptimiz );

		cdEXIT(EXIT_FAILURE);
	}

	for( i=0; i < optimize.nvary; i++ )
	{
		optimize.varmax[i] = (float)MAX2(optimize.varmax[i],optimize.vpused[i]);
		optimize.varmin[i] = (float)MIN2(optimize.varmin[i],optimize.vpused[i]);
	}

	lgBAD = cloudy();

	/* Do this stuff once.  Malloc some arrays, and get the energies. */
	if( FirstRun )
	{
		grid.numEnergies = rfield.nupper-2;

		if( (grid.Energies = (float*)MALLOC(sizeof(float)*(unsigned)(grid.numEnergies) ) )==NULL )
			BadMalloc();

		if( (grid.Spectra = (float**)MALLOC(sizeof(float*)*(unsigned)(grid.totNumModels) ) )==NULL )
			BadMalloc();

		if( (grid.ExpMinusTau = (float**)MALLOC(sizeof(float*)*(unsigned)(grid.totNumModels) ) )==NULL )
			BadMalloc();
			
		for( i=0; i<grid.totNumModels; i++ )
		{
			if( (grid.Spectra[i] = (float*)MALLOC(sizeof(float)*(unsigned)(grid.numEnergies) ) )==NULL )
				BadMalloc();

			if( (grid.ExpMinusTau[i] = (float*)MALLOC(sizeof(float)*(unsigned)(grid.numEnergies) ) )==NULL )
				BadMalloc();
		}

		for( i=0; i<grid.numEnergies; i++ )
		{
			grid.Energies[i] = rfield.AnuOrg[i];
		}

		FirstRun = FALSE;
	}

	/* Grab spectrum for xspec printout */
	cdSPEC2( option, grid.numEnergies, grid.Spectra[optimize.nOptimiz-1]);
	/* e^-tau for this run, at each energy bin */
	cdSPEC2( 8, grid.numEnergies, grid.ExpMinusTau[optimize.nOptimiz-1]);
	
	if( lgBAD )
	{
		fprintf( ioQQQ, " Cloudy returned error condition - what happened?\n" );
	}

	if( LineSave.ipNormWavL < 0 )
	{
		fprintf( ioQQQ, 
			" Normalization line array index is bad.  What has gone wrong?\n" );
		puts( "[Stop in grid_func]" );
		cdEXIT(EXIT_FAILURE);
	}
	snorm = LineSv[LineSave.ipNormWavL].sumlin;

	if( snorm == 0. )
	{
		fprintf( ioQQQ, " Normalization line has zero intensity.  What has gone wrong?\n" );
		fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
		puts( "[Stop in grid_func]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* first print all warnings */
	cdWarnings(ioQQQ);

#	ifdef DEBUG_FUN
	fputs( " <->grid_func()\n", debug_fp );
#	endif
	return;
}

void GridGather( void )
{
#	ifdef DEBUG_FUN
	fputs( " <+>GridGather()\n", debug_fp );
#	endif

	if( optimize.nOptimiz==grid.totNumModels && grid.lgGrid )
	{
		ASSERT( optimize.nOptimiz >= 1 );
		/* Grab spectrum for xspec printout */
		cdSPEC2( option, grid.numEnergies, grid.Spectra[optimize.nOptimiz-1]);
		/* e^-tau for this run, at each energy bin */
		cdSPEC2( 8, grid.numEnergies, grid.ExpMinusTau[optimize.nOptimiz-1]);
	}

#	ifdef DEBUG_FUN
	fputs( " <->GridGather()\n", debug_fp );
#	endif

	return;
}
