source: trunk/MagicSoft/Cosy/main/MMoonPointing.cc @ 9551

Last change on this file since 9551 was 9551, checked in by tbretz, 11 years ago
*** empty log message ***
File size: 3.7 KB
Line 
1#include "MMoonPointing.h"
2
3#include <TMath.h>
4#include <TFile.h>
5
6ClassImp(MMoonPointing);
7
8using namespace std;
9
10class MPositionOffsetCalc
11{
12private:
13    Double_t fCosOffset;
14    Double_t fSinOffset;
15
16    Double_t fCosAngle;
17    Double_t fSinAngle;
18
19public:
20    MPositionOffsetCalc(Double_t offset) :
21        fCosOffset(offset<=0?0:TMath::Cos(offset*TMath::DegToRad())),
22        fSinOffset(offset<=0?0:TMath::Sin(offset*TMath::DegToRad())),
23        fCosAngle(0), fSinAngle(0)
24    {
25    }
26
27    MPositionOffsetCalc(Double_t offset, Double_t angle) :
28        fCosOffset(offset<=0?0:TMath::Cos(offset*TMath::DegToRad())),
29        fSinOffset(offset<=0?0:TMath::Sin(offset*TMath::DegToRad())),
30        fCosAngle(TMath::Cos(angle*TMath::DegToRad())),
31        fSinAngle(TMath::Sin(angle*TMath::DegToRad()))
32    {
33    }
34
35    MPositionOffsetCalc(Double_t offset, Double_t cos, Double_t sin) :
36        fCosOffset(offset<=0?0:TMath::Cos(offset*TMath::DegToRad())),
37        fSinOffset(offset<=0?0:TMath::Sin(offset*TMath::DegToRad())),
38        fCosAngle(cos), fSinAngle(sin)
39    {
40    }
41
42    void SetAngle(Double_t c, Double_t s) { fCosAngle=c; fSinAngle=s; }
43    void SetAngle(Double_t a) { fCosAngle=TMath::Cos(a*TMath::DegToRad()); fSinAngle=TMath::Sin(a*TMath::DegToRad()); }
44
45    ZdAz Calc(const ZdAz &pos, Double_t *sn=0, Double_t *cs=0) const
46    {
47        if (fCosOffset==0 && fSinOffset==0)
48        {
49            if (sn)
50                *sn = fSinAngle;
51            if (cs)
52                *cs = fCosAngle;
53            return pos;
54        }
55
56        const Double_t coszd = TMath::Cos(pos.Zd());
57        const Double_t sinzd = TMath::Sin(pos.Zd());
58
59        if (sinzd<=0)
60            return ZdAz();
61
62        const Double_t costheta = coszd*fCosOffset + sinzd*fSinOffset*fCosAngle;
63        if (costheta*costheta >= 1)
64            return ZdAz();
65
66        const Double_t sintheta = TMath::Sqrt(1 - costheta*costheta);
67
68        const Double_t cosdeltaaz = (fCosOffset - coszd*costheta)/(sinzd*sintheta);
69        const Double_t sindeltaaz = fSinAngle*fSinOffset/sintheta;
70
71        if (sn)
72            *sn = fSinAngle*coszd*sindeltaaz - cosdeltaaz*fCosAngle;
73        if (cs)
74            *cs = fSinAngle*sinzd/sintheta;
75
76        return ZdAz(TMath::ACos(costheta), TMath::ATan2(sindeltaaz, cosdeltaaz) + pos.Az());
77    }
78
79    ZdAz GetOffset(const ZdAz &pos, Double_t *sn=0, Double_t *cs=0) const
80    {
81        return Calc(pos, sn, cs) - pos;
82    }
83};
84
85MMoonPointing::MMoonPointing(const char *filename) : fOffsetShadow(-1), fOffsetWobble(-1)
86{
87    if (TString(filename).IsNull())
88    {
89        MakeZombie();
90        return;
91    }
92
93    TFile f(filename);
94    if (f.IsZombie())
95    {
96        MakeZombie();
97        return;
98    }
99
100    if (fCos.Read("cosAngle")<=0)
101    {
102        MakeZombie();
103        return;
104    }
105
106    if (fSin.Read("sinAngle")<=0)
107    {
108        MakeZombie();
109        return;
110    }
111
112    fCos.SetDirectory(0);
113    fSin.SetDirectory(0);
114}
115
116// position of the moon [rad]
117// position of the Moon Shadow [rad]
118// position of telescope pointing [rad]
119Bool_t MMoonPointing::CalcPosition(const ZdAz &moon, ZdAz &srcpos, ZdAz &pointpos) const
120{
121    const Double_t fcos = const_cast<TGraph2D&>(fCos).Interpolate(moon.Az(), moon.Zd());
122    const Double_t fsin = const_cast<TGraph2D&>(fSin).Interpolate(moon.Az(), moon.Zd());
123    if (fcos*fcos + fsin*fsin<0.98)
124        return kFALSE;
125
126    // Calculate Moon Shadow Position:
127    Double_t sn=fsin;
128    Double_t cs=fcos;
129
130    const MPositionOffsetCalc calc1(fOffsetShadow, fcos, fsin);
131    srcpos = calc1.Calc(moon, &sn, &cs);
132    if (srcpos.Zd()==0)
133        return kFALSE;
134
135    // Calaculate Telescope pointing position:
136    const MPositionOffsetCalc calc2(fOffsetWobble, cs, sn);
137    pointpos = calc2.Calc(srcpos);
138    if (pointpos.Zd()==0)
139        return kFALSE;
140
141    return kTRUE;
142}
143
144
Note: See TracBrowser for help on using the repository browser.