
<!--Define JavaScript functions.-->

var M = 50; // Global variable, MAX number of known data pairs, (x, y).
var N = 1;  // Global variable, MAX number of points at which interpolating values to be calculated.




xi = new Array(25, 50, 75, 100, 200, 500);
yi = new Array(141, 167, 191, 213, 299, 530);

function csSolve(xValue)
{

 numRows = xi.length;

// var dataFormElements = dataForm.elements; // Reference to the form elements array.

 var textareaIndex = 0;  // The form array index for the textarea element
 var textareaString = "";
 var rawArray = new Array(); // Array for holding raw data entry from textarea box
 var rawArrayLength = 0; // The length of the raw array entered from textarea box

// var numRows = 0; //Indicates the number of rows for which valid data has been entered

 var dVec = new Array(M);  // Array of derivatives at each of the xi points
 var xInt = new Array(N);        // Array of x values at which interpolating y-values will be calculated
 var xk = new Array(N);        // Array of indices for the x values
 var yInt = new Array(N);        // Array of y values, calculated for each of the input xInt values

 var ierr = 0; // Error Code

 var wk = new Array(2); // wk is a 2*M work array, used in the calculation of the derivatives at each of the entered x-values
 wk[0] = new Array(M);
 wk[1] = new Array(M);

 var next = new Array(2);

 var ir = 1, j, jfirst = 0, k, nj;              //Array indices
 var xval_kount = 0; // Counter for number of points at which interpolating values will be calculated
 var dummy, tempx, tempy; //Dummy variables

 var c2, c2t2, c3, c3t3, delta, del1, del2, h, xma, xmi; // Real variables used in the calculation of the interpolant

 // First, confirm that some evaluation points have been entered.
 // This check is relatively quick and easy. And if no valid evaluation points have been entered,
 // it is better to find out before going through the string manipulation/clean-up processing.

 k = 0; // The form index for the "x1" field
  for (var i = 0; i < N; i++) { // Input the points at which an interpolating value will be calculated
 //  tempx = parseFloat(dataFormElements[k].value);
    tempx = parseFloat(xValue);
   if (!isNaN(tempx)){ // The field contains a valid number; otherwise ignore
    xk[i] = xval_kount;
    xInt[xval_kount] = tempx;
    xval_kount++;
   } // End  if (!isNaN(tempx))
   else
    xk[i] = -1; // Keep track of the fact that the field did NOT contain a valid number
   k++;
   k++;
  } // End for i




 // At this point, xi should contain the x-values entered and yi should contain the y-values entered
 // numRows should indicate the number of valid data pairs that were entered

 // Check that pairs are entered in order of increasing x-value
 for (var i = 1; i < numRows; i++){
  if (xi[i] <= xi[i-1]){
   alert("The data pairs have NOT been entered in order of increasing x-value. No further action taken.");
   return;
  }  // End if (xi[i] <= xi[i-1])
 } // End for i

 // Differences between the x-values are now stored in wk[0][. . .], starting with wk[0][1]
 // Divided differences between y-values are stored in wk[1][. . .], starting with wk[1][1]
 // Note that wk[0][0] and wk[1][0] are not filled in this loop; they are presently left unassigned

 j= 1;
 for (var i = 0; i < (numRows - 1); i++){
  wk[0][j] = xi[j] - xi[i];
  wk[1][j] = (yi[j] - yi[i])/wk[0][j];
  j++;
 }//End for i

 if (numRows == 2){
  wk[1][0] = wk[0][0] = 1.0;
  dVec[0] = 2.0*wk[1][1];
 }// End if (numRows == 2)
 else { // else numRows > 2
  tempx = dummy = wk[0][1];
  wk[1][0] = wk[0][2];
  wk[0][0] = tempx + wk[0][2];
  dummy *= dummy*wk[1][2];
  dVec[0] = ((tempx + 2.0*wk[0][0])*wk[1][1]*wk[0][2] + dummy)/wk[0][0];
 } // End else numRows > 2

 nj = numRows - 1;
 for (var i = 1; i < nj; i++){
  tempx = -(wk[0][i+1]/wk[1][i-1]);
  dVec[i] = tempx*dVec[i-1] + 3.0*(wk[0][i]*wk[1][i+1] + wk[0][i+1]*wk[1][i]);
  wk[1][i] = tempx*wk[0][i-1] + 2.0*(wk[0][i] + wk[0][i+1]);
 }//End for i

 if (numRows == 2){
   dVec[1] = wk[1][1];
 } // End if numRows == 2
 else { //else numRows != 2
   if (numRows == 3){
    dVec[2] = 2.0*wk[1][2];
    wk[1][2] = 1.0;
    tempx = -(1.0/wk[1][1]);
   }// End if (numRows == 3)
   else {
    tempx = wk[0][numRows-2] + wk[0][numRows-1];
    dummy = wk[0][numRows-1]*wk[0][numRows-1]*(yi[numRows-2] - yi[numRows-3]);
	dummy /= wk[0][numRows-2];
    dVec[numRows-1] = ((wk[0][numRows-1] + 2.0*tempx)*wk[1][numRows-1]*wk[0][numRows-2] + dummy)/tempx;
    tempx = -(tempx/wk[1][numRows-2]);
    wk[1][numRows-1] = wk[0][numRows-2];
   }//End else

   // Complete forward pass of Gauss Elimination

   wk[1][numRows-1] = tempx*wk[0][numRows-2] + wk[1][numRows-1];
   dVec[numRows-1] = (tempx*dVec[numRows-2] + dVec[numRows-1])/wk[1][numRows-1];

 } // End else numRows != 2

 //Carry out back substitution

 for (var i= numRows-2; i >= 0; i--){
  dVec[i] = (dVec[i] - wk[0][i]*dVec[i+1])/wk[1][i];
 }//End for i

 // End of PCHEZ

 // Begin PCHEV

 // Main loop. Go through and calculate interpolant at each xInt value

 while (jfirst < xval_kount){

  // Locate all points in interval

  for (k = jfirst; k < xval_kount; k++){
   if (xInt[k] >= xi[ir]) break;
  } // End for k loop

  if (k < xval_kount){
   if (ir == (numRows-1))
    k = xval_kount;
  } // End if (k < xval_kount)

  nj = k - jfirst;

  // Skip evaluation if no points in interval

  if (nj > 0){

  // Evaluate Cubic at xInt[k], j = jfirst (1) to k-1
  // =========================================================
  // Begin CHFDV

   xma = h = xi[ir] - xi[ir - 1];

   next[1] = next[0] = 0;
   xmi = 0.0;

   // Compute Cubic Coefficients (expanded about x1)

  delta = (yi[ir] - yi[ir - 1])/h;
  del1 = (dVec[ir - 1] - delta)/h;
  del2 = (dVec[ir] - delta)/h;

  //delta is no longer needed

  c2 = -(del1 + del1 + del2);
  c2t2 = c2 + c2;
  c3 = (del1 + del2)/h;

  // h, del1, and del2 are no longer needed

  c3t3 = c3 + c3 + c3;

  // Evaluation loop

  for (j = 0; j < nj; j++){
   dummy = xInt[jfirst + j] - xi[ir - 1];
   yInt[jfirst + j] = yi[ir - 1] + dummy*(dVec[ir - 1] + dummy*(c2 + dummy*c3));
   if (dummy < xmi) next[0] = next[0] + 1;
   if (dummy > xma) next[1] = next[1] + 1;
   // Note the redundancy: if either condition is true, other is false
  } // End for j loop

  // End CHFDV

  // ========================================================

  if ((next[1] > 0) && (ir != (numRows - 1))) {
    alert("Error Code 5005 for next[1]. No further action taken.");  // This option should never happen.
    return;
  } // End if ((next[1] > 0) && (ir != (numRows - 1)))

  if ((next[0] > 0)  && (ir != 1)) {
    // xInt is not ordered relative to xi, so must adjust evaluation interval
    // First, locate first point to left of xi[ir - 1]
    for (j = jfirst; j < k; j++){
	 if (xInt[j] < xi[ir - 1])
	  break;
	} //End for j
    if (j == k){
     alert("Error Code 5005 for next[0]. No further action taken.");  // This option should never happen.
     return;}
    k = j; //Reset k. This will be the first new jfirst

	// Now find out how far to back up in the xi array

	for (j = 0; j < ir; j++){
	 if (xInt[k] < xi[j])
	  break;
	} // End for j

	// The above loop should NEVER run to completion because xInt[k] < xi[ir - 1]

    // At this point, either xInt[k] < xi[0] or
	//                       xi[j-1] <= xInt[k] < xi[j]
	// Reset ir, recognizing that it will be incremented before cycling

	ir = (((j-1) > 0) ? (j-1) : 0);

  } // End if ((next[0] > 0)  && (ir != 1))

  jfirst = k;

  } // End if (nj > 0)

  ir++;
  if (ir >= numRows) break;

  }// End while (jfirst < xval_kount)

  // End PCHEV



retval="";

 // Now output the interpolating values into the "yInt" field
  k = 2; // The form index for the "y1" field
  for (var i = 0; i < N; i++) {
   if (xk[i] == -1)
   {
    //dataFormElements[k].value = " ";
   }
   else {
   	retval = (yInt[xk[i]]);
    
    if ((xInt[xk[i]] > xi[numRows - 1]) || (xInt[xk[i]] < xi[0]))    ierr++;  // Count the number of points which were extrapolated
   } // End else
   k++;
   k++;
 } // End for i



 return retval;
}  //End of csSolve