source: trunk/Cosy/main/MMoonPointing.cc@ 19333

Last change on this file since 19333 was 10025, checked in by tbretz, 14 years ago
Increased the valid range for cos^2+sin^2 in MMoonPointing
File size: 5.0 KB
Line 
1#include "MMoonPointing.h"
2
3#include <TMath.h>
4#include <TFile.h>
5
6#include "MLog.h"
7#include "MLogManip.h"
8
9ClassImp(MMoonPointing);
10
11using namespace std;
12
13class MPositionOffsetCalc
14{
15private:
16 Double_t fCosOffset;
17 Double_t fSinOffset;
18
19 Double_t fCosAngle;
20 Double_t fSinAngle;
21
22public:
23 MPositionOffsetCalc(Double_t offset) :
24 fCosOffset(offset==0 ? 0 : TMath::Cos(offset)),
25 fSinOffset(offset==0 ? 0 : TMath::Sin(offset)),
26 fCosAngle(0), fSinAngle(0)
27 {
28 }
29
30 MPositionOffsetCalc(Double_t offset, Double_t angle) :
31 fCosOffset(offset==0 ? 0 : TMath::Cos(offset)),
32 fSinOffset(offset==0 ? 0 : TMath::Sin(offset)),
33 fCosAngle(TMath::Cos(angle*TMath::DegToRad())),
34 fSinAngle(TMath::Sin(angle*TMath::DegToRad()))
35 {
36 }
37
38 MPositionOffsetCalc(Double_t offset, Double_t cos, Double_t sin) :
39 fCosOffset(offset==0 ? 0 : TMath::Cos(offset)),
40 fSinOffset(offset==0 ? 0 : TMath::Sin(offset)),
41 fCosAngle(cos), fSinAngle(sin)
42 {
43 }
44
45 void SetAngle(Double_t c, Double_t s) { fCosAngle=c; fSinAngle=s; }
46 void SetAngle(Double_t a) { fCosAngle=TMath::Cos(a*TMath::DegToRad()); fSinAngle=TMath::Sin(a*TMath::DegToRad()); }
47
48 ZdAz Calc(const ZdAz &pos, Double_t *sn=0, Double_t *cs=0) const
49 {
50 // This is to avoid unnecessary calculations
51 // in case the offset is zero
52 if (fCosOffset==0 && fSinOffset==0)
53 {
54 if (sn)
55 *sn = -fCosAngle;
56 if (cs)
57 *cs = fSinAngle;
58 return pos;
59 }
60
61 const Double_t coszd = TMath::Cos(pos.Zd());
62 const Double_t sinzd = TMath::Sin(pos.Zd());
63
64 if (sinzd<=0)
65 {
66 gLog << warn << "WARNING - |Zd|>=90°" << endl;
67 return ZdAz();
68 }
69
70 const Double_t costheta = coszd*fCosOffset + sinzd*fSinOffset*fCosAngle;
71 if (costheta >= 1)
72 {
73 gLog << warn << "WARNING - cos(Zd) > 1." << endl;
74 return ZdAz();
75 }
76
77 const Double_t sintheta = TMath::Sqrt(1 - costheta*costheta);
78
79 const Double_t cosdeltaaz = (fCosOffset - coszd*costheta)/(sinzd*sintheta);
80 const Double_t sindeltaaz = fSinAngle*fSinOffset/sintheta;
81
82 if (sn)
83 *sn = fSinAngle*coszd*sindeltaaz - cosdeltaaz*fCosAngle;
84 if (cs)
85 *cs = fSinAngle*sinzd/sintheta;
86
87 return ZdAz(TMath::ACos(costheta), TMath::ATan2(sindeltaaz, cosdeltaaz) + pos.Az());
88 }
89
90 ZdAz GetOffset(const ZdAz &pos, Double_t *sn=0, Double_t *cs=0) const
91 {
92 return Calc(pos, sn, cs) - pos;
93 }
94};
95
96MMoonPointing::MMoonPointing(const char *filename) : fOffsetShadow(-1), fOffsetWobble(-1)
97{
98 if (TString(filename).IsNull())
99 {
100 MakeZombie();
101 return;
102 }
103
104 TFile f(filename);
105 if (f.IsZombie())
106 {
107 MakeZombie();
108 gLog << warn << "WARNING - Reading of " << filename << " failed." << endl;
109 return;
110 }
111
112 if (fCos.Read("cosAngle")<=0)
113 {
114 MakeZombie();
115 gLog << warn << "WARNING - Reading of TGraph2D 'cosAngle' from " << filename << " failed." << endl;
116 return;
117 }
118
119 if (fSin.Read("sinAngle")<=0)
120 {
121 MakeZombie();
122 gLog << warn << "WARNING - Reading of TGraph2D 'sinAngle' from " << filename << " failed." << endl;
123 return;
124 }
125
126 fCos.SetDirectory(0);
127 fSin.SetDirectory(0);
128}
129
130// position of the moon [rad]
131// position of the Moon Shadow [rad]
132// position of telescope pointing [rad]
133Bool_t MMoonPointing::CalcPosition(const ZdAz &moon, ZdAz &srcpos, ZdAz &pointpos) const
134{
135 const Double_t az = fmod(moon.Az()+TMath::TwoPi(), TMath::TwoPi());
136
137 const Double_t fcos = const_cast<TGraph2D&>(fCos).Interpolate(az, moon.Zd());
138 const Double_t fsin = const_cast<TGraph2D&>(fSin).Interpolate(az, moon.Zd());
139
140 // This is a sanity check for the case the Interpolation failed
141 const Double_t norm = TMath::Hypot(fcos, fsin);
142 if (norm<0.996 || norm>1.001)
143 {
144 gLog << warn << "WARNING - Local moon position Zd/Az=";
145 gLog << moon.Zd()*TMath::RadToDeg() << "/" << moon.Az()*TMath::RadToDeg();
146 gLog << " out of range of table for interpolation." << endl;
147 gLog << " fcos=" << fcos << endl;
148 gLog << " fsin=" << fsin << endl;
149 gLog << " norm=" << norm << endl;
150 return kFALSE;
151 }
152
153 // Apply a small numerical improvement
154 Double_t sn=fsin/norm;
155 Double_t cs=fcos/norm;
156
157 // Calculate moon shadow position
158 const MPositionOffsetCalc calc1(fOffsetShadow, fcos, fsin);
159 srcpos = calc1.Calc(moon, &sn, &cs);
160 if (srcpos.Zd()==0)
161 {
162 gLog << warn << "WARNING - Calculation of moon shadow position failed." << endl;
163 return kFALSE;
164 }
165
166 // Calaculate Telescope pointing position:
167 const MPositionOffsetCalc calc2(fOffsetWobble, cs, sn);
168 pointpos = calc2.Calc(srcpos);
169 if (pointpos.Zd()==0)
170 {
171 gLog << warn << "WARNING - Calculation of moon shadow's wobble position failed." << endl;
172 return kFALSE;
173 }
174
175 return kTRUE;
176}
177
178
Note: See TracBrowser for help on using the repository browser.