#include "Engine.h"
#include <Math.hpp>
//#include "viewmv.h"
//#include "viewmv.h"
//fix for pow

HINSTANCE LogHandle_Engine;
AnsiString MessageCache_Engine;
setengloghangle(HINSTANCE h) {LogHandle_Engine = h;}

void useotherfit(bfit *fit){fit->fitl(1); }
void errormessage_Engine(AnsiString S){MessageCache_Engine = S; SendMessage(LogHandle_Engine,WM_USER + 301,332,1);}
AnsiString messageis_Engine(){ AnsiString r = MessageCache_Engine; MessageCache_Engine = ""; return r;}

double poww(double x,double n){

        if (n==0) return 1;
        if (x==0) return 0;
        return pow(x,n);

}


int fact(int k) { int ans = 1;
for (int i = 1;i <=k; i++) {ans *= i;}
return ans; }

int cocombin(int n, int k) { return fact(n) / ( fact(k) * fact(n-k) ); }

double logg(double ll){
        if (ll > 0) { return log10(ll);}
        return -1* log10(-1 * ll);

}

int WhereisLastCharPos(const AnsiString &S,const char c) {
        AnsiString T = S;
        int ans=0;
        while (int ii = T.Pos(c)) {  T=T.SubString(ii+1,S.Length()-ans);  ans += ii;};
        return ans+1;
}

AnsiString FromLastCharToEnd(const AnsiString &S,const char c){
        return S.SubString(WhereisLastCharPos(S,c),S.Length());
        }

//Vector Funtions----------------------------------------------------------------
vec& vec::assign(vec b){

newlength(b.no());

for (int k=0;k<b.no();k++)
{
(*this)[k] = b[k];
}

return *this;


}              // */


vec vec::operator*(double b) {
   	vec ans(n);
        for (int c=0;c<n;c++) {
      	        ans[c]=p[c]*b ; }
      return ans;
      }
vec vec::operator+(vec b) {
   	vec ans(n);
   	for (int c=0;c<n;c++) { ans[c]=p[c]+b[c] ; }
      return ans;
      }
vec   vec::operator-(vec b) {
   	vec ans(n);
   	for (int c=0;c<n;c++) { ans[c]=p[c]-b[c] ; }
      return ans;
      }
vec   vec::operator^(double b) {
   	vec ans(n);
        for (int c=0;c<n;c++) { ans[c]=poww (p[c], b) ; }
      return ans;
      }
/*vec& vec::operator= (vec b){
//delete this;
newlength(b.no());
//vec *r = new vec(b.no());
for (int k=0;k<b.no();k++)
{
(*this)[k] = b[k];
}
;
return *this;


}              // */

double vec::maxx() {
      double maxi=p[0];
   	for (int c=0;c<n-1;c++) {
      	if (((p[c]>maxi) && (p[c]>0)) || ((p[c]<0) && (p[c]<maxi))) {maxi=p[c];}
         }
      return maxi;
   }
double vec::abssum(){
        double sumv =0;
        for (int k=0;k < no() ;k++) {
        sumv +=sqrt( poww((*this)[k],2));
        }
        return sumv;
}

vec vec::increment(){
vec ans(no());
for (int k=0;k<no();k++) {
        ans[k] = (1+4e-8) * (*this)[k];
}
return ans;
}

bool operator<(vec &a,vec &b) {   //
     if (a.no()!=b.no()) return 0;
     int t=0,r=0;
     for (int c=0;c<(a.no());c++) { if((!t) & (a[c]<b[c])) {t=0; r=0;} }
     return r; }

bool operator<=(vec &a,vec &b) {   //
     if (a.no()!=b.no()) return 0;
     bool r=true;
     for (int c=0;c<a.no();c++) {
        if (a[c]*a[c] < b[c]*b[c]) {r = true; break;  }
        if (a[c]*a[c] > b[c]*b[c]) {r = false; break; } }
      return r; }


// Matrix functions--------------------------------------------------------------
matr::matr(int r,int c):nrows(r),ncols(c) {p=new vec*[nrows];
		for (int i=0;i<nrows;i++) {p[i]=new vec(ncols);}}


double matr::twobytwo(matr cm) { return 0; }


double matr::cofactorw(int j,int k,matr &m) {
        if ((j >= nrows) || (nrows != ncols)) return 0;
        //int a1 = m[0][0];
        //int a2 = m[0][1];
        matr res(m.nrows-1,m.ncols-1); // = *result;
        if (m.nrows==2) {
                return m[!j][!k];}
        bool ofsetr = 0;
        //int l = poww(-1,k+j+2)  ;
        for (int a=0;a<res.nrows ;a++)
        {
                if (j==a) {ofsetr = 1;} // a++; }else{
                bool ofsetc = 0;
                for (int b=0;b<res.ncols;b++)
                {
                        if (k==b) {ofsetc = 1; }//b++; //}else{
                        res[a][b] = m[a+ofsetr][b+ofsetc];
        }  }

        if (res.nrows == 2) {
                double  answer = ((res[0][0] * res[1][1]) - (res[0][1] * res[1][0]));
                return answer;
        }else{
                return detm(&res);
                }






   }
double matr::detm(matr *xx){
        if (xx == 0) {
                xx = this;}
        double det =0;
        if (xx->nrows == 2) {
                return (*xx)[1][1] * (*xx)[0][0] - (*xx)[1][0] * (*xx)[0][1];  }
        for (int k=0; k < xx->ncols;k++)
        {
                det +=  poww(-1,k+2) * cofactorw(0,k,*xx) * (*xx)[0][k]  ;
        }
   return det;
   }

matr matr::tran() {
        matr answer(ncols,nrows);
                for (int a=0;a < ncols;a++) {
                        for (int b=0;b<nrows;b++) {
                        (answer)[a][b] = (*this)[b][a];
                        }
                }
        return answer;
        }       //  */
matr matr::operator/(double &d) {
         matr answer(nrows,ncols);
        for (int a=0;a < ncols;a++) {
                for (int b=0;b<nrows;b++) {
                        (answer)[a][b] = (*this)[a][b] / d;  }
                        }
     return answer;
                }
matr matr::operator/(float d) {
         matr answer(nrows,ncols);
        for (int a=0;a < ncols;a++) {
                for (int b=0;b<nrows;b++) {
                        (answer)[a][b] = (*this)[a][b] / d;  }
                        }
        return answer;
                }

vec matr::operator* (vec &d) {

          if (d.no() != nrows) {return 0;}
          vec answer(nrows);
          for (int j=0;j < nrows;j++) {
                (answer)[j] = 0;
                for (int l=0;l < nrows;l++)  {
                        (answer)[j] += (*this)[l][j] * d[l];
                }
                //float a1 = answer[j];
          }
     return answer;

     }

vec matr::gaussianelim(vec aa) {  // guassian after sorting
        struct vandm{
        vec *v;
        matr *m;};
        vandm ans;
        ans.v = new vec(aa.no());
        ans.m = new matr(nrows,ncols);
        vec &vv = (*ans.v);
        matr &mm = *ans.m;
        for (int kk = 0;kk < nrows;kk++) {
        for (int j=0;j < ncols;j++) {
        mm[kk][j] = (*this)[kk][j]; }
        vv[kk] = aa[kk]; }

        for (int r =0;r<ncols-1;r++){
                mm.sort(vv);
        for (int q = r+1;q<nrows;q++) {
                double piv = mm[q][r]  / mm[r][r];
                for (int c=r+1;c<ncols;c++)
                {
                mm[q][c] -= piv * mm[r][c];

                }
                mm[q][r] = 0;
                vv[q] -= piv * vv[r];

        }                            }
        vec ss(nrows);
        for (int c=nrows-1;c>=0;c--)
        {
        ss[c] = vv[c] ;
        for (int k=c+1;k<nrows;k++) {
                //double adjust = ss[k] * mm[c][k];
                ss[c] -= ss[k] * mm[c][k];  }
        ss[c]/= mm[c][c] ;
        }
        return ss;
}



/*vec matr::operator% (vec &d) {
        if (d.no() != nrows) {return 0;}
        vec answer(nrows);
        for (int j=0;j < nrows;j++) {
                (answer)[j] = 0;
                for (int l=0;l < nrows;l++)  {
                        (answer)[j] += pow ((*this)[l][j] * d[l], 2) ;
                }
                float a1 = answer[j]; }
     return answer^0.5;
     }       */

   matr matr::inv(){ //matr answer = (
        matr cof(ncols,nrows);
        double dettest = det();

        {
        for (int a=0;a < ncols;a++) {
                for (int b=0;b<nrows;b++) {
                        cof[a][b] =poww(-1,a+b+2) * cofactorw(a,b,*this)  ;
                }}    }
                cof = cof / dettest; return cof;}

void matr::sort(vec &vv){
        while(!sorted()) {
                for (int a=0;a<nrows-1;a++) {
                if ((*this)[a+1]<=(*this)[a]) {continue;
                }else {
                        swap((*this)[a+1],(*this)[a]);
                        swap(vv[a+1],vv[a]);
                }}  }

}

bool matr::sorted(){
        //while((*this)[a]<=(*this)[a+1])
        bool r = true;
        for (int a=0;a<nrows-1;a++) {
                if ((*this)[a+1]<=(*this)[a]) {continue;
                }else {r  = false; break;}}
        return r;

}


matr matr::operator+ (matr mm){
        if ((mm.nrows != nrows) || (mm.ncols != ncols)) {
                throw (622869); }
        matr ans(ncols,nrows);
        for (int c=0; c < ncols;c++) {
                for (int r=0; r < nrows;r++) {
                ans[c][r] = mm[c][r] + (*this)[c][r];
                }}
                return ans;
                }
matr matr::operator- (matr &mm){
        if ((mm.nrows != nrows) || (mm.ncols != ncols)) {throw (622869); }
        matr ans(ncols,nrows);
        for (int c=0; c < ncols;c++) {
                for (int r=0; r < nrows;r++) {
                ans[c][r] = (*this)[c][r] - mm[c][r];
                }}
                return ans;}
matr matr::operator- (double &dd){
        matr ans(ncols,nrows);
        for (int c=0; c < ncols;c++) {
                for (int r=0; r < nrows;r++) {
                ans[c][r] = (*this)[c][r] - dd;
                }}
                return ans;
        }
matr matr::operator+ (double &dd){
        matr ans(ncols,nrows);
        for (int c=0; c < ncols;c++) {
                for (int r=0; r < nrows;r++) {
                ans[c][r] = dd + (*this)[c][r];
                }}
                return ans;

        }


matr jacgen(vec(*fx)(vec *sol,dataset *d),dataset *d,vec *cursol, vec *oldsol=0){  // generates a jacobian matrix via secent method
vec oldsol2(cursol->no());
if (oldsol==0) {
vec oldsol2 = cursol->increment();
oldsol = &oldsol2;}


matr jac(cursol->no(),cursol->no());
vec parm(cursol->no());
jac.setzero();
vec fx1 = fx(cursol,d);
for (int k=0;k<cursol->no();k++)
{
        for (int h=0;h<cursol->no();h++)
        {
        parm.assign(*cursol);
        parm[h] = (*oldsol)[h]; //dh
        vec fx2 = fx(&parm,d); //feval(FUNC,DATA,MODSOL) ;
        //if ((*oldsol)[h]-(*cursol)[h] < 1e-150) { jac[h][k] = 0;
        //}else {
        jac[k][h] = (fx2[k]-fx1[k])/((*oldsol)[h]-(*cursol)[h]);
        //}
        }
}
return jac;
}


//Dataset Functions--------------------------------------------------------------
dataset dataset::copy(){
        dataset ans;
        ans.x = x->xcopy();
        ans.y = y->xcopy();
        ans.erry = erry->xcopy();
        ans.w = w->xcopy();
        return ans;
        }


bool dataset::bettermethod() {
if (onminerrdone) return onminerr;
bool ans = false;
double xwid = dmaxx() - dminx();
int xworder = 0;
int xminorder = 0;
int xmaxorder = 0;
if (xwid) xworder = floor(logg(xwid));
if (dminx()) xminorder = floor(logg(dminx()));
if (dmaxx()) xmaxorder = floor(logg(dmaxx()));
if (xworder) {
        if ((float)xminorder / (float)xworder < 0.75){

                ans = true; onminerrdone = true;
        }
}
double ywid = dmaxy() - dminy();
int yworder = 0;
int yminorder = 0;
int ymaxorder = 0;
if (ywid) yworder = floor(logg(ywid));
if (dminy()) yminorder = floor(logg(dminy()));
if (dmaxy()) ymaxorder = floor(logg(dmaxy()));
if (yworder) {
        if ((float)yminorder / (float)yworder < 0.75){

                ans = true; onminerrdone = true;
        }
}



return ans;
}

double dataset::dmeanx(){ double ans = 0;
        for (int k=0;k< x->no() ; k++ ) { ans += (*x)[k];}
        ans /= x->no(); return ans; }
double dataset::dmeany(){
        double ans = 0;
        for (int k=0;k< y->no() ; k++ ) { ans += (*y)[k];}
        ans /= y->no(); return ans;  }

double dataset::dminx(){
        double ans = (*x)[0];
        for (int k=0;k< x->no() ; k++ ) { if (ans > (*x)[k]) {ans = (*x)[k];}}
        return ans; }
double dataset::dminy(){        double ans = (*y)[0];
        for (int k=0;k< y->no() ; k++ ) { if (ans > (*y)[k]) {ans = (*y)[k];}}
        return ans; }
double dataset::dmaxx(){        double ans = (*x)[0];
        for (int k=0;k< x->no() ; k++ ) { if (ans < (*x)[k]) {ans = (*x)[k];}}
        return ans; }
double dataset::dmaxy(){        double ans = (*y)[0];
        for (int k=0;k< y->no() ; k++ ) { if (ans < (*y)[k]) {ans = (*y)[k];}}
        return ans; }    //*/
double dataset::dmaxywerr(){        double ans = (*y)[0]+(*erry)[0];
        for (int k=0;k< y->no() ; k++ ) { if (ans < (*y)[k] + (*erry)[k]) {ans = (*y)[k] + (*erry)[k];}}
        return ans; }    //*/
double dataset::dminywerr(){        double ans = (*y)[0]-(*erry)[0];
        for (int k=0;k< y->no() ; k++ ) { if (ans > (*y)[k]-(*erry)[k]) {ans = (*y)[k]-(*erry)[k];}}
        return ans; }
void __fastcall dataset::calw() {
        if (!erry) { throw(63035);  }
        if (w) delete w;
        w = new row<double>(erry->no()) ;
        for (int k=0;k< erry->no() ; k++ ) {
        (*w)[k] = poww((*erry)[k],-2) ;        }

}
long double dataset::wxiN(int i,unsigned int no) {
        if ((w->no() != x->no()) || ( w->no() != y->no() ) ) {throw(63044); }
        long double answer;
        //if (no==0) {answer = 1 /  (*erry)[i] ;}else {
        answer = poww((*x)[i] , no) * (*w)[i] / w->no(); //}
        return answer;
}


long double dataset::EwxN(unsigned int no) {
        if ((w->no() != x->no()) || ( w->no() != y->no() ) ) {throw(63044); }
        long double answer = 0;
        for (int k=0;k< w->no() ; k++ ) {

                if (high) {
                double ii;
                if (no==0) {ii = (*w)[k]; }else{ii = (*w)[k] * poww((*x)[k] , no);}
                        if (ii) {answer = ii*(answer / ii + 1); }
                        //if (k == 0) { answer = ii;
                        //}

                }else {
                if (no==0) {answer += (*w)[k];}else{answer +=(*w)[k] * poww((*x)[k], no); }

        }                 }
        return answer;
}
long double dataset::EwyxN(unsigned int no) {
        if ((w->no() != x->no()) || ( w->no() != y->no() ) ) {throw(63044); }
        long double answer = 0;
        for (int k=0;k< w->no() ; k++ ) {

                if (high) {
                        double ii;
                        if (no==0) {ii = (*w)[k] * (*y)[k] ;}else{ii = (*w)[k] * (*y)[k] * poww((*x)[k] , no);}
                        if (ii) {answer = ii*(answer / ii + 1); }
                        //if (k == 0) { answer = ii;}
                }else {
                if (no==0) {answer += (*w)[k]* (*y)[k];}else{answer +=(*w)[k] * (*y)[k] * poww((*x)[k], no); }  }
          }

        return answer;
}
long double dataset::EwerrxN(unsigned int no) {
        if ((w->no() != x->no()) || ( w->no() != erry->no() ) ) {throw(63044); }
        long double answer = 0;
        for (int k=0;k< w->no() ; k++ ) {
                if (high) {double ii = (*w)[k] * (*erry)[k] * poww((*x)[k] , no);
                        if (answer) {answer = ii*(answer / ii + 1); }
                        //if (k == 0) { answer = ii; }
                }else {answer +=(*w)[k] * (*erry)[k] * poww((*x)[k], no); }

        }
        return answer;
}

long double dataset::Ew() {
        if ((w->no() != x->no()) || ( w->no() != y->no() ) ) {throw(63044); }
        long double answer = 0;
        for (int k=0;k< w->no() ; k++ ) {
                if (high) {double ii = (*w)[k] ;
                        if (answer) {answer = ii*(answer / ii + 1); }
                        if (k == 0) { answer = ii; }

                }else {answer += (*w)[k] ; }

        }
        return answer;
}
long double dataset::Ewy() {
        if ((w->no() != x->no()) || ( w->no() != y->no() ) ) {throw(63044); }
        long double answer = 0;
        for (int k=0;k < w->no() ; k++ ) {
                if (high) {double ii = (*w)[k] * (*y)[k];
                        if (answer) {answer = ii*(answer / ii + 1);}
                        if (k == 0) { answer = ii; }
                }else {answer +=(*w)[k] * (*y)[k]; }

        }
        return answer;
}

long double dataset::EIwxerrI2N(unsigned int no){
                if ((w->no() != x->no()) || ( w->no() != y->no() ) ) {throw(63044); }
        long double answer = 0;
        for (int k=0;k< w->no() ; k++ ) {
                if (high) {double ii = poww ((*w)[k] * (*erry)[k] * poww((*x)[k] , no),2);
                        if (answer) {answer = ii*(answer / ii + 1); }
                        if (k == 0) { answer = ii; }
                }else {answer += poww((*w)[k] * (*y)[k] * poww((*x)[k], no),2); }

        }
        return sqrt(answer)/w->no();

}
//best line fit functions ------------------------------------------------------


int bfit::goodchi(){ // good chi 0, suspect chi -1 or 1, bad chi 2 , -2;
        //int ans = 0 ;
        float chi05;
        float chi20;
        float chi80;
        float chi95;
        int dofree  = dofis() ;
        if (dofree<30) {
                 chi05 = chi05table[dof];
                 chi20 = chi20table[dof];
                 chi80 = chi80table[dof];
                 chi95 = chi95table[dof];
        } else { // estimate degree of freedom
                float xx2 =  2.0/9.0/(float)dofree;
                float xx = sqrt(xx2);
                xx2 = 1 - xx2;
                chi95 = dofree * poww( xx2 + xx * 1.64485 , 3);
                chi80 = dofree * poww( xx2 + xx * 0.84162 , 3);
                chi20 = dofree * poww( xx2 + xx * -0.84162 , 3);
                chi05 = dofree * poww( xx2 + xx * -1.64485 , 3);
        }
        double chis = chisqu();
        if (chis < chi05) return -2;
        if (chis > chi95) return 2;
        if (chis < chi20) return -1;
        if (chis > chi80) return 1;
        return 0;}

AnsiString bfit::stringchi(){
        switch(goodchi()){
        case 0 : return "good";
        case 1 : return "suspicious" ;
        case 2 : return "bad";
          }
}


void bfit::go(TChart *chart,TChart *rchart,int ii,bool mmethod) {


downer->calw();



fitl(mmethod);
plot(chart,rchart,ii);
}
void __fastcall bfit::plot(TChart *chart,TChart *rchart,int ii) {

bline = new TFastLineSeries(chart);
resid = new TErrPointSeries(rchart);
lnull = new TFastLineSeries(rchart);
bline->SeriesColor = (*downer->colour);
chart->AddSeries(bline);
rchart->AddSeries(resid);
rchart->AddSeries(lnull);
}

double bfit::chisqu(){
        double answer =0;
        for (int k = 0;k<downer->x->no();k++) {
                answer += (*downer->w)[k] * poww( (*downer->y)[k] - equ((*downer->x)[k]) , 2)  ;
}
return answer;
}


bool dataset::ToChart(TChart* chart){ return 0;
        }

// Fitting methods
void lineronlym::fitl(bool high) {
        lconst[0] = downer->EwyxN() / downer->EwxN(2);
        errconst[0] = sqrt(1 / downer->EwxN(2));
        }

double linermc::equ(double x){ return (x*lconst[0] + lconst[1]);}

void linermc::fitl(bool high) {
        if (!downer) throw(61301)  ;
        double ew = downer->Ew();
        double ewx = downer->EwxN() / ew;
        double ewy = downer->Ewy() / ew ;
        double ewyx = downer->EwyxN() /ew  ;
        double ewx2 = downer->EwxN(2) /ew ;
        if (!high) { lconst[0] = (ewyx -  (ewy * ewx))  / ( ewx2 - poww( ewx , 2) ) ;
                lconst[1] = ( ewy - lconst[0] * ewx ) ;

        }else{}
}

void bfit::newtonsmethod(vec(*fx)(vec *sol,dataset *d),vec(*fxerri)(vec *sol,double x,double y,double w)){
        vec fun(np+1);
        vec ans(np+1);
        vec oldans(np+1);
        vec stopvec(np+1);
        for (unsigned int c=0; c <= np ; c++ ) { ans[c] = lconst[c]; }
        oldans = ans.increment();
        matr jacobin(np+1,np+1);
        int q =0;
        try{
        bool loopwhile = true;
        list<vec> deltas;
        while (loopwhile) {
                fun.setzero();
                jacobin.setzero();
                fun =  fx(&ans,downer);
                jacobin = jacgen(fx,downer,&ans);
                vec delta = jacobin % fun;
                vec *d = new vec(delta.no());
                *d = delta;
                deltas.grow(d);
                if (deltas.len >= 10) { //test to see if converged or not
                                if ( (np+1) *  1.605e-8*ans.abssum() >= deltas.tail()->i->abssum()){loopwhile = false; continue;}
                                if (1000 * deltas.root->i->abssum() >= deltas.tail()->i->abssum() ){throw(0);}
                                deltas.remove(0);}
                oldans = ans;
                ans = ans - delta;
                }
        if (IsNan(ans[0])) { throw(0); }
        vec err(2);  err.setzero();

        for (int k=0;k< downer->w->no() ;k++)
        {
                double &xx = (*downer->x)[k];
                double &yy = (*downer->y)[k];
                double &ww = (*downer->w)[k];
                vec vv = fxerri(&ans,xx,yy,ww);
                vec anserri = jacobin % vv;
                for (unsigned int c=0; c <= np ; c++ ) { err[c] += poww( anserri[c],2) * poww((*downer->erry)[k],2);}
        }
        for (unsigned int c=0; c <= np ; c++ ) { lconst[c] = ans[c];
                errconst[c] = sqrt( err[c] ) ; }     //  */



        }catch(...) { errormessage_Engine("Fit failed to converge try improving estimates"); throw(0); return;}
        }



double linermcPC1181D::equ(double x){ return (x*lconst[0] + lconst[1]);}

void linermcPC1181D::fitl(bool high) {
        if (!downer) throw(61301)  ;
        double IwI = downer->Ew() / downer->x->no();
        double IwxI = downer->EwxN() / downer->x->no();
        double IwyI = downer->Ewy() / downer->x->no();
        double IwxyI = downer->EwyxN() / downer->x->no()  ;
        double Iwx2I = downer->EwxN(2) / downer->x->no();
        if (!high) {
        //m =
        lconst[0] = ( ( IwI * IwxyI ) - ( IwxI * IwyI ) ) /
                   ( ( IwI * Iwx2I )  - ( IwxI * IwxI ) );
        // c =
        lconst[1] = ( ( IwyI * Iwx2I ) -  ( IwxyI * IwxI ) ) /
                     ( (IwI * Iwx2I ) - ( IwxI * IwxI ) ) ;
        }else{}
}

double lineronlym::equ(double x){ return (x*lconst[0]);}

void lineronlym::genstrs(TStrings *S,bool on) {
        double val = lconst[0];
        double err = errconst[0];
        if (errconst[0]) { int sig = -1* floor(log10(errconst[0]));
        val = floor(val * poww(10,sig+2)) ;
        val /=poww(10,sig+2);
        err = floor(err * poww(10,sig+3))/poww(10,sig+3); }
        AnsiString m = FloatToStrF(val,ffGeneral,8,8) ;
        AnsiString errm = FloatToStrF(err,ffGeneral,8,8) ;
        if (on) S->Clear();
        S->Add("Fit of the form y=mx with");
        switch (goodchi()) {
                case 0 : break;
                case -2 : {S->Add(chi05msg); return;}
                case -1 : {S->Add(chi20msg); break;}
                case 1 : {S->Add(chi80msg); break;}
                case 2 : {S->Add(chi95msg); return;}
        }
        S->Add("m = " + m + " +/- " + errm);
}
void linermc::genstrs(TStrings *S,bool on) {
        AnsiString m = FloatToStrF(lconst[0],ffGeneral,8,8) ;
        AnsiString c = FloatToStrF(lconst[1],ffGeneral,8,8) ;
        if (on) S->Clear();
        S->Add("Fit of the form y=mx+c with");
        S->Add("m = " + m);
        S->Add("c = " + c);


        }
AnsiString linermc::func(int ii){
AnsiString fxx;
fxx += "y";//
fxx += ii ;
fxx += "(x) = ";
fxx += FloatToStrF(lconst[0],ffGeneral,8,8) + " *x";
return fxx;
}


AnsiString lineronlym::func(int ii){
AnsiString fxx;
fxx += "y";//
fxx += ii ;
fxx += "(x) = ";
fxx += FloatToStrF(lconst[0],ffGeneral,8,8) + " *x";
return fxx;
}

AnsiString linermcPC1181D::func(int ii){
AnsiString fxx;
fxx += "y";//
fxx += ii ;
fxx += "(x) = ";
fxx += FloatToStrF(lconst[0],ffGeneral,8,8) + " *x";
fxx += " + " + FloatToStrF(lconst[1],ffGeneral,8,8);
return fxx;
}

void linermcPC1181D::genstrs(TStrings *S,bool on) {
        AnsiString m = FloatToStrF(lconst[0],ffGeneral,8,8) ;
        AnsiString c = FloatToStrF(lconst[1],ffGeneral,8,8) ;
        if (on) S->Clear();
        S->Add("PC1181D Fit of the form y=mx+c with ");
        S->Add("m = " + m);
        S->Add("c = " + c);
        }
void polynom::fitl(bool high) {
        //np=2;
        //if (high) {
                double meanxx = -1 * downer->dmeanx();
                double meanyy = -1 * downer->dmeany();

                for (int j=0;j < downer->x->no();j++) {
                        (*downer->x)[j] += meanxx;
                        (*downer->y)[j] += meanyy;
                }
                int pindex = floor(logg(downer->dmaxx() - downer->dminx()));
                int pindey = floor(logg(downer->dmaxy() - downer->dminy()));
                for (int j=0;j < downer->x->no();j++) {
                        (*downer->x)[j] /= poww(10,pindex);
                        (*downer->y)[j] /= poww(10,pindey);
                }
        matr polym(np+1,np+1);
        vec polyv(np+1);
        vec polyerr(np+1);
        //vec anserr (np)

        for (int c=0; c <= np ; c++ )          {

                polyv[c] = downer->EwyxN(np-c)/downer->w->no();
                polyerr[c] = 0;//downer->EIwxerrI2N(np-c); //downer->EwerrxN(np-c);
                for (int r=0; r <= np ; r++ )          {
                        polym[r][c] = downer->EwxN(2*np-c-r)/downer->w->no();
        }        }

if (!high) {
        matr polyinv = polym.inv();
        for (int k=0;k< downer->w->no() ;k++)
        {
                vec polyerri(np+1);
                for ( int c=0; c <= np ; c++ ) { polyerri[c] = downer->wxiN(k,np-c);}
                vec anserri = polyinv * polyerri;
                for (int c=0; c <= np ; c++ ) { polyerr[c] += poww( anserri[c],2) * poww((*downer->erry)[k],2);}
        }
        vec ans = polyinv * polyv ;
        for (int c=0; c <= np ; c++ ) { lconst[c] = ans[c];
                errconst[c] = sqrt( polyerr[c] ) ;  }

}else{
        for (int k=0;k< downer->w->no() ;k++)
        {
                vec polyerri(np+1);
                for ( int c=0; c <= np ; c++ ) { polyerri[c] = downer->wxiN(k,np-c);}
                vec anserri = polym % polyerri;
                for ( int c=0; c <= np ; c++ ) { polyerr[c] += poww( anserri[c],2) * poww((*downer->erry)[k],2);}
        }
        vec ans = polym % polyv ;
        for (int c=0; c <= np ; c++ ) { lconst[c] = ans[c];
                errconst[c] = sqrt( polyerr[c] ) ; }
        }
                for (int j=0;j < downer->x->no();j++) {
                        (*downer->x)[j] *= poww(10,pindex);
                        (*downer->y)[j] *= poww(10,pindey);
                }
                for (int c = 0;c <= np;c++) {lconst[c] *= poww(10, pindex*(c-np)) * poww(10,pindey); }

                double par[10];
                for (int c = 0;c <= np;c++) {par[c] = 0; }
                for (int g = 0;g < np;g++) {

                        for (int h = g+1; h <= np ;h++)
                        {
                                //float coff = lconst[g] * cocombin(np-g,h-g) * poww (meanxx,h-g);
                                par[h] += lconst[g] * cocombin(np-g,h-g) * poww (meanxx,h-g);
                }     }
                lconst[np] -= meanyy;
                for (int c = 0;c <= np;c++) { lconst[c] += par[c]; }
                for (int j=0;j < downer->x->no();j++) {
                        (*downer->x)[j] -= meanxx;
                        (*downer->y)[j] -= meanyy;   }

                        return;
                                }
double polynom::equ(double x){
        double ans = 0;
        for (int k=0;k <= np; k++ ) { ans += lconst[k] * poww( x , np - k ); }
        return ans;}

AnsiString polynom::func(int ii){
AnsiString fxx;
fxx += "y";//
fxx += ii ;
fxx += "(x) = ";
for (int k=0;k <= np;k++) { if (k) fxx += " + ";
fxx += FloatToStrF(lconst[k],ffGeneral,8,8) + " *x**";
fxx += IntToStr(np-k);
}
return fxx;
}


void polynom::genstrs(TStrings *S,bool on) {
        row<AnsiString> aj(np+1);
        row<AnsiString> erraj(np+1);
        if (on) S->Clear();

        switch (np) {
        case 1 : {  S->Add("Fit of the form y=ax+b with "); break; }
        case 2 : {  S->Add("Fit of the form y=ax^2+bx+c with "); break; }
        case 3 : {  S->Add("Fit of the form y=ax^3+bx^2+cx+d with "); break; }
        default : {  S->Add("Fit of the form y= SUM{i=0:N}( a{i}x^(N-i) ) "); break; }
        }
                switch (goodchi()) {
                case 0 : break;
                case -2 : {S->Add(chi05msg); break;}
                case -1 : {S->Add(chi20msg); break;}
                case 1 : {S->Add(chi80msg); break;}
                case 2 : {S->Add(chi95msg); break;}
        }
        char coeq = 'a';
        for (int k=0;k <= np;k++) {
                        double val = lconst[k];
                        double err = errconst[k];
                        if (errconst[k]) { int sig = -1* floor(log10(errconst[k]));
                        val = floor(val * poww(10,sig+2)) ;
                        val /=poww(10,sig+2);
                        err = floor(err * poww(10,sig+2))/poww(10,sig+2);
                        }

                        aj[k] = FloatToStrF(val,ffGeneral,8,8) ;
                        erraj[k] = FloatToStrF(err,ffGeneral,8,8) ;
        AnsiString  coeqA = coeq;
        S->Add(  coeqA +" = " + aj[k] + " +/- " + erraj[k]);
        coeq++;
        }}

void __fastcall Mlinear::set(int amount) { for (int k=0; k < amount; k++) {  bfit *ntmp  = new polynom(1); iline->grow(ntmp); } };

void Mlinearsamem::fitl(TChart* MChart,TChart* resid,bool high){
                matrix = new matr(dat.len + 1,dat.len + 1 );
                equals = new vec(dat.len + 1 );
                vec ans(dat.len + 1);
                double sumx2 = 0;
                double sumxy = 0;
                for (unsigned int k=0;k < dat.len;k++) {
                        dat.get(k)->i->calw();
                        sumx2 += dat.get(k)->i->EwxN(2);
                        sumxy += dat.get(k)->i->EwyxN();
                        }
                        (*matrix)[0][0] = sumx2;
                        (*equals)[0] = sumxy;

                for (unsigned int c=1; c < dat.len + 1 ; c++){ // setting up of the matrix to solve
                        (*matrix)[0][c] = dat.get(c-1)->i->EwxN();
                                // dont get confused here doing 1st row not col
                        (*matrix)[c][0] = dat.get(c-1)->i->EwxN();
                        (*equals)[c] = dat.get(c-1)->i->Ewy();
                        for (unsigned int r=1; r < dat.len + 1 ; r++){
                                (*matrix)[r][c] = 0;
                                if  (c==r) {
                                        (*matrix)[r][c] = dat.get(r-1)->i->Ew();
                                } else { (*matrix)[r][c] = 0; }
                }   }
                vec polyerr(dat.len + 1 );
                polyerr.setzero();

// error analysis starts here

                if (high) { ans = (*matrix) % (*equals);  // Gaussian Elimination
        for (unsigned int ii=0;ii<dat.len;ii++) {
        dataset *downer = dat.get(ii)->i;
        for (int k=0;k< downer->x->no() ;k++)
        {
                vec polyerri(dat.len + 1 );
                polyerri.setzero();
                polyerri[0] = downer->wxiN(k,1);
                polyerri[ii+1] =(*downer->w)[k];
                vec anserri = (*matrix) % polyerri;
                polyerr[0] += poww( anserri[0],2) * poww((*downer->erry)[k],2);
                polyerr[ii+1] += poww( anserri[ii+1],2) * poww((*downer->erry)[k],2);
                }
        }
                }
                else {                           // Cramers Mehtod
                matr polyinv = matrix->inv();
        for (unsigned int ii=0;ii<dat.len;ii++) {
        dataset *downer = dat.get(ii)->i;
        for (int k=0;k< downer->x->no() ;k++)
        {
                vec polyerri(dat.len + 1 );
                polyerri.setzero();
                polyerri[0] = downer->wxiN(k,1);
                polyerri[ii+1] =(*downer->w)[k];
                vec anserri = polyinv * polyerri;
                polyerr[0] += poww( anserri[0],2) * poww((*downer->erry)[k],2);
                polyerr[ii+1] += poww( anserri[ii+1],2) * poww((*downer->erry)[k],2);
                }
        }
                ans = (*matrix).inv() *  (*equals); }

//ends here

                for (int k = 1; k< ans.no(); k++) {
                        iline->get(k-1)->i->errconst[0] = polyerr[0];
                        iline->get(k-1)->i->errconst[1] = polyerr[k];
                        iline->get(k-1)->i->lconst[0] = ans[0];
                        iline->get(k-1)->i->lconst[1] = ans[k];
                        iline->get(k-1)->i->plot(MChart,resid);
                }
}
void Mlinearsamec::fitl(TChart* MChart,TChart* resid,bool high){
                                   matrix = new matr(dat.len + 1,dat.len + 1 );
                equals = new vec(dat.len + 1 );

                double sumw = 0;
                double sumy = 0;
                for (unsigned int k=0;k < dat.len;k++) {
                        dat.get(k)->i->calw();
                        sumw += dat.get(k)->i->Ew();
                        sumy += dat.get(k)->i->Ewy();
                        }
                        (*matrix)[0][0] = sumw;
                        (*equals)[0] = sumy;

                for (unsigned int c=1; c < dat.len + 1 ; c++){ // setting up of the matrix to solve
                        (*matrix)[0][c] = dat.get(c-1)->i->EwxN();
                                // dont get confused here doing 1st row not col
                        (*matrix)[c][0] = dat.get(c-1)->i->EwxN();
                        (*equals)[c] = dat.get(c-1)->i->EwyxN();
                        for (unsigned int r=1; r < dat.len + 1 ; r++){
                                (*matrix)[r][c] = 0;
                                if  (c==r) {
                                        (*matrix)[r][c] = dat.get(r-1)->i->EwxN(2);
                                } else { (*matrix)[r][c] = 0; }
                }   }
                vec ans(dat.len + 1);
                vec polyerr(dat.len + 1 );
                polyerr.setzero();
// error analysis starts here

                if (high) { ans = (*matrix) % (*equals);  // Gaussian Elimination
        for (unsigned int ii=0;ii<dat.len;ii++) {
        dataset *downer = dat.get(ii)->i;
        for (int k=0;k< downer->x->no() ;k++)
        {
                vec polyerri(dat.len + 1 );
                polyerri.setzero();
                polyerri[0] = (*downer->w)[k];
                polyerri[ii+1] = downer->wxiN(k,1);
                vec anserri = (*matrix) % polyerri;
                polyerr[0] += poww( anserri[0],2) * poww((*downer->erry)[k],2);
                polyerr[ii+1] += poww( anserri[ii+1],2) * poww((*downer->erry)[k],2);
                }
        }
                }
                else {                           // Cramers Mehtod
                matr polyinv = matrix->inv();
        for (unsigned int ii=0;ii<dat.len;ii++) {
        dataset *downer = dat.get(ii)->i;
        for (int k=0;k< downer->x->no() ;k++)
        {
                vec polyerri(dat.len + 1 );
                polyerri.setzero();
                polyerri[0] = (*downer->w)[k];
                polyerri[ii+1] = downer->wxiN(k,1);
                vec anserri = polyinv * polyerri;
                polyerr[0] += poww( anserri[0],2) * poww((*downer->erry)[k],2);
                polyerr[ii+1] += poww( anserri[ii+1],2) * poww((*downer->erry)[k],2);
                }
        }
                ans = (*matrix).inv() *  (*equals); }

//ends here
                for (int k = 1; k< ans.no(); k++) {
                        iline->get(k-1)->i->errconst[1] = sqrt(polyerr[0]);
                        iline->get(k-1)->i->errconst[0] = sqrt(polyerr[k]);
                        iline->get(k-1)->i->lconst[1] = ans[0];
                        iline->get(k-1)->i->lconst[0] = ans[k];
                        iline->get(k-1)->i->plot(MChart,resid);
                }


}
//sample : power

double powersq::equ(double x){ return (lconst[0] * poww(x*x,(lconst[1]/2)));}

void powersq::fitl(bool high) { //fit routine
        dataset tmpd(0);   // create new dataset
        tmpd.x = new row<double>(downer->x->no()); // create a new row with length of the x component of the orginal dataset
        tmpd.y = new row<double>(downer->y->no()); // same but y
        tmpd.erry = new row<double>(downer->erry->no()); // same but error in y
        for (int a = 0; a < downer->x->no(); a++) { // populate x,y,erry for new dataset
                if (((*downer->x)[a]==0) || ((*downer->y)[a]==0) || ((*downer->erry)[a]==0)) {
                        errormessage_Engine("This fit routine can not have 0 value in the data for x or y"); throw(0); return; }
                (*tmpd.x)[a] = logg((*downer->x)[a] * (*downer->x)[a]); // new x=log(x^2)
                (*tmpd.y)[a] = logg((*downer->y)[a]);    // new y = log(y)
                (*tmpd.erry)[a] = (*downer->erry)[a] / (*downer->y)[a]; } // new err = err/y
        tmpd.calw(); // calculate w
        polynom tmpf(1);   // create a instance of a fit polynomial of the form y=ax+b lconst[0] = a, lconst[1] = b
        tmpf.downer = &tmpd; // set tmpd as the data to analysis
        useotherfit(&tmpf);  // perform the fit log(y)=log(A) + (B/2)*log(x^2)
        lconst[0] = poww(10 , tmpf.lconst[1]); // determine A from fit
        lconst[1] = tmpf.lconst[0]*2;          // determine B from fit
        errconst[0] = poww(10, tmpf.lconst[1]-1) * tmpf.lconst[1] * tmpf.errconst[1];  // determine error in A
        errconst[1] = tmpf.errconst[0];} // determine error in B

void  powersq::genstrs(TStrings *S,bool on) { // displays results
        //AnsiString m = FloatToStrF(lconst[0],ffGeneral,8,8) ;
        //AnsiString c = FloatToStrF(2*lconst[1],ffGeneral,8,8) ;
        //AnsiString em = FloatToStrF(errconst[0],ffGeneral,8,8) ;
        //AnsiString ec = FloatToStrF(2*errconst[1],ffGeneral,8,8) ;
        if (on) S->Clear();
        S->Add("Fit of the form y=A*(x^2)^(B/2) with ");
        switch (goodchi()) { // check if good chi
                case 0 : break;
                case -2 : {S->Add(chi05msg); break;}
                case -1 : {S->Add(chi20msg); break;}
                case 1 : {S->Add(chi80msg); break;}
                case 2 : {S->Add(chi95msg); break;}   }
        row<AnsiString> aj(np+1);
        row<AnsiString> erraj(np+1);
        char coeq = 'A';
        for (int k=0;k <= np;k++) {
                        double val = lconst[k];
                        double err = errconst[k];
                        if (errconst[k]) { int sig = -1* floor(log10(errconst[k]));
                        val = floor(val * poww(10,sig+2)) ;
                        val /=poww(10,sig+2);
                        err = floor(err * poww(10,sig+2))/poww(10,sig+2);
                        }

                        aj[k] = FloatToStrF(val,ffGeneral,8,8) ;
                        erraj[k] = FloatToStrF(err,ffGeneral,8,8) ;
        AnsiString  coeqA = coeq;
        //S->Add(  coeqA +" = " + aj[k] + " +/- " + erraj[k]);
        coeq++;
        }
        S->Add("A = " + aj[0]  +" +/- " + erraj[0] );
        S->Add("2*B = " + aj[1]  +" +/- " + erraj[1] );         }

AnsiString powersq::func(int ii){ // different way of repersenting the equation : when exporting to gnuplot.
        AnsiString fxx;
        fxx += "y"; fxx += ii ; fxx += "(x) = ";
        fxx += FloatToStrF(lconst[0],ffGeneral,8,8) + " *(x**2)**" + FloatToStrF(lconst[1]/2,ffGeneral,8,8) ;
        return fxx; }

//sample fit for resonance

double resonance::equ(double x){ return lconst[0] / ( pow(lconst[1] - x,2) + lconst[2] ) ;}

void resonance::fitl(bool high) { //fit routine
        try{
        dataset tmpd(0);   // create new dataset
        tmpd.x = new row<double>(downer->x->no()); // create a new row with length of the x component of the orginal dataset
        tmpd.y = new row<double>(downer->y->no()); // same but y
        tmpd.erry = new row<double>(downer->erry->no()); // same but error in y
        for (int a = 0; a < downer->x->no(); a++) { // populate x,y,erry for new dataset
                (*tmpd.x)[a] = (*downer->x)[a]; // new x=x
                (*tmpd.y)[a] = 1 / (*downer->y)[a];    // new y = 1/y
                (*tmpd.erry)[a] = 0.5 * (*downer->erry)[a] / poww((*downer->y)[a],2); } // new err = 0.5 * err/y^2
        tmpd.calw(); // calculate w
        polynom tmpf(2);   // create a instance of a fit polynomial of the form y=ax+b lconst[0] = a, lconst[1] = b
        tmpf.downer = &tmpd; // set tmpd as the data to analysis
        useotherfit(&tmpf);  // perform the fit 1/y=x/A -2(B*x) + (C+B^2) / A ;

        lconst[0] = 1 / tmpf.lconst[0];  // determine A from fit
        lconst[1] = lconst[0] * tmpf.lconst[1] / -2.0 ;  // determine B from fit
        lconst[2] = lconst[0] * tmpf.lconst[2] - pow(lconst[1],2); // determine C from fit
        errconst[0] = sqrt( poww( 0.5 * tmpf.errconst[0] / ( tmpf.lconst[0] * tmpf.lconst[0]),2)) ;
        errconst[1] = sqrt( poww( 0.5 * tmpf.errconst[1] / ( tmpf.lconst[0] ),2)
                + poww( 0.25 * tmpf.errconst[0] * tmpf.lconst[1] / ( tmpf.lconst[0] * tmpf.lconst[0] ),2) ) ;
        errconst[2] = sqrt( poww(-2.0 * lconst[1] * errconst[1] ,2)
                + poww( tmpf.errconst[2] / tmpf.lconst[0] , 2)
                + poww( 0.5 * tmpf.errconst[0]  / ( tmpf.lconst[0] * tmpf.lconst[0]),2)) ;
        //poww(10, tmpf.lconst[1]-1) * tmpf.lconst[1] * tmpf.errconst[1];  // determine error in A
        //errconst[1] = tmpf.errconst[0]; // determine error in B

        }catch(...) { errormessage_Engine("There is an error with your fit routine"); throw(0); return;}
        }
void resonance::genstrs(TStrings *S,bool on) { // displays results
        //AnsiString AA = FloatToStrF(lconst[0],ffGeneral,8,8) ;
        //AnsiString BB = FloatToStrF(lconst[1],ffGeneral,8,8) ;
        //AnsiString CC = FloatToStrF(lconst[2],ffGeneral,8,8) ;
        //AnsiString errAA = FloatToStrF(errconst[0],ffGeneral,8,8) ;
        //AnsiString errBB = FloatToStrF(errconst[1],ffGeneral,8,8) ;
        //AnsiString errCC = FloatToStrF(errconst[2],ffGeneral,8,8) ;
        if (on) S->Clear();
        S->Add("Fit of the form y=A*(x^2)^(B/2) with ");
        switch (goodchi()) { // check if good chi
                case 0 : break;
                case -2 : {S->Add(chi05msg); break;}
                case -1 : {S->Add(chi20msg); break;}
                case 1 : {S->Add(chi80msg); break;}
                case 2 : {S->Add(chi95msg); break;}   }
        char coeq = 'A';
        row<AnsiString> aj(np+1);
        row<AnsiString> erraj(np+1);
        for (int k=0;k <= np;k++) {
                        double val = lconst[k];
                        double err = errconst[k];
                        if (errconst[k]) { int sig = -1* floor(log10(errconst[k]));
                        val = floor(val * poww(10,sig+2)) ;
                        val /=poww(10,sig+2);
                        err = floor(err * poww(10,sig+2))/poww(10,sig+2);
                        }

                        aj[k] = FloatToStrF(val,ffGeneral,8,8) ;
                        erraj[k] = FloatToStrF(err,ffGeneral,8,8) ;
        AnsiString  coeqA = coeq;
        S->Add(  coeqA +" = " + aj[k] + " +/- " + erraj[k]);
        coeq++;
        }
        //S->Add("A = " + AA +" +/- " + errAA);
        //S->Add("B = " + BB + " +/- " + errBB );
        //S->Add("C = " + CC +" +/- " + errCC);
        }

AnsiString resonance::func(int ii){ // different way of repersenting the equation : when exporting to gnuplot.
        AnsiString fxx;
        fxx += "y"; fxx += ii ; fxx += "(x) = ";
        fxx += FloatToStrF(lconst[0],ffGeneral,8,8) + " / ( (" + FloatToStrF(lconst[1],ffGeneral,8,8) + " + x)**2 + " + FloatToStrF(lconst[2],ffGeneral,8,8) + " )";
        return fxx; }

// sample : perodicdecay


double periodicdecay::equ(double x){ return (lconst[0] * exp(lconst[1] * x) * cos(lconst[2] *x + lconst[3])) ;}

vec dfunperodic_dy(double A,double B,double c,double d,double x,double y,double w){
vec v(4);
v[0] = 2.0*w*exp(B*x)*cos(c*x+d);
v[1] = 2.0*w*A*x*exp(B*x)*cos(c*x+d);
v[2] = -2.0*w*A*exp(B*x)*sin(c*x+d)*x;
v[3] = -2.0*w*A*exp(B*x)*sin(c*x+d);
return v; }


vec funperodicdecay(vec *sol,dataset *dat){
vec v(4);
v.setzero();
double &A = (*sol)[0];
double &B = (*sol)[1];
double &c = (*sol)[2];
double &d = (*sol)[3];
try{
for (int h=0;h<dat->x->no();h++) {
double &x = (*dat->x)[h];
double &y = (*dat->y)[h];
double &w = (*dat->w)[h];
v[0] += w*(A*exp(B*x)*cos(c*x+d)-y)*exp(B*x)*cos(c*x+d);
v[1] += w*(A*exp(B*x)*cos(c*x+d)-y)*A*x*exp(B*x)*cos(c*x+d);
v[2] += -1.0*w*(A*exp(B*x)*cos(c*x+d)-y)*A*exp(B*x)*sin(c*x+d)*x;
v[3] += -1.0*w*(A*exp(B*x)*cos(c*x+d)-y)*A*exp(B*x)*sin(c*x+d);  }
}catch(...) {throw(90); }
return v; }

vec errperodicdecay(vec *sol,double x,double y,double w){
vec v(4);
double &A = (*sol)[0];
double &B = (*sol)[1];
double &c = (*sol)[2];
double &d = (*sol)[3];
try{
v[0] = w*exp(B*x)*cos(c*x+d);
v[1] = w*A*x*exp(B*x)*cos(c*x+d);
v[2] = -1.0*w*A*exp(B*x)*sin(c*x+d)*x;
v[3] = -1.0*w*A*exp(B*x)*sin(c*x+d);
}catch(...) {throw(90); }
return v; }

void periodicdecay::fitl(bool high) { //fit routine
        newtonsmethod(funperodicdecay,errperodicdecay);
        }
void periodicdecay::genstrs(TStrings *S,bool on) { // displays results
        if (on) S->Clear();
        S->Add("Fit of the form y=A*exp(B*x)*cos(C*x+D) with ");
        switch (goodchi()) { // check if good chi
                case 0 : break;
                case -2 : {S->Add(chi05msg); break;}
                case -1 : {S->Add(chi20msg); break;}
                case 1 : {S->Add(chi80msg); break;}
                case 2 : {S->Add(chi95msg); break;}   }
        char coeq = 'A';
        row<AnsiString> aj(np+1);
        row<AnsiString> erraj(np+1);
        for (int k=0;k <= np;k++) {
                        double val = lconst[k];
                        double err = errconst[k];
                        if (errconst[k]) { int sig = -1* floor(log10(errconst[k]));
                        val = floor(val * poww(10,sig+2)) ;
                        val /=poww(10,sig+2);
                        err = floor(err * poww(10,sig+2))/poww(10,sig+2);
                        }

                        aj[k] = FloatToStrF(val,ffGeneral,8,8) ;
                        erraj[k] = FloatToStrF(err,ffGeneral,8,8) ;
        AnsiString  coeqA = coeq;
        S->Add(  coeqA +" = " + aj[k] + " +/- " + erraj[k]);
        coeq++;
        }
        }

AnsiString periodicdecay::func(int ii){ // different way of repersenting the equation : when exporting to gnuplot.
        AnsiString fxx;
        fxx += "y"; fxx += ii ; fxx += "(x) = ";
        fxx += FloatToStrF(lconst[0],ffGeneral,8,8) + "*exp(x * "  + FloatToStrF(lconst[1],ffGeneral,8,8) +" )*cos(x * " + FloatToStrF(lconst[2],ffGeneral,8,8) + "  + "
                + FloatToStrF(lconst[3],ffGeneral,8,8) + " ) " ;
        return fxx; }

//-----------------------------------------------------------------------------

double periodic::equ(double x){ return (lconst[0] * cos(lconst[1]*x)) ;}
vec funperodic(vec *sol,dataset *d){
vec v(2);
v.setzero();
double &A = (*sol)[0];
double &B = (*sol)[1];
try{
for (int h=0;h<d->x->no();h++) {
double &x = (*d->x)[h];
double &y = (*d->y)[h];
double &w = (*d->w)[h];
v[0] -= w*(y-A*cos(B*x))*cos(B*x);
v[1] += w*(y-A*cos(B*x))*A*x*sin(B*x); }
}catch(...) {throw(90); }
return v; }

vec errperodics(vec *sol,double x,double y,double w){
vec v(2);
double &A = (*sol)[0];
double &B = (*sol)[1];
try{
v[0] = w*cos(B*x);
v[1] = w*A*x*sin(B*x);
}catch(...) {throw(90); }
return v; }


matr jacperodics(vec *sol,double x,double y,double w){  // never used
double &A = (*sol)[0];
double &B = (*sol)[1];
matr m(2,2);
try{

    m[0][0] = poww(w*cos(B*x),2);
    m[0][1] = w*A*sin(B*x)*x*cos(B*x)+w*(y-A*cos(B*x))*sin(B*x)*x;
    m[1][0] = w*A*sin(B*x)*x*cos(B*x)+w*(y-A*cos(B*x))*sin(B*x)*x;
    m[1][1] = w*poww(A,2)*poww(sin(B*x),2)*x*x+w*(y-A*cos(B*x))*A*cos(B*x)*x*x;
}catch(...) {throw(90); }
return m;}




void periodic::fitl(bool high) { //fit routine
        newtonsmethod(funperodic,errperodics);
        }
void periodic::genstrs(TStrings *S,bool on) { // displays results
        //AnsiString AA = FloatToStrF(lconst[0],ffGeneral,8,8) ;
        //AnsiString BB = FloatToStrF(lconst[1],ffGeneral,8,8) ;
        //AnsiString errAA = FloatToStrF(errconst[0],ffGeneral,8,8) ;
        //AnsiString errBB = FloatToStrF(errconst[1],ffGeneral,8,8) ;
        if (on) S->Clear();
        S->Add("Fit of the form y=A*cos(B*x) with ");
        switch (goodchi()) { // check if good chi
                case 0 : break;
                case -2 : {S->Add(chi05msg); break;}
                case -1 : {S->Add(chi20msg); break;}
                case 1 : {S->Add(chi80msg); break;}
                case 2 : {S->Add(chi95msg); break;}   }
        char coeq = 'A';
        row<AnsiString> aj(np+1);
        row<AnsiString> erraj(np+1);
        for (int k=0;k <= np;k++) {
                        double val = lconst[k];
                        double err = errconst[k];
                        if (errconst[k]) { int sig = -1* floor(log10(errconst[k]));
                        val = floor(val * poww(10,sig+2)) ;
                        val /=poww(10,sig+2);
                        err = floor(err * poww(10,sig+2))/poww(10,sig+2);
                        }

                        aj[k] = FloatToStrF(val,ffGeneral,8,8) ;
                        erraj[k] = FloatToStrF(err,ffGeneral,8,8) ;
        AnsiString  coeqA = coeq;
        S->Add(  coeqA +" = " + aj[k] + " +/- " + erraj[k]);
        coeq++;
        }
        }

AnsiString periodic::func(int ii){ // different way of repersenting the equation : when exporting to gnuplot.
        AnsiString fxx;
        fxx += "y"; fxx += ii ; fxx += "(x) = ";
        fxx += FloatToStrF(lconst[0],ffGeneral,8,8) + " *cos(" + FloatToStrF(lconst[1],ffGeneral,8,8) + " *x)";
        return fxx; }