#include <fstream.h>
#include <math.h>

void main(void)
	{
   char output_file[50];
   char input_file[50];
   int n, outpnts,full_ave, n_total;

   double *a,*b,*c,*d,*x,*y,*dydx,*d2ydx2,*xnu,*ynu, *all_x, *all_y;
   void create_input(int n, double *x, double *y);
   void create_output(int n, int n_out, double *x, double *y, double *dydx,
                   double *d2ydx2, double *a, double *b, double *c, double *d,
                   double *xnu, double *ynu, double xo, double xf);
   void read_datafile(char *filename, int n, double *x, double *y);
   void write_results(char *filename, int n, double *x, double *y,
                   double *dydx, double *d2ydx2);
   void spline(double *a, double *b, double *c, double *d, double *x,
                   double *y, int n);
   void f(double x_val, double &y_val, double &dydx_val, double &d2ydx2_val,
                       double *a, double *b, double *c, double *d, double *x);
   void f(double x_val, double &y_val, double &dydx_val, double &d2ydx2_val,
                 double *a, double *b, double *c, double *d, double *x, int n);
   void create_subset(int sub_num, int n_total, double *all_x, double *all_y,
                 int n, double *x, double *y);

   cout << "Enter input filename: ";
   cin >> input_file;
   cout << "Enter output filename: ";
   cin >> output_file;
   cout << "Enter the number of data points to read: ";
   cin >> n_total;
   cout << "Enter the number of points to output: ";
   cin >> outpnts;
   cout << "Enter number of full set averages: " ;
   cin >> full_ave;

   n = n_total/full_ave; // n is number of points during each "spline" call
   n_total = n*full_ave; // n_total is number of points to read from input

   all_x = new double[n_total];  // all the data read from input
   all_y = new double[n_total];
   a = new double[n];   // Coefficients of interpolation cubics.
   b = new double[n];   // Only n-1 a,b,c,d are really needed.
   c = new double[n];
   d = new double[n];
	x = new double[n];
   y = new double[n];
   dydx = new double[outpnts];   // Variables to output
   d2ydx2 = new double[outpnts];
   xnu = new double[outpnts];
   ynu = new double[outpnts];

   for (int i=0; i<outpnts; i++)  // Initialize to zero so averaging works
   	{
      ynu[i] = dydx[i] = d2ydx2[i] = 0.0;
      }

   read_datafile(input_file,n_total,all_x,all_y); // Reads input data file

   for (int i = 0;i<full_ave;i++)
   	{
	   create_subset(i, n_total, all_x, all_y, n, x, y);
      spline(a,b,c,d,x,y,n);

      // "create_output" sums values in ynu, dydx, d2ydx2
      create_output(n,outpnts,x,y,dydx,d2ydx2,a,b,c,d,
      	xnu,ynu,all_x[0],all_x[n_total-1]);
      }

   for (int i=0; i<outpnts; i++)    // Finish the averaging.
   	{
      ynu[i] /= (double)full_ave;
      dydx[i] /= (double)full_ave;
      d2ydx2[i] /= (double)full_ave;
      }

   write_results(output_file, outpnts, xnu,ynu,dydx,d2ydx2);

   delete []a;
   delete []b;
   delete []c;
   delete []d;
   delete []x;
   delete []y;
   delete []dydx;
   delete []d2ydx2;
   delete []xnu;
   delete []ynu;
   delete []all_x;
   delete []all_y;

   // Windows closes windows and error messages may never be read, so I pause.
   cout << "Program completely done.  Enter number to exit." << endl;
   int wait; cin >> wait;

   }

 /*
  Takes n number of (x,y) data points and returns the coefficients of the
  cubic polynomials fit in each zone (n-1 zones) after completing a cubic
  spline interpolation. Hence the inpolation function in one zone is
  y = a + bx + cx^2 + dx^3.
 */
void spline(double *a, double *b, double *c, double *d, double *x, double *y, int n)
	{
   double set_konst(double *x, double *y, int index_of_d2y_to_return);

   double *aa,*bb;
   aa = new double[n]; // n-1 needed.
   bb = new double[n];

   // aa and bb are constants so that the following is true...
   // d2ym[i] = aa[i] dym[i] + bb[i]

   aa[0] = 0.0;
   bb[0] = set_konst(x,y,0);  // This sets d2ydx2 at left-most boundary
   									// to a value computed by fitting a cubic
                              // through the first 4 points of the data.

   for (int i=1; i<n-1; i++)  // loop to calculate aa and bb constants.
   	{
      double xm,xp,ym,yp;
      xm = x[i-1];
      xp = x[i];
      ym = y[i-1];
      yp = y[i];

      aa[i] = 3.0/(-xm + xp) + aa[i-1]/(4.0 - aa[i-1]*xm + aa[i-1]*xp);
      bb[i] = (2.0* ( bb[i-1]*(xm - xp)*(xm - xp) + 3.0*(-2.0 + aa[i-1]*
              (xm - xp))*(ym - yp)))
              /( (xm - xp)*(xm - xp)*(-4.0 + aa[i-1]*xm - aa[i-1]*xp) );
      }

   double known_dpyp;  // Originally set by fitting a cubic through last 4
   						  // data points of the set.
	known_dpyp = set_konst(&x[n-4],&y[n-4], 3);

   for (int i=n-2; i>=0; i--)  // loop to compute coefficients
   	{
      double xm,xp,ym,yp,AA,BB;
      xm = x[i];
      xp = x[i+1];
      ym = y[i];
      yp = y[i+1];
      AA = aa[i];
      BB = bb[i];

      a[i] = (-2.0*BB*xm*(xm - 2.0 * xp)*(xm - xp)* xp -
        4.0 * xp*xp*(3.0 + AA*xp)*ym + 4.0*AA* pow(xm,3)*yp +
        12.0*xm*xp*(ym + AA*xp*ym + yp) -
        6.0*xm*xm * (2.0*yp + AA* xp*(ym + yp)) +
        xm*(xm - xp)* xp*(xm*(-4.0 + AA*(xm - xp)) + 2.0*xp)* known_dpyp)
        / (4.0*(xm - xp)*(xm - xp)*(-3.0 + AA* xm - AA* xp));

      b[i] = (2.0*(BB*(xm*xm*xm - 3.0*xm*xm*xp + 2.0*xp*xp*xp) +
              3.0*(xm*(-2.0+ AA*(xm - 2.0* xp)) + 2.0* xp)*
                (ym - yp)) +(-xm*xm*xm * (-4.0 + AA*xm) +
              3.0* xm*(-2.0 + AA* xm)*xp*xp -
              2.0*(-1 + AA* xm)*xp*xp*xp)* known_dpyp)
              / (4.0*(xm - xp)*(xm - xp)*(-3.0 + AA* xm - AA* xp));

      c[i] = (6.0*xp*(BB*(xm - xp) + AA*(ym - yp)) +
        (2.0* AA* xm*xm*xm + 6.0* xm* xp + AA*xp*xp*xp - 3.0*xm*xm*
        (2.0 + AA*xp))* known_dpyp)
        /(4.0*(xm - xp)*(xm - xp)*(-3.0 + AA* xm - AA* xp));

      d[i] = -(2.0 *BB*xm-2.0*BB*xp + 2.0*AA*ym-2.0*AA*yp+
        (-2.0+AA*(xm-xp))* (xm-xp)* known_dpyp)
        /(4.0*(xm - xp)*(xm - xp)*(-3.0 + AA* xm - AA* xp));

      known_dpyp = -(6.0*BB*xm - 6.0*BB*xp + 6.0*AA*ym - 6.0*AA*yp +
          AA*(xm - xp)*(xm - xp)*known_dpyp)
          /(2.0*(xm - xp)*(-3.0 + AA*xm - AA*xp));

      }

   delete []aa;
   delete []bb;
   }

   /*
    Opens the text file named "filename" and reads n number of x,y data points
    where the input is tab delimited. If the end of the file is reached it
    will not read the full amount of data and display an error message.
   */
void read_datafile(char *filename, int n, double *x, double *y)
	{
   fstream input_file;
   double xpoint,ypoint;
   input_file.open(filename, ios::in);

   int i;
   for(i=0; ( i<n ) && ( input_file >> xpoint)  ; i++)
   	{

      input_file >> ypoint;
      x[i] = xpoint;
      y[i] = ypoint;
      }
   if (i != n) cerr << "Full amount of data was not read." << endl;

   input_file.close();

	}

   /*
    Opens text file named "filename" and outputs n datapoints in the format:
    x tab y tab dydx tab d2ydx2 endline.
   */
void write_results(char *filename, int n, double *x, double *y, double *dydx,
                   double *d2ydx2)
	{
   fstream output_file;
   output_file.open(filename, ios::out);
   if(!output_file) {cerr << "Error opening output file."; return;}
   for(int i=0; i<n; i++)
   	{
      output_file << x[i] << '\t' << y[i] << '\t' << dydx[i] << '\t'
                  << d2ydx2[i] << endl;
      }
   output_file.close();

	}

	/*
    n = number of points sampled, number zones is n-1.
    n_out = number of points to sample
    x and y are original data
    a,b,c,d are cubic coefficients as returned by spline function
    xnu, ynu, dydx, d2ydx2 are computed and returned. Note the values are
    	added to ynu, dydx, d2ydx2 so averaging may occur.
    The xnu values are evenly spaced from xo to xf, with points on the
       boundaries.
   */

void create_output(int n, int n_out, double *x, double *y, double *dydx,
                   double *d2ydx2, double *a, double *b, double *c, double *d,
                   double *xnu, double *ynu, double xo, double xf)
	{
   double dx;
   dx = (xf - xo) / (double)(n_out-1);

   	int z; // index for which coefficients to use
	   for (int i=0; i<n_out; i++)
   		{
   	   xnu[i] = xo + (double)i*dx;

         z = 0;
         while(xnu[i] > x[z]) z++;
         if(z>0) z--;
         if(z >= n-1) z = n-2;

      	ynu[i] += a[z] + xnu[i]*( b[z] + xnu[i]*(c[z] + xnu[i]*d[z]));
   	   dydx[i] += b[z] + xnu[i]*(2.0*c[z] + xnu[i]*3.0*d[z]);
	      d2ydx2[i] += 2.0*c[z] + 6.0*d[z]*xnu[i];
      	}

   }

   /*
   Function to create a dataset for debugging purposes.
   n points are sampled.  In this case the sin(x) is computed from 0 to 2pi.
   */

void create_input(int n, double *x, double *y)
	{
   double dx, xo = 0.0, xf = 6.28;
   dx = (xf-xo)/(double)(n-1);
   for (int i=0;i<n;i++)
   	{
      x[i] = xo+dx*(double)i;
      y[i] = sin(x[i]);
      }
   }

   /* Function fits a cubic polynomial to four data points:
      x[0], x[1], x[2], x[3]. It returns the
      second derivative of the polynomial evaluated at x[n]
      where (obviously) 0 <= n <= 3.
   */

double set_konst(double *x, double *y, int n)
	{
   if(n < 0 || n > 3)
   	{
      cerr << "Index of point for calculating d2ydx2 is out of range (0-3).";
      return NULL;
      }
   else
   	{
      double xa,xb,xc,xd,ya,yb,yc,yd;
	   xa = x[(n+3)%4]; xb = x[(n+2)%4]; xc = x[(n+1)%4]; xd = x[n];
   	ya = y[(n+3)%4]; yb = y[(n+2)%4]; yc = y[(n+1)%4]; yd = y[n];

	   return ( (2.0*(-xa *(xa - 2.0*xd) * (xa - xd)*xd*(yb - yc) +
            3.0*xc*xc*xd*(xd*ya + xa*yb - xd*yb - xa*yd) +
            xc*xc*xc*(-xd*ya - xa*yb + xd*yb + xa*yd) +
            3.0*xb*xb*xd*(-xd*ya - xa*yc + xd*yc + xc*(ya - yd) + xa*yd)
            + xc*(-2.0*xd*xd*xd*ya + xa*xa*xa*yb - 3.0*xa*xa*xd*yb +
            2.0*xd*xd*xd*yb - xa*xa*(xa - 3.0*xd)*yd) +
            xb*xb*xb*(xd*ya + xa*yc - xd*yc - xa*yd + xc*(-ya+yd))
            + xb*(2.0*xd*xd*xd * ya - xa*xa*xa * yc + 3.0*xa*xa*xd*yc
            - 2.0 *xd*xd*xd*yc + xc*xc*xc*(ya - yd) + xa*xa*xa*yd
            - 3.0*xa*xa*xd*yd + 3.0*xc*xc*xd*(-ya + yd))))
            /((xa - xb)*(xa - xc)*(xb - xc)*(xa - xd)*(xb - xd)*(xc - xd)));
      }
   }

/* This function is not actually used, but may be useful.  It returns the
   value (y_val) of the interpolated function at x_val as well as the values of
   the derivatives in dydx_val and d2ydx2_val. It is supplied the coefficients
   of the interpolation cubics (a,b,c,d) and the number of datapoints, n.  Note
   this means n-1 zones exist so the highest index for the coefficients is n-2.
*/

void f(double x_val, double &y_val, double &dydx_val, double &d2ydx2_val,
            double *a, double *b, double *c, double *d, double *x, int n)
	{
   int z = 0;
   while(x_val > x[z]) z++;
   if(z>0) z--;
   if(z >= n-1) z = n-2;
   y_val = a[z] + x_val*( b[z] + x_val*(c[z] + x_val*d[z]));
   dydx_val = b[z] + x_val*(2.0*c[z] + x_val*3.0*d[z]);
   d2ydx2_val = 2.0*c[z] + 6.0*d[z]*x_val;
   }

   /*
    This function takes the entire data set which is n points in the arrays
    all_x and all_y. It returns n points in x and y corresponding to every
    (n_total / n_th point.  The sub_num refers to an initial offset.  This
    function is useful for averaging.  For example, one may put every odd
    indexed point in another array and every even indexed point in an array
    with two calls.
   */

void create_subset(int sub_num, int n_total, double *all_x, double *all_y,
                   int n, double *x, double *y)
	{
   int di;
   di = n_total / n;
   for (int i=0;i<n;i++)
   	{
      x[i] = all_x[sub_num+di*i];
      y[i] = all_y[sub_num+di*i];
      }
   }