/**
 * 
 */
package pacemaker;

/**
 * Useful <code>static</code> code related to
 * <b>Phase Response Curve</b> (PRC) related to 
 * {@link Firefly#PAR_MODEL "MAS" model} (Mirollo and Strogatz 1990).<p> 
 * 
 * PRC is a function in the x,y-plane (in the square (0,0)-(1,1)) 
 * that shows the relation between:<ol>
 * 		<li>the phase (the moment during the period)
 * 		<li>the amount of eccitation accumulated (eccitation increases during the period,
 * 		when it reach the threshold, the firefly emits a flash) 
 * </ol>
 * <p>
 * (threshold is at value 1,1)
 * <p>
 * The function will be the classic<p> 
 * <code>f(x)=1/<b>b</b>*log(1+x*(exp(<b>b</b>)-1))</code> (Mirollo and Strogatz 1990, Fig. 1)<p>
 * that is a monotonic increasing, concave-down function moving from (0,0) to (1,1).<p>
 *
 * The parameter named '<b>b</b>' influences the shape of the PRC in the plane. 
 * (see Mirollo and Strogatz 1990, Fig. 3).<p>
 * <b>b</b> must be in the range ]0,+inf] for the definition of f(x).
 * When <b>b</b> approaches 0, f(x) resembles the identity function (y=x). 
 * But when <b>b</b> exactly =0, the function shape degenerates cause of mathematical 
 * infinites, although reasonably she should fit exacly the identity. 
 * We overcome it in this implementation: when you supply b=0 we will
 * substitute the above f(x) with identity f(x)=x. 
 * @author Tony Binci
 */
public class MaSModel {
	
	/**
	 * Evaluates the PRC for the given phase (x).<p>
	 * <ol>
	 * <li>f(x)=x if <b>b</b>=0.<p>
	 * <li><code>f(x)=1/<b>b</b>*log(1+x*(exp(<b>b</b>)-1))</code> if <b>b</b>=0.<p>
	 * </ol>  
	 * @param phase - the phase (x value) proportional to period. In range [0,1]. 
	 * If longer than the actual period, is beyond the
	 * f(x) definition. In this case, always 1 is returned (and in implementation, 
	 * this will probably trigger the execution of a flash immediately) 
	 * @return PRC evaluation (y value): in range [0,1].
	 */
	public static double eval_PRC(double phase, int b) {
		//double x = phase/period;
		if ( phase < 0.0d ) return 0.0d;
		if ( phase > 1.0d ) return 1.0d;
		if ( b == 0 ) return phase;
		return 1.0d/b * Math.log( 1.0d + phase*( Math.exp(b)-1 ) );
	}
	
	/**
	 * Evaluate the inverse function <code>f<sup>-1</sup></code> of the PRC:
	 * given the y value, you get the x value (the phase).<p>
	 * eval_prc(double, int)
	 * @param y - y value, in range [0,1]
	 * @return x value - the phase. In range [0,1].
	 */
	public static double eval_inverse_PRC(double y, int b) {
		if ( y < 0.0d ) return 0.0d;
		if ( y > 1.0d ) return 1.0d;
		if ( b == 0 ) return y;
		return ( Math.exp( b*y ) - 1.0d)/( Math.exp(b) - 1.0d);
	}

	/** For testing purposes */
	public static void main(String[] args) {
		System.out.println("Testing PRC...");
		for (int b=1; b<4; b++)
			for (double x=0.0d; x<2.0d; x+=0.1d)
				System.out.println("b="+b+" x="+x+" PRC="+eval_PRC(x, b));
		
		System.out.println("x=0.991 bug:"+eval_PRC(0.991, 4));
	}
}
