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

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