package integrators;

import main.*;

/** Simple predictor-corrector with simple-minded adaptive stepsizing...
probably unstable and rather inefficient, and shouldn't be used for
any real work, but maybe for igniting a better multistep method... or
maybe with methods that use extrapolation to zero stepsize.  Basically
it makes an estimate of y(t+h) by first making such an estimate with
the forward Euler step y(t+h) = y(t) + hf(y(t)), then using that estimate
to calculate the backward Euler step y(t+h) = y(t) + hf(y(t+h)), then
averaging the two estimates.  Because of the averaging, and because no
iteration is performed on the correction step, this is not really
an implicit method.
 *
 * @deprecated */

/* INTERNAL TEST METHOD.  DO NOT USE*/
public class HuenIntegrator extends Integrator	{
	
	final String		type = "Huen";
	
	NodeValues		[]	yT;
	NodeValues		[]	fT;
	NodeValues		[]	yP;
	NodeValues		[]	fP;
	NodeValues		[]	yC;
	
	boolean				adaptStepSize=true;
	final float			maxcut=10f, maxboost=4f;
	float				minTimeStep=0.00001f;

	boolean				gotCurrentInfo=false, finalsSet=false;	// used in recursion
	boolean				warned = false;

	/** no-argument constructor allows instantiating a HuenIntegrator
	generically from its name alone, with default minimum step size of
	10E-5 and with step size adaptation turned ON. */
	public HuenIntegrator() {}
	
	/** Other than allocating the necessary storage arrays, HuenIntegrator.init()
	also sets the error tolerance to 10E-4 and minimum step size to 10E-5.  To override
	these defaults, follow init() with a call to Integrator.setErrorTolerance(float),
	then call HuenIntegrator.setMinStepSize(float).  The minimum step size should
	probably be about one tenth of the error tolerance. */
	public void init(Model model)	{
		super.init(model);
		
		errorTolerance = 0.0001f;
		minTimeStep = 0.00001f;
		
		yT = new NodeValues[num_cells * num_nodes];
		fT = new NodeValues[num_cells * num_nodes];
		yP = new NodeValues[num_cells * num_nodes];
		fP = new NodeValues[num_cells * num_nodes];
		yC = new NodeValues[num_cells * num_nodes];
		
		for(int i = 0; i < num_cells; i++) {
			for(int j = 0; j < num_nodes; j++) {
				int k = (i * num_nodes) + j;
				yT[k] = new NodeValues();
				yT[k].numSides = model.cells[0].nodes[j].numSides;
				fT[k] = new NodeValues();
				fT[k].numSides = model.cells[0].nodes[j].numSides;
				yP[k] = new NodeValues();
				yP[k].numSides = model.cells[0].nodes[j].numSides;
				fP[k] = new NodeValues();
				fP[k].numSides = model.cells[0].nodes[j].numSides;
				yC[k] = new NodeValues();
				yC[k].numSides = model.cells[0].nodes[j].numSides;
			}
		}
	}
	
	/** Resets timestep to default, and cleans up various internal flags;
	call before starting a new run.  Does NOT affect the error estimate
	(which is overwritten each step anyway), the default time step size,
	the minimum step size, or whether the method adapts the step size.*/
	public void reset()	{
		timestep = defaultStepsize;
		gotCurrentInfo = finalsSet = warned = false;
	}
	
	public float IntegrateOneStep()	{
		
		int i, j, k, l;
		float yPi, yCi, errest, maxerr = 0.0f, stepTaken = timestep;
		finalsSet = false;
		
		// Euler Step... crude prediction of y at t+timestep
		if(!gotCurrentInfo)	{
			getCurrentValues(yT);
			getDerivatives(fT);
			gotCurrentInfo = true;
		}
		for(i = 0; i < num_cells; i++) {
			for(j = 0; j < num_nodes; j++) {
				k = (i*num_nodes) + j;
				for(l = 0; l < yP[k].numSides; l++)	{
					yP[k].values[l] = yT[k].values[l] + (timestep * fT[k].values[l]);
				}
			}
		}
		setIntegrationValues(yP);
		
		// Backward Euler Step, averaged with the forward one...
		getDerivatives(fP);
		for(i = 0; i < num_cells; i++) {
			for(j = 0; j < num_nodes; j++) {
				k = (i*num_nodes) + j;
				for(l = 0; l < yP[k].numSides; l++)	{
					yC[k].values[l] = 0.5f * (yP[k].values[l] + (yT[k].values[l] + (timestep * fP[k].values[l])));
					yCi = yC[k].values[l]; yPi = yP[k].values[l];
					errest = Math.abs(yCi-yPi);
					if(errest > maxerr) maxerr = errest;
				}
			}
		}
		setIntegrationValues(yC);
		
		eEst = maxerr;		// crude estimate of the local error per timestep; since eEst is set in here,
							// and maxerr is a local variable, subsequent recursions may improve eEst.

		if(adaptStepSize)	{
			// if the error is too great and we haven't hit the mimimum stepsize,
			// start over with a smaller timestep...
			maxerr /= errorTolerance;
			if(maxerr > 1.0f)	{
				if(timestep > minTimeStep)	{	// NOTE: if the timestep is exactly minTimeStep, leave it alone and don't recurse again
					timestep = Math.max((timestep / (1.0f + maxerr)),(timestep / maxcut));
					if(timestep < minTimeStep)	{
						timestep = minTimeStep;
					}
					stepTaken = this.IntegrateOneStep();
				}
				else if(timestep == minTimeStep)	{
					if(!warned) {	// only issue the warning once, and don't recurse again because we've already done the calculation with the minimum step size
						System.out.println("Warning from HuenIntegrator: using minimum timestep!");
						warned = true;
					}
				}
			}
			// ...but if there is plenty of headroom, cautiously increase the timestep
			else if(maxerr < 0.25f)	{
				timestep /= ((1/maxboost) + maxerr);
				warned = false;	// must've gotten off the floor somehow, but warn next time we hit it
			}
		}
		
		if(!finalsSet)	{
			setFinalValues(yC);
			finalsSet = true;
		}
		gotCurrentInfo = false;
		return stepTaken;

	}
	
	/** Turn on or off the recursion this simple-minded method uses to try and
	make the error estimate acceptable.  Turn it off if this method is being used
	as a subsidiary for another method, or if exact step sizes are required. */
	public void setAdaptiveStepsizing(boolean b)	{	adaptStepSize = b;	}
	
	/** Set the minimum timestep.  If the method has to use this minimum size,
	it will belch a warning and just let the errors accumulate, rather than recurse
	further with a smaller step.  The warning is issued to System.out only once each
	time the floor is hit (I hope), so multiple warnings mean the Integrator is
	bouncing off the floor set by the minimum time step.*/
	public void setMinStepSize(float ts)	{	minTimeStep = ts;	}
	
	/** HuenIntegrator overrides Integrator.setStepSize(float) to enforce the minimum
	time step.  Use setMinStepSize(float) first if you need to. */
	public void setStepSize(float ts)	{
		if(ts > minTimeStep) timestep = ts;
		else timestep = minTimeStep;
	}
	
}