/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Sebastian Raducci 01/2004  <mailto:raducci@fisica.uniud.it>
!
!   Copyright: MAGIC Software Development, 2001-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          // 
//  Cubic Spline Interpolation                                              //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MCubicSpline.h"

#include "MCubicCoeff.h"

#include "TMath.h"

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MCubicSpline);

using namespace std;

//---------------------------------------------------------------------------
//
// Contructor 
//
//
MCubicSpline::MCubicSpline(Byte_t *y, Byte_t *x, Bool_t areAllEq,
			   Int_t n, Double_t begSD, Double_t endSD):
    fN(n)
{
    Init(y,x,areAllEq,n,begSD,endSD);
}

//---------------------------------------------------------------------------
//
// Constructor for FADC slice (only the FADC counts are needed)
//
//
MCubicSpline::MCubicSpline(Byte_t *y)
{
    Byte_t x[]={0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0A,0x0B,0x0C,0x0D,0x0E};
    fN = 15;
    Init(y,x,kTRUE,15,0.0,0.0);
}

//---------------------------------------------------------------------------
//
// Constructors common part
//
//
void MCubicSpline::Init(Byte_t *y, Byte_t *x, Bool_t areAllEq,
			   Int_t n, Double_t begSD, Double_t endSD)

{
    Double_t *temp = new Double_t[fN-1];
    Double_t *ysd = new Double_t[fN-1];
    Double_t p,h;
    MCubicCoeff *tempCoeff;
    fCoeff = new TObjArray(fN-1,0);
    
    if (areAllEq)
    {
	h = x[1]-x[0];
	ysd[0]=temp[0]=begSD;
	ysd[n-1]=endSD;
	for(Int_t i = 1; i < n-1; i++)
	{
	    p = 0.5*ysd[i-1]+2.0;
	    ysd[i] = (-0.5)/p;
	    temp[i] = (y[i+1]-2*y[i]+y[i-1])/h;
	    temp[i] = (6.0*temp[i]/h-0.5*temp[i-1])/p;
	}
	for(Int_t i = n-2; i > 0; i--)
	    ysd[i]=ysd[i]*ysd[i+1]+temp[i];
	for(Int_t i = 0; i < n-1; i++)
	{
	    tempCoeff = new MCubicCoeff(x[i],x[i+1],y[i],y[i+1],(ysd[i+1]-ysd[i])/(6*h),
					ysd[i]/2.0,(y[i+1]-y[i])/h-(h*(ysd[i+1]+2*ysd[i]))/6);
	    fCoeff->AddAt(tempCoeff,i);
	}
	delete [] temp;
	delete [] ysd;
    }
    else
    {
	Double_t sig;
	ysd[0]=temp[0]=begSD;
	ysd[n-1]=endSD;
	for(Int_t i = 1; i < n-1; i++)
	{
	    sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
	    p = sig*ysd[i-1]+2.0;
	    ysd[i] = (sig-1.0)/p;
	    temp[i] = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
	    temp[i] = (6.0*temp[i]/(x[i+1]-x[i-1])-sig*temp[i-1])/p;
	}
	for(Int_t i = n-2; i > 0; i--)
	    ysd[i]=ysd[i]*ysd[i+1]+temp[i];
	for(Int_t i = 0; i < n-1; i++)
	{
	    h = x[i+1]-x[i];
	    tempCoeff = new MCubicCoeff(x[i],x[i+1],y[i],y[i+1],(ysd[i+1]-ysd[i])/(6*h),
					ysd[i]/2.0,(y[i+1]-y[i])/h-(h*(ysd[i+1]+2*ysd[i]))/6);
	    fCoeff->AddAt(tempCoeff,i);
	}
	delete [] temp;
	delete [] ysd;
    }
}

MCubicSpline::~MCubicSpline()
{
    fCoeff->Delete();
}

//---------------------------------------------------------------------------
//
// Evaluate the spline at a given point
//
Double_t MCubicSpline :: Eval(Double_t x)
{
    for (Int_t i = 0; i < fN-1; i++)
	if (((MCubicCoeff*)fCoeff->At(i))->IsIn(x))
	    return ((MCubicCoeff*)fCoeff->At(i))->Eval(x);
    gLog << warn << "Cannot evaluate Spline at " << x << "; returning 0";
    return 0.0;
}

//----------------------------------------------------------------------------
//
// Search for max
//
Double_t MCubicSpline :: EvalMax()
{
    Double_t temp_max;
    Double_t max = ((MCubicCoeff*)fCoeff->At(0))->GetMax();
    for (Int_t i = 1; i < fN-1; i++)
    {
	temp_max = ((MCubicCoeff*)fCoeff->At(i))->GetMax();
        if (temp_max > max)
	    max = temp_max;
    }
    return max;
}

//----------------------------------------------------------------------------
//
// Search for min
//
Double_t MCubicSpline :: EvalMin()
{
    Double_t temp_min;
    Double_t min = ((MCubicCoeff*)fCoeff->At(0))->GetMin();
    for (Int_t i = 1; i < fN-1; i++)
    {
	temp_min = ((MCubicCoeff*)fCoeff->At(i))->GetMin();
        if (temp_min < min)
	    min = temp_min;
    }
    return min;
}

//----------------------------------------------------------------------------
//
// Search for abscissa of the max
//
Double_t MCubicSpline :: EvalAbMax()
{
    Double_t temp_max;
    Double_t abMax = ((MCubicCoeff*)fCoeff->At(0))->GetAbMax();
    Double_t max = ((MCubicCoeff*)fCoeff->At(0))->GetMax();
    for (Int_t i = 1; i < fN-1; i++)
    {
	temp_max = ((MCubicCoeff*)fCoeff->At(i))->GetMax();
        if (temp_max > max)
	{
	    max = temp_max;
	    abMax = ((MCubicCoeff*)fCoeff->At(i))->GetAbMax();
	}
    }
    return abMax;
}

//----------------------------------------------------------------------------
//
// Search for abscissa of the min
//
Double_t MCubicSpline :: EvalAbMin()
{
    Double_t temp_min;
    Double_t abMin = ((MCubicCoeff*)fCoeff->At(0))->GetAbMin();
    Double_t min = ((MCubicCoeff*)fCoeff->At(0))->GetMin();
    for (Int_t i = 1; i < fN-1; i++)
    {
	temp_min = ((MCubicCoeff*)fCoeff->At(i))->GetMin();
        if (temp_min < min)
	{
	    min = temp_min;
	    abMin = ((MCubicCoeff*)fCoeff->At(i))->GetAbMin();
	}
    }
    return abMin;
}

//----------------------------------------------------------------------------
//
// Finds the abscissa where the spline reaches y starting from x0 going in
// direction direction
// You have to give as input a starting point and a direction ("l" or "r")
//
Double_t MCubicSpline :: FindVal(Double_t y, Double_t x0, Char_t direction = 'l')
{
    Short_t whichRoot;
    Double_t tempRoot;
    Double_t *roots = new Double_t[3];

    for (Int_t i = 0; i < fN-1; i++)
    {
	if(((MCubicCoeff*)fCoeff->At(i))->IsIn(x0))
	{
	    if(direction == 'l')
	    {
		for (Int_t j = i; j >= 0; j--)
		{
		    whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
		    if (whichRoot >= 0 )
		    {
			tempRoot = roots[whichRoot];
			delete [] roots;
			return tempRoot;
		    }
		}
	    }
	    if(direction == 'r')
	    {
		for (Int_t j = i; j < fN-1; j++)
		{
		    whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
		    if (whichRoot >= 0)
		    {
			tempRoot = roots[whichRoot];
			delete [] roots;
			return tempRoot;
		    }
		}
	    }
	}
    }
    gLog << warn << "Nothing found calling MCubicSpline :: FindVal(), returning 0" << endl;
    return 0.0;
}
