source: trunk/Mars/mstarcam/MStarCamTrans.cc@ 19593

Last change on this file since 19593 was 4966, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 35.8 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): Wolfgang Wittek, 08/2004 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MStarCamTrans
28// ---------------
29//
30// this is a collection of transformations between
31//
32// a) celestial (declination, hourangle) = ( dec, h); units (deg,hour)
33// b) local (zenith angle, azimuth angle) = (theta, phi); units (deg, deg)
34// c) camera (x-coordinate, y-coordinate) = ( X, Y); units ( mm, mm)
35//
36// coordinates. As to the definition of the coordinates, the conventions
37// of TDAS 00-11 are used.
38//
39// The transformations use
40// - cos(Lat) and sin(Lat) from an MObservatory container,
41// where 'Lat' is the geographical latitude of the telescope.
42// - fDistCam from an MGeomCam container,
43// which is the distance of the camera from the reflector center
44//
45// in order to apply one of the transformations the
46// 'MGeomCam' and 'MObservatory' containers have to be present;
47// the corresponding objects 'mgeom' and 'mobserv' have to be used as
48// arguments of the constructor
49//
50// MStarCamTrans(mgeom, mobserv);
51//
52// The tranformations ignore effects like nutation, precession, refraction, ...
53// This is not a real problem as long as relative positions are considered, as
54// in all the transformations, except in CelToLoc and LocToCel.
55//
56// The camera coordinates (X, Y) are assumed to be the coordinates for an ideal
57// imaging, without imaging errors. The x-axis is assumed to be horizonthal,
58// the y-axis is pointing upwards.
59//
60// the azimuthal angle is defined as : 0 degrees = north
61// 90 degrees = east
62// 180 degrees = south
63// 270 degrees = west
64//
65// Examples for the use of the transformations :
66//
67// - if the local coordinates (theta1, phi1) of a point (X1, Y1) in the camera
68// are known one may caculate the local coordinates (theta2, phi2) of any
69// other point (X2, Y2) in the camera :
70//
71// LocCamCamToLoc(theta1, phi1, X1, Y1, X2, Y2, theta2, phi2);
72//
73// if (X1, Y1) = (0, 0) one may use
74// Loc0CamToLoc(theta1, phi1, X2, Y2, theta2, phi2);
75//
76
77// - if the local coordinates (theta1, phi1) of a point (X1, Y1) in the camera
78// are known one may caculate the position (X2, Y2) for another direction
79// (theta2, phi2) :
80//
81// LocCamLocToCam(theta1, phi1, X1, Y1, theta2, phi2, X2, Y2);
82//
83// if (X1, Y1) = (0, 0) one may use
84// Loc0LocToCam(theta1, phi1, X2, Y2, theta2, phi2);
85//
86// - if the celestial coordinates (dec1, h1) of a point (X1, Y1) in the camera
87// are known one may caculate the celestial coordinates (dec2, h2) of any
88// other point (X2, Y2) in the camera :
89//
90// CelCamCamToCel(dec1, h1, X1, Y1, X2, Y2, dec2, h2);
91//
92// if (X1, Y1) = (0, 0) one may use
93// Cel0CamToCel(dec1, h1, X2, Y2, dec2, h2);
94//
95//
96// - if the celestial coordinates (dec1, h1) of a point (X1, Y1) in the camera
97// are known one may caculate the position (X2, Y2) for any other direction :
98//
99// CelCamCelToCam(dec1, h1, X1, Y1, dec2, h2, X2, Y2);
100//
101// if (X1, Y1) = (0, 0) one may use
102// Cel0CelToCam(dec1, h1, dec2, h2, X2, Y2);
103//
104//
105//
106////////////////////////////////////////////////////////////////////////////
107#include "MStarCamTrans.h"
108
109#include <math.h>
110
111#include <TH1.h>
112#include <TH2.h>
113#include <TPad.h>
114#include <TStyle.h>
115#include <TGraph.h>
116#include <TLatex.h>
117#include <TCanvas.h>
118
119//#include "MH.h"
120#include "MLog.h"
121#include "MLogManip.h"
122#include "MGeomCam.h"
123#include "MObservatory.h"
124
125using namespace std;
126
127ClassImp(MStarCamTrans);
128
129
130// --------------------------------------------------------------------------
131//
132// Constructor
133//
134// get distance between camera and reflector from an MGeomCam container
135// get cos(Lat) and sin(Lat) from an MObservatory container
136//
137MStarCamTrans::MStarCamTrans(const MGeomCam &cam, const MObservatory &obs)
138{
139 if (&cam == NULL)
140 {
141 gLog << err << "MStarCamTrans::MStarCamTrans; MGeomCam container is not available"
142 << endl;
143 return;
144 }
145
146 if (&obs == NULL)
147 {
148 gLog << err << "MStarCamTrans::MStarCamTrans; MObservatory container is not available"
149 << endl;
150 return;
151 }
152
153 fDistCam = cam.GetCameraDist() * 1000.0; // [mm]
154 fCLat = obs.GetCosPhi();
155 fSLat = obs.GetSinPhi();
156
157 if (fDistCam*fCLat*fSLat == 0.0)
158 {
159 gLog << "MStarCamTrans::TransCelLocCam; one of 'fDistCam, fCLat, fSLat' is zero !!!$$$$$$$$$$$$$$$$$ : "
160 << fDistCam << ", " << fCLat << ", " << fSLat << endl;
161 }
162
163
164 fGridBinning = 0.50; // degrees
165 fGridFineBin = 0.01; // degrees
166
167}
168
169// --------------------------------------------------------------------------
170//
171// Loc0CamToLoc
172//
173// Input : (theta0, phi0) direction for the position (0,0) in the camera
174// ( X, Y) a position in the camera
175//
176// Output : ( theta, phi) direction for the position (X,Y) in the camera
177//
178Bool_t MStarCamTrans::Loc0CamToLoc(Double_t theta0deg, Double_t phi0deg,
179 Double_t X, Double_t Y,
180 Double_t &thetadeg, Double_t &phideg)
181{
182 Double_t theta0 = theta0deg / kRad2Deg;
183 Double_t phi0 = phi0deg / kRad2Deg;
184
185 Double_t XC = X / fDistCam;
186 Double_t YC = Y / fDistCam;
187
188 Double_t theta;
189 Double_t phiminphi0;
190
191 //--------------------------------------------
192 Double_t sip = -XC;
193 Double_t cop = sin(theta0) + YC*cos(theta0);
194
195 Double_t sit = sqrt(cop*cop + XC*XC);
196 Double_t cot = cos(theta0) - YC*sin(theta0);
197
198 // there is an ambiguity in the sign of cos(theta)
199 // - if theta0 is close to 0 or Pi
200 // choose that theta which is closer to theta0,
201 // i.e. require the same sign for cos(theta) and cos(theta0)
202 // - otherwise choose that theta which is compatible with a small
203 // difference |phi-phi0|, i.e. cos(phi-phi0) > 0
204
205 if( fabs(theta0deg - 90.0) > 45.0 )
206 theta = TMath::ATan2( sit, TMath::Sign(cot, cos(theta0)) );
207 else
208 theta = TMath::ATan2( sit, TMath::Sign(cot, cop/cot) );
209
210 Double_t sig = TMath::Sign(1.0, cot*tan(theta));
211 phiminphi0 = TMath::ATan2(sig*sip, sig*cop);
212
213 //--------------------------------------------
214 phideg = (phi0 + phiminphi0) * kRad2Deg;
215 thetadeg = theta * kRad2Deg;
216
217 //gLog << "MStarCamTrans::Loc0Cam2Log; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
218 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
219 // << thetadeg << ", " << phideg << endl;
220
221 return kTRUE;
222}
223
224// --------------------------------------------------------------------------
225//
226// Loc0LocToCam
227//
228// Input : (theta0, phi0) direction for the position (0,0) in the camera
229// ( theta, phi) some other direction
230//
231// Output : (X, Y) position in the camera corresponding to (theta, phi)
232//
233// (see also MAstro::GetDistOnPlain())
234//
235Bool_t MStarCamTrans::Loc0LocToCam(Double_t theta0deg, Double_t phi0deg,
236 Double_t thetadeg, Double_t phideg,
237 Double_t &X, Double_t &Y)
238{
239 Double_t theta0 = theta0deg / kRad2Deg;
240 Double_t phi0 = phi0deg / kRad2Deg;
241
242 Double_t theta = thetadeg / kRad2Deg;
243 Double_t phi = phideg / kRad2Deg;
244
245 //--------------------------------------------
246 Double_t YC = (cos(theta0)*tan(theta)*cos(phi-phi0) - sin(theta0))
247 / (cos(theta0) + sin(theta0)*tan(theta));
248 Double_t XC = -sin(phi-phi0) * (cos(theta0) - YC*sin(theta0)) * tan(theta);
249
250 //--------------------------------------------
251 X = XC * fDistCam;
252 Y = YC * fDistCam;
253
254
255 //gLog << "MStarCamTrans::Loc0LocToCam; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
256 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
257 // << thetadeg << ", " << phideg << endl;
258
259 return kTRUE;
260}
261
262// --------------------------------------------------------------------------
263//
264// LocCamToLoc0
265//
266// Input : ( X, Y) a position in the camera
267// (theta, phi) direction for the position (X,Y) in the camera
268//
269// Output : ( theta0, phi0) direction for the position (0,0) in the camera
270//
271Bool_t MStarCamTrans::LocCamToLoc0(Double_t thetadeg, Double_t phideg,
272 Double_t X, Double_t Y,
273 Double_t &theta0deg, Double_t &phi0deg)
274{
275 Double_t theta = thetadeg / kRad2Deg;
276 Double_t phi = phideg / kRad2Deg;
277
278 Double_t XC = X / fDistCam;
279 Double_t YC = Y / fDistCam;
280
281 //--------------------------------------------
282
283 Double_t theta0;
284 Double_t phiminphi0;
285
286 Double_t tant1 = 0.0;
287 Double_t tant2 = 0.0;
288
289 Double_t ip;
290
291 // calculate tan(theta0)
292 // dummy loop to avoid a 'go to'
293 for (Int_t i=0; i<1; i++)
294 {
295 if ( tan(theta) == 0.0 )
296 {
297 if (XC != 0.0)
298 {
299 // this can never occur
300 gLog << "MStarCamTrans::LocCam2Loc0; thetadeg, XC = " << thetadeg
301 << ", " << XC << endl;
302 return kFALSE;
303 }
304 tant1 = -YC;
305
306 theta0 = TMath::Pi() *
307 modf( (atan(tant1)+TMath::Pi()) / TMath::Pi(), &ip );
308 phiminphi0 = 0.0;
309
310 break;
311 }
312
313 Double_t a = 1.0 + XC*XC - YC*YC*tan(theta)*tan(theta);
314 Double_t b = 2.0 * YC * (1.0 + tan(theta)*tan(theta));
315 Double_t c = XC*XC + YC*YC - tan(theta)*tan(theta);
316
317 if (a == 0.0)
318 {
319 if (b == 0.0)
320 {
321 // this can never occur
322 gLog << "MStarCamTrans::LocCam2Loc0; a, b = " << a << ", "
323 << b << endl;
324 return kFALSE;
325 }
326 tant1 = -c/b;
327 }
328 else
329 {
330 Double_t discr = b*b - 4.0*a*c;
331 if (discr < 0.0)
332 {
333 gLog << "MStarCamTrans::LocCam2Loc0; discr = " << discr << endl;
334 return kFALSE;
335 }
336
337 // two solutions for tan(theta0)
338 tant1 = (-b + sqrt(discr)) / (2.0*a);
339 tant2 = (-b - sqrt(discr)) / (2.0*a);
340
341 if (fabs(tant1-tant2) < 1.e-5)
342 tant2 = 0.0;
343 }
344
345 // define the sign of tan(theta0) and
346 // reject the solution with the wrong sign
347
348 if ( fabs(thetadeg - 90.0) > 45.0 )
349 {
350 Double_t sig = TMath::Sign(1.0, cos(theta));
351 if ( tant1 != 0.0 && TMath::Sign(1.0, tant1) != sig )
352 tant1 = 0.0;
353
354 if ( tant2 != 0.0 && TMath::Sign(1.0, tant2) != sig )
355 tant2 = 0.0;
356 }
357 else
358 {
359 // require sign of cos(phi-phi0) to be > 0
360 if ( tant1 != 0.0 &&
361 TMath::Sign( 1.0, (tant1+YC)/ ((1.0-YC*tant1)*tan(theta)) ) != 1.0 )
362 tant1 = 0.0;
363
364 if ( tant2 != 0.0 &&
365 TMath::Sign( 1.0, (tant2+YC)/ ((1.0-YC*tant2)*tan(theta)) ) != 1.0 )
366 tant2 = 0.0;
367 }
368
369 if (tant1 != 0.0 && tant2 != 0.0)
370 {
371 gLog << "MStarCamTrans::LocCam2Loc0; there are 2 solutions for tan(theta0), tant1, tant2 = "
372 << tant1 << ", " << tant2 << endl;
373 }
374
375 if (tant1 != 0.0)
376 theta0 = TMath::Pi() *
377 modf( (atan(tant1)+TMath::Pi()) / TMath::Pi(), &ip );
378 else if (tant2 != 0.0)
379 theta0 = TMath::Pi() *
380 modf( (atan(tant2)+TMath::Pi()) / TMath::Pi(), &ip );
381 else
382 {
383 theta0 = theta;
384 gLog << "MStarCamTrans::LocCam2Loc0; there is no solution for tan(theta0)"
385 << endl;
386 }
387
388 Double_t sip = -XC;
389 Double_t cop = sin(theta0) + YC*cos(theta0);
390 Double_t cot = cos(theta0) - YC*sin(theta0);
391 Double_t sig = TMath::Sign(1.0, cot*tan(theta));
392
393 phiminphi0 = TMath::ATan2(sig*sip, sig*cop);
394 } // end of dummy loop
395
396 //--------------------------------------------
397 phi0deg = (phi - phiminphi0) * kRad2Deg;
398 theta0deg = theta0 * kRad2Deg;
399
400 //gLog << "MStarCamTrans::LocCamToLoc0; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
401 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
402 // << thetadeg << ", " << phideg << endl;
403
404
405 return kTRUE;
406}
407
408// --------------------------------------------------------------------------
409//
410// Cel0CamToCel
411//
412// Input : (dec0, h0) direction for the position (0,0) in the camera
413// ( X, Y) a position in the camera
414//
415// Output : (dec, h) direction for the position (X,Y) in the camera
416//
417Bool_t MStarCamTrans::Cel0CamToCel(Double_t dec0deg, Double_t h0hour,
418 Double_t X, Double_t Y,
419 Double_t &decdeg, Double_t &hhour)
420{
421 //--------------------------------------------
422
423 // transform celestial coordinates ( dec0, h0)
424 // into local coordinates (theta0, phi0)
425 Double_t theta0deg;
426 Double_t phi0deg;
427 CelToLoc(dec0deg, h0hour, theta0deg, phi0deg);
428
429 // get the local coordinates (theta, phi)
430 // for the position (X, Y) in the camera
431 Double_t thetadeg;
432 Double_t phideg;
433 Loc0CamToLoc(theta0deg, phi0deg, X, Y, thetadeg, phideg);
434
435 // transform local coordinates (theta, phi)
436 // into celestial coordinates ( dec, h)
437 LocToCel(thetadeg, phideg, decdeg, hhour);
438
439 //--------------------------------------------
440
441
442 //gLog << "MStarCamTrans::Cel0CamToCel; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
443 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
444 // << thetadeg << ", " << phideg << endl;
445
446 return kTRUE;
447}
448
449// --------------------------------------------------------------------------
450//
451// Cel0CelToCam
452//
453// Input : (dec0, h0) direction for the position (0,0) in the camera
454// (dec, h) some other direction
455//
456// Output : (X, Y) position in the camera corresponding to (dec, h)
457//
458Bool_t MStarCamTrans::Cel0CelToCam(Double_t dec0deg, Double_t h0hour,
459 Double_t decdeg, Double_t hhour,
460 Double_t &X, Double_t &Y)
461{
462 //--------------------------------------------
463
464 // transform celestial coordinates ( dec0, h0)
465 // into local coordinates (theta0, phi0)
466 Double_t theta0deg;
467 Double_t phi0deg;
468 CelToLoc(dec0deg, h0hour, theta0deg, phi0deg);
469
470 // transform celestial coordinates ( dec, h)
471 // into local coordinates (theta, phi)
472 Double_t thetadeg;
473 Double_t phideg;
474 CelToLoc(decdeg, hhour, thetadeg, phideg);
475
476 // get the position (X, Y) in the camera
477 // from the local coordinates (theta, phi)
478 Loc0LocToCam(theta0deg, phi0deg, thetadeg, phideg, X, Y);
479
480 //--------------------------------------------
481
482 //gLog << "MStarCamTrans::Cel0CelToCam; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
483 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
484 // << thetadeg << ", " << phideg << endl;
485
486 return kTRUE;
487}
488
489// --------------------------------------------------------------------------
490//
491// CelCamToCel0
492//
493// Input : ( X, Y) a position in the camera
494// (dec, h) direction for the position (X,Y) in the camera
495//
496// Output : (dec0, h0) direction for the position (0,0) in the camera
497//
498
499Bool_t MStarCamTrans::CelCamToCel0(Double_t decdeg, Double_t hhour,
500 Double_t X, Double_t Y,
501 Double_t &dec0deg, Double_t &h0hour)
502{
503 //--------------------------------------------
504
505 // transform celestial coordinates ( dec, h)
506 // into local coordinates (theta, phi)
507 Double_t thetadeg;
508 Double_t phideg;
509 CelToLoc(decdeg, hhour, thetadeg, phideg);
510
511 // get the local coordinates (theta, phi)
512 // for the position (0, 0) in the camera
513 Double_t theta0deg;
514 Double_t phi0deg;
515 LocCamToLoc0(thetadeg, phideg, X, Y, theta0deg, phi0deg);
516
517 // transform local coordinates (theta0, phi0)
518 // into celestial coordinates ( dec0, h0)
519 LocToCel(theta0deg, phi0deg, dec0deg, h0hour);
520
521 //--------------------------------------------
522
523
524 //gLog << "MStarCamTrans::CelCamToCel0; theta0deg, phi0deg, X, Y, thetadeg, phideg = "
525 // << theta0deg << ", " << phi0deg << ", " << X << ", " << Y << ", "
526 // << thetadeg << ", " << phideg << endl;
527
528 return kTRUE;
529}
530
531// --------------------------------------------------------------------------
532//
533// LocCamCamToLoc
534//
535// Input : (theta1, phi1) direction for the position (X1,Y1) in the camera
536// ( X1, Y1) a position in the camera
537// ( X2, Y2) another position in the camera
538//
539// Output : (theta2, phi2) direction for the position (X2,Y2) in the camera
540//
541Bool_t MStarCamTrans::LocCamCamToLoc(Double_t theta1deg, Double_t phi1deg,
542 Double_t X1, Double_t Y1,
543 Double_t X2, Double_t Y2,
544 Double_t &theta2deg, Double_t &phi2deg)
545{
546 if (X1 == 0.0 && Y1 == 0.0)
547 {
548 Loc0CamToLoc(theta1deg, phi1deg, X2, Y2, theta2deg, phi2deg);
549 return kTRUE;
550 }
551
552 if (X2 == 0.0 && Y2 == 0.0)
553 {
554 LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta2deg, phi2deg);
555 return kTRUE;
556 }
557
558 Double_t theta0deg;
559 Double_t phi0deg;
560 LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta0deg, phi0deg);
561
562 Loc0CamToLoc(theta0deg, phi0deg, X2, Y2, theta2deg, phi2deg);
563
564 return kTRUE;
565}
566
567// --------------------------------------------------------------------------
568//
569// LocCamLocToCam
570//
571// Input : (theta1, phi1) direction for the position (X1,Y1) in the camera
572// ( X1, Y1) a position in the camera
573// (theta2, phi2) another direction
574//
575// Output : ( X2, Y2) position corresponding to (theta2, phi2)
576//
577Bool_t MStarCamTrans::LocCamLocToCam(Double_t theta1deg, Double_t phi1deg,
578 Double_t X1, Double_t Y1,
579 Double_t theta2deg, Double_t phi2deg,
580 Double_t &X2, Double_t &Y2)
581
582{
583 if (X1 == 0.0 && Y1 == 0.0)
584 {
585 Loc0LocToCam(theta1deg, phi1deg, theta2deg, phi2deg, X2, Y2);
586
587 return kTRUE;
588 }
589
590
591 Double_t theta0deg;
592 Double_t phi0deg;
593 LocCamToLoc0(theta1deg, phi1deg, X1, Y1, theta0deg, phi0deg);
594
595 Loc0LocToCam(theta0deg, phi0deg, theta2deg, phi2deg, X2, Y2);
596
597 return kTRUE;
598}
599
600// --------------------------------------------------------------------------
601//
602// CelCamCamToCel
603//
604// Input : (dec1, h1) direction for the position (X1,Y1) in the camera
605// ( X1, Y1) a position in the camera
606// ( X2, Y2) another position in the camera
607//
608// Output : (dec2, h2) direction for the position (X2,Y2) in the camera
609//
610Bool_t MStarCamTrans::CelCamCamToCel(Double_t dec1deg, Double_t h1deg,
611 Double_t X1, Double_t Y1,
612 Double_t X2, Double_t Y2,
613 Double_t &dec2deg, Double_t &h2deg)
614{
615 if (X1 == 0.0 && Y1 == 0.0)
616 {
617 Cel0CamToCel(dec1deg, h1deg, X2, Y2, dec2deg, h2deg);
618 return kTRUE;
619 }
620
621 if (X2 == 0.0 && Y2 == 0.0)
622 {
623 CelCamToCel0(dec1deg, h1deg, X1, Y1, dec2deg, h2deg);
624 return kTRUE;
625 }
626
627 Double_t dec0deg;
628 Double_t h0deg;
629 CelCamToCel0(dec1deg, h1deg, X1, Y1, dec0deg, h0deg);
630
631 Cel0CamToCel(dec0deg, h0deg, X2, Y2, dec2deg, h2deg);
632
633 return kTRUE;
634}
635
636// --------------------------------------------------------------------------
637//
638// CelCamCelToCam
639//
640// Input : (dec1, h1) direction for the position (X1,Y1) in the camera
641// ( X1, Y1) a position in the camera
642// (dec2, h2) another direction
643//
644// Output : ( X2, Y2) position corresponding to (dec2, h2)
645//
646Bool_t MStarCamTrans::CelCamCelToCam(Double_t dec1deg, Double_t h1deg,
647 Double_t X1, Double_t Y1,
648 Double_t dec2deg, Double_t h2deg,
649 Double_t &X2, Double_t &Y2)
650
651{
652 if (X1 == 0.0 && Y1 == 0.0)
653 {
654 Cel0CelToCam(dec1deg, h1deg, dec2deg, h2deg, X2, Y2);
655
656 return kTRUE;
657 }
658
659
660 Double_t dec0deg;
661 Double_t h0deg;
662 CelCamToCel0(dec1deg, h1deg, X1, Y1, dec0deg, h0deg);
663
664 Cel0CelToCam(dec0deg, h0deg, dec2deg, h2deg, X2, Y2);
665
666 return kTRUE;
667}
668
669// --------------------------------------------------------------------------
670//
671// CelToLoc
672//
673// Input : (dec, h) celestial coordinates
674//
675// Output : (theta, phi) corresponding local coordinates
676//
677// (see also MAstroCatalog and MAstroSky2Local)
678//
679Bool_t MStarCamTrans::CelToLoc(Double_t decdeg, Double_t hhour,
680 Double_t &thetadeg, Double_t &phideg)
681{
682 Double_t a1 = fCLat;
683 Double_t a3 = -fSLat;
684
685 Double_t dec = decdeg / kRad2Deg;
686 Double_t h = hhour / 24.0 * TMath::TwoPi();
687
688 // transform celestial coordinates ( dec, h)
689 // into local coordinates (theta, phi)
690 Double_t xB = cos(dec) * cos(h);
691 Double_t yB = cos(dec) * sin(h);
692 Double_t zB = - sin(dec);
693
694 Double_t xA = a3*xB - a1*zB;
695 Double_t yA = - yB;
696 Double_t zA = -a1*xB - a3*zB;
697
698 thetadeg = acos(-zA) * kRad2Deg;
699 phideg = TMath::ATan2(yA, xA) * kRad2Deg;
700
701 return kTRUE;
702}
703
704// --------------------------------------------------------------------------
705//
706// LocToCel
707//
708// Input : (theta, phi) local coordinates
709//
710// Output : (dec, h) corresponding celestial coordinates
711//
712// (see also MAstroCatalog and MAstroSky2Local)
713//
714Bool_t MStarCamTrans::LocToCel(Double_t thetadeg, Double_t phideg,
715 Double_t &decdeg, Double_t &hhour)
716{
717 Double_t a1 = fCLat;
718 Double_t a3 = -fSLat;
719
720 Double_t theta = thetadeg / kRad2Deg;
721 Double_t phi = phideg / kRad2Deg;
722
723 //--------------------------------------------
724
725 // transform local coordinates (theta, phi)
726 // into celestial coordinates ( dec, h)
727 Double_t xA = sin(theta) * cos(phi);
728 Double_t yA = sin(theta) * sin(phi);
729 Double_t zA = - cos(theta);
730
731 Double_t xB = a3*xA - a1*zA;
732 Double_t yB = - yA;
733 Double_t zB = -a1*xA - a3*zA;
734
735 Double_t dec = asin(-zB);
736 Double_t h = TMath::ATan2(yB, xB);
737
738 //--------------------------------------------
739 decdeg = dec * kRad2Deg;
740 hhour = h * 24.0 / TMath::TwoPi();
741
742
743 return kTRUE;
744}
745
746// --------------------------------------------------------------------------
747//
748// PlotGridAtDec0H0
749//
750// set the celestial coordinates of the camera center
751// and plot the lines of constant (Theta, Phi)
752// and the lines of constant (Dec, H )
753//
754// (see also MAstroCatalog::Draw and MAstroCamera::Draw)
755//
756// Warning: Leaks Memory!
757//
758Bool_t MStarCamTrans::PlotGridAtDec0H0(TString name,
759 Double_t dec0deg, Double_t h0hour)
760{
761 Double_t theta0deg;
762 Double_t phi0deg;
763 CelToLoc(dec0deg, h0hour, theta0deg, phi0deg);
764
765 PlotGrid(name, theta0deg, phi0deg, dec0deg, h0hour);
766
767 return kTRUE;
768}
769
770// --------------------------------------------------------------------------
771//
772// PlotGridAtTheta0Phi0
773//
774// set the local coordinates of the camera center
775// and plot the lines of constant (Theta, Phi)
776// and the lines of constant (Dec, H )
777//
778// (see also MAstroCatalog::Draw and MAstroCamera::Draw)
779//
780// Warning: Leaks Memory!
781//
782Bool_t MStarCamTrans::PlotGridAtTheta0Phi0(TString name,
783 Double_t theta0deg, Double_t phi0deg)
784{
785 Double_t dec0deg;
786 Double_t h0hour;
787 LocToCel(theta0deg, phi0deg, dec0deg, h0hour);
788
789 PlotGrid(name, theta0deg, phi0deg, dec0deg, h0hour);
790
791 return kTRUE;
792}
793
794// --------------------------------------------------------------------------
795//
796// SetGridParameters
797//
798// set the binning for the grid (fGridBinning)
799// set the binning along the lines (fGridFineBin)
800//
801Bool_t MStarCamTrans::SetGridParameters(
802 Double_t gridbinning, Double_t gridfinebin)
803{
804 fGridBinning = gridbinning;
805 fGridFineBin = gridfinebin;
806
807 return kTRUE;
808}
809
810// --------------------------------------------------------------------------
811//
812// PlotGrid
813//
814// - plot the lines of constant (Theta, Phi)
815// with the camera center at (Theta0, Phi0)
816// - and the lines of constant (Dec, H )
817// with the camera center at (Dec0, H0 )
818//
819// (see also MAstroCatalog::Draw and MAstroCamera::Draw)
820//
821// Warning: Leaks Memory!
822//
823Bool_t MStarCamTrans::PlotGrid(TString name,
824 Double_t theta0deg, Double_t phi0deg,
825 Double_t dec0deg, Double_t h0hour)
826{
827 Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;
828 //gLog << "mmtodeg = " << mmtodeg << endl;
829
830 // aberr is the ratio r_optaberr/r_ideal between
831 // the distance from the camera center with optical aberration and
832 // the distance from the camera center with an ideal imaging
833 // the value 1.07 is valid if the expected position (with aberration)
834 // in the camera is calculated as the average of the positions obtained
835 // from the reflections at the individual mirrors
836 Double_t aberr = 1.07;
837
838 //--------------------------------------------------------------------
839
840 TCanvas *c1 = new TCanvas(name, name, 0, 0, 300, 600);
841
842 //gROOT->Reset();
843 gStyle->SetCanvasColor(0);
844 gStyle->SetCanvasBorderMode(0);
845 gStyle->SetPadBorderMode(0);
846 gStyle->SetFrameBorderMode(0);
847 gStyle->SetOptStat(0000000);
848 gStyle->SetPalette(1);
849
850 c1->Divide(1,2);
851 c1->SetBorderMode(0);
852
853 Double_t xlo = -534.0 * mmtodeg;
854 Double_t xup = 534.0 * mmtodeg;
855
856 Double_t ylo = -534.0 * mmtodeg;
857 Double_t yup = 534.0 * mmtodeg;
858
859 TString nam1 = name;
860 nam1 += "_:_Theta-Phi";
861 TString tit1 = name;
862 tit1 += "_:_Theta-Phi-lines";
863 TH2D *plot1 = new TH2D(nam1, tit1, 1, xlo, xup, 1, ylo, yup);
864 plot1->GetXaxis()->SetTitle("x [deg]");
865 plot1->GetYaxis()->SetTitle("y [deg]");
866
867
868 TString nam2 = name;
869 nam2 += "_:_Dec-Hour";
870 TString tit2 = name;
871 tit2 += "_:_Dec-Hour-lines";
872 TH2D *plot2 = new TH2D(nam2, tit2, 1, xlo, xup, 1, ylo, yup);
873 plot2->GetXaxis()->SetTitle("x [deg]");
874 plot2->GetYaxis()->SetTitle("y [deg]");
875
876
877 c1->cd(1);
878 plot1->Draw();
879 //delete plot1;
880
881 c1->cd(2);
882 plot2->Draw();
883 //delete plot2;
884
885 TGraph *graph1;
886 TLatex *pix;
887
888 Char_t tit[100];
889 Double_t xtxt;
890 Double_t ytxt;
891
892 Double_t xx;
893 Double_t yy;
894
895 Double_t gridbinning = fGridBinning;
896 Double_t gridfinebin = fGridFineBin;
897 Int_t numgridlines = (Int_t)(4.0/gridbinning);
898
899
900 //--------------------------------------------------------------------
901 // Theta-Phi lines
902
903 // direction for the center of the camera
904 Double_t theta0 = theta0deg;
905 Double_t phi0 = phi0deg;
906
907
908 //----- lines for fixed theta ------------------------------------
909
910 const Int_t Ntheta = numgridlines;
911 Double_t theta[Ntheta];
912 Double_t dtheta = gridbinning;
913
914 Int_t nphi = (Int_t)(4.0/gridfinebin);
915 const Int_t Nphi = nphi+1;
916 Double_t phi[Nphi];
917 Double_t dphi = gridfinebin/sin(theta0/180.0*3.1415926);
918 if ( dphi > 360.0/((Double_t)(Nphi-1)) )
919 dphi = 360.0/((Double_t)(Nphi-1));
920
921 for (Int_t j=0; j<Ntheta; j++)
922 {
923 theta[j] = theta0 + ((Double_t)j)*dtheta
924 - ((Double_t)(Ntheta/2)-1.0)*dtheta;
925 }
926
927 for (Int_t k=0; k<Nphi; k++)
928 {
929 phi[k] = phi0 + ((Double_t)k)*dphi - ((Double_t)(Nphi/2)-1.0)*dphi;
930 }
931
932
933 Double_t x[Nphi];
934 Double_t y[Nphi];
935
936
937
938 for (Int_t j=0; j<Ntheta; j++)
939 {
940 if (theta[j] < 0.0) continue;
941
942 for (Int_t k=0; k<Nphi; k++)
943 {
944 Loc0LocToCam(theta0, phi0, theta[j], phi[k],
945 xx, yy);
946 x[k] = xx * mmtodeg * aberr;
947 y[k] = yy * mmtodeg * aberr;
948
949 //gLog << "theta0, phi0 = " << theta0 << ", " << phi0
950 // << " : theta, phi, x, y = " << theta[j] << ", "
951 // << phi[k] << "; " << x[k] << ", " << y[k] << endl;
952 }
953
954 c1->cd(1);
955 graph1 = new TGraph(Nphi,x,y);
956 graph1->SetLineColor(j+1);
957 graph1->SetLineStyle(1);
958 graph1->SetMarkerColor(j+1);
959 graph1->SetMarkerSize(.2);
960 graph1->SetMarkerStyle(20);
961 graph1->Draw("PL");
962 //delete graph1;
963
964 sprintf(tit,"Theta = %6.2f", theta[j]);
965 xtxt = xlo + (xup-xlo)*0.1;
966 ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0;
967 pix = new TLatex(xtxt, ytxt, tit);
968 pix->SetTextColor(j+1);
969 pix->SetTextSize(.03);
970 pix->Draw("");
971 //delete pix;
972
973 }
974
975 //----- lines for fixed phi ------------------------------------
976
977 Int_t ntheta1 = (Int_t)(4.0/gridfinebin);
978 const Int_t Ntheta1 = ntheta1;
979 Double_t theta1[Ntheta1];
980 Double_t dtheta1 = gridfinebin;
981
982 const Int_t Nphi1 = numgridlines;
983 Double_t phi1[Nphi1];
984 Double_t dphi1 = gridbinning/sin(theta0/180.0*3.1415926);
985 if ( dphi1 > 360.0/((Double_t)Nphi1) )
986 dphi1 = 360.0/((Double_t)Nphi1);
987
988 for (Int_t j=0; j<Ntheta1; j++)
989 {
990 theta1[j] = theta0 + ((Double_t)j)*dtheta1
991 - ((Double_t)(Ntheta1/2)-1.0)*dtheta1;
992 }
993
994 for (Int_t k=0; k<Nphi1; k++)
995 {
996 phi1[k] = phi0 + ((Double_t)k)*dphi1 - ((Double_t)(Nphi1/2)-1.0)*dphi1;
997 }
998
999
1000 Double_t x1[Ntheta1];
1001 Double_t y1[Ntheta1];
1002
1003 for (Int_t k=0; k<Nphi1; k++)
1004 {
1005 Int_t count = 0;
1006 for (Int_t j=0; j<Ntheta1; j++)
1007 {
1008 if (theta1[j] < 0.0) continue;
1009
1010 Loc0LocToCam(theta0, phi0, theta1[j], phi1[k],
1011 xx, yy);
1012 x1[count] = xx * mmtodeg * aberr;
1013 y1[count] = yy * mmtodeg * aberr;
1014
1015 //gLog << "theta0, phi0 = " << theta0 << ", " << phi0
1016 // << " : theta1, phi1, x1, y1 = " << theta1[j] << ", "
1017 // << phi1[k] << "; " << x1[count] << ", " << y1[count] << endl;
1018 count++;
1019 }
1020
1021 c1->cd(1);
1022 graph1 = new TGraph(count,x1,y1);
1023 graph1->SetLineColor(k+1);
1024 graph1->SetLineStyle(2);
1025 graph1->SetMarkerColor(k+1);
1026 graph1->SetMarkerSize(.2);
1027 graph1->SetMarkerStyle(20);
1028 graph1->Draw("PL");
1029 //delete graph1;
1030
1031 sprintf(tit,"Phi = %6.2f", phi1[k]);
1032 Double_t xtxt = xlo + (xup-xlo)*0.75;
1033 Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0;
1034 pix = new TLatex(xtxt, ytxt, tit);
1035 pix->SetTextColor(k+1);
1036 pix->SetTextSize(.03);
1037 pix->Draw("");
1038 //delete pix;
1039
1040 }
1041
1042 c1->cd(1);
1043 sprintf(tit,"Theta0 = %6.2f Phi0 = %6.2f", theta0, phi0);
1044 xtxt = xlo + (xup-xlo)*0.05;
1045 ytxt = ylo + (yup-ylo)*0.92;
1046 pix = new TLatex(xtxt, ytxt, tit);
1047 pix->SetTextColor(1);
1048 pix->SetTextSize(.06);
1049 pix->Draw("");
1050 //delete pix;
1051
1052 sprintf(tit," [deg] [deg]");
1053 xtxt = xlo + (xup-xlo)*0.05;
1054 ytxt = ylo + (yup-ylo)*0.86;
1055 pix = new TLatex(xtxt, ytxt, tit);
1056 pix->SetTextColor(1);
1057 pix->SetTextSize(.06);
1058 pix->Draw("");
1059 //delete pix;
1060
1061
1062 //--------------------------------------------------------------------
1063 // Dec-Hour lines
1064
1065 // direction for the center of the camera
1066 Double_t dec0 = dec0deg;
1067 Double_t h0 = h0hour;
1068 //trans.LocToCel(theta0, phi0, dec0, h0);
1069
1070
1071 //----- lines for fixed dec ------------------------------------
1072
1073 const Int_t Ndec = numgridlines;
1074 Double_t dec[Ndec];
1075 Double_t ddec = gridbinning;
1076
1077 Int_t nh = (Int_t)(4.0/gridfinebin);
1078 const Int_t Nh = nh+1;
1079 Double_t h[Nh];
1080 Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
1081 if ( dh > 360.0/((Double_t)(Nh-1)) )
1082 dh = 360.0/((Double_t)(Nh-1));
1083
1084
1085 for (Int_t j=0; j<Ndec; j++)
1086 {
1087 dec[j] = dec0 + ((Double_t)j)*ddec
1088 - ((Double_t)(Ndec/2)-1.0)*ddec;
1089 }
1090
1091 for (Int_t k=0; k<Nh; k++)
1092 {
1093 h[k] = h0 + ((Double_t)k)*dh - ((Double_t)(Nh/2)-1.0)*dh;
1094 }
1095
1096
1097 Double_t xh[Nh];
1098 Double_t yh[Nh];
1099
1100
1101
1102 for (Int_t j=0; j<Ndec; j++)
1103 {
1104 if (fabs(dec[j]) > 90.0) continue;
1105
1106 for (Int_t k=0; k<Nh; k++)
1107 {
1108 Double_t hh0 = h0 *24.0/360.0;
1109 Double_t hx = h[k]*24.0/360.0;
1110 Cel0CelToCam(dec0, hh0, dec[j], hx,
1111 xx, yy);
1112 xh[k] = xx * mmtodeg * aberr;
1113 yh[k] = yy * mmtodeg * aberr;
1114
1115 //gLog << "dec0, h0 = " << dec0 << ", " << h0
1116 // << " : dec, h, xh, yh = " << dec[j] << ", "
1117 // << h[k] << "; " << xh[k] << ", " << yh[k] << endl;
1118 }
1119
1120 c1->cd(2);
1121 graph1 = new TGraph(Nh,xh,yh);
1122 graph1->SetLineColor(j+1);
1123 graph1->SetLineStyle(1);
1124 graph1->SetMarkerColor(j+1);
1125 graph1->SetMarkerSize(.2);
1126 graph1->SetMarkerStyle(20);
1127 graph1->Draw("PL");
1128 //delete graph1;
1129
1130 sprintf(tit,"Dec = %6.2f", dec[j]);
1131 xtxt = xlo + (xup-xlo)*0.1;
1132 ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0;
1133 pix = new TLatex(xtxt, ytxt, tit);
1134 pix->SetTextColor(j+1);
1135 pix->SetTextSize(.03);
1136 pix->Draw("");
1137 //delete pix;
1138
1139 }
1140
1141 //----- lines for fixed hour angle ------------------------------------
1142
1143 Int_t ndec1 = (Int_t)(4.0/gridfinebin);
1144 const Int_t Ndec1 = ndec1;
1145 Double_t dec1[Ndec1];
1146 Double_t ddec1 = gridfinebin;
1147
1148 const Int_t Nh1 = numgridlines;
1149 Double_t h1[Nh1];
1150 Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
1151 if ( dh1 > 360.0/((Double_t)Nh1) )
1152 dh1 = 360.0/((Double_t)Nh1);
1153
1154
1155 for (Int_t j=0; j<Ndec1; j++)
1156 {
1157 dec1[j] = dec0 + ((Double_t)j)*ddec1
1158 - ((Double_t)(Ndec1/2)-1.0)*ddec1;
1159 }
1160
1161 for (Int_t k=0; k<Nh1; k++)
1162 {
1163 h1[k] = h0 + ((Double_t)k)*dh1 - ((Double_t)(Nh1/2)-1.0)*dh1;
1164 }
1165
1166
1167 Double_t xd[Ndec1];
1168 Double_t yd[Ndec1];
1169
1170 for (Int_t k=0; k<Nh1; k++)
1171 {
1172 Int_t count = 0;
1173 for (Int_t j=0; j<Ndec1; j++)
1174 {
1175 if (fabs(dec1[j]) > 90.0) continue;
1176
1177 Double_t hh0 = h0 *24.0/360.0;
1178 Double_t hhx = h1[k]*24.0/360.0;
1179 Cel0CelToCam(dec0, hh0, dec1[j], hhx,
1180 xx, yy);
1181 xd[count] = xx * mmtodeg * aberr;
1182 yd[count] = yy * mmtodeg * aberr;
1183
1184 //gLog << "dec0, h0 = " << dec0 << ", " << h0
1185 // << " : dec1, h1, xd, yd = " << dec1[j] << ", "
1186 // << h1[k] << "; " << xd[count] << ", " << yd[count] << endl;
1187
1188 count++;
1189 }
1190
1191 c1->cd(2);
1192 graph1 = new TGraph(count,xd,yd);
1193 graph1->SetLineColor(k+1);
1194 graph1->SetLineStyle(2);
1195 graph1->SetMarkerColor(k+1);
1196 graph1->SetMarkerSize(.2);
1197 graph1->SetMarkerStyle(20);
1198 graph1->Draw("PL");
1199 //delete graph1;
1200
1201 sprintf(tit,"H = %6.2f", h1[k]);
1202 Double_t xtxt = xlo + (xup-xlo)*0.75;
1203 Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0;
1204 pix = new TLatex(xtxt, ytxt, tit);
1205 pix->SetTextColor(k+1);
1206 pix->SetTextSize(.03);
1207 pix->Draw("");
1208 //delete pix;
1209
1210 }
1211
1212 c1->cd(2);
1213 sprintf(tit,"Dec0 = %6.2f H0 = %6.2f", dec0, h0);
1214 xtxt = xlo + (xup-xlo)*0.05;
1215 ytxt = ylo + (yup-ylo)*0.92;
1216 pix = new TLatex(xtxt, ytxt, tit);
1217 pix->SetTextColor(1);
1218 pix->SetTextSize(.06);
1219 pix->Draw("");
1220 //delete pix;
1221
1222 sprintf(tit," [deg] [deg]");
1223 xtxt = xlo + (xup-xlo)*0.05;
1224 ytxt = ylo + (yup-ylo)*0.86;
1225 pix = new TLatex(xtxt, ytxt, tit);
1226 pix->SetTextColor(1);
1227 pix->SetTextSize(.06);
1228 pix->Draw("");
1229 //delete pix;
1230
1231 return kTRUE;
1232}
1233
Note: See TracBrowser for help on using the repository browser.