#include "MBending.h"

#include <fstream.h>

#include <TMinuit.h>


ClassImp(MBending);

void MBending::Load(const char *name)
{
    /*
     ! MMT 1987 July 8
     ! T   36   7.3622   41.448  -0.0481
     !   IA        -37.5465    20.80602
     !   IE        -13.9180     1.25217
     !   NPAE       +7.0751    26.44763
     !   CA         -6.9149    32.05358
     !   AN         +0.5053     1.40956
     !   AW         -2.2016     1.37480
     ! END
     */

    ifstream fin(name);
    if (!fin)
    {
        cout << "Error: Cannot open file '" << name << "'" << endl;
        return;
    }

    char c;
    while (fin && fin.get()!='\n');
    fin >> c;

    if (c!='S' && c!='s')
    {
        cout << "Error: This in not a model correcting the star position (" << c << ")" << endl;
        return;
    }

    Clear();

    cout << endl;

    Double_t val;
    fin >> val;
    cout << "Number of observed stars: " << val << endl;
    fin >> val;
    cout << "Sky RMS: " << val << "\"" << endl;
    fin >> val;
    cout << "Refraction Constant A: " << val << "\"" << endl;
    fin >> val;
    cout << "Refraction Constant B: " << val << "\"" << endl;

    cout << endl;

    cout << "  & = Name      Value           Sigma" << endl;
    cout << "------------------------------------------" << endl;

    while (fin)
    {
        TString str;
        fin >> str;

        if (str=="END")
            break;

        if (str[0]=='&')
        {
            cout << " & ";
            str.Remove(0);
        }
        else
            cout << "   ";

        if (str[1]=='=')
        {
            cout << "=  ";
            str.Remove(0);
        }
        else
            cout << "   ";

        fin >> val;
        cout << str << "\t" << val << "°      \t";
        val *= kDeg2Rad;

        if (str=="IA")   fIa   = val;
        if (str=="IE")   fIe   = val;
        if (str=="CA")   fCa   = val;
        if (str=="AN")   fAn   = val;
        if (str=="AW")   fAw   = val;
        if (str=="NRX")  fNrx  = val;
        if (str=="NRY")  fNry  = val;
        if (str=="CRX")  fCrx  = val;
        if (str=="CRY")  fCry  = val;
        if (str=="NPAE") fNpae = val;
        if (str=="ECES") fEces = val;
        if (str=="ACES") fAces = val;
        if (str=="ECEC") fEcec = val;
        if (str=="ACEC") fAcec = val;

        fin >> val;
        cout << val*kDeg2Rad << endl;
    }
    cout << endl;
}

void MBending::Save(const char *name) {}

AltAz MBending::Correct(const AltAz &aa) const
{
    // Correct [rad]
    // zdaz    [rad]
    AltAz p = aa;

    const AltAz CES(-fEces*kDeg2Rad*sin(p.Alt()), -fAces*kDeg2Rad*sin(p.Az()));
    const AltAz CEC(-fEcec*kDeg2Rad*cos(p.Alt()), -fAcec*kDeg2Rad*cos(p.Az()));
    p += CES;
    p += CEC;

    const AltAz CRX(-fCrx*kDeg2Rad*sin(p.Az()-p.Alt()),  fCrx*kDeg2Rad*cos(p.Az()-p.Alt())/cos(p.Alt()));
    const AltAz CRY(-fCry*kDeg2Rad*cos(p.Az()-p.Alt()), -fCry*kDeg2Rad*sin(p.Az()-p.Alt())/cos(p.Alt()));
    p += CRX;
    p += CRY;

    const AltAz NRX(fNrx*kDeg2Rad*sin(p.Alt()), -fNrx*kDeg2Rad);
    const AltAz NRY(fNry*kDeg2Rad*cos(p.Alt()), -fNry*kDeg2Rad*tan(p.Alt()));
    p += NRX;
    p += NRY;

    const AltAz AW( fAw*kDeg2Rad*sin(p.Az()), -fAw*kDeg2Rad*cos(p.Az())*tan(p.Alt()));
    const AltAz AN(-fAn*kDeg2Rad*cos(p.Az()), -fAn*kDeg2Rad*sin(p.Az())*tan(p.Alt()));
    p += AW;
    p += AN;

    const AltAz CA(0, -fCa*kDeg2Rad/cos(p.Alt()));
    p += CA;

    const AltAz NPAE(0, -fNpae*kDeg2Rad*tan(p.Alt()));
    p += NPAE;

    const AltAz I(fIe*kDeg2Rad, fIa*kDeg2Rad);
    p += I;

    return p;
}

ZdAz MBending::Correct(const ZdAz &zdaz) const
{
    // Correct [rad]
    // zdaz    [rad]
    AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
    AltAz c = Correct(p);
    return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
}

void MBending::SetParameters(const Double_t *par, Int_t n=14)
{
    Clear();

    switch (n)
    {
    case 14:
        fAcec=par[13] ; // Azimuth Centering Error (cos)
    case 13:
        fEcec=par[12] ; // Elevation Centering Error (cos)
    case 12:
        fAces=par[11] ; // Azimuth Centering Error (sin)
    case 11:
        fEces=par[10] ; // Elevation Centering Error (sin)
    case 10:
        fCry =par[9] ; // Alt/Az Coude Displacement (E-W)
    case 9:
        fCrx =par[8] ; // Alt/Az Coude Displacement (N-S)
    case 8:
        fNry =par[7] ; // Nasmyth rotator displacement, vertical
    case 7:
        fNrx =par[6] ; // Nasmyth rotator displacement, horizontan
    case 6:
        fAw  =par[5] ; // Azimuth Axis Misalignment (E-W)
    case 5:
        fAn  =par[4] ; // Azimuth Axis Misalignment (N-S)
    case 4:
        fCa  =par[3] ; // Left-Right Collimation Error
    case 3:
        fNpae=par[2] ; // Az-El Nonperpendicularity
    case 2:
        fIe  =par[1] ; // Index Error in Elevation
    case 1:
        fIa  =par[0] ; // Index Error in Azimuth
    }
}

void MBending::GetParameters(Double_t *par, Int_t n=14) const
{
    switch (n)
    {
    case 14:
        par[13]=fAcec; // Azimuth Centering Error (cos)
    case 13:
        par[12]=fEcec; // Elevation Centering Error (cos)
    case 12:
        par[11]=fAces ; // Azimuth Centering Error (sin)
    case 11:
        par[10]=fEces ; // Elevation Centering Error (sin)
    case 10:
        par[9]=fCry ; // Alt/Az Coude Displacement (E-W)
    case 9:
        par[8]=fCrx ; // Alt/Az Coude Displacement (N-S)
    case 8:
        par[7]=fNry ; // Nasmyth rotator displacement, vertical
    case 7:
        par[6]=fNrx ; // Nasmyth rotator displacement, horizontan
    case 6:
        par[5]=fAw ; // Azimuth Axis Misalignment (E-W)
    case 5:
        par[4]=fAn ; // Azimuth Axis Misalignment (N-S)
    case 4:
        par[3]=fCa ; // Left-Right Collimation Error
    case 3:
        par[2]=fNpae ; // Az-El Nonperpendicularity
    case 2:
        par[1]=fIe ; // Index Error in Elevation
    case 1:
        par[0]=fIa ; // Index Error in Azimuth
    }
}

void MBending::SetMinuitParameters(TMinuit &m, Int_t n=-1) const
{
    if (n<0)
        n = m.GetNumPars();

    Int_t ierflg = 0;

    switch (n)
    {
    case 14:
    case 13:
    case 12:
    case 11:
    case 10:
    case 9:
    case 8:
    case 7:
    case 5:
    case 4:
        m.mnparm(3, "CA",     fCa,        1, -360, 360, ierflg);
        // cout << "Init 3 CA:    " << fCa << endl;
    case 3:
        m.mnparm(2, "NPAE",   fNpae,      1, -360, 360, ierflg);
        // cout << "Init 2 NPAE:  " << fNpae << endl;
    case 2:
        m.mnparm(1, "IE",     fIe,        1, -360, 360, ierflg);
        // cout << "Init 1 IE:    " << fIe << endl;
    case 1:
        m.mnparm(0, "IA",     fIa,        1, -360, 360, ierflg);
        // cout << "Init 0 IA:    " << fIa << endl;
    }
}

void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1)
{
    if (n<0)
        n = m.GetNumPars();

    Double_t err;

    switch (n)
    {
    case 14:
    case 13:
    case 12:
    case 11:
    case 10:
    case 9:
    case 8:
    case 7:
    case 6:
    case 5:
    case 4:
        m.GetParameter(3, fCa, err);
    case 3:
        m.GetParameter(2, fNpae, err);
    case 2:
        m.GetParameter(1, fIe, err);
    case 1:
        m.GetParameter(0, fIa, err);
    }
}

void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const
{
    if (n<0)
        n = m.GetNumPars();

    Double_t par, err;

    switch (n)
    {
    case 14:
    case 13:
    case 12:
    case 11:
    case 10:
    case 9:
    case 8:
    case 7:
    case 6:
    case 5:
    case 4:
        m.GetParameter(3, par, err);
        cout << " 3 CA:  " << par << " +- " << err << endl;
    case 3:
        m.GetParameter(2, par, err);
        cout << " 2 NPAE:  " << par << " +- " << err << endl;
    case 2:
        m.GetParameter(1, par, err);
        cout << " 1 IE:  " << par << " +- " << err << endl;
    case 1:
        m.GetParameter(0, par, err);
        cout << " 0 IA:  " << par << " +- " << err << endl;
    }
}