/* 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 "rfield.h"
#include "stars.h"

static float RebinSingleCell(float,float,float[],float[],float[],long);
static long RebinFind(float[],long,float);

/*====================================================================== */
/*RebinAtmosphere: generic routine for rebinning atmospheres onto Cloudy grid */
void RebinAtmosphere(
  /* the number of points in the incident continuum*/
  long nCont , 
  /* the energy grid for the stellar continuum in Ryd*/
  /*@in@*/ float StarEner[] , /* StarEner[nCont] */
  /* the predicted continuum flux on the stellar continuum scale */
  /*@in@*/ float StarFlux[] , /* StarFlux[nCont] */
  /* what we want, the continuum on the cloudy grid */
  /*@out@*/ float CloudyFlux[] , /* CloudyFlux[NC_ELL] */
  /* the number of edges in AbsorbEdge */
  long nEdge ,
  /* the list of absorption edges that require special attention */
  /*@in@*/ float AbsorbEdge[] ) /* AbsorbEdge[nEdge] */
{
	int lgDone;
	long int ind,
	  j,
	  k;
	/* >>chng 00 jun 02, demoted next two to float, PvH */
	float BinHigh, 
	  BinLow,
	  BinMid,
	  BinNext,
	  *EdgeLow=NULL,
	  *EdgeHigh=NULL,
	  *StarPower;

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

	if( nEdge > 0 )
	{
		if( (EdgeLow = (float *)MALLOC( sizeof(float)*(unsigned)nEdge )) == NULL )
		{
			fprintf( ioQQQ, " not enough memory to allocate EdgeLow in RebinAtmosphere\n" );
			puts( "[Stop in RebinAtmosphere]" );
			cdEXIT(EXIT_FAILURE);
		}
		if( (EdgeHigh = (float *)MALLOC( sizeof(float)*(unsigned)nEdge )) == NULL )
		{
			fprintf( ioQQQ, " not enough memory to allocate EdgeHigh in RebinAtmosphere\n" );
			puts( "[Stop in RebinAtmosphere]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* this loop should be before the next loop, otherwise models with a
	 * very strong He II edge (i.e. no flux beyond that edge) will fail */
	for( j=0; j < nEdge; j++ )
	{
		ind = RebinFind(StarEner,nCont,AbsorbEdge[j]);

		/* sanity check */
		ASSERT( ind+1 < nCont );

		EdgeLow[j] = StarEner[ind];
		EdgeHigh[j] = StarEner[ind+1];
	}

	/* cut off that part of the Wien tail that evaluated to zero */
	for( j=nCont-1; j >= 0; j-- )
	{
		if( StarFlux[j] == 0.f )
		{
			nCont = j;
		}
	}

	ASSERT( nCont>0);
	if( (StarPower = (float *)MALLOC( sizeof(float)*(unsigned)(nCont-1) )) == NULL )
	{
		fprintf( ioQQQ, " not enough memory to allocate StarPower in RebinAtmosphere\n" );
		puts( "[Stop in RebinAtmosphere]" );
		cdEXIT(EXIT_FAILURE);
	}

	for( j=0; j < nCont-1; j++ )
	{
		StarPower[j] = (float)(log(StarFlux[j+1]/StarFlux[j])/log(StarEner[j+1]/StarEner[j]));
			  
	}

	for( j=0; j < rfield.nupper; j++ )
	{
		/* BinHigh is upper bound of this continuum cell */
		BinHigh = rfield.anu[j] + 0.5f*rfield.widflx[j];

		/* BinLow is lower bound of this continuum cell */
		BinLow = rfield.anu[j] - 0.5f*rfield.widflx[j];

		/* BinNext is upper bound of next continuum cell */
		BinNext = ( j+1 < rfield.nupper ) ? rfield.anu[j+1] + 0.5f*rfield.widflx[j+1] : BinHigh;

		lgDone = FALSE;

		/* >>chng 00 aug 14, take special care not to interpolate over major edges,
		 * the region inbetween EdgeLow and EdgeHigh should be avoided,
		 * the spectrum is extremely steep there, leading to significant roundoff error, PvH */
		for( k=0; k < nEdge; k++ )
		{
			if( BinLow < EdgeLow[k] && BinNext > EdgeHigh[k] )
			{
				BinMid = 0.99999f*EdgeLow[k];
				CloudyFlux[j] = RebinSingleCell(BinLow,BinMid,StarEner,StarFlux,StarPower,nCont);
				j++;

				/* sanity check */
				ASSERT( j < rfield.nupper );

				BinMid = 1.00001f*EdgeHigh[k];
				CloudyFlux[j] = RebinSingleCell(BinMid,BinNext,StarEner,StarFlux,StarPower,nCont);
				lgDone = TRUE;
				break;
			}
		}

		/* default case when we are not close to an edge */
		if( !lgDone )
		{
			CloudyFlux[j] = RebinSingleCell(BinLow,BinHigh,StarEner,StarFlux,StarPower,nCont);
		}
	}

	free(StarPower);
	if( nEdge > 0 )
	{
		free(EdgeHigh);
		free(EdgeLow);
	}

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

static float RebinSingleCell(float BinLow,
			float BinHigh,
			/*@in@*/ float StarEner[],  /* StarEner[nCont] */
			/*@in@*/ float StarFlux[],  /* StarFlux[nCont] */
			/*@in@*/ float StarPower[], /* StarPower[nCont-1] */
			long nCont)
{
	long int i, 
	  ipHi, 
	  ipLo;
	double anu,
	  retval,
	  widflx;
	double power, 
	  slope2, 
	  sum, 
	  v1, 
	  val, 
	  x1, 
	  x2;

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

	anu = (BinLow+BinHigh)/2.;
	widflx = BinHigh-BinLow;

	/* this is the Rayleigh-Jeans slope in F_nu, for low frequency extrapolation */
	slope2 = 2.0e00;

	if( BinLow < StarEner[0] )
	{
		/* this is case where cloudy's continuum is below stellar continuum,
		 * (at least for part of the cell), so we do Rayleigh Jeans extrapolation */
		retval = (float)(StarFlux[0]*(pow(anu/StarEner[0],slope2)));
	}
	else if( BinLow > StarEner[nCont-1] )
	{
		/* case where cloudy continuum is entirely above highest stellar point */
		retval = 0.0e00;
	}
	else
	{
		/* now go through stellar continuum to find bins corresponding to
		 * this cloudy cell, stellar continuum defined through nCont cells */
		ipLo = RebinFind(StarEner,nCont,BinLow);
		ipHi = RebinFind(StarEner,nCont,BinHigh);
		/* sanity check */
		ASSERT( ipLo >= 0 && ipLo < nCont-1 && ipHi >= ipLo );

		if( ipHi == ipLo )
		{
			/* Do the case where the cloudy cell and its edges are between
			 * two adjacent stellar model points: do power-law interpolation  */

			power = StarPower[ipLo];

			val = power*log(anu/StarEner[ipLo]);

			retval = (float)(StarFlux[ipLo]*exp(val));
		}
		else
		{
			/* Do the case where the cloudy cell and its edges span two or more
			 * stellar model cells:  add segments with power-law interpolation up to
			 * do the averaging.*/

			sum = 0.0e00;

			/* ipHi points to stellar point at high end of cloudy continuum cell,
			 * if the Cloudy cell extends beyond the stellar grid, ipHi == nCont-1
			 * and the MIN2(ipHi,nCont-2) prevents access beyond allocated memory
			 * ipLo points to low end, above we asserted that 0 <= ipLo < nCont-1 */
			for( i=ipLo; i <= MIN2(ipHi,nCont-2); i++ )
			{
				power = StarPower[i];

				if( i == ipLo )
				{
					x1 = BinLow;
					x2 = StarEner[i+1];
					v1 = StarFlux[i]*exp(power*log(x1/StarEner[i]));
					/*v2 = StarFlux[i+1];*/
				}

				else if( i == ipHi )
				{
					x2 = BinHigh;
					x1 = StarEner[i];
					/*v2 = StarFlux[i]*exp(power*log(x2/StarEner[i]));*/
					v1 = StarFlux[i];
				}

				/*if( i > ipLo && i < ipHi )*/
				else
				{
					x1 = StarEner[i];
					x2 = StarEner[i+1];
					v1 = StarFlux[i];
					/*v2 = StarFlux[i+1];*/
				}

				if( fabs(power+1.0) < 0.001 )
				{
					val = x1*v1*log(x2/x1);
				}
				else
				{
					val = pow(x2/x1,power + 1.e00) - 1.e00;
					val = val*x1*v1/(power + 1.e00);
				}
				sum += val;
			}

			retval = sum/widflx;
		}
	}

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

#define SIGN3(X) ((X) < 0. ? -1L : ((X) == 0. ? 0L : 1L))

static long RebinFind(/*@in@*/ float array[], /* array[nArr] */
		long nArr,
		float val)
{
	long i1,
	  i2,
	  i3,
	  ind = -2,
	  sgn;

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

	/* sanity check */
	ASSERT( nArr > 1 );

	/* return ind(val) such that array[ind] <= val <= array[ind+1],
	 *
	 * NB NB: this routine assumes that array[] increases monotonically !
	 *
	 * the first two clauses indicate out-of-bounds conditions and
	 * guarantee that when val1 <= val2, also ind(val1) <= ind(val2) */

	if( val < array[0] )
	{
		ind = -1;
	}
	else if( val > array[nArr-1] )
	{
		ind = nArr-1;
	}
	else
	{
		/* do a binary search for ind */
		i1 = 0;
		i3 = nArr-1;
		while( i3-i1 > 1 )
		{
			i2 = (i1+i3)/2;
			sgn = SIGN3(val-array[i2]);

			switch(sgn)
			{
			case -1:
				i3 = i2;
				break;
			case 0:
				ind = i2;

#				ifdef DEBUG_FUN
				fputs( " <->RebinFind()\n", debug_fp );
#				endif
				return( ind );
			case 1:
				i1 = i2;
				break;
			}
		}
		ind = i1;
	}

	/* sanity check */
	ASSERT( ind > -2 );

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

	return( ind );
}
