source: trunk/FACT++/src/calcsource.cc@ 19054

Last change on this file since 19054 was 19053, checked in by tbretz, 7 years ago
This is a draft to fill in source positions to the database.
File size: 28.3 KB
Line 
1#include <boost/algorithm/string/join.hpp>
2
3#include "Database.h"
4
5#include "pal.h"
6#include "nova.h"
7#include "Time.h"
8#include "Configuration.h"
9
10#include <TROOT.h>
11#include <TSystem.h>
12#include <TFile.h>
13#include <TTree.h>
14#include <TLeaf.h>
15#include <TError.h>
16#include <TVector3.h>
17#include <TRotation.h>
18
19using namespace std;
20
21// ------------------------------------------------------------------------
22
23struct Map : pair<string, string>
24{
25 Map() { }
26};
27
28std::istream &operator>>(std::istream &in, Map &m)
29{
30 string txt;
31 in >> txt;
32
33 const auto p = txt.find_first_of('/');
34 if (p==string::npos)
35 throw runtime_error("Missing / in map argument.");
36 if (p!=txt.find_last_of('/'))
37 throw runtime_error("Only one / allowed in map argument.");
38
39 m.first = txt.substr(0, p);
40 m.second = txt.substr(p+1);
41
42 return in;
43}
44
45
46void SetupConfiguration(Configuration &conf)
47{
48 po::options_description control("Root to Database");
49 control.add_options()
50 ("uri,u", var<string>()->required(), "Database link as in\n\tuser:password@server[:port]/database.")
51 ("ra", var<double>()/*5.575539)*/, "")
52 ("dec", var<double>()/*22.014500)*/, "")
53 ("focal-dist", var<double>(4889.), "")
54 //("source-key", var<uint32_t>(5), "")
55 //("source-name", var<string>(""), "")
56 ("file", var<uint32_t>(171011115)->required(), "")
57 ("drop", po_switch(), "Drop the table (implied create)")
58 ("list-files", po_switch(), "")
59 ("tree,t", var<string>("Events"), "Name of the root tree to convert")
60 ("table", var<string>("Events"), "")
61 ("update", var<string>(""), "")
62 ("create", po_switch(), "")
63 ("verbose,v", var<uint16_t>(1), "Verbosity (0: quiet, 1: default, 2: more, 3, ...)")
64 ("no-insert", po_switch(), "Does not insert any data to the table")
65 //("print-insert", po_switch(), "Print the INSERT query (note that it contains all data)")
66 //("print-create", po_switch(), "Print the CREATE query")
67 ;
68
69 po::positional_options_description p;
70 p.add("file", 1); // The 1st positional options (n=1)
71
72 conf.AddOptions(control);
73 conf.SetArgumentPositions(p);
74}
75
76void PrintUsage()
77{
78 cout <<
79 "calcsource - Fill a table with source positions\n"
80// "\n"
81// "Usage: rootifysql [rootify.sql [rootify.root]] [-u URI] [-q query|-f file] [-i] [-o out] [-f] [-cN] [-t tree] [-vN]\n"
82 "\n"
83 ;
84 cout << endl;
85}
86
87enum BasicType_t
88{
89 kNone = 0,
90 kFloat,
91 kDouble,
92 kInt16,
93 kUInt16,
94 kInt32,
95 kUInt32,
96 kInt64,
97 kUInt64,
98};
99
100static const map<string, pair<BasicType_t, string>> ConvRoot =
101{
102 { "Float_t", { kFloat, "FLOAT" } },
103 { "Double_t", { kDouble, "DOUBLE" } },
104 { "ULong64_t", { kUInt64, "BIGINT UNSIGNED" } },
105 { "Long64_t", { kInt64, "BIGINT" } },
106 { "UInt_t", { kUInt32, "INT UNSIGNED" } },
107 { "Int_t", { kInt32, "INT" } },
108 { "UShort_t", { kUInt16, "SMALLINT UNSIGNED" } },
109 { "Short_t", { kInt16, "SMALLINT" } },
110};
111
112struct Container
113{
114 string branch; // branch name
115 string column; // column name
116 BasicType_t type;
117 void *ptr;
118
119 Container(const string &b, const string &c, const BasicType_t &t) : branch(b), column(c), type(t), ptr(0)
120 {
121 switch (t)
122 {
123 case kFloat: ptr = new Float_t; break;
124 case kDouble: ptr = new Double_t; break;
125 case kInt16: ptr = new Short_t; break;
126 case kUInt16: ptr = new UShort_t; break;
127 case kInt32: ptr = new Int_t; break;
128 case kUInt32: ptr = new UInt_t; break;
129 case kInt64: ptr = new Long64_t; break;
130 case kUInt64: ptr = new ULong64_t; break;
131 case kNone:
132 break;
133 }
134 }
135 ~Container()
136 {
137 //::operator delete(ptr); // It seems root is deleting it already
138 }
139
140 string fmt() const
141 {
142 ostringstream str;
143
144 switch (type)
145 {
146 case kFloat: str << setprecision(8) << *reinterpret_cast<Float_t*>(ptr); break;
147 case kDouble: str << setprecision(16) << *reinterpret_cast<Double_t*>(ptr); break;
148 case kInt16: str << *reinterpret_cast<Short_t*>(ptr); break;
149 case kUInt16: str << *reinterpret_cast<UShort_t*>(ptr); break;
150 case kInt32: str << *reinterpret_cast<Int_t*>(ptr); break;
151 case kUInt32: str << *reinterpret_cast<UInt_t*>(ptr); break;
152 case kInt64: str << *reinterpret_cast<Long64_t*>(ptr); break;
153 case kUInt64: str << *reinterpret_cast<ULong64_t*>(ptr); break;
154 case kNone:
155 break;
156 }
157
158 return str.str();
159 }
160};
161
162class MRotation : public TRotation
163{
164public:
165 MRotation() : TRotation(1, 0, 0, 0, -1, 0, 0, 0, 1)
166 {
167 }
168};
169
170
171void ErrorHandlerAll(Int_t level, Bool_t abort, const char *location, const char *msg)
172{
173 if (string(msg).substr(0,24)=="no dictionary for class ")
174 return;
175
176 DefaultErrorHandler(level, abort, location, msg);
177}
178
179int main(int argc, const char* argv[])
180{
181 Time start;
182
183 gROOT->SetBatch();
184 SetErrorHandler(ErrorHandlerAll);
185
186 Configuration conf(argv[0]);
187 conf.SetPrintUsage(PrintUsage);
188 SetupConfiguration(conf);
189
190 if (!conf.DoParse(argc, argv))
191 return 127;
192
193 // ----------------------------- Evaluate options --------------------------
194 const bool has_radec = conf.Has("ra") && conf.Has("dec");
195
196 const string uri = conf.Get<string>("uri");
197 const uint32_t file = conf.Get<uint32_t>("file");
198 const string tree = conf.Get<string>("tree");
199 const string table = conf.Get<string>("table");
200 //string source_name = conf.Get<string>("source-name");
201 //uint32_t source_key = conf.Has("source-key") ? conf.Get<uint32_t>("source-key") : 0;
202 double source_ra = conf.Has("ra") ? conf.Get<double>("ra") : 0;
203 double source_dec = conf.Has("dec") ? conf.Get<double>("dec") : 0;
204 double focal_dist = conf.Get<double>("focal-dist");
205 //const bool print_create = conf.Get<bool>("print-create");
206 //const bool print_insert = conf.Get<bool>("print-insert");
207 const bool drop = conf.Get<bool>("drop");
208 const bool create = conf.Get<bool>("create") || drop;
209 const string update = conf.Get<string>("update");
210 const bool list_files = conf.Get<bool>("list-files");
211 const uint16_t verbose = conf.Get<uint16_t>("verbose");
212
213 // -------------------------------------------------------------------------
214
215 if (list_files)
216 {
217 const string query =
218 "SELECT FileId FROM Events GROUP BY FileId";
219
220 cout << query << endl;
221
222 const mysqlpp::StoreQueryResult res =
223 Database(uri).query(query).store();
224
225 for (size_t i=0; i<res.num_rows(); i++)
226 cout << "calcsource " << res[i][0] << '\n';
227 cout << endl;
228
229 return 0;
230 }
231
232 if (verbose>0)
233 cout << "\n------------------------- Evaluating source ------------------------" << endl;
234
235 cout << "Requesting coordinates from RunInfo/Source table for 20" << file/1000 << "/" << file%1000 << endl;
236
237 if (!has_radec)
238 {
239 const string query =
240 "SELECT Source.fRightAscension, Source.fDeclination, Source.fSourceName"
241 " FROM RunInfo"
242 " LEFT JOIN Source"
243 " USING (fSourceKey)"
244 " WHERE fNight=20"+to_string(file/1000)+
245 " AND fRunID="+to_string(file%1000);
246
247 cout << query << endl;
248
249 const mysqlpp::StoreQueryResult res =
250 Database(uri).query(query).store();
251
252 if (res.num_rows()!=1)
253 {
254 cerr << "No coordinates from RunInfo for " << file << endl;
255 return 1;
256 }
257
258 source_ra = res[0][0];
259 source_dec = res[0][1];
260
261 cout << "Using coordinates " << source_ra << "h / " << source_dec << " deg for " << res[0][2] << endl;
262 }
263 else
264 cout << "Using coordinates " << source_ra << "h / " << source_dec << " deg from resources." << endl;
265
266/*
267 if (!source_name.empty())
268 {
269 cout << "Requesting coordinates for '" << source_name << "'" << endl;
270
271 const mysqlpp::StoreQueryResult res =
272 Database(uri).query("SELECT `Ra`, `Dec` WHERE fSourceName='"+source_name+"'").store();
273
274 if (res.num_rows()!=1)
275 {
276 cerr << "No " << (res.num_rows()>1?"unique ":"") << "coordinates found for '" << source_name << "'" << endl;
277 return 1;
278 }
279
280 source_ra = res[0][0];
281 source_dec = res[0][1];
282 }
283*/
284
285 if (verbose>0)
286 cout << "\n-------------------------- Evaluating file ------------------------" << endl;
287
288
289 Database connection(uri); // Keep alive while fetching rows
290
291 const string query =
292 "SELECT `Ra`, `Dec`, MJD, MilliSec, NanoSec, Zd, Az, EvtNumber"
293 " FROM `"+table+"`"
294 " WHERE FileId="+to_string(file);
295
296 cout << query << endl;
297
298 const mysqlpp::UseQueryResult res1 =
299 connection.query(query).use();
300
301 Nova::RaDecPosn source(source_ra, source_dec);
302
303 source_ra *= M_PI/12;
304 source_dec *= M_PI/180;
305
306 auto obs = Nova::kORM;
307
308 obs.lng *= M_PI/180;
309 obs.lat *= M_PI/180;
310
311 //const double mm2deg = 1.17193246260285378e-02;
312
313 ostringstream ins;
314 ins << setprecision(16);
315
316 ostringstream upd;
317 upd << setprecision(16);
318
319 size_t count = 0;
320 while (auto row=res1.fetch_row())
321 {
322 count++;
323
324 double point_ra = row[0];
325 double point_dec = row[1];
326 uint32_t mjd = row[2];
327 int64_t millisec = row[3];
328 //uint32_t nanosec = row[4];
329 //double zd = row[5];
330 //double az = row[6];
331 uint32_t event = row[7];
332
333 Nova::RaDecPosn point(point_ra, point_dec);
334
335 point_ra *= M_PI/12;
336 point_dec *= M_PI/180;
337
338 /*
339 // ============================ Mars ================================
340
341 TVector3 pos; // pos: source position
342 TVector3 pos0; // pos: source position
343
344 pos.SetMagThetaPhi(1, M_PI/2-source_dec, source_ra);
345 pos0.SetMagThetaPhi(1, M_PI/2-point_dec, point_ra);
346
347 const double ut = (nanosec/1e6+millisec)/(24*3600000);
348
349 // Julian centuries since J2000.
350 const double t = (ut -(51544.5-mjd)) / 36525.0;
351
352 // GMST at this UT1
353 const double r1 = 24110.54841+(8640184.812866+(0.093104-6.2e-6*t)*t)*t;
354 const double r2 = 86400.0*ut;
355
356 const double sum = (r1+r2)/(3600*24);
357
358 double gmst = fmod(sum, 1) * 2*M_PI;
359
360 MRotation conv;
361 conv.RotateZ(gmst + obs.lng);
362 conv.RotateY(obs.lat-M_PI/2);
363 conv.RotateZ(M_PI);
364
365 pos *= conv;
366 pos0 *= conv;
367
368 pos.RotateZ(-pos0.Phi());
369 pos.RotateY(-pos0.Theta());
370 pos.RotateZ(-M_PI/2); // exchange x and y
371 pos *= -focal_dist/pos.Z();
372
373 TVector2 v = pos.XYvector();
374
375 //if (fDeviation)
376 // v -= fDeviation->GetDevXY()/fGeom->GetConvMm2Deg();
377
378 //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n';
379 */
380
381 // ================================= Nova ===============================
382
383 Nova::ZdAzPosn ppos = Nova::GetHrzFromEqu(source, 2400000.5+mjd+millisec/1000./3600/24);
384 Nova::ZdAzPosn ppos0 = Nova::GetHrzFromEqu(point, 2400000.5+mjd+millisec/1000./3600/24);
385
386 TVector3 pos;
387 TVector3 pos0;
388 pos.SetMagThetaPhi( 1, ppos.zd *M_PI/180, ppos.az *M_PI/180);
389 pos0.SetMagThetaPhi(1, ppos0.zd*M_PI/180, ppos0.az*M_PI/180);
390
391 pos.RotateZ(-pos0.Phi());
392 pos.RotateY(-pos0.Theta());
393 pos.RotateZ(-M_PI/2); // exchange x and y
394 pos *= -focal_dist/pos.Z();
395
396 TVector2 v = pos.XYvector();
397
398 //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n';
399
400 /*
401 // =============================== Slalib ==========================
402 const double height = 2200;
403 const double temp = 10;
404 const double hum = 0.25;
405 const double press = 780;
406
407 const double _mjd = mjd+millisec/1000./3600/24;
408
409 const double dtt = palDtt(_mjd); // 32.184 + 35
410
411 const double tdb = _mjd + dtt/3600/24;
412 const double dut = 0;
413
414 // prepare calculation: Mean Place to geocentric apperent
415 // (UTC would also do, except for the moon?)
416 double fAmprms[21];
417 palMappa(2000.0, tdb, fAmprms); // Epoche, TDB
418
419 // prepare: Apperent to observed place
420 double fAoprms[14];
421 palAoppa(_mjd, dut, // mjd, Delta UT=UT1-UTC
422 obs.lng, obs.lat, height, // long, lat, height
423 0, 0, // polar motion x, y-coordinate (radians)
424 273.155+temp, press, hum, // temp, pressure, humidity
425 0.40, 0.0065, // wavelength, tropo lapse rate
426 fAoprms);
427
428 // ---- Mean to apparent ----
429 double r=0, d=0;
430 palMapqkz(point_ra, point_dec, fAmprms, &r, &d);
431
432 double _zd, _az, ha, ra, dec;
433 // -- apparent to observed --
434 palAopqk(r, d, fAoprms,
435 &_az, // observed azimuth (radians: N=0,E=90) [-pi, pi]
436 &_zd, // observed zenith distance (radians) [-pi/2, pi/2]
437 &ha, // observed hour angle (radians)
438 &dec, // observed declination (radians)
439 &ra); // observed right ascension (radians)
440
441 //cout << _zd*180/M_PI << " " << _az*180/M_PI << endl;
442
443 pos0.SetMagThetaPhi(1, _zd, _az);
444
445 r=0, d=0;
446 palMapqkz(source_ra, source_dec, fAmprms, &r, &d);
447
448 // -- apparent to observed --
449 palAopqk(r, d, fAoprms,
450 &_az, // observed azimuth (radians: N=0,E=90) [-pi, pi]
451 &_zd, // observed zenith distance (radians) [-pi/2, pi/2]
452 &ha, // observed hour angle (radians)
453 &dec, // observed declination (radians)
454 &ra); // observed right ascension (radians)
455
456 pos.SetMagThetaPhi(1, _zd, _az);
457
458 //cout << _zd*180/M_PI << " " << _az*180/M_PI << endl;
459
460 pos.RotateZ(-pos0.Phi());
461 pos.RotateY(-pos0.Theta());
462 pos.RotateZ(-M_PI/2); // exchange x and y
463 pos *= -focal_dist/pos.Z();
464
465 v = pos.XYvector();
466
467 //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n';
468 */
469
470 if (1/*insert*/)
471 ins << "( " << file << ", " << event << ", " << v.X() << ", " << v.Y() << " ),\n";
472
473 if (!update.empty())
474 upd << "UPDATE " << update << " SET X=" << v.X() << ", Y=" << v.Y() << " WHERE FileId=" << file << " AND EvtNumber=" << event <<";\n";
475 };
476
477 if (connection.errnum())
478 {
479 cerr << "SQL error fetching row: " << connection.error() << endl;
480 return 7;
481 }
482
483
484 if (drop)
485 {
486 cout << "Drop table Position." << endl;
487 const mysqlpp::SimpleResult res2 =
488 Database(uri).query("DROP TABLE Position").execute();
489
490 //cout << res.rows() << " rows affected." << res.info() << endl;
491 }
492
493
494 if (create)
495 {
496 cout << "Create table Position." << endl;
497 const mysqlpp::SimpleResult res3 =
498 Database(uri).query("CREATE TABLE IF NOT EXISTS Position\n"
499 "(\n"
500 " FileId INT UNSIGNED NOT NULL,\n"
501 " EvtNumber INT UNSIGNED NOT NULL,\n"
502 " X FLOAT NOT NULL,\n"
503 " Y FLOAT NOT NULL,\n"
504 " PRIMARY KEY (FileId, EvtNumber)\n"
505 ")\n"
506 "DEFAULT CHARSET=latin1 COLLATE=latin1_general_ci\n"
507 "ENGINE=MyISAM\n"
508 "COMMENT='created by root2db'\n").execute();
509
510 //cout << res.rows() << " rows affected." << res.info() << endl;
511
512 }
513 else
514 {
515 // FIXME: Only if entries exist?
516 cout << "Delete old entries from Position." << endl;
517 const mysqlpp::SimpleResult res3 =
518 Database(uri).query("DELETE FROM Position WHERE FileId="+to_string(file)).execute();
519
520 cout << res3.rows() << " rows affected." << endl;
521 }
522
523 if (1/*insert*/)
524 {
525 cout << "Insert data into table Position." << endl;
526 const string query3 =
527 "INSERT `Position` (FileId, EvtNumber, X, Y) VALUES\n"+
528 ins.str().substr(0, ins.str().size()-2)+
529 "\n";
530
531 if (0/*print_insert*/)
532 cout << query3 << endl;
533
534 const mysqlpp::SimpleResult res3 =
535 Database(uri).query(query3).execute();
536
537 cout << res3.info() << endl;
538 }
539
540 if (!update.empty())
541 {
542 cout << upd.str() << endl;
543 const mysqlpp::StoreQueryResult res3 =
544 connection.query(upd.str()).store();
545
546 /*
547 mysqlpp::Connection con(db, server, user, pass);
548 // Set option to allow multiple queries to be issued.
549 mysqlpp::MultiStatementsOption* opt = new mysqlpp::MultiStatementsOption(true);
550 con.set_option(opt);
551 // Build either single queries or multiple queries strung together.
552 mysqlpp::Query query = con.query();
553 ....
554 // Issue query
555 query.store();
556 // Do NOT delete 'opt'; it will be destroyed by the 'con' object when it
557 // falls out of scope.
558 */
559 }
560
561 //cout << ins.str().substr(0, ins.str().size()-2) << endl;
562 //cout << count << endl;
563
564 return 0;
565}
566
567/*
568 TRotation & TRotation::RotateX(Double_t a) {
569 //rotate around x
570 Double_t c = TMath::Cos(a);
571 Double_t s = TMath::Sin(a);
572 Double_t x = fyx, y = fyy, z = fyz;
573 fyx = c*x - s*fzx;
574 fyy = c*y - s*fzy;
575 fyz = c*z - s*fzz;
576 fzx = s*x + c*fzx;
577 fzy = s*y + c*fzy;
578 fzz = s*z + c*fzz;
579 return *this;
580}
581
582TRotation & TRotation::RotateY(Double_t a){
583 //rotate around y
584 Double_t c = TMath::Cos(a);
585 Double_t s = TMath::Sin(a);
586 Double_t x = fzx, y = fzy, z = fzz;
587 fzx = c*x - s*fxx;
588 fzy = c*y - s*fxy;
589 fzz = c*z - s*fxz;
590 fxx = s*x + c*fxx;
591 fxy = s*y + c*fxy;
592 fxz = s*z + c*fxz;
593 return *this;
594}
595
596TRotation & TRotation::RotateZ(Double_t a) {
597 //rotate around z
598 Double_t c = TMath::Cos(a);
599 Double_t s = TMath::Sin(a);
600 Double_t x = fxx, y = fxy, z = fxz;
601 fxx = c*x - s*fyx;
602 fxy = c*y - s*fyy;
603 fxz = c*z - s*fyz;
604 fyx = s*x + c*fyx;
605 fyy = s*y + c*fyy;
606 fyz = s*z + c*fyz;
607 return *this;
608}
609
610 */
611
612
613// Reuired:
614// pointing_ra[8]
615// pointing_dec[8]
616// fNanoSec[4], fTime[fMiliSec,8], fMjd[4]
617// source_ra -> rc
618// source_dec -> rc
619// fLong, fLat -> rc
620// fCameraDist -> rc
621// 32 byte per row, 1 Monat = 4mio rows => 120 MB
622
623/*
624
625void TVector3::SetThetaPhi(Double_t theta, Double_t phi)
626{
627 //setter with mag, theta, phi
628 fX = TMath::Sin(theta) * TMath::Cos(phi);
629 fY = TMath::Sin(theta) * TMath::Sin(phi);
630 fZ = TMath::Cos(theta);
631}
632
633
634 // Set Sky coordinates of source, taken from container "MSourcePos"
635 // of type MPointingPos. The sky coordinates must be J2000, as the
636 // sky coordinates of the camera center that we get from the container
637 // "MPointingPos" filled by the Drive.
638 //MVector3 pos; // pos: source position
639 //pos.SetRaDec(fSourcePos->GetRaRad(), fSourcePos->GetDecRad());
640
641 TVector3 pos; // pos: source position
642 pos.SetMagThetaPhi(1, TMath::Pi()/2-source_dec, source_ra);
643
644 // Set sky coordinates of camera center in pos0 (for details see below)
645 //MVector3 pos0; // pos0: camera center
646 //pos0.SetRaDec(fPointPos->GetRaRad(), fPointPos->GetDecRad());
647
648 TVector3 pos0; // pos: source position
649 pos0.SetMagThetaPhi(1, TMath::Pi()/2-pointing_dec, pointing_ra);
650
651 // Convert sky coordinates of source to local coordinates. Warning! These are not the "true" local
652 // coordinates, since this transformation ignores precession and nutation effects.
653
654 //const MAstroSky2Local conv(*fTime, *fObservatory);
655 //pos *= conv;
656
657 const Double_t ut = (Double_t)(fNanoSec/1e6+(Long_t)fTime)/kDay;
658
659 // Julian centuries since J2000.
660 const Double_t t = (ut -(51544.5-fMjd)) / 36525.0;
661
662 // GMST at this UT1
663 const Double_t r1 = 24110.54841+(8640184.812866+(0.093104-6.2e-6*t)*t)*t;
664 const Double_t r2 = 86400.0*ut;
665
666 const Double_t sum = (r1+r2)/kDaySec;
667
668 double gmst = fmod(sum, 1)*TMath::TwoPi();
669
670 // TRotation::TRotation(Double_t mxx, Double_t mxy, Double_t mxz,
671 // Double_t myx, Double_t myy, Double_t myz,
672 // Double_t mzx, Double_t mzy, Double_t mzz)
673 // : fxx(mxx), fxy(mxy), fxz(mxz),
674 // fyx(myx), fyy(myy), fyz(myz),
675 // fzx(mzx), fzy(mzy), fzz(mzz) {}
676
677 TRotation conv(1, 0, 0, 0, -1, 0, 0, 0, 1);
678 conv.RotateZ(gmst + obs.GetElong());
679 conv.RotateY(obs.GetPhi()-TMath::Pi()/2);
680 conv.RotateZ(TMath::Pi());
681
682 // TVector3 & TVector3::operator *= (const TRotation & m){
683 // return *this = m * (*this);
684 // }
685
686 // inline TVector3 TRotation::operator * (const TVector3 & p) const {
687 // return TVector3(fxx*p.X() + fxy*p.Y() + fxz*p.Z(),
688 // fyx*p.X() + fyy*p.Y() + fyz*p.Z(),
689 // fzx*p.X() + fzy*p.Y() + fzz*p.Z());
690 // }
691
692 // Convert sky coordinates of camera center convert to local.
693 // Same comment as above. These coordinates differ from the true
694 // local coordinates of the camera center that one could get from
695 // "MPointingPos", calculated by the Drive: the reason is that the Drive
696 // takes into account precession and nutation corrections, while
697 // MAstroSky2Local (as of Jan 27 2005 at least) does not. Since we just
698 // want to get the source position on the camera from the local
699 // coordinates of the center and the source, it does not matter that
700 // the coordinates contained in pos and pos0 ignore precession and
701 // nutation... since the shift would be the same in both cases. What
702 // would be wrong is to set in pos0 directly the local coordinates
703 // found in MPointingPos!
704 pos0 *= conv;
705
706 if (fDeviation)
707 {
708 // Position at which the starguider camera is pointing in real:
709 // pointing position = nominal position - dev
710 //
711 //vx = CalcXYinCamera(pos0, pos)*fGeom->GetCameraDist()*1000;
712 pos0.SetZdAz(pos0.Theta()-fDeviation->GetDevZdRad(),
713 pos0.Phi() -fDeviation->GetDevAzRad());
714 }
715
716 // Calculate source position in camera, and convert to mm:
717 TVector2 v = MAstro::GetDistOnPlain(pos0, pos, -fGeom->GetCameraDist()*1000);
718
719 pos.RotateZ(-pos0.Phi());
720 pos.RotateY(-pos0.Theta());
721 pos.RotateZ(-TMath::Pi()/2); // exchange x and y
722 pos *= -fGeom->GetCameraDist()*1000/pos.Z();
723
724 TVector2 v = pos.XYvector();
725
726 if (fDeviation)
727 v -= fDeviation->GetDevXY()/fGeom->GetConvMm2Deg();
728
729 SetSrcPos(v);
730
731 // ================================================
732
733 if (fMode==kWobble)
734 {
735 // The trick here is that the anti-source position in case
736 // of the off-source regions is always the on-source positon
737 // thus a proper anti-source cut is possible.
738 fSrcPosAnti->SetXY(v);
739 if (fCallback)
740 v = Rotate(v, fCallback->GetNumPass(), fCallback->GetNumPasses());
741 else
742 v *= -1;
743 fSrcPosCam->SetXY(v);
744 }
745 else
746 {
747 // Because we don't process this three times like in the
748 // wobble case we have to find another way to determine which
749 // off-source region is the right anti-source position
750 // We do it randomly.
751 fSrcPosCam->SetXY(v);
752 if (fNumRandomOffPositions>1)
753 v = Rotate(v, gRandom->Integer(fNumRandomOffPositions), fNumRandomOffPositions);
754 else
755 v *= -1;
756 fSrcPosAnti->SetXY(v);
757 }
758
759
760
761 */
762
763
764
765
766/*
767 PointingData CalcPointingPos(const PointingSetup &setup, double _mjd, const Weather &weather, uint16_t timeout, bool tpoint=false)
768 {
769 PointingData out;
770 out.mjd = _mjd;
771
772 const double elong = Nova::kORM.lng * M_PI/180;
773 const double lat = Nova::kORM.lat * M_PI/180;
774 const double height = 2200;
775
776 const bool valid = weather.time+boost::posix_time::seconds(timeout) > Time();
777
778 const double temp = valid ? weather.temp : 10;
779 const double hum = valid ? weather.hum : 0.25;
780 const double press = valid ? weather.press : 780;
781
782 const double dtt = palDtt(_mjd); // 32.184 + 35
783
784 const double tdb = _mjd + dtt/3600/24;
785 const double dut = 0;
786
787 // prepare calculation: Mean Place to geocentric apperent
788 // (UTC would also do, except for the moon?)
789 double fAmprms[21];
790 palMappa(2000.0, tdb, fAmprms); // Epoche, TDB
791
792 // prepare: Apperent to observed place
793 double fAoprms[14];
794 palAoppa(_mjd, dut, // mjd, Delta UT=UT1-UTC
795 elong, lat, height, // long, lat, height
796 0, 0, // polar motion x, y-coordinate (radians)
797 273.155+temp, press, hum, // temp, pressure, humidity
798 0.40, 0.0065, // wavelength, tropo lapse rate
799 fAoprms);
800
801 out.source.ra = setup.source.ra * M_PI/ 12;
802 out.source.dec = setup.source.dec * M_PI/180;
803
804 if (setup.planet!=kENone)
805 {
806 // coordinates of planet: topocentric, equatorial, J2000
807 // One can use TT instead of TDB for all planets (except the moon?)
808 double ra, dec, diam;
809 palRdplan(tdb, setup.planet, elong, lat, &ra, &dec, &diam);
810
811 // ---- apparent to mean ----
812 palAmpqk(ra, dec, fAmprms, &out.source.ra, &out.source.dec);
813 }
814
815 if (setup.wobble_offset<=0 || tpoint)
816 {
817 out.pointing.dec = out.source.dec;
818 out.pointing.ra = out.source.ra;
819 }
820 else
821 {
822 const double dphi =
823 setup.orbit_period==0 ? 0 : 2*M_PI*(_mjd-setup.start)/setup.orbit_period;
824
825 const double phi = setup.wobble_angle + dphi;
826
827 const double cosdir = cos(phi);
828 const double sindir = sin(phi);
829 const double cosoff = cos(setup.wobble_offset);
830 const double sinoff = sin(setup.wobble_offset);
831 const double cosdec = cos(out.source.dec);
832 const double sindec = sin(out.source.dec);
833
834 const double sintheta = sindec*cosoff + cosdec*sinoff*cosdir;
835
836 const double costheta = sintheta>1 ? 0 : sqrt(1 - sintheta*sintheta);
837
838 const double cosdeltara = (cosoff - sindec*sintheta)/(cosdec*costheta);
839 const double sindeltara = sindir*sinoff/costheta;
840
841 out.pointing.dec = asin(sintheta);
842 out.pointing.ra = atan2(sindeltara, cosdeltara) + out.source.ra;
843 }
844
845 // ---- Mean to apparent ----
846 double r=0, d=0;
847 palMapqkz(out.pointing.ra, out.pointing.dec, fAmprms, &r, &d);
848
849 //
850 // Doesn't work - don't know why
851 //
852 // slaMapqk (radec.Ra(), radec.Dec(), rdpm.Ra(), rdpm.Dec(),
853 // 0, 0, (double*)fAmprms, &r, &d);
854 //
855
856 // -- apparent to observed --
857 palAopqk(r, d, fAoprms,
858 &out.sky.az, // observed azimuth (radians: N=0,E=90) [-pi, pi]
859 &out.sky.zd, // observed zenith distance (radians) [-pi/2, pi/2]
860 &out.apparent.ha, // observed hour angle (radians)
861 &out.apparent.dec, // observed declination (radians)
862 &out.apparent.ra); // observed right ascension (radians)
863
864 // ----- fix ambiguity -----
865 if (out.sky.zd<0)
866 {
867 out.sky.zd = -out.sky.zd;
868 out.sky.az += out.sky.az<0 ? M_PI : -M_PI;
869 }
870
871 // Star culminating behind zenith and Az between ~90 and ~180deg
872 if (out.source.dec<lat && out.sky.az>0)
873 out.sky.az -= 2*M_PI;
874
875 out.mount = SkyToMount(out.sky);
876
877 return out;
878 }
879};
880
881*/
882
Note: See TracBrowser for help on using the repository browser.