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

static const double Tref = 273.15;
static const double Tc = 513.0;
using namespace std;

double f
(
    const double *coeffs,
    const double& t,
    const int& deg
)
{
    double s = 0.0;
    double p = t/Tc - 1.0;
    for(int i=0; i<=deg; i++)
    {
      s += coeffs[i]*pow(p, i);
    }

    return s;
}

double df
(
 const double& t,
 const int& i
)
{
    double p = t/Tc -1.0;

    return pow(p, i);
}

double s
(
    const double *t,
    const double *coeffs,
    const double *y,
    const int& deg,
    const int& N
)
{
    double sum = 0.0;
    for(int i=0; i<N; i++)
    {
        sum += pow(f(coeffs, t[i], deg) - y[i], 2.0);
    }

    return sum;
}

double ds
(
    const double *t,
    const double *coeffs,
    const double *y,
    const int& j,
    const int& deg,
    const int& N
)
{
    double sum = 0.0;
    for(int i=0; i<N; i++)
    {
      sum += (f(coeffs, t[i], deg) - y[i])*df(t[i], j);
    }

    return sum;
}

int main()
{
    ifstream dataFile;
    dataFile.open("cp.dat");

    int N = 130;
    int deg = 6;
    double coeffs[deg+1];
    int i=1;
    double T[N], cp1[N], cp2[N];
    double err = 1.0e+15;
    double tol = 1.0;
    int Nmax = 1000000;
    /*
    // cp1 coeffs
    coeffs[0] = -53683.29;
    coeffs[1] = 102907.404;
    coeffs[2] = 22864.594;
    coeffs[3] = 7491.1549;
    coeffs[4] = 1536.495;
    */

    // cp2 coeffs
    /*
    coeffs[0] = -40475.5;
    coeffs[1] = 102108.04;
    coeffs[2] = -30209.237;
    coeffs[3] = -5900.68;
    coeffs[4] = -138.03;
    coeffs[5] = 329.51;
    coeffs[6] = 68.86;
    */
    coeffs[0] = -40469.5;
    coeffs[1] = -102091;
    coeffs[2] = -30316.1;
    coeffs[3] = 5986;
    coeffs[4] = 125;
    coeffs[5] = -561;
    coeffs[6] = 126.7;

    //while(!dataFile.eof())
    for(int i=0; i<N; i++)
    {
        dataFile >> T[i] >> cp1[i] >> cp2[i];
	
	cout << "T = " << T[i]
	     << ", cp1 = " << cp1[i]
	     << ", cp2 = " << cp2[i]
	     << endl;
	//t[i] = Tc/T[i];
	//i++;
    }
    dataFile.close();
    int n=0;
    int k=0;
    double dsda[deg];
    while((err > tol) && (k < Nmax))
    {
        double sc1 = 0.0;
	for(int i=0; i<=deg; i++)
	{
	  sc1 += coeffs[i]*coeffs[i];
	}
	sc1 = sqrt(sc1);

        err = s(T, coeffs, cp2, deg, N);
	for(int i=0; i<=deg; i++)
	{
	    dsda[i] = ds(T, coeffs, cp2, i, deg, N);
	    coeffs[i] -= 1.0e-4*dsda[i];
	}

        double sc2 = 0.0;
	for(int i=0; i<=deg; i++)
	{
	  sc2 += dsda[i]*dsda[i];
	}
	sc2 = sqrt(sc2);

	for(int i=0; i<=deg; i++)
	{
	  //coeffs[i] -= 1.0e-6*(1.0/((sc2+1.0)))*dsda[i];
	}

	if (n>100)
	{
	    cout << "err = " << err << endl;
            for(int i=0; i<=deg; i++)
	    {
	        cout << coeffs[i] << " ";
	    }
	    cout << endl;
            for(int i=0; i<=deg; i++)
	    {
	        cout << dsda[i] << " ";
	    }
	    cout << endl;
	    n = 0;
	}
	k++;
	n++;
    }

    cout.precision(10);
    for(int i=0; i<=deg; i++)
    {
      cout << "coeff[" << i << "] = " << coeffs[i] << endl;
    }
    //cout << "number of lines = " << i << endl;

    ofstream fitFile;
    fitFile.open("cpFit.dat");
    for(int i=0; i<N; i++)
    {
      fitFile << T[i]
	      << " " << f(coeffs, T[i], deg) - cp2[i]
	      << " " << f(coeffs, T[i], deg)
	      << " " << cp2[i]
	      << endl;
    }
    fitFile.close();
    return 0;
}