source: trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc@ 3022

/* ======================================================================== *\
!
! *
! * 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 "MCubicCoeff.h"

#include "TMath.h"

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

ClassImp(MCubicCoeff);

using namespace std;

//----------------------------------------------------------------------------
//
// Constructor (The spline is: fA(x-fX)3+fB(x-fX)2+fC(x-fX)+fY
// where x is the independent variable, 2 and 3 are exponents
//

MCubicCoeff::MCubicCoeff(Double_t x, Double_t xNext, Double_t y, Double_t yNext, 
                         Double_t a, Double_t b, Double_t c) : 
    fX(x), fXNext(xNext), fA(a), fB(b), fC(c), fY(y), fYNext(yNext)
{
    fH = fXNext - fX;
    if(!EvalMinMax())
    {
        gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl;
        fMin = 0;
        fMax = 0;
    }
}

//----------------------------------------------------------------------------
//
// Evaluate the spline at a given point
//

Double_t MCubicCoeff::Eval(Double_t x)
{
    Double_t dx = x - fX;
    return (fY+dx*(fC+dx*(fB+dx*fA)));
}

//----------------------------------------------------------------------------
//
// Find min and max using derivatives. The min and max could be at the begin
// or at the end of the interval or somewhere inside the interval (in this case
// a comparison between the first derivative and zero is made)
// The first derivative coefficients are obviously: 3*fA, 2*fB, fC
//
Bool_t MCubicCoeff::EvalMinMax()
{
    fMin = fMax = fY;
    fAbMin = fAbMax = fX;
    if (fYNext < fMin)
    {
        fMin = fYNext;
        fAbMin = fXNext;
    }
    if (fYNext > fMax)
    {
        fMax = fYNext;
        fAbMax = fXNext;
    }
    Double_t tempMinMax;
    Double_t delta = 4.0*fB*fB - 12.0*fA*fC;
        if (delta >= 0.0 && fA != 0.0)
        {
            Double_t sqrtDelta = sqrt(delta);
            Double_t xPlus  = (-2.0*fB + sqrtDelta)/(6.0*fA);
            Double_t xMinus = (-2.0*fB - sqrtDelta)/(6.0*fA);
            if (xPlus >= 0.0 && xPlus <= fH)
            {
                tempMinMax = this->Eval(fX+xPlus);
                if (tempMinMax < fMin)
                {
                    fMin = tempMinMax;
                    fAbMin = fX + xPlus;
                }
                if (tempMinMax > fMax)
                {
                    fMax = tempMinMax;
                    fAbMax = fX + xPlus;
                }
            }
            if (xMinus >= 0.0 && xMinus <= fH)
            {
                tempMinMax = this->Eval(fX+xMinus);
                if (tempMinMax < fMin)
                {
                    fMin = tempMinMax;
                    fAbMin = fX + xMinus;
                }
                if (tempMinMax > fMax)
                {
                    fMax = tempMinMax;
                    fAbMax = fX + xMinus;
                }
            }
        }
/* if fA is zero then we have only one possible solution */
        else if (fA == 0.0 && fB != 0.0)
        {
            Double_t xSolo = - (fC/(2.0*fB));
            if (xSolo >= 0.0 && xSolo <= fH)
            {
                tempMinMax = this->Eval(fX+xSolo);
                if (tempMinMax < fMin)
                {
                    fMin = tempMinMax;
                    fAbMin = fX + xSolo;
                }
                if (tempMinMax > fMax)
                {
                    fMax = tempMinMax;
                    fAbMax = fX + xSolo;
                }
            }
        }
    return kTRUE;
}
//-------------------------------------------------------------------------
//
// Given y finds x using the cubic (cardan) formula.
//
// we consider the following form: x3 + ax2 + bx + c = 0 where
//   a = fB/fA, b = fC/fA, c = (fY - y)/fA  
//
// There could be three or one real solutions
//

Short_t MCubicCoeff::FindCardanRoot(Double_t y, Double_t *x)
{

    Short_t whichRoot = -1;

    Double_t a = fB/fA;
    Double_t b = fC/fA;
    Double_t c = (fY - y)/fA;

    Double_t q = (a*a-3.0*b)/9.0;
    Double_t r = (2.0*a*a*a-9.0*a*b+27.0*c)/54.0;

    Double_t aOver3 = a/3.0;
    Double_t r2 = r*r;
    Double_t q3 = q*q*q;
    
    if (r2 < q3) //3 real sol
    {
        Double_t sqrtQ = TMath::Sqrt(q);
        Double_t min2SqQ = -2.0*sqrtQ;
        Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ));

        x[0] = min2SqQ * TMath::Cos(theta/3.0) - aOver3;
        x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3.0) - aOver3;
        x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3.0) - aOver3;
        for (Int_t i = 0; i < 3; i++)
            if (x[i] >= 0.0 && x[i] <= fH)
            {
                x[i] = x[i] + fX;
                whichRoot = i;
                break;
            }
    }
    else //only 1 real sol
    {
        Double_t s;
        if (r == 0.0)
            s = 0.0;
        else if (r < 0.0)
            s = TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0);
        else // r > 0.0
            s = - TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0);
        if (s == 0.0)
            x[0] = - aOver3;
        else
            x[0] = (s + q/s) - aOver3;
        if (x[0] >= 0.0 && x[0] <= fH)
        {
            x[0] = x[0] + fX;
            whichRoot = 0;
        }
    }
    return whichRoot;
}

//------------------------------------------------------------------------------------
//
// return true if x is in this interval
//

Bool_t MCubicCoeff :: IsIn(Double_t x)
{
    if (x >= fX && x <= fXNext)
        return kTRUE;
    else
        return kFALSE;
}