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

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