package integrators;

import main.*;

/** Modified Midpoint Method as described in Numerical Recipes in C;
this is NOT a translation of their code (which I didn't even look at),
but merely my implementation of the recipe outlined in the text.  This
should be used as a subsidiary of an Bulirsch-Stoer-type Integrator, but
can also be used as an integrator in it own right, although since I have
a hard time imagining why anyone would do so, the error control implemented
here is rudimentary and inefficient.  The goal of this approach is to cross
a large interval H with a series of n substeps h.  The first substep
is a forward Euler step.  Subsequent steps are midpoint estimates with
y(t+h) = y(t-h) + 2hf(y(t)).  The final step averages the first estimate
of y(t+nh) with a backward Euler step from y(t+(n-1)h).  Calling
MidpointIntegrator.IntegrateOneStep(int n) causes the Integrator to cross
the interval set by its timestep using n substeps.  This is the call to
make if the MidpointIntegrater is being used as a susidiary method.  The
generic IntegrateOneStep() method, on the other hand, allows one to use
this method in its own right according to a formula in NRC (see comment
in MidpointIntegrator.IntegrateOneStep()). */

/*  MAY NOT BE FUNCTIONAL */
public class MidpointIntegrator extends Integrator	{
	
	final String		type = "Midpoint";
	
	private int			nDiv = 8;			// number of substeps to take
	private int			n = 0;				// step counter
	private float		substep;			// size of substep, = timestep / nDiv
	
	NodeValues		[]	yT=null;			// initial state
	NodeValues		[]	estY=null;			// final estimate
	NodeValues		[]	Fn=null;			// storage for derivative at current point
	NodeValues		[]	yPrev,yCurr,yNext;	// this method keeps track of only three steps at a time
	
	/** no-argument constructor allows instantiating a MidpointIntegrator
	generically from its name alone, with default number of substeps = 8. */
	public MidpointIntegrator()	{}
	
	/** init() allocates the necessary storage, nothing else. */
	public void init(Model model)	{
		super.init(model);
		
		timestep = defaultStepsize = 0.1f;
		
		yT = new NodeValues[num_cells * num_nodes];
		estY = new NodeValues[num_cells * num_nodes];
		Fn = new NodeValues[num_cells * num_nodes];
		yPrev = new NodeValues[num_cells * num_nodes];
		yCurr = new NodeValues[num_cells * num_nodes];
		yNext = 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;
				estY[k] = new NodeValues();
				estY[k].numSides = model.cells[0].nodes[j].numSides;
				Fn[k] = new NodeValues();
				Fn[k].numSides = model.cells[0].nodes[j].numSides;
				yPrev[k] = new NodeValues();
				yPrev[k].numSides = model.cells[0].nodes[j].numSides;
				yCurr[k] = new NodeValues();
				yCurr[k].numSides = model.cells[0].nodes[j].numSides;
				yNext[k] = new NodeValues();
				yNext[k].numSides = model.cells[0].nodes[j].numSides;
			}
		}
	}
	
	/** Call IntegrateOneStep() to use the Modified Midpoint Method by itself;
	this call uses an estimate give in NRC based on a weighted sum of results
	from crossing the interval with n/2 substeps and then with n substeps...
	use MidpointIntegrator.setNumberOfSubSteps(int) to set n; yEst = (4/3)yEst(n)
	- (1/3)yEst(n/2).  The error is the maximum one-dimensional difference
	between the coarse estimate and finer one with all n substeps.  Probably
	not the most useful estimate in the world, but I haven't thought of a
	better one, and I probably won't because I don't intend to use this method
	on its own.  The rudimentary error control here seems very inefficient. */
	public float IntegrateOneStep()	{
		float maxerr = 0.0f;
		
		// cross the interval with n/2 substeps...
		float temp = CrossInterval(nDiv/2);
		// ...then get the coarse part of the weighted sum...
		for(int i = 0; i < num_cells; i++)	{
			for(int j = 0; j < num_nodes; j++)	{
				int k = (i * num_nodes) + j;
				for(int l = 0; l < yCurr[k].numSides; l++)	{
					estY[k].values[l] = yNext[k].values[l];
				}
			}
		}
		// ...reset model to y(t) to take a finer-resolution estimate...
		setIntegrationValues(yT);
		setFinalValues(yT);
		temp = CrossInterval(nDiv);
		// ...then get the fine part of the weighted sum...
		for(int i = 0; i < num_cells; i++)	{
			for(int j = 0; j < num_nodes; j++)	{
				int k = (i * num_nodes) + j;
				for(int l = 0; l < yCurr[k].numSides; l++)	{
					maxerr = Math.max(Math.abs(estY[k].values[l] - yNext[k].values[l]),maxerr);
					estY[k].values[l] = ((4.0f * yNext[k].values[l]) - estY[k].values[l]) / 3.0f;
				}
			}
		}
		
		eEst = maxerr;
		maxerr /= errorTolerance;			// pretty lame, but what else to do?
		if(maxerr > 1.0f)	{				// ...if the error is intolerable,...
			timestep = Math.max(timestep/10,timestep/(1.0f + maxerr));	// ...reduce the timestep...
			setIntegrationValues(yT);		// ...reset everything nicely...
			setFinalValues(yT);
			temp = IntegrateOneStep();		// ...and try again.
		}
		else if(maxerr < 0.25f)	{			// ...but if we've got room,...
			timestep /= 0.25f + maxerr;		// ...boldly raise it for the next round.
		}
		
		// ...finally, set the model to the weighted sum rather than the fine estimate.
		setIntegrationValues(estY);
		setFinalValues(estY);
		
		return temp;
	}
	
	/** Call IntegrateOneStep(int) with the desired number of substeps as the
	argument when using the Modified Midpoint Method as part of another method,
	like the BulirschStoerIntegrator.  Keep in mind that the result is actually
	copied into the Model's Nodes, so copy the current state before calling this
	method to get an estimate of the next state, then restore it from the copy
	before trying to get another estimate!  Also note that this call does NOT
	make any estimate of the error per timestep, so calls to getErrorEstimate()
	will return a meaningless value */
	public float IntegrateOneStep(int num_substeps)	{
		return CrossInterval(num_substeps);
	}
	
	private float CrossInterval(int num_substeps)	{
		n = 0;
		substep = timestep / num_substeps;
		TakeInitialStep();
		setIntegrationValues(yCurr);
		while(n < num_substeps - 1)	{
			TakeInternalStep();
			NodeValues [] temp = yPrev;		// rotate the three storage arrays
			yPrev = yCurr;
			yCurr = yNext;
			yNext = temp;
			setIntegrationValues(yCurr);
		}
		TakeFinalStep();
		setIntegrationValues(yNext);
		setFinalValues(yNext);
		return timestep;
	}
	
	private void TakeInitialStep()	{
		n++;
		getCurrentValues(yT);
		copyNVarray(yT, yPrev);
		getDerivatives(Fn);
		for(int i = 0; i < num_cells; i++)	{
			for(int j = 0; j < num_nodes; j++)	{
				int k = (i * num_nodes) + j;
				for(int l = 0; l < yCurr[k].numSides; l++)	{
					yCurr[k].values[l] = yPrev[k].values[l] + substep * Fn[k].values[l];
				}
			}
		}
	}
	
	private void TakeInternalStep()	{
		n++;
		getDerivatives(Fn);
		for(int i = 0; i < num_cells; i++)	{
			for(int j = 0; j < num_nodes; j++)	{
				int k = (i * num_nodes) + j;
				for(int l = 0; l < yCurr[k].numSides; l++)	{
					yNext[k].values[l] = yPrev[k].values[l] + 2 * substep * Fn[k].values[l];
				}
			}
		}
	}
	
	private void TakeFinalStep()	{
		n++;
		getDerivatives(Fn);
		for(int i = 0; i < num_cells; i++)	{
			for(int j = 0; j < num_nodes; j++)	{
				int k = (i * num_nodes) + j;
				for(int l = 0; l < yCurr[k].numSides; l++)	{
					yNext[k].values[l] = 0.5f * (yCurr[k].values[l] + yPrev[k].values[l] + substep * Fn[k].values[l]);
				}
			}
		}
	}
	
	/** If one is going to use this object as an integrator in its own right, rather
	than as a subsidiary of a more sophisticate one, then one has the option of
	setting how many substeps it uses to cross an interval.  The default is 8. */
	public void setNumberOfSubsteps(int n)	{	nDiv = n;	}
	
}