Index: trunk/MagicSoft/Cosy/Changelog
===================================================================
--- trunk/MagicSoft/Cosy/Changelog	(revision 1392)
+++ trunk/MagicSoft/Cosy/Changelog	(revision 1393)
@@ -1,3 +1,13 @@
                                                                   -*-*- END -*-*-
+
+ 2002/07/10 - Thomas Bretz:
+
+   * base/coord.h:
+     - added some member functions
+
+   * main/MBending.[h,cc]:
+     - added
+
+
 
  2002/06/03 - Thomas Bretz:
Index: trunk/MagicSoft/Cosy/base/coord.h
===================================================================
--- trunk/MagicSoft/Cosy/base/coord.h	(revision 1392)
+++ trunk/MagicSoft/Cosy/base/coord.h	(revision 1393)
@@ -3,4 +3,5 @@
 
 #include <math.h>          // floor
+#include <fstream.h>
 
 /* pi/180:  degrees to radians */
@@ -34,4 +35,7 @@
 class XY
 {
+    friend ifstream& operator>>(ifstream &in, XY &xy);
+    friend ofstream& operator<<(ofstream &in, XY &xy);
+
 protected:
     double fX;
@@ -46,4 +50,7 @@
     double X() const { return fX; }
     double Y() const { return fY; }
+
+    void X(double x) { fX=x; }
+    void Y(double y) { fY=y; }
 
     void operator/=(double c) { fX/=c; fY/=c; }
@@ -64,4 +71,7 @@
 };
 
+inline ifstream& operator>>(ifstream &in,  XY &xy) { in  >> xy.fX; in  >> xy.fY; return in; }
+inline ofstream& operator<<(ofstream &out, XY &xy) { out << xy.fX << " " << xy.fY; return out; }
+
 class AltAz : public XY
 {
@@ -72,8 +82,12 @@
     double Az()  const { return fY; }
 
+    void operator*=(double c) { fX*=c; fY*=c; }
+    void operator/=(double c) { fX*=c; fY*=c; }
+
     void Alt(double d) { fX=d; }
     void Az(double d)  { fY=d; }
     void operator*=(const XY &c)    { fX*=c.X(); fY*=c.Y(); }
     void operator-=(const AltAz &c) { fX-=c.fX; fY-=c.fY; }
+    void operator+=(const AltAz &c) { fX+=c.fX; fY+=c.fY; }
 
     AltAz operator/(double c) const { return AltAz(fX/c, fY/c); }
@@ -91,4 +105,7 @@
     ZdAz(double zd=0, double az=0) : XY(zd, az) {}
     ZdAz(const ZdAz &c) : XY(c) {}
+
+    void operator*=(double c) { fX*=c; fY*=c; }
+    void operator/=(double c) { fX*=c; fY*=c; }
 
     double Zd() const { return fX; }
@@ -118,4 +135,7 @@
     double Dec() const { return fY; }
 
+    void operator*=(double c) { fX*=c; fY*=c; }
+    void operator/=(double c) { fX*=c; fY*=c; }
+
     void Ra(double x)  { fX = x; }
     void Dec(double y) { fY = y; }
Index: trunk/MagicSoft/Cosy/main/MBending.cc
===================================================================
--- trunk/MagicSoft/Cosy/main/MBending.cc	(revision 1393)
+++ trunk/MagicSoft/Cosy/main/MBending.cc	(revision 1393)
@@ -0,0 +1,325 @@
+#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;
+    }
+}
Index: trunk/MagicSoft/Cosy/main/MBending.h
===================================================================
--- trunk/MagicSoft/Cosy/main/MBending.h	(revision 1393)
+++ trunk/MagicSoft/Cosy/main/MBending.h	(revision 1393)
@@ -0,0 +1,106 @@
+#ifndef COSY_MBending
+#define COSY_MBending
+
+#include <TROOT.h>
+
+#include <fstream.h>
+
+#include "coord.h"
+
+class TMinuit;
+
+class MBending
+{
+private:
+    Double_t fIe   ; // [rad] Index Error in Elevation
+    Double_t fIa   ; // [rad] Index Error in Azimuth
+    Double_t fCa   ; // [rad] Left-Right Collimation Error
+    Double_t fAn   ; // [rad] Azimuth Axis Misalignment (N-S)
+    Double_t fAw   ; // [rad] Azimuth Axis Misalignment (E-W)
+    Double_t fNpae ; // [rad] Az-El Nonperpendicularity
+    Double_t fNrx  ; // [rad] Nasmyth rotator displacement, horizontan
+    Double_t fNry  ; // [rad] Nasmyth rotator displacement, vertical
+    Double_t fCrx  ; // [rad] Alt/Az Coude Displacement (N-S)
+    Double_t fCry  ; // [rad] Alt/Az Coude Displacement (E-W)
+    Double_t fEces ; // [rad] Elevation Centering Error (sin)
+    Double_t fAces ; // [rad] Azimuth Centering Error (sin)
+    Double_t fEcec ; // [rad] Elevation Centering Error (cos)
+    Double_t fAcec ; // [rad] Azimuth Centering Error (cos)
+
+    void Clear()
+    {
+        fIe  =0 ; // Index Error in Elevation
+        fIa  =0 ; // Index Error in Azimuth
+        fCa  =0 ; // Left-Right Collimation Error
+        fAn  =0 ; // Azimuth Axis Misalignment (N-S)
+        fAw  =0 ; // Azimuth Axis Misalignment (E-W)
+        fNpae=0 ; // Az-El Nonperpendicularity
+        fNrx =0 ; // Nasmyth rotator displacement, horizontan
+        fNry =0 ; // Nasmyth rotator displacement, vertical
+        fCrx =0 ; // Alt/Az Coude Displacement (N-S)
+        fCry =0 ; // Alt/Az Coude Displacement (E-W)
+        fEces=0 ; // Elevation Centering Error (sin)
+        fAces=0 ; // Azimuth Centering Error (sin)
+        fEcec=0 ; // Elevation Centering Error (cos)
+        fAcec=0 ; // Azimuth Centering Error (cos)
+    }
+
+public:
+    MBending() { Clear(); }
+
+    void Load(const char *name);
+    void Save(const char *name);
+
+    ZdAz  Correct(const ZdAz &zdaz) const;
+    AltAz Correct(const AltAz &aaz) const;
+
+    ZdAz operator()(const ZdAz &zdaz) const { return Correct(zdaz); }
+    AltAz operator()(const AltAz &aaz) const { return Correct(aaz); }
+
+    ZdAz operator()(const ZdAz &zdaz, void (*fcn)(ZdAz &zdaz, Double_t *par)) const
+    {
+        Double_t par[14];
+        GetParameters(par);
+        ZdAz za = zdaz;
+        fcn(za, par);
+        return za;
+    }
+
+    AltAz operator()(const AltAz &aaz, void (*fcn)(AltAz &aaz, Double_t *par)) const
+    {
+        Double_t par[14];
+        GetParameters(par);
+        AltAz aa = aaz;
+        fcn(aa, par);
+        return aa;
+    }
+                                 
+
+    void SetParameters(const Double_t *par, Int_t n=14);
+    void GetParameters(Double_t *par, Int_t n=14) const;
+
+    void SetMinuitParameters(TMinuit &m, Int_t n=-1) const;
+    void GetMinuitParameters(TMinuit &m, Int_t n=-1);
+    void PrintMinuitParameters(TMinuit &m, Int_t n=-1) const;
+
+    /*
+     Double_t GetIe() const { return fIe; }
+     Double_t GetIa() const { return fIa; }
+     Double_t GetCa() const { return fCa; }
+     Double_t GetAn() const { return fAn; }
+     Double_t GetAw() const { return fAw; }
+     Double_t GetNrx() const { return fNrx; }
+     Double_t GetNry() const { return fNry; }
+     Double_t GetCrx() const { return fNrx; }
+     Double_t GetCry() const { return fNry; }
+     Double_t GetEces() const { return fEces; }
+     Double_t GetEcec() const { return fEcec; }
+     Double_t GetAces() const { return fAces; }
+     Double_t GetAcec() const { return fAcec; }
+     Double_t GetNpae() const { return fNpae; }
+     */
+
+    ClassDef(MBending, 1)
+};
+
+#endif
Index: trunk/MagicSoft/Cosy/main/MainLinkDef.h
===================================================================
--- trunk/MagicSoft/Cosy/main/MainLinkDef.h	(revision 1392)
+++ trunk/MagicSoft/Cosy/main/MainLinkDef.h	(revision 1393)
@@ -7,3 +7,5 @@
 #pragma link C++ class MCosy;
 
+#pragma link C++ class MBending;
+
 #endif
Index: trunk/MagicSoft/Cosy/main/Makefile
===================================================================
--- trunk/MagicSoft/Cosy/main/Makefile	(revision 1392)
+++ trunk/MagicSoft/Cosy/main/Makefile	(revision 1393)
@@ -33,4 +33,5 @@
 
 SRCFILES = MCosy.cc \
+           MBending.cc
            MStarguider.cc