source: trunk/MagicSoft/Mars/mstarcam/MStarCamTrans.cc@ 4692

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