#include "MBending.h"

#include <fstream.h>
#include <iomanip.h>

#include <TMinuit.h>

#include "timer.h"

ClassImp(MBending);

void MBending::Reset()
{
    Clear();
}

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;
        //if (str=="MAGIC1") fMagic1 = val;
        //if (str=="MAGIC2") fMagic2 = val;

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

void MBending::Save(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
     */

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

    Timer t;
    t.Now();

    fout << "MAGIC1 " << t.GetTimeStr() << endl;
    fout << "S   00   000000   000000  0000000" << endl;
    fout << setprecision(8);
    fout << " IA     " << kRad2Deg*fIa   << " -1" << endl;
    fout << " IE     " << kRad2Deg*fIe   << " -1" << endl;
    fout << " CA     " << kRad2Deg*fNpae << " -1" << endl;
    fout << " NPAE   " << kRad2Deg*fCa   << " -1" << endl;
    fout << " AN     " << kRad2Deg*fAn   << " -1" << endl;
    fout << " AW     " << kRad2Deg*fAw   << " -1" << endl;
    fout << " NRX    " << kRad2Deg*fNrx  << " -1" << endl;
    fout << " NRY    " << kRad2Deg*fNry  << " -1" << endl;
    fout << " CRX    " << kRad2Deg*fCrx  << " -1" << endl;
    fout << " CRY    " << kRad2Deg*fCry  << " -1" << endl;
    fout << " ECES   " << kRad2Deg*fEces << " -1" << endl;
    fout << " ACES   " << kRad2Deg*fAces << " -1" << endl;
    fout << " ECEC   " << kRad2Deg*fEcec << " -1" << endl;
    fout << " ACEC   " << kRad2Deg*fAcec << " -1" << endl;
    //fout << " MAGIC1 " << -kRad2Deg*fMagic1 << " -1" << endl;
    //fout << " MAGIC2 " << -kRad2Deg*fMagic2 << " -1" << endl;
    fout << "END" << endl;
}

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

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

    //    const AltAz MAGIC(0, -fMagic1*tan(p.Alt())-fMagic2);
    //    p += MAGIC;

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

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

    //    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
    //    p += MAGIC;

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

    //    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
    //    p += MAGIC;

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

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

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

    return p;
}

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

    const AltAz I(fIe, fIa);
    p -= I;

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

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

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

    //    const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);
    //    p -= MAGIC;

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

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

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

    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());
}

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

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

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

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

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

    Int_t ierflg = 0;

    switch (n)
    {
    case 16:
        m.mnparm(15,"MAGIC2", fMagic2*kRad2Deg,  1, -360, 360, ierflg);
    case 15:
        m.mnparm(14,"MAGIC1", fMagic1*kRad2Deg,  1, -360, 360, ierflg);
    case 14:
        m.mnparm(13,"ACEC", fAcec*kRad2Deg,  1, -360, 360, ierflg);
    case 13:
        m.mnparm(12,"ECEC", fEcec*kRad2Deg,  1, -360, 360, ierflg);
    case 12:
        m.mnparm(11,"ACES", fAcec*kRad2Deg,  1, -360, 360, ierflg);
    case 11:
        m.mnparm(10,"ECES", fEcec*kRad2Deg,  1, -360, 360, ierflg);
    case 10:
        m.mnparm(9, "CRY",  fCry*kRad2Deg,  1, -360, 360, ierflg);
    case 9:
        m.mnparm(8, "CRX",  fCrx*kRad2Deg,  1, -360, 360, ierflg);
    case 8:
        m.mnparm(7, "NRY",  fNry*kRad2Deg,  1, -360, 360, ierflg);
    case 7:
        m.mnparm(6, "NRX",  fNrx*kRad2Deg,  1, -360, 360, ierflg);
    case 6:
        m.mnparm(5, "AW",   fAw*kRad2Deg,   1, -360, 360, ierflg);
    case 5:
        m.mnparm(4, "AN",   fAn*kRad2Deg,   1, -360, 360, ierflg);
    case 4:
        m.mnparm(3, "CA",   fCa*kRad2Deg,   1, -360, 360, ierflg);
    case 3:
        m.mnparm(2, "NPAE", fNpae*kRad2Deg, 1, -360, 360, ierflg);
    case 2:
        m.mnparm(1, "IE",   fIe*kRad2Deg,   1,  -90,  90, ierflg);
    case 1:
        m.mnparm(0, "IA",   fIa*kRad2Deg,   1, -360, 360, ierflg);
    }
}

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

    Double_t err;

    switch (n)
    {
    case 16:
        m.GetParameter(15, fMagic2, err);
        fMagic2 /= kRad2Deg;
    case 15:
        m.GetParameter(14, fMagic1, err);
        fMagic1 /= kRad2Deg;
    case 14:
        m.GetParameter(13, fAcec, err);
        fAcec /= kRad2Deg;
    case 13:
        m.GetParameter(12, fEcec, err);
        fEcec /= kRad2Deg;
    case 12:
        m.GetParameter(11, fAces, err);
        fAces /= kRad2Deg;
    case 11:
        m.GetParameter(10, fEces, err);
        fEces /= kRad2Deg;
    case 10:
        m.GetParameter(9, fCry, err);
        fCry /= kRad2Deg;
    case 9:
        m.GetParameter(8, fCrx, err);
        fCrx /= kRad2Deg;
    case 8:
        m.GetParameter(7, fNry, err);
        fNry /= kRad2Deg;
    case 7:
        m.GetParameter(6, fNrx, err);
        fNrx /= kRad2Deg;
    case 6:
        m.GetParameter(5, fAw, err);
        fAw /= kRad2Deg;
    case 5:
        m.GetParameter(4, fAn, err);
        fAn /= kRad2Deg;
    case 4:
        m.GetParameter(3, fCa, err);
        fCa /= kRad2Deg;
    case 3:
        m.GetParameter(2, fNpae, err);
        fNpae /= kRad2Deg;
    case 2:
        m.GetParameter(1, fIe, err);
        fIe /= kRad2Deg;
    case 1:
        m.GetParameter(0, fIa, err);
        fIa /= kRad2Deg;
    }
}
/*
void FormatPar(TMinuit &m, Int_t n)
{
    Double_t par, err;
    m.GetParameter(n, par, err);

    int expp = (int)log10(par);
    int expe = (int)log10(err);

    if (err<2*pow(10, expe))
        expe--;

    Int_t exp = expe>expp ? expp : expe;

    par = (int)(par/pow(10, exp)) * pow(10, exp);
    err = (int)(err/pow(10, exp)) * pow(10, exp);

    cout << par << " +- " << err << flush;
}
*/
void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const
{
    if (n<0)
        n = m.GetNumPars();

    cout << setprecision(3);

    Double_t par, err;

    switch (n)
    {
    case 16:
        m.GetParameter(15, par, err);
        cout << " 15 MAGIC2: " << setw(8) << par << " +- " << setw(4) <<  err << endl;
    case 15:
        m.GetParameter(14, par, err);
        cout << " 14 MAGIC1: " << setw(8) << par << " +- " << setw(4) <<  err << endl;
    case 14:
        m.GetParameter(13, par, err);
        cout << " 13 ACEC:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Azimuth Centering Error (cos)" << endl;
    case 13:
        m.GetParameter(12, par, err);
        cout << " 12 ECEC:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Elevation Centering Error (cos)" << endl;
    case 12:
        m.GetParameter(11, par, err);
        cout << " 11 ACES:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Azimuth Centering Error (sin)" << endl;
    case 11:
        m.GetParameter(10, par, err);
        cout << " 10 ECES:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Elevation Centering Error (sin)" << endl;
    case 10:
        m.GetParameter(9, par, err);
        cout << "  9 CRY:    " << setw(8) << par << " +- " << setw(4) <<  err << "  Alt/Az Coude Displacement (E-W)" << endl;
    case 9:
        m.GetParameter(8, par, err);
        cout << "  8 CRX:    " << setw(8) << par << " +- " << setw(4) <<  err << "  Alt/Az Coude Displacement (N-S)" << endl;
    case 8:
        m.GetParameter(7, par, err);
        cout << "  7 NRY:    " << setw(8) << par << " +- " << setw(4) <<  err << "  Nasmyth rotator displacement, vertical" << endl;
    case 7:
        m.GetParameter(6, par, err);
        cout << "  6 NRX:    " << setw(8) << par << " +- " << setw(4) <<  err << "  Nasmyth rotator displacement, horizontan" << endl;
    case 6:
        m.GetParameter(5, par, err);
        cout << "  5 AW:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Azimuth Axis Misalignment (E-W)" << endl;
    case 5:
        m.GetParameter(4, par, err);
        cout << "  4 AN:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Azimuth Axis Misalignment (N-S)" << endl;
    case 4:
        m.GetParameter(3, par, err);
        cout << "  3 CA:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Left-Right Collimation Error" << endl;
    case 3:
        m.GetParameter(2, par, err);
        cout << "  2 NPAE:   " << setw(8) << par << " +- " << setw(4) <<  err << "  Az-El Nonperpendicularity" << endl;
    case 2:
        m.GetParameter(1, par, err);
        cout << "  1 IE:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Index Error Elevation (Offset)" << endl;
    case 1:
        m.GetParameter(0, par, err);
        cout << "  0 IA:     " << setw(8) << par << " +- " << setw(4) <<  err << "  Index Error Azimuth (Offset)" << endl;
    }
}