source: trunk/MagicSoft/Cosy/main/MBending.cc@ 1714

Last change on this file since 1714 was 1699, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 9.3 KB
Line 
1#include "MBending.h"
2
3#include <fstream.h>
4
5#include <TMinuit.h>
6
7
8ClassImp(MBending);
9
10void MBending::Reset()
11{
12 fIa = 0;
13 fIe = 0;
14 fCa = 0;
15 fAn = 0;
16 fAw = 0;
17 fNrx = 0;
18 fNry = 0;
19 fCrx = 0;
20 fCry = 0;
21 fNpae = 0;
22 fEces = 0;
23 fAces = 0;
24 fEcec = 0;
25 fAcec = 0;
26}
27
28void MBending::Load(const char *name)
29{
30 /*
31 ! MMT 1987 July 8
32 ! T 36 7.3622 41.448 -0.0481
33 ! IA -37.5465 20.80602
34 ! IE -13.9180 1.25217
35 ! NPAE +7.0751 26.44763
36 ! CA -6.9149 32.05358
37 ! AN +0.5053 1.40956
38 ! AW -2.2016 1.37480
39 ! END
40 */
41
42 ifstream fin(name);
43 if (!fin)
44 {
45 cout << "Error: Cannot open file '" << name << "'" << endl;
46 return;
47 }
48
49 char c;
50 while (fin && fin.get()!='\n');
51 fin >> c;
52
53 if (c!='S' && c!='s')
54 {
55 cout << "Error: This in not a model correcting the star position (" << c << ")" << endl;
56 return;
57 }
58
59 Clear();
60
61 cout << endl;
62
63 Double_t val;
64 fin >> val;
65 cout << "Number of observed stars: " << val << endl;
66 fin >> val;
67 cout << "Sky RMS: " << val << "\"" << endl;
68 fin >> val;
69 cout << "Refraction Constant A: " << val << "\"" << endl;
70 fin >> val;
71 cout << "Refraction Constant B: " << val << "\"" << endl;
72
73 cout << endl;
74
75 cout << " & = Name Value Sigma" << endl;
76 cout << "------------------------------------------" << endl;
77
78 while (fin)
79 {
80 TString str;
81 fin >> str;
82
83 if (str=="END")
84 break;
85
86 if (str[0]=='&')
87 {
88 cout << " & ";
89 str.Remove(0);
90 }
91 else
92 cout << " ";
93
94 if (str[1]=='=')
95 {
96 cout << "= ";
97 str.Remove(0);
98 }
99 else
100 cout << " ";
101
102 fin >> val;
103 cout << str << "\t" << val << "° \t";
104 val *= kDeg2Rad;
105
106 if (str=="IA") fIa = val;
107 if (str=="IE") fIe = val;
108 if (str=="CA") fCa = val;
109 if (str=="AN") fAn = val;
110 if (str=="AW") fAw = val;
111 if (str=="NRX") fNrx = val;
112 if (str=="NRY") fNry = val;
113 if (str=="CRX") fCrx = val;
114 if (str=="CRY") fCry = val;
115 if (str=="NPAE") fNpae = val;
116 if (str=="ECES") fEces = val;
117 if (str=="ACES") fAces = val;
118 if (str=="ECEC") fEcec = val;
119 if (str=="ACEC") fAcec = val;
120
121 fin >> val;
122 cout << val*kDeg2Rad << endl;
123 }
124 cout << endl;
125}
126
127void MBending::Save(const char *name) {}
128
129AltAz MBending::Correct(const AltAz &aa) const
130{
131 // Correct [rad]
132 // zdaz [rad]
133 AltAz p = aa;
134
135 const AltAz CES(-fEces*kDeg2Rad*sin(p.Alt()), -fAces*kDeg2Rad*sin(p.Az()));
136 const AltAz CEC(-fEcec*kDeg2Rad*cos(p.Alt()), -fAcec*kDeg2Rad*cos(p.Az()));
137 p += CES;
138 p += CEC;
139
140 const AltAz CRX(-fCrx*kDeg2Rad*sin(p.Az()-p.Alt()), fCrx*kDeg2Rad*cos(p.Az()-p.Alt())/cos(p.Alt()));
141 const AltAz CRY(-fCry*kDeg2Rad*cos(p.Az()-p.Alt()), -fCry*kDeg2Rad*sin(p.Az()-p.Alt())/cos(p.Alt()));
142 p += CRX;
143 p += CRY;
144
145 const AltAz NRX(fNrx*kDeg2Rad*sin(p.Alt()), -fNrx*kDeg2Rad);
146 const AltAz NRY(fNry*kDeg2Rad*cos(p.Alt()), -fNry*kDeg2Rad*tan(p.Alt()));
147 p += NRX;
148 p += NRY;
149
150 const AltAz AW( fAw*kDeg2Rad*sin(p.Az()), -fAw*kDeg2Rad*cos(p.Az())*tan(p.Alt()));
151 const AltAz AN(-fAn*kDeg2Rad*cos(p.Az()), -fAn*kDeg2Rad*sin(p.Az())*tan(p.Alt()));
152 p += AW;
153 p += AN;
154
155 const AltAz CA(0, -fCa*kDeg2Rad/cos(p.Alt()));
156 p += CA;
157
158 const AltAz NPAE(0, -fNpae*kDeg2Rad*tan(p.Alt()));
159 p += NPAE;
160
161 const AltAz I(fIe/**kDeg2Rad*/, fIa/**kDeg2Rad*/);
162 p += I;
163
164 return p;
165}
166
167AltAz MBending::CorrectBack(const AltAz &aa) const
168{
169 // Correct [rad]
170 // zdaz [rad]
171 AltAz p = aa;
172
173 const AltAz CES(-fEces*kDeg2Rad*sin(p.Alt()), -fAces*kDeg2Rad*sin(p.Az()));
174 const AltAz CEC(-fEcec*kDeg2Rad*cos(p.Alt()), -fAcec*kDeg2Rad*cos(p.Az()));
175 p -= CES;
176 p -= CEC;
177
178 const AltAz CRX(-fCrx*kDeg2Rad*sin(p.Az()-p.Alt()), fCrx*kDeg2Rad*cos(p.Az()-p.Alt())/cos(p.Alt()));
179 const AltAz CRY(-fCry*kDeg2Rad*cos(p.Az()-p.Alt()), -fCry*kDeg2Rad*sin(p.Az()-p.Alt())/cos(p.Alt()));
180 p -= CRX;
181 p -= CRY;
182
183 const AltAz NRX(fNrx*kDeg2Rad*sin(p.Alt()), -fNrx*kDeg2Rad);
184 const AltAz NRY(fNry*kDeg2Rad*cos(p.Alt()), -fNry*kDeg2Rad*tan(p.Alt()));
185 p -= NRX;
186 p -= NRY;
187
188 const AltAz AW( fAw*kDeg2Rad*sin(p.Az()), -fAw*kDeg2Rad*cos(p.Az())*tan(p.Alt()));
189 const AltAz AN(-fAn*kDeg2Rad*cos(p.Az()), -fAn*kDeg2Rad*sin(p.Az())*tan(p.Alt()));
190 p -= AW;
191 p -= AN;
192
193 const AltAz CA(0, -fCa*kDeg2Rad/cos(p.Alt()));
194 p -= CA;
195
196 const AltAz NPAE(0, -fNpae*kDeg2Rad*tan(p.Alt()));
197 p -= NPAE;
198
199 const AltAz I(fIe/**kDeg2Rad*/, fIa/**kDeg2Rad*/);
200 p -= I;
201
202 return p;
203}
204
205ZdAz MBending::Correct(const ZdAz &zdaz) const
206{
207 // Correct [rad]
208 // zdaz [rad]
209 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
210 AltAz c = Correct(p);
211 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
212}
213
214ZdAz MBending::CorrectBack(const ZdAz &zdaz) const
215{
216 // Correct [rad]
217 // zdaz [rad]
218 AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
219 AltAz c = CorrectBack(p);
220 return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
221}
222
223void MBending::SetParameters(const Double_t *par, Int_t n=14)
224{
225 Clear();
226
227 switch (n)
228 {
229 case 14:
230 fAcec=par[13] ; // Azimuth Centering Error (cos)
231 case 13:
232 fEcec=par[12] ; // Elevation Centering Error (cos)
233 case 12:
234 fAces=par[11] ; // Azimuth Centering Error (sin)
235 case 11:
236 fEces=par[10] ; // Elevation Centering Error (sin)
237 case 10:
238 fCry =par[9] ; // Alt/Az Coude Displacement (E-W)
239 case 9:
240 fCrx =par[8] ; // Alt/Az Coude Displacement (N-S)
241 case 8:
242 fNry =par[7] ; // Nasmyth rotator displacement, vertical
243 case 7:
244 fNrx =par[6] ; // Nasmyth rotator displacement, horizontan
245 case 6:
246 fAw =par[5] ; // Azimuth Axis Misalignment (E-W)
247 case 5:
248 fAn =par[4] ; // Azimuth Axis Misalignment (N-S)
249 case 4:
250 fCa =par[3] ; // Left-Right Collimation Error
251 case 3:
252 fNpae=par[2] ; // Az-El Nonperpendicularity
253 case 2:
254 fIe =par[1] ; // Index Error in Elevation
255 case 1:
256 fIa =par[0] ; // Index Error in Azimuth
257 }
258}
259
260void MBending::GetParameters(Double_t *par, Int_t n=14) const
261{
262 switch (n)
263 {
264 case 14:
265 par[13]=fAcec; // Azimuth Centering Error (cos)
266 case 13:
267 par[12]=fEcec; // Elevation Centering Error (cos)
268 case 12:
269 par[11]=fAces ; // Azimuth Centering Error (sin)
270 case 11:
271 par[10]=fEces ; // Elevation Centering Error (sin)
272 case 10:
273 par[9]=fCry ; // Alt/Az Coude Displacement (E-W)
274 case 9:
275 par[8]=fCrx ; // Alt/Az Coude Displacement (N-S)
276 case 8:
277 par[7]=fNry ; // Nasmyth rotator displacement, vertical
278 case 7:
279 par[6]=fNrx ; // Nasmyth rotator displacement, horizontan
280 case 6:
281 par[5]=fAw ; // Azimuth Axis Misalignment (E-W)
282 case 5:
283 par[4]=fAn ; // Azimuth Axis Misalignment (N-S)
284 case 4:
285 par[3]=fCa ; // Left-Right Collimation Error
286 case 3:
287 par[2]=fNpae ; // Az-El Nonperpendicularity
288 case 2:
289 par[1]=fIe ; // Index Error in Elevation
290 case 1:
291 par[0]=fIa ; // Index Error in Azimuth
292 }
293}
294
295void MBending::SetMinuitParameters(TMinuit &m, Int_t n=-1) const
296{
297 if (n<0)
298 n = m.GetNumPars();
299
300 Int_t ierflg = 0;
301
302 switch (n)
303 {
304 case 14:
305 case 13:
306 case 12:
307 case 11:
308 case 10:
309 case 9:
310 case 8:
311 case 7:
312 case 5:
313 case 4:
314 m.mnparm(3, "CA", fCa, 1, -360, 360, ierflg);
315 // cout << "Init 3 CA: " << fCa << endl;
316 case 3:
317 m.mnparm(2, "NPAE", fNpae, 1, -360, 360, ierflg);
318 // cout << "Init 2 NPAE: " << fNpae << endl;
319 case 2:
320 m.mnparm(1, "IE", fIe, 1, -360, 360, ierflg);
321 // cout << "Init 1 IE: " << fIe << endl;
322 case 1:
323 m.mnparm(0, "IA", fIa, 1, -360, 360, ierflg);
324 // cout << "Init 0 IA: " << fIa << endl;
325 }
326}
327
328void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1)
329{
330 if (n<0)
331 n = m.GetNumPars();
332
333 Double_t err;
334
335 switch (n)
336 {
337 case 14:
338 case 13:
339 case 12:
340 case 11:
341 case 10:
342 case 9:
343 case 8:
344 case 7:
345 case 6:
346 case 5:
347 case 4:
348 m.GetParameter(3, fCa, err);
349 case 3:
350 m.GetParameter(2, fNpae, err);
351 case 2:
352 m.GetParameter(1, fIe, err);
353 case 1:
354 m.GetParameter(0, fIa, err);
355 }
356}
357
358void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const
359{
360 if (n<0)
361 n = m.GetNumPars();
362
363 Double_t par, err;
364
365 switch (n)
366 {
367 case 14:
368 case 13:
369 case 12:
370 case 11:
371 case 10:
372 case 9:
373 case 8:
374 case 7:
375 case 6:
376 case 5:
377 case 4:
378 m.GetParameter(3, par, err);
379 cout << " 3 CA: " << par << " +- " << err << endl;
380 case 3:
381 m.GetParameter(2, par, err);
382 cout << " 2 NPAE: " << par << " +- " << err << endl;
383 case 2:
384 m.GetParameter(1, par, err);
385 cout << " 1 IE: " << par << " +- " << err << endl;
386 case 1:
387 m.GetParameter(0, par, err);
388 cout << " 0 IA: " << par << " +- " << err << endl;
389 }
390}
Note: See TracBrowser for help on using the repository browser.