source: trunk/MagicSoft/Mars/mastro/MAstroCatalog.cc@ 3440

Last change on this file since 3440 was 3415, checked in by tonello, 21 years ago
*** empty log message ***
File size: 11.1 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 03/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2002-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MAstroCatalog
28//
29// THIS IMPLEMENTATION IS PRELIMINARY AND WILL BE MERGED WITH
30// SOME PARTS OF THE DRIVE SOFTWARE SOON!
31//
32//////////////////////////////////////////////////////////////////////////////
33#include "MAstroCatalog.h"
34
35#include <fstream>
36
37#include <TLine.h>
38#include <TMarker.h>
39#include <TArrayI.h>
40#include <TRotation.h>
41#include <TStopwatch.h>
42#include <TVirtualPad.h>
43
44#include <stdlib.h>
45
46#include "MLog.h"
47#include "MLogManip.h"
48
49#include "MAstro.h"
50#include "MTime.h"
51#include "MObservatory.h"
52
53ClassImp(MVector3);
54ClassImp(MAstroCatalog);
55
56using namespace std;
57
58MVector3 MVector3::GetZdAz(const MObservatory &obs, Double_t gmst) const
59{
60 if (!fType==kIsRaDec)
61 return MVector3();
62
63 const Double_t alpha = gmst + obs.GetElong();
64
65 MVector3 zdaz;
66 zdaz.SetZdAz(Theta(), alpha-Phi(), Mag());
67 zdaz.RotateY(obs.GetPhi()-TMath::Pi()/2);
68 return zdaz;
69
70 /*
71 // ------ The same using slalib, tested in the drive system -------
72 const Double_t alpha = slaGmst(mjd) + obs.GetElong();
73 Double_t el;
74 slaDe2h(fAlpha-ra, dec, obs.GetPhi(), &az, &el);
75 zd = TMath::Pi()/2-el;
76 return;
77 */
78}
79
80MVector3 MVector3::GetZdAz(const MTime &time, MObservatory &obs) const
81{
82 return GetZdAz(obs, time.GetGmst());
83}
84
85TString MAstroCatalog::FindToken(TString &line, Char_t tok)
86{
87 Ssiz_t token = line.First(tok);
88 if (token<0)
89 {
90 const TString copy(line);
91 line = "";
92 return copy;
93 }
94
95 const TString res = line(0, token);
96 line.Remove(0, token+1);
97 return res;
98}
99
100Int_t MAstroCatalog::atoi(const TSubString &sub)
101{
102 return atoi(TString(sub));
103}
104
105Float_t MAstroCatalog::atof(const TSubString &sub)
106{
107 return atof(TString(sub));
108}
109
110Int_t MAstroCatalog::atoi(const TString &s)
111{
112 return std::atoi(s);
113}
114
115Float_t MAstroCatalog::atof(const TString &s)
116{
117 return std::atof(s);
118}
119
120Int_t MAstroCatalog::ReadXephem(TString catalog)
121{
122 gLog << inf << "Reading Xephem catalog: " << catalog << endl;
123
124 ifstream fin(catalog);
125
126 Int_t add =0;
127 Int_t cnt =0;
128 Int_t line=0;
129
130 Double_t maxmag=0;
131
132 while (1)
133 {
134 TString row;
135 row.ReadLine(fin);
136 if (!fin)
137 break;
138
139 line++;
140
141 if (row[0]=='#')
142 continue;
143
144 TString line(row);
145
146 TString name = FindToken(line);
147 TString dummy = FindToken(line);
148 TString r = FindToken(line);
149 TString d = FindToken(line);
150 TString m = FindToken(line);
151 TString epoch = FindToken(line);
152
153 if (name.IsNull() || r.IsNull() || d.IsNull() || m.IsNull() || epoch.IsNull())
154 {
155 gLog << warn << "Invalid Entry Line #" << line << ": " << row << endl;
156 continue;
157 }
158
159 cnt++;
160
161 const Double_t mag = atof(m);
162
163 maxmag = TMath::Max(maxmag, mag);
164
165 if (mag>fLimMag)
166 continue;
167
168 if (epoch!="2000")
169 {
170 gLog << warn << "Epoch != 2000... skipped." << endl;
171 continue;
172 }
173
174 Double_t ra0, dec0;
175 MAstro::Coordinate2Angle(r, ra0);
176 MAstro::Coordinate2Angle(d, dec0);
177
178 ra0 *= TMath::Pi()/12;
179 dec0 *= TMath::Pi()/180;
180
181 MVector3 *star=new MVector3;
182 star->SetRaDec(ra0, dec0, mag);
183 if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV)
184 {
185 delete star;
186 continue;
187 }
188
189 fList.Add(star);
190 add++;
191 }
192 gLog << inf << "Read " << add << " out of " << cnt << " (Total max mag=" << maxmag << ")" << endl;
193
194 return add;
195}
196
197Int_t MAstroCatalog::ReadNGC2000(TString catalog)
198{
199 gLog << inf << "Reading NGC2000 catalog: " << catalog << endl;
200
201 ifstream fin(catalog);
202
203 Int_t add=0;
204 Int_t cnt=0;
205 Int_t n =0;
206
207 Double_t maxmag=0;
208
209 while (1)
210 {
211 TString row;
212 row.ReadLine(fin);
213 if (!fin)
214 break;
215
216 cnt++;
217
218 const Int_t rah = atoi(row(13, 2));
219 const Float_t ram = atof(row(16, 4));
220 const Char_t decs = row(22);
221 const Int_t decd = atoi(row(23, 2));
222 const Int_t decm = atoi(row(26, 2));
223 const TString m = row(43, 4);
224
225 if (m.Strip().IsNull())
226 continue;
227
228 n++;
229
230 const Double_t mag = atof(m);
231
232 maxmag = TMath::Max(maxmag, mag);
233
234 if (mag>fLimMag)
235 continue;
236
237 const Double_t ra = MAstro::Hms2Rad(rah, (int)ram, fmod(ram, 1)*60);
238 const Double_t dec = MAstro::Dms2Rad(decd, decm, 0, decs);
239
240 MVector3 *star=new MVector3;
241 star->SetRaDec(ra, dec, mag);
242 if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV)
243 {
244 delete star;
245 continue;
246 }
247
248 fList.Add(star);
249 add++;
250 }
251
252 gLog << inf << "Read " << add << " out of " << n << " (Total max mag=" << maxmag << ")" << endl;
253
254 return add;
255}
256
257Int_t MAstroCatalog::ReadBSC(TString catalog)
258{
259 gLog << inf << "Reading Bright Star Catalog (BSC5) catalog: " << catalog << endl;
260
261 ifstream fin(catalog);
262
263 Int_t add=0;
264 Int_t cnt=0;
265 Int_t n =0;
266
267 Double_t maxmag=0;
268
269 while (1)
270 {
271 TString row;
272 row.ReadLine(fin);
273 if (!fin)
274 break;
275
276 cnt++;
277
278 const Int_t rah = atoi(row(75, 2));
279 const Int_t ram = atoi(row(77, 2));
280 const Float_t ras = atof(row(79, 4));
281 const Char_t decsgn = row(83);
282 const Int_t decd = atoi(row(84, 2));
283 const Int_t decm = atoi(row(86, 2));
284 const Int_t decs = atoi(row(88, 2));
285 const TString m = row(102, 5);
286
287 if (m.Strip().IsNull())
288 continue;
289
290 n++;
291
292 const Double_t mag = atof(m.Data());
293
294 maxmag = TMath::Max(maxmag, mag);
295
296 if (mag>fLimMag)
297 continue;
298
299 const Double_t ra = MAstro::Hms2Rad(rah, ram, ras);
300 const Double_t dec = MAstro::Dms2Rad(decd, decm, decs, decsgn);
301
302 MVector3 *star=new MVector3;
303 star->SetRaDec(ra, dec, mag);
304 if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV)
305 {
306 delete star;
307 continue;
308 }
309
310 fList.Add(star);
311 add++;
312 }
313
314 gLog << inf << "Read " << add << " out of " << n << " (Total max mag=" << maxmag << ")" << endl;
315
316 return add;
317}
318
319Bool_t MAstroCatalog::Convert(const TRotation &rot, TVector2 &v, Int_t type)
320{
321 MVector3 w;
322 w.SetRaDec(v.Y(), v.X()-fRaDec.Phi(), v.Y());
323 w *= rot;
324
325 v.Set(w.Phi(), w.Theta());
326
327 return w.Angle(TVector3(1, 0, 0))*TMath::RadToDeg()<=fRadiusFOV;
328}
329
330Bool_t MAstroCatalog::PaintLine(const TVector2 &v, Double_t dx, Double_t dy, const TRotation &rot, Int_t type)
331{
332 const TVector2 add(dx*TMath::DegToRad(), dy*TMath::DegToRad());
333
334 TVector2 v0 = v;
335 TVector2 v1 = v+add;
336
337 const Bool_t rc0 = Convert(rot, v0, type);
338 const Bool_t rc1 = Convert(rot, v1, type);
339 if (!rc0 && !rc1)
340 return kFALSE;
341
342 TLine line;
343 line.SetLineColor(kGreen);
344 line.PaintLine(v0.X(), TMath::Pi()/2-v0.Y(),
345 v1.X(), TMath::Pi()/2-v1.Y());
346
347 return kTRUE;
348}
349
350void MAstroCatalog::Paint(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Byte_t type)
351{
352 Int_t idx[] = {1, 1, 1, 1};
353
354 Int_t dirs[4][2] = { {0, 1}, {1, 0}, {0, -1}, {-1, 0} };
355
356 for (int i=0; i<dx.GetSize(); i++)
357 {
358 for (int j=0; j<4; j++)
359 if (dx[i]==dx[0]+dirs[j][0] && dy[i]==dy[0]+dirs[j][1])
360 idx[j] = 0;
361 }
362
363 TVector2 v1 = v0 + TVector2(dx[0]*TMath::DegToRad(), dy[0]*TMath::DegToRad());
364
365 for (int i=0; i<4; i++)
366 if (idx[i])
367 {
368 dx.Set(dx.GetSize()+1);
369 dy.Set(dy.GetSize()+1);
370
371 dx[dx.GetSize()-1] = dx[0];
372 dy[dy.GetSize()-1] = dy[0];
373
374 dx[0] += dirs[i][0];
375 dy[0] += dirs[i][1];
376
377 if (PaintLine(v1, dirs[i][0], dirs[i][1], rot, type))
378 Paint(v0, rot, dx, dy, type);
379
380 dx[0] -= dirs[i][0];
381 dy[0] -= dirs[i][1];
382 }
383}
384
385void MAstroCatalog::PaintNet(const TVector2 &v0, const TRotation &rot, Int_t type)
386{
387 //const Double_t step = TMath::DegToRad();
388
389 TArrayI dx(1);
390 TArrayI dy(1);
391
392 TVector2 v = v0*TMath::RadToDeg();
393 v.Set(TMath::Floor(v.X()), TMath::Floor(v.Y()));
394 v *= TMath::DegToRad();
395
396 Paint(v, rot, dx, dy, type);
397}
398
399void MAstroCatalog::Paint(Option_t *o)
400{
401 Double_t ra = fRaDec.Phi();
402 Double_t dec = TMath::Pi()/2-fRaDec.Theta();
403
404 TIter Next(&fList);
405 TVector3 *v;
406
407 Double_t minra=360, maxra=0, mindec=360, maxdec=0;
408
409 while ((v=(TVector3*)Next()))
410 {
411 minra = TMath::Min(minra, v->Phi());
412 maxra = TMath::Max(maxra, v->Phi());
413
414 mindec = TMath::Min(mindec, TMath::Pi()/2-v->Theta());
415 maxdec = TMath::Max(maxdec, TMath::Pi()/2-v->Theta());
416 }
417
418 cout << gPad << endl;
419
420 cout << "Minra: " << (minra-ra)*TMath::RadToDeg() << endl;
421 cout << "Maxra: " << (maxra-ra)*TMath::RadToDeg() << endl;
422
423 cout << "Mindec: " << (mindec-dec)*TMath::RadToDeg() << endl;
424 cout << "Maxdec: " << (maxdec-dec)*TMath::RadToDeg() << endl;
425
426 //gPad->Range(minra-ra, mindec-dec, maxra-ra, maxdec-dec);
427 gPad->Range(-fRadiusFOV*TMath::DegToRad()/TMath::Sqrt(2.), -fRadiusFOV*TMath::DegToRad()/TMath::Sqrt(2.),
428 fRadiusFOV*TMath::DegToRad()/TMath::Sqrt(2.), fRadiusFOV*TMath::DegToRad()/TMath::Sqrt(2.));
429
430 Next.Reset();
431
432 cout << "Dec: " << dec*TMath::RadToDeg() << endl;
433 cout << "Ra: " << ra*TMath::RadToDeg() << endl;
434
435 // Precalc Sin/Cos...
436 TRotation trans;
437 trans.Rotate(dec, TVector3(0, 1, 0));
438
439 TStopwatch clk;
440 clk.Start();
441
442 TMarker mark;
443 mark.SetMarkerColor(kBlack);
444 mark.SetMarkerStyle(kCircle);
445 while ((v=(TVector3*)Next()))
446 {
447 MVector3 v0;
448 v0.SetRaDec(v->Phi()-ra, TMath::Pi()/2-v->Theta(), 1);
449 v0 *= trans;
450
451 mark.SetMarkerSize((fLimMag-TMath::Log(v->Mag()))/4);
452 mark.PaintMarker(v0.Phi(), TMath::Pi()/2-v0.Theta());
453 }
454
455 TMarker m;
456 m.SetMarkerStyle(kCross);
457 m.PaintMarker(0, 0);
458
459 m.SetMarkerColor(kRed);
460 m.SetMarkerStyle(kFullDotSmall);
461
462 TVector2 v0(ra, dec);
463 PaintNet(v0, trans);
464
465 clk.Stop();
466 clk.Print();
467}
Note: See TracBrowser for help on using the repository browser.