////////////////////////////////////////////////////////////////////////////
//
//                       Copyright (c) 2006-2008
//                       Future Team Aps 
//                       Denmark
//
//                       All Rights Reserved
//
//   This source file is subject to the terms and conditions of the
//   Future Team Software License Agreement which restricts the manner
//   in which it may be used.
//   Mail: hve@hvks.com
//
/////////////////////////////////////////////////////////////////////////////
//
// Module name     :   interpolation.js
// Module ID Nbr   :   
// Description     :   Javascript for cubic spline & polynomial interpolation
// --------------------------------------------------------------------------
// Change Record   :   
//
// Version	Author/Date		Description of changes
// -------  -------------	----------------------
// 01.01	HVE/06-Dec-18	Initial release
// 01.02 	HVE/2007-Apr-8	Added cubic splines
// 01.03	HVE/2008-Jul-15	Added more verbose information to polynomial curve fitting
// End of Change Record
//
/////////////////////////////////////////////////////////////////////////////

// Function argument checkking 
function check(args,func_name)
 {var actual=args.length; var expected=args.callee.length;
  if(actual!=expected) throw Error("Wrong number of arguments in function: "+func_name+"; Expected: "+expected+"; Actually passed "+actual); }
  
 // Print the Interpolation points
function print_interpolation_points(form,Ptrs)
	{ 
	form.eqField.value+="\nPoints Pair to Interpolate: "+Ptrs.length+"\n";
	for(var i=0;i<Ptrs.length;i++) form.eqField.value+="X["+i+"]="+Ptrs[i][0]+" Y["+i+"]="+Ptrs[i][1]+"\n";
	}

// Print the polynominal to be solved
function print_polynomial(form,coeff)
	{form.eqField.value+="\nInterpolation Polynomial:\n"+coeff.toString()+"\n";}


// Graph the solutions
function graph_solution(jg,Ptrs,degree,y2,form)
	{
	var i, co;
	var max_x, max_y, min_x, min_y, dx, dy;
	var xptrs=new Array;
	var yptrs=new Array;
	var zptrs=new Array;

	max_x=max_y=-Infinity;
	min_x=min_y=+Infinity;
	for(i=0;i<Ptrs.length;i++)
		{ 
		if(max_x<Ptrs[i][0]) max_x=Ptrs[i][0];
		if(min_x>Ptrs[i][0]) min_x=Ptrs[i][0];
		if(max_y<Ptrs[i][1]) max_y=Ptrs[i][1];
		if(min_y>Ptrs[i][1]) min_y=Ptrs[i][1];
		}
	max_x++; min_x--;max_y++; min_y--;
	max=max_x; mix=min_x;
	max_x=Math.max(max_x,max_y); min_x=Math.min(min_x,min_y);
	max_y=max_x; min_y=min_x;		
	max_x=Math.ceil(max_x); min_x=Math.floor(min_x);
	max_y=Math.ceil(max_y); min_y=Math.floor(min_y);
	dx=max_x-min_x; dy=max_y-min_y;
	if(dx==0&&dy==0) {max_x=+1, min_x=-1; max_y=1; min_y=-1; dx=2; dy=2 }
	//alert("("+min_x+","+min_y+");("+max_x+","+max_y+") dx="+dx+" dy="+dy);
	co=new Coordinate(jg,400,400,Number(min_x), Number(min_y), Number(max_x), Number(max_y), Math.ceil(dx/10), Math.ceil(dy/10),false,false);
	var ss="";
	if(form.Cubic.checked==true) ss+="Natural Cubic Spline ";
	if(form.Cubic.checked==true&&form.Polynomial.checked==true) ss+="& ";
	if(form.Polynomial.checked==true) ss+="Polynomial";
	ss+="Interpolation" 
	co.addTitle(ss);
	
	for(i=0;i<Ptrs.length;i++)
		{ 
		jg.setStroke(2);
		co.drawPoint(Number(Ptrs[i][0]),Number(Ptrs[i][1]));
		jg.setStroke(1);
		}
	if(form.Cubic.checked==true) co.addLegend("#ff0000", "Cubic Spline");
	if(form.Polynomial.checked==true) co.addLegend("#00ff00", "Polynomial ("+degree+" deg.)");
	
	for(j=0,i=mix;i<=max;i+=(max-mix)/50)
		{
		if(form.Cubic.checked) {
		var val2=cubic_splint(Ptrs,y2,i);
		if(val2>=min_y&&val2<=max_y) zptrs[j]=val2;
		}
		if(form.Polynomial.checked) {
		var val=value(i,degree);
		if(val>=min_y&&val<=max_y) yptrs[j]=val; 
		}
		xptrs[j]=i; 
//		alert("i="+i+" F(x)="+yptrs[j]);
		j++;
		} 
	
	if(form.Cubic.checked) {
	jg.setColor("#ff0000");
	jg.setStroke(2);
	co.drawPolyline(xptrs,zptrs); 
	jg.setStroke(1);jg.setColor("#000000");}
	if(form.Polynomial.checked) {
	jg.setColor("#00ff00");
	jg.setStroke(2);
	co.drawPolyline(xptrs,yptrs); 
	jg.setStroke(1);jg.setColor("#000000");}
	
	co.paint();
	}

// Print the Interpolation points
function print_cubic_spline_points(form,Ptrs)
	{ 
	form.eqField.value+="\nPoints Pair for Cubic SPlines: "+Ptrs.length+"\n";
	for(var i=0;i<Ptrs.length;i++) form.eqField.value+="X["+i+"]="+Ptrs[i][0]+" Y["+i+"]="+Ptrs[i][1]+"\n";
	}

// find the natrural cubic spline and return the 2nd derrivative of the cubic spline 
function cubic_spline(Ptrs,verbose)
	{var p;
	var u=new Array; var y2=new Array;
	y2[0]=u[0]=0;
	for(var i=1;i<Ptrs.length-1;i++)
		{
		var sig=(Ptrs[i][0]-Ptrs[i-1][0])/(Ptrs[i+1][0]-Ptrs[i-1][0]);
		p=sig*y2[i-1]+2.0;
		y2[i]=(sig-1.0)/p;
		u[i]=(Ptrs[i+1][1]-Ptrs[i][1])/(Ptrs[i+1][0]-Ptrs[i][0])-(Ptrs[i][1]-Ptrs[i-1][1])/(Ptrs[i][0]-Ptrs[i-1][0]);
		u[i]=(6.0*u[i]/(Ptrs[i+1][0]-Ptrs[i-1][0])-sig*u[i-1])/p;
		}
	y2[Ptrs.length-1]=(0.0-0.0*u[Ptrs.length-2])/(0.0*y2[Ptrs.length-2]+1.0);
	for(var k=Ptrs.length-2;k>=0;k--)
		y2[k]=y2[k]*y2[k+1]+u[k];
	return y2;
	}

// Return 	the y value of the input x
function cubic_splint(Ptrs,y2,x)
	{
	var y, h, k;
	var klo=1, khi=Ptrs.length-1;
	while(khi-klo>1)
		{
		k=(khi+klo)>>1;
		if(Ptrs[k][0]>x) khi=k; else klo=k;
		}
	h=Ptrs[khi][0]-Ptrs[klo][0];
	if(h==0) alert("X points must be disctinct");
	a=(Ptrs[khi][0]-x)/h;
	b=(x-Ptrs[klo][0])/h;
	y=a*Ptrs[klo][1]+b*Ptrs[khi][1]+((a*a*a-a)*y2[klo]+(b*b*b-b)*y2[khi])*(h*h)/6.0;
	return y;
	}	
	
	
var jg = new jsGraphics("theCanvas");

// Miscellaneous parsePolynomial() function in it's various format: (float [+-]i float) or (float) or ([+-]i float)
// It also allows input without paranthese e.g. float [+-]i float, float, [+-]i float
function parseInterpolation(eq)
	{
	// Remove all whitespace
	function removeWhiteSpace(eq){var ws=/[\s]+/;while(eq.search(ws)!= -1){eq=eq.replace(ws, "" );}return eq;}
	// Test for leading sign. If omitted add a + sign
	function leadingSign(eq){var lsign=/^[+-]/;if( lsign.test(eq)==false) return "+" + eq;else return eq;}
	// Remove empties for the split function (Firefox only)
	function removeEmpties(arIn)
	{var arOut = new Array;
    var dstIdx=0;
    for(var srcIdx=0;srcIdx<arIn.length;srcIdx++){if(arIn[srcIdx].length!=0){arOut[dstIdx++]=arIn[srcIdx];}	  }
    return arOut;
	}//for(var srcIdx=0;srcIdx<arIn.length;){if(arIn[srcIdx].length==0) arIn.splice(srcIdx,1) else srcIdx++}retyrb arIn;}
	// Convert req to rep
	function converting(eq,reg,rep){if(eq.search(reg)!=-1) eq=eq.replace(reg,rep);return eq;}
	// Convert the equation to an un-ambigous form
// Read and interpret equation to solve
var re=/^[+-]?(([\d]*([\.][\d]+)?)|([\d]+[\.][\d]*))?([eE][+-]?[\d]+)?$/
var rec=/^([+-]\([+-]?(([\d]*([\.][\d]+)?)|([\d]+[\.][\d]*))?([eE][+-]?[\d]+)?[+-]?([iI](([\d]*([\.][\d]+)?)|([\d]+[\.][\d]*))?([eE][+-]?[\d]+)?)?\)([xX]([\^][\d]+)?)?)+$/
var repp=/([+-]{2,})|([+-]$)/

var i, j,k;
var x,y;
var pairs, args;
var Ptrs=new Array;
	
 	eq=removeWhiteSpace(eq);
	//eq=leadingSign(eq);
	pairs=removeEmpties(eq.split(/;/));		// Separate the pairs
//	alert("Pairs="+pairs);
	for(i=0,k=0;i<pairs.length;i++)
		{
		args=removeEmpties(pairs[i].split(/,/));		// Separate the arguments
		//alert("Args="+args);
//		if(args.length!=2) alert("Pair incorrect");
		x=args[0];
		y=args[1];
//		alert("parsing x,y="+x+","+y);
		if(re.test(x)&&re.test(y))
			{
			xp=parseFloat(x);
			yp=parseFloat(y);
			Ptrs[k++]=new Array(xp,yp);
//			alert("Float="+xp+","+yp+" Ptrs="+Ptrs+" k="+k);
//			alert("Pair"+String(k)+": X["+Number(k-1)+"]="+Ptrs[k-1][0]+" Y["+Number(k-1)+"]="+Ptrs[k-1][1]);
			}
		else
			{
			alert("Invalid numbers in pair:"+Number(k+1)+" ("+x+","+y+")");
			}
		}
	return Ptrs;
	}

var alpha=new Array;
var beta=new Array;
var c=new Array;

function value(x,degree)
	{var d2,d1,d0;
	d2=d1=0;
	for(var j=degree;j>=0;j--)
		{
		d0=c[j]+d1*(x-(j>=degree?0:alpha[j]))-d2*(j>=degree-1?0:beta[j+1]);
		d2=d1;
		d1=d0;
		}
//	alert("Value return d0="+d0+" x="+x);		
	return d0;
	}

function findit(Ptrs,degree,verbose)
	{
	var n=Ptrs.length;
	var s, s0;
	var pk1=new Array;
	var pk2=new Array;
	var p=new Array;
	var x,y;

	degree=Math.min(Ptrs.length-1,degree);
	alpha.length=0; beta.length=0; c.length=0;	
	s0=n; alpha[0]=0; beta[0]=0; c[0]=0;
	for(var j=0;j<n;j++)
		{
		c[0]+=Ptrs[j][1];
		alpha[0]+=Ptrs[j][0];
		pk2[j]=0;
		pk1[j]=1;
		}
	c[0]/=s0; alpha[0]/=s0;
	for(var k=1;k<=degree;k++)
		{
		var pk=new Array;
		
		for(j=0,s=0;j<n;j++) 
			{
			x=Ptrs[j][0];
			pk[j]=(x-alpha[k-1])*pk1[j]-beta[k-1]*pk2[j];
			s+=pk[j]*pk[j];
			}
		for(j=0,c[k]=0,alpha[k]=0;j<n;j++)
			{
			x=Ptrs[j][0]; y=Ptrs[j][1];
			c[k]+=y*pk[j];
			alpha[k]+=x*pk[j]*pk[j];
			}
		c[k]/=s;
		alpha[k]/=s;
		beta[k]=s/s0;
		s0=s;pk2=pk1;pk1=pk;
//		alert("iter="+k+":"+pk);	
		}
	
	}
	
// Evaluate and parse the equation. If OK then soilve the equation and print the result
function solveit(form)
	{
	var Ptrs=new Array;
	var cubic, poly, y2; 
	
	/*Ptrs[0]=new Array(-2,1);
	Ptrs[1]=new Array(-1.5,2);
	Ptrs[2]=new Array(-1,2);
	Ptrs[3]=new Array(-0.5,1.5);
	Ptrs[4]=new Array(0,1);	
	Ptrs[5]=new Array(0.5,1.5);
	Ptrs[6]=new Array(1,2);	
	Ptrs[7]=new Array(1.5,3);	
	Ptrs[8]=new Array(2,5);
	*/
	jg.clear();
	eq=form.equation.value;
	verbose=form.Verbose.checked;
	cubic=form.Cubic.checked;
	poly=form.Polynomial.checked;
	degree=form.Degree.options[form.Degree.selectedIndex ].value;
 	Ptrs=parseInterpolation(eq);
	if(Ptrs!=undefined&&Ptrs.length>=2)
		{
		var p;
		//print_interpolation_points(form,Ptrs); 	
		//p=Polynomial.interpolation(Ptrs,degree);
		//print_polynomial(form,p)
		if(cubic==true)
			{y2=cubic_spline(Ptrs,verbose);
			if(verbose)
				{var max_x, Max_y, min_x, min_y;
				var i, j;
				max_x=max_y=-Infinity;
				min_x=min_y=+Infinity;
				for(i=0;i<Ptrs.length;i++)
					{ 
					if(max_x<Ptrs[i][0]) max_x=Ptrs[i][0];
					if(min_x>Ptrs[i][0]) min_x=Ptrs[i][0];
					if(max_y<Ptrs[i][1]) max_y=Ptrs[i][1];
					if(min_y>Ptrs[i][1]) min_y=Ptrs[i][1];
					}
				form.eqField.value+="X		Y\n";
				for(j=0,i=min_x;i<=max_x;i+=(max_x-min_x)/35)
					{
					var val2=cubic_splint(Ptrs,y2,i);
					form.eqField.value+=i+", "+Math.round(val2)+"\n";
					}
				}
			}
		if(poly==true)
			{p=new Array; 
			if(degree>Ptrs.length-1){degree=Ptrs.length-1;if(verbose) form.eqField.value+="Degree reduced to "+Number(Ptrs.length-1)+"\n";}
			findit(Ptrs,degree,verbose);
			if(verbose)	{
			form.eqField.value+="C-length="+c.length+"\n";
			form.eqField.value+="C-coefficients="+c+"\n";
			form.eqField.value+="Alpha-length="+alpha.length+"\n";
			form.eqField.value+="Alpha-coefficients="+alpha+"\n";
			form.eqField.value+="Beta-length="+beta.length+"\n";
			form.eqField.value+="Beta-coefficients="+beta+"\n";
			p[1]=new Polynomial(1); p[0]=new Polynomial(0);
			//alert("Before");
			for(var i=0;i<=degree-2;i++)
				{var tmp, a, b;
				a=new Polynomial(1,-alpha[i]);
				b=new Polynomial(beta[i]);
				p[i+2]=Polynomial.mul(a,p[i+1]);
				tmp=Polynomial.mul(b,p[i]);
				p[i+2]=Polynomial.sub(p[i+2],tmp);
			//	alert("p["+(i+2)+"]="+p[i+2].toString());
				}
			form.eqField.value+="Interpolation Polynomial=";
			for(i=degree-1;i>=0;i--)
				form.eqField.value+=(c[i+1]<0?"-":"+")+Math.abs(c[i+1])+"*P("+i+")";
			form.eqField.value+="\nWhere:\n";
			for(i=0;i<=degree-1;i++)
				form.eqField.value+="P("+i+")="+p[i].toString()+"\n";
			}
			}
		graph_solution(jg,Ptrs,degree,y2,form);
		}
	}
	

function reset_textarea(form)
	{
	form.eqField.value="";
	jg.clear();
	}
function save_textarea(form)
	{
//	form.eqField.value.print();
	}

function print_textarea(form)
	{
	var printing=form.eqField.value;
	newWindow=null;
	newWindow=window.open("","newWin","toolbar=yes,location=yes,menubar,status,toolbar,scrollbars=yes,width=400,height=400,resizable=yes");
    newWindow.document.write('<html>\r\n\r\n');
    newWindow.document.write('<head>\r\n');
    newWindow.document.write('<meta http-equiv="Content-Language" content="en-us">\r\n');
	newWindow.document.write('<meta http-equiv="Content-Type" content="text/html; charset=windows-1252">\r\n');    
    newWindow.document.write('<title>Web based Polynomial root finder result</title>\r\n');
    newWindow.document.write('</head>\r\n\r\n');
    newWindow.document.write('<body>\r\n');
//    newWindow.document.write('<form><input type="button" value="Print This Page" onClick="javascript:window.print();" /></form>');
    newWindow.document.write('<a href="javascript:window.print();">Print This Page</a>');
    newWindow.document.write('\t');
    newWindow.document.write('<a href="javascript:window.close()">\tClose This Page </a>');
    newWindow.document.write('<pre>\r\n');
    newWindow.document.write(printing);
    newWindow.document.write('</pre>\r\n');
    newWindow.document.write('</body>\r\n\r\n');
    newWindow.document.write('</html>');
	}