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

Last change on this file since 19108 was 19107, checked in by tbretz, 7 years ago
Implmented ignore-errors option.
File size: 33.8 KB
Line 
1#include "Database.h"
2
3#include "pal.h"
4#include "nova.h"
5#include "Time.h"
6#include "Configuration.h"
7
8#include <TROOT.h>
9#include <TVector3.h>
10#include <TRotation.h>
11
12using namespace std;
13
14// ------------------------------------------------------------------------
15
16void SetupConfiguration(Configuration &conf)
17{
18 po::options_description control("Calcsource options");
19 control.add_options()
20 ("uri,u", var<string>()
21#if BOOST_VERSION >= 104200
22 ->required()
23#endif
24 , "Database link as in\n\tuser:password@server[:port]/database.")
25 //("source-key", var<uint32_t>(5), "")
26 //("source-name", var<string>(""), "")
27 ("file", var<uint32_t>(uint32_t(0)),"FileId (YYMMDDXXX), defined the events to be processed (if omitted (=0), a list of possible numbers is printed)")
28 ("table.events", var<string>("Events"), "Name of the table where the events are stored")
29 ("table.runinfo", var<string>("RunInfo"), "Name of the table where the run-info data is stored")
30 ("table.source", var<string>("Source"), "Name of the table where the sources are stored")
31 ("table.position", var<string>("Position"), "Name of the table where the calculated posiiton will be stored")
32 ("ignore-errors", po_switch(), "Adds the IGNORE keyword to the INSERT query (turns errors into warnings, ignores rows with errors)")
33 ("force", po_switch(), "Force processing even if there is no database connection")
34 ("drop", po_switch(), "Drop the table (implies create)")
35 ("create", po_switch(), "Create the table if it does not yet exist")
36 ("update", po_switch(), "Update the table.position instead of inserting it (disables drop, create and delete)")
37 ("ra", var<double>(), "Right ascension of the source (use together with --dec)")
38 ("dec", var<double>(), "Declination of the source (use together with --ra)")
39 ("focal-dist", var<double>(4889.), "Focal distance of the camera in millimeter")
40 ;
41
42 po::options_description debug("Debug options");
43 debug.add_options()
44 ("no-insert", po_switch(), "Does not insert or update any data to any table")
45 ("dry-run", po_switch(), "Skip any query which changes the databse (might result in consecutive failures)")
46 ("print-meta", po_switch(), "Print meta-queries (DROP, CREATE, DELETE, SELECT)")
47 ("print-insert", po_switch(), "Print the INSERT/UPDATE queries")
48 ("verbose,v", var<uint16_t>(1), "Verbosity (0: quiet, 1: default, 2: more, 3, ...)")
49 ;
50
51 po::positional_options_description p;
52 p.add("file", 1); // The 1st positional options (n=1)
53
54 conf.AddOptions(control);
55 conf.AddOptions(debug);
56 conf.SetArgumentPositions(p);
57}
58
59void PrintUsage()
60{
61 cout <<
62 "calcsource - Fill a table with source positions\n"
63 "\n"
64 "This tool is to calculate the source position in the camera and fill "
65 "it into a database table for all events which come from a single run. "
66 "The run is idetified by its FileId (YYMMDDXXX), where YYMMDD is "
67 "its date and XXX is the run-number. For this run the corresponding "
68 "source key `fSourceKey` is obtained from the RunInfo table "
69 "(--table.runinfo) and the corresponding fRightAscension and fDeclination "
70 "from the Source table (--table.source). If --ra and --dec was specified "
71 "by the user, this step is skipped and these two values are used instead."
72 "\n\n"
73 "Then for each event with the given FileId, its pointing position "
74 "(`Ra`, `Dec`), its time (`MJD`, `MilliSec`) and its event number "
75 "(`EvtNumber`) are requested from the table given by --table.events. "
76 "Based on this information and the focal distance (--focal-distance) the "
77 "`X` and `Y` coordinate of the source in the camera plane is calculated. "
78 "The result can then be filled into a database."
79 "\n\n"
80 "The table to be filled or updated should contain the following columns:\n"
81 " - FileId INT UNSIGNED NOT NULL\n"
82 " - EvtNumber INT UNSIGNED NOT NULL\n"
83 " - X FLOAT NOT NULL\n"
84 " - Y FLOAT NOT NULL\n"
85 "\n"
86 "If the table does not exist, it can be created automatically specifying "
87 "the --crate option. If a table already exists and it should be dropped "
88 "before creation, the --drop option can be used:\n"
89 "\n"
90 " calcsource 170909009 --create\n"
91 "\n"
92 "Process the data from the run 009 from 09/09/2017. If no table exists "
93 "a table with the name given by --table.position is created.\n"
94 "\n"
95 " calcsource 170909009 --create --drop\n"
96 "\n"
97 "Same as before, but if a table with the name given by --table.position "
98 "exists, it is dropped before.\n"
99 "\n"
100 "For each event, a new row is inserted. If existing rows should be updated, "
101 "use:\n"
102 "\n"
103 " calcsource 170909009 --updated\n"
104 "\n"
105 "The --create option is compatible with that. The --drop option is ignored.\n"
106 "\n"
107 "To avoid failure in case an entry does already exist, you can add the IGNORE "
108 "keyword to the INSERT query by --ignore-errors, which essentially ignores "
109 "all errors and turns them into warnings which are printed after the query "
110 "succeeded.\n"
111 "\n"
112 "\n"
113 "For debugging purposes several print options and options to avoid irreversible "
114 "changes to the database exist."
115 "\n\n"
116 "Usage: calcsource YYMMDDXXX [-u URI] [options]\n"
117 "\n"
118 ;
119 cout << endl;
120}
121
122class MRotation : public TRotation
123{
124public:
125 MRotation() : TRotation(1, 0, 0, 0, -1, 0, 0, 0, 1)
126 {
127 }
128};
129
130
131int main(int argc, const char* argv[])
132{
133 Time start;
134
135 gROOT->SetBatch();
136
137 Configuration conf(argv[0]);
138 conf.SetPrintUsage(PrintUsage);
139 SetupConfiguration(conf);
140
141 if (!conf.DoParse(argc, argv))
142 return 127;
143
144 // ----------------------------- Evaluate options --------------------------
145 if (conf.Has("ra")^conf.Has("dec"))
146 throw runtime_error("--ra and --dec can only be used together");
147
148 const bool has_radec = conf.Has("ra") && conf.Has("dec");
149
150 double source_ra = conf.Has("ra") ? conf.Get<double>("ra") : 0;
151 double source_dec = conf.Has("dec") ? conf.Get<double>("dec") : 0;
152
153 //string source_name = conf.Get<string>("source-name");
154 //uint32_t source_key = conf.Has("source-key") ? conf.Get<uint32_t>("source-key") : 0;
155
156 const string uri = conf.Get<string>("uri");
157 const string tab_events = conf.Get<string>("table.events");
158 const string tab_runinfo = conf.Get<string>("table.runinfo");
159 const string tab_source = conf.Get<string>("table.source");
160 const string tab_position = conf.Get<string>("table.position");
161
162 const bool ignore_errors = conf.Get<bool>("ignore-errors");
163
164 const uint32_t file = conf.Get<uint32_t>("file");
165
166 const double focal_dist = conf.Get<double>("focal-dist");
167
168 const bool print_meta = conf.Get<bool>("print-meta");
169 const bool print_insert = conf.Get<bool>("print-insert");
170
171 const bool force = conf.Get<bool>("force");
172 const bool drop = conf.Get<bool>("drop");
173 const bool create = conf.Get<bool>("create") || drop;
174 const bool update = conf.Get<bool>("update");
175 const bool noinsert = conf.Get<bool>("no-insert");
176 const bool dry_run = conf.Get<bool>("dry-run");
177 const uint16_t verbose = conf.Get<uint16_t>("verbose");
178
179 // -------------------------------------------------------------------------
180 // Checking for database connection
181
182 Database connection(uri); // Keep alive while fetching rows
183
184 try
185 {
186 if (!force)
187 connection.connected();
188 }
189 catch (const exception &e)
190 {
191 cerr << "SQL connection failed: " << e.what() << endl;
192 return 1;
193 }
194
195 // -------------------------------------------------------------------------
196 // list file-ids in the events table
197
198 if (file==0)
199 {
200 const string query =
201 "SELECT FileId FROM `"+tab_events+"` GROUP BY FileId";
202
203 if (print_meta)
204 cout << query << endl;
205
206 const mysqlpp::StoreQueryResult res =
207 connection.query(query).store();
208
209 for (size_t i=0; i<res.num_rows(); i++)
210 cout << "calcsource " << res[i][0] << '\n';
211 cout << endl;
212
213 return 0;
214 }
215
216 // -------------------------------------------------------------------------
217 // retrieve source from the database
218
219 if (verbose>0)
220 {
221 cout << "\n------------------------- Evaluating source ------------------------" << endl;
222 cout << "Requesting coordinates from " << tab_runinfo << "/" << tab_source << " table for 20" << file/1000 << "/" << file%1000 << endl;
223 }
224
225 if (!has_radec)
226 {
227 const string query =
228 "SELECT "+tab_source+".fRightAscension, "+tab_source+".fDeclination, "+tab_source+".fSourceName"
229 " FROM `"+tab_runinfo+"`"+
230 " LEFT JOIN `"+tab_source+"`"+
231 " USING (fSourceKey)"
232 " WHERE fNight=20"+to_string(file/1000)+
233 " AND fRunID="+to_string(file%1000);
234
235 if (print_meta)
236 cout << query << endl;
237
238 try
239 {
240 const mysqlpp::StoreQueryResult res =
241 connection.query(query).store();
242
243 if (res.num_rows()!=1)
244 {
245 cerr << "No coordinates from " << tab_runinfo << " for " << file << endl;
246 return 2;
247 }
248
249 source_ra = res[0][0];
250 source_dec = res[0][1];
251
252 if (verbose>0)
253 cout << "Using coordinates " << source_ra << "h / " << source_dec << " deg for '" << res[0][2] << "'" << endl;
254 }
255 catch (const exception &e)
256 {
257 cerr << query << "\n";
258 cerr << "SQL query failed:\n" << e.what() << endl;
259 return 3;
260 }
261 }
262 else
263 if (verbose>0)
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 connection.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 // -------------------------------------------------------------------------
286 // create INSERT/UPDATE query (calculate positions)
287
288 if (verbose>0)
289 {
290 cout << "\n-------------------------- Evaluating file ------------------------";
291 cout << "\nRequesting data from table `" << tab_events << "`" << endl;
292 }
293
294 const string query =
295 "SELECT `Ra`, `Dec`, MJD, MilliSec, NanoSec, Zd, Az, EvtNumber"
296 " FROM `"+tab_events+"`"
297 " WHERE FileId="+to_string(file);
298
299 if (print_meta)
300 cout << query << endl;
301
302 const mysqlpp::UseQueryResult res1 =
303 connection.query(query).use();
304
305 Nova::RaDecPosn source(source_ra, source_dec);
306
307 source_ra *= M_PI/12;
308 source_dec *= M_PI/180;
309
310 auto obs = Nova::kORM;
311
312 obs.lng *= M_PI/180;
313 obs.lat *= M_PI/180;
314
315 //const double mm2deg = 1.17193246260285378e-02;
316
317 ostringstream ins;
318 ins << setprecision(16);
319
320 vector<string> upd;
321
322 size_t count = 0;
323 while (auto row=res1.fetch_row())
324 {
325 count++;
326
327 double point_ra = row[0];
328 double point_dec = row[1];
329 uint32_t mjd = row[2];
330 int64_t millisec = row[3];
331 //uint32_t nanosec = row[4];
332 //double zd = row[5];
333 //double az = row[6];
334 uint32_t event = row[7];
335
336 Nova::RaDecPosn point(point_ra, point_dec);
337
338 point_ra *= M_PI/12;
339 point_dec *= M_PI/180;
340
341 /*
342 // ============================ Mars ================================
343
344 TVector3 pos; // pos: source position
345 TVector3 pos0; // pos: source position
346
347 pos.SetMagThetaPhi(1, M_PI/2-source_dec, source_ra);
348 pos0.SetMagThetaPhi(1, M_PI/2-point_dec, point_ra);
349
350 const double ut = (nanosec/1e6+millisec)/(24*3600000);
351
352 // Julian centuries since J2000.
353 const double t = (ut -(51544.5-mjd)) / 36525.0;
354
355 // GMST at this UT1
356 const double r1 = 24110.54841+(8640184.812866+(0.093104-6.2e-6*t)*t)*t;
357 const double r2 = 86400.0*ut;
358
359 const double sum = (r1+r2)/(3600*24);
360
361 double gmst = fmod(sum, 1) * 2*M_PI;
362
363 MRotation conv;
364 conv.RotateZ(gmst + obs.lng);
365 conv.RotateY(obs.lat-M_PI/2);
366 conv.RotateZ(M_PI);
367
368 pos *= conv;
369 pos0 *= conv;
370
371 pos.RotateZ(-pos0.Phi());
372 pos.RotateY(-pos0.Theta());
373 pos.RotateZ(-M_PI/2); // exchange x and y
374 pos *= -focal_dist/pos.Z();
375
376 TVector2 v = pos.XYvector();
377
378 //if (fDeviation)
379 // v -= fDeviation->GetDevXY()/fGeom->GetConvMm2Deg();
380
381 //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n';
382 */
383
384 // ================================= Nova ===============================
385
386 Nova::ZdAzPosn ppos = Nova::GetHrzFromEqu(source, 2400000.5+mjd+millisec/1000./3600/24);
387 Nova::ZdAzPosn ppos0 = Nova::GetHrzFromEqu(point, 2400000.5+mjd+millisec/1000./3600/24);
388
389 TVector3 pos;
390 TVector3 pos0;
391 pos.SetMagThetaPhi( 1, ppos.zd *M_PI/180, ppos.az *M_PI/180);
392 pos0.SetMagThetaPhi(1, ppos0.zd*M_PI/180, ppos0.az*M_PI/180);
393
394 pos.RotateZ(-pos0.Phi());
395 pos.RotateY(-pos0.Theta());
396 pos.RotateZ(-M_PI/2); // exchange x and y
397 pos *= -focal_dist/pos.Z();
398
399 TVector2 v = pos.XYvector();
400
401 //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n';
402
403 /*
404 // =============================== Slalib ==========================
405 const double height = 2200;
406 const double temp = 10;
407 const double hum = 0.25;
408 const double press = 780;
409
410 const double _mjd = mjd+millisec/1000./3600/24;
411
412 const double dtt = palDtt(_mjd); // 32.184 + 35
413
414 const double tdb = _mjd + dtt/3600/24;
415 const double dut = 0;
416
417 // prepare calculation: Mean Place to geocentric apperent
418 // (UTC would also do, except for the moon?)
419 double fAmprms[21];
420 palMappa(2000.0, tdb, fAmprms); // Epoche, TDB
421
422 // prepare: Apperent to observed place
423 double fAoprms[14];
424 palAoppa(_mjd, dut, // mjd, Delta UT=UT1-UTC
425 obs.lng, obs.lat, height, // long, lat, height
426 0, 0, // polar motion x, y-coordinate (radians)
427 273.155+temp, press, hum, // temp, pressure, humidity
428 0.40, 0.0065, // wavelength, tropo lapse rate
429 fAoprms);
430
431 // ---- Mean to apparent ----
432 double r=0, d=0;
433 palMapqkz(point_ra, point_dec, fAmprms, &r, &d);
434
435 double _zd, _az, ha, ra, dec;
436 // -- apparent to observed --
437 palAopqk(r, d, fAoprms,
438 &_az, // observed azimuth (radians: N=0,E=90) [-pi, pi]
439 &_zd, // observed zenith distance (radians) [-pi/2, pi/2]
440 &ha, // observed hour angle (radians)
441 &dec, // observed declination (radians)
442 &ra); // observed right ascension (radians)
443
444 //cout << _zd*180/M_PI << " " << _az*180/M_PI << endl;
445
446 pos0.SetMagThetaPhi(1, _zd, _az);
447
448 r=0, d=0;
449 palMapqkz(source_ra, source_dec, fAmprms, &r, &d);
450
451 // -- apparent to observed --
452 palAopqk(r, d, fAoprms,
453 &_az, // observed azimuth (radians: N=0,E=90) [-pi, pi]
454 &_zd, // observed zenith distance (radians) [-pi/2, pi/2]
455 &ha, // observed hour angle (radians)
456 &dec, // observed declination (radians)
457 &ra); // observed right ascension (radians)
458
459 pos.SetMagThetaPhi(1, _zd, _az);
460
461 //cout << _zd*180/M_PI << " " << _az*180/M_PI << endl;
462
463 pos.RotateZ(-pos0.Phi());
464 pos.RotateY(-pos0.Theta());
465 pos.RotateZ(-M_PI/2); // exchange x and y
466 pos *= -focal_dist/pos.Z();
467
468 v = pos.XYvector();
469
470 //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n';
471 */
472
473 if (!update)
474 ins << "( " << file << ", " << event << ", " << v.X() << ", " << v.Y() << " ),\n";
475 else
476 {
477 ostringstream out;
478 out << setprecision(16);
479 out << "UPDATE `" << tab_position << "` SET X=" << v.X() << ", Y=" << v.Y() << " WHERE FileId=" << file << " AND EvtNumber=" << event;
480 upd.emplace_back(out.str());
481 }
482 }
483
484 if (connection.errnum())
485 {
486 cerr << "SQL error fetching row: " << connection.error() << endl;
487 return 4;
488 }
489
490 if (verbose>0)
491 cout << "Processed " << count << " events.\n" << endl;
492
493 if (count==0)
494 {
495 if (verbose>0)
496 cout << "Total execution time: " << Time().UnixTime()-start.UnixTime() << "s\n" << endl;
497 return 0;
498 }
499
500 // -------------------------------------------------------------------------
501 // drop table if requested
502
503 if (drop && !update)
504 {
505 try
506 {
507 if (verbose>0)
508 cout << "Dropping table `" << tab_position << "`" << endl;
509
510 if (!dry_run)
511 connection.query("DROP TABLE `"+tab_position+"`").execute();
512
513 if (verbose>0)
514 cout << "Table `" << tab_position << "` dropped.\n" << endl;
515 }
516 catch (const exception &e)
517 {
518 cerr << "DROP TABLE `" << tab_position << "`\n\n";
519 cerr << "SQL query failed:\n" << e.what() << endl;
520 return 5;
521 }
522 }
523
524 // -------------------------------------------------------------------------
525 // crate table if requested
526
527 if (create && !update)
528 {
529 if (verbose>0)
530 cout << "Creating table `" << tab_position << "`" << endl;
531
532 const string query2 =
533 "CREATE TABLE IF NOT EXISTS `"+tab_position+"`\n"
534 "(\n"
535 " FileId INT UNSIGNED NOT NULL,\n"
536 " EvtNumber INT UNSIGNED NOT NULL,\n"
537 " X FLOAT NOT NULL,\n"
538 " Y FLOAT NOT NULL,\n"
539 " PRIMARY KEY (FileId, EvtNumber)\n"
540 ")\n"
541 "DEFAULT CHARSET=latin1 COLLATE=latin1_general_ci\n"
542 "ENGINE=MyISAM\n"
543 "COMMENT='created by "+conf.GetName()+"'\n";
544
545 try
546 {
547 if (!dry_run)
548 connection.query(query2).execute();
549 }
550 catch (const exception &e)
551 {
552 cerr << query2 << "\n\n";
553 cerr << "SQL query failed:\n" << e.what() << endl;
554 return 6;
555 }
556
557 if (print_meta)
558 cout << query2 << endl;
559
560 if (verbose>0)
561 cout << "Table `" << tab_position << "` created.\n" << endl;
562 }
563
564 // -------------------------------------------------------------------------
565 // delete old entries from table
566
567 if (!drop && !update)
568 {
569 if (verbose>0)
570 cout << "Deleting old entries from table `" << tab_position << "`" << endl;
571
572 const string query2 =
573 "DELETE FROM `"+tab_position+"` WHERE FileId="+to_string(file);
574
575 try
576 {
577 if (!dry_run)
578 {
579 const mysqlpp::SimpleResult res =
580 connection.query(query2).execute();
581
582 if (verbose>0)
583 cout << res.rows() << " row(s) affected.\n" << endl;
584 }
585 }
586 catch (const exception &e)
587 {
588 cerr << query2 << "\n\n";
589 cerr << "SQL query failed:\n" << e.what() << endl;
590 return 7;
591 }
592
593 if (print_meta)
594 cout << query2 << endl;
595 }
596
597 // -------------------------------------------------------------------------
598 // insert data into table
599
600 if (!update)
601 {
602 if (verbose>0)
603 cout << "Inserting data into table " << tab_position << "." << endl;
604
605 string query2 = "INSERT ";
606 if (ignore_errors)
607 query2 += "IGNORE ";
608 query2 += "`"+tab_position+"` (FileId, EvtNumber, X, Y) VALUES\n"+
609 ins.str().substr(0, ins.str().size()-2)+
610 "\n";
611
612 try
613 {
614 if (!noinsert)
615 {
616 const mysqlpp::SimpleResult res =
617 connection.query(query2).execute();
618
619 if (verbose>0)
620 cout << res.info() << '\n' << endl;
621 }
622 }
623 catch (const exception &e)
624 {
625 cerr << query2 << "\n\n";
626 cerr << "SQL query (" << query2.length() << " bytes) failed:\n" << e.what() << endl;
627 return 8;
628 }
629
630 if (print_insert)
631 cout << query2 << endl;
632 }
633
634
635
636 // -------------------------------------------------------------------------
637 // update data in table
638
639 if (update)
640 {
641 if (verbose>0)
642 cout << "Updating data in table " << tab_position << "." << endl;
643
644 size_t cnt = 0;
645 for (const auto &str : upd)
646 {
647 try
648 {
649 if (!noinsert && !dry_run)
650 {
651 const mysqlpp::SimpleResult res =
652 connection.query(str).execute();
653
654 cnt += res.rows();
655 }
656 }
657 catch (const exception &e)
658 {
659 cerr << str << "\n\n";
660 cerr << "SQL query failed:\n" << e.what() << endl;
661 return 9;
662 }
663 }
664
665 if (verbose>0)
666 cout << cnt << " row(s) out of " << upd.size() << " affected.\n" << endl;
667
668 if (print_insert)
669 {
670 for (const auto &str : upd)
671 cout << str << '\n';
672 cout << endl;
673 }
674 }
675
676 if (verbose>0)
677 {
678 cout << "Total execution time: " << Time().UnixTime()-start.UnixTime() << "s\n" << endl;
679
680 try
681 {
682 const auto resw =
683 connection.query("SHOW WARNINGS").store();
684
685 for (size_t i=0; i<resw.num_rows(); i++)
686 {
687 const mysqlpp::Row &roww = resw[i];
688
689 cout << roww["Level"] << '[' << roww["Code"] << "]: ";
690 cout << roww["Message"] << '\n';
691 }
692 cout << endl;
693
694 }
695 catch (const exception &e)
696 {
697 cerr << "\nSHOW WARNINGS\n\n";
698 cerr << "SQL query failed:\n" << e.what() << '\n' <<endl;
699 return 8;
700 }
701 }
702
703 return 0;
704}
705
706/*
707 TRotation & TRotation::RotateX(Double_t a) {
708 //rotate around x
709 Double_t c = TMath::Cos(a);
710 Double_t s = TMath::Sin(a);
711 Double_t x = fyx, y = fyy, z = fyz;
712 fyx = c*x - s*fzx;
713 fyy = c*y - s*fzy;
714 fyz = c*z - s*fzz;
715 fzx = s*x + c*fzx;
716 fzy = s*y + c*fzy;
717 fzz = s*z + c*fzz;
718 return *this;
719}
720
721TRotation & TRotation::RotateY(Double_t a){
722 //rotate around y
723 Double_t c = TMath::Cos(a);
724 Double_t s = TMath::Sin(a);
725 Double_t x = fzx, y = fzy, z = fzz;
726 fzx = c*x - s*fxx;
727 fzy = c*y - s*fxy;
728 fzz = c*z - s*fxz;
729 fxx = s*x + c*fxx;
730 fxy = s*y + c*fxy;
731 fxz = s*z + c*fxz;
732 return *this;
733}
734
735TRotation & TRotation::RotateZ(Double_t a) {
736 //rotate around z
737 Double_t c = TMath::Cos(a);
738 Double_t s = TMath::Sin(a);
739 Double_t x = fxx, y = fxy, z = fxz;
740 fxx = c*x - s*fyx;
741 fxy = c*y - s*fyy;
742 fxz = c*z - s*fyz;
743 fyx = s*x + c*fyx;
744 fyy = s*y + c*fyy;
745 fyz = s*z + c*fyz;
746 return *this;
747}
748
749 */
750
751
752// Reuired:
753// pointing_ra[8]
754// pointing_dec[8]
755// fNanoSec[4], fTime[fMiliSec,8], fMjd[4]
756// source_ra -> rc
757// source_dec -> rc
758// fLong, fLat -> rc
759// fCameraDist -> rc
760// 32 byte per row, 1 Monat = 4mio rows => 120 MB
761
762/*
763
764void TVector3::SetThetaPhi(Double_t theta, Double_t phi)
765{
766 //setter with mag, theta, phi
767 fX = TMath::Sin(theta) * TMath::Cos(phi);
768 fY = TMath::Sin(theta) * TMath::Sin(phi);
769 fZ = TMath::Cos(theta);
770}
771
772
773 // Set Sky coordinates of source, taken from container "MSourcePos"
774 // of type MPointingPos. The sky coordinates must be J2000, as the
775 // sky coordinates of the camera center that we get from the container
776 // "MPointingPos" filled by the Drive.
777 //MVector3 pos; // pos: source position
778 //pos.SetRaDec(fSourcePos->GetRaRad(), fSourcePos->GetDecRad());
779
780 TVector3 pos; // pos: source position
781 pos.SetMagThetaPhi(1, TMath::Pi()/2-source_dec, source_ra);
782
783 // Set sky coordinates of camera center in pos0 (for details see below)
784 //MVector3 pos0; // pos0: camera center
785 //pos0.SetRaDec(fPointPos->GetRaRad(), fPointPos->GetDecRad());
786
787 TVector3 pos0; // pos: source position
788 pos0.SetMagThetaPhi(1, TMath::Pi()/2-pointing_dec, pointing_ra);
789
790 // Convert sky coordinates of source to local coordinates. Warning! These are not the "true" local
791 // coordinates, since this transformation ignores precession and nutation effects.
792
793 //const MAstroSky2Local conv(*fTime, *fObservatory);
794 //pos *= conv;
795
796 const Double_t ut = (Double_t)(fNanoSec/1e6+(Long_t)fTime)/kDay;
797
798 // Julian centuries since J2000.
799 const Double_t t = (ut -(51544.5-fMjd)) / 36525.0;
800
801 // GMST at this UT1
802 const Double_t r1 = 24110.54841+(8640184.812866+(0.093104-6.2e-6*t)*t)*t;
803 const Double_t r2 = 86400.0*ut;
804
805 const Double_t sum = (r1+r2)/kDaySec;
806
807 double gmst = fmod(sum, 1)*TMath::TwoPi();
808
809 // TRotation::TRotation(Double_t mxx, Double_t mxy, Double_t mxz,
810 // Double_t myx, Double_t myy, Double_t myz,
811 // Double_t mzx, Double_t mzy, Double_t mzz)
812 // : fxx(mxx), fxy(mxy), fxz(mxz),
813 // fyx(myx), fyy(myy), fyz(myz),
814 // fzx(mzx), fzy(mzy), fzz(mzz) {}
815
816 TRotation conv(1, 0, 0, 0, -1, 0, 0, 0, 1);
817 conv.RotateZ(gmst + obs.GetElong());
818 conv.RotateY(obs.GetPhi()-TMath::Pi()/2);
819 conv.RotateZ(TMath::Pi());
820
821 // TVector3 & TVector3::operator *= (const TRotation & m){
822 // return *this = m * (*this);
823 // }
824
825 // inline TVector3 TRotation::operator * (const TVector3 & p) const {
826 // return TVector3(fxx*p.X() + fxy*p.Y() + fxz*p.Z(),
827 // fyx*p.X() + fyy*p.Y() + fyz*p.Z(),
828 // fzx*p.X() + fzy*p.Y() + fzz*p.Z());
829 // }
830
831 // Convert sky coordinates of camera center convert to local.
832 // Same comment as above. These coordinates differ from the true
833 // local coordinates of the camera center that one could get from
834 // "MPointingPos", calculated by the Drive: the reason is that the Drive
835 // takes into account precession and nutation corrections, while
836 // MAstroSky2Local (as of Jan 27 2005 at least) does not. Since we just
837 // want to get the source position on the camera from the local
838 // coordinates of the center and the source, it does not matter that
839 // the coordinates contained in pos and pos0 ignore precession and
840 // nutation... since the shift would be the same in both cases. What
841 // would be wrong is to set in pos0 directly the local coordinates
842 // found in MPointingPos!
843 pos0 *= conv;
844
845 if (fDeviation)
846 {
847 // Position at which the starguider camera is pointing in real:
848 // pointing position = nominal position - dev
849 //
850 //vx = CalcXYinCamera(pos0, pos)*fGeom->GetCameraDist()*1000;
851 pos0.SetZdAz(pos0.Theta()-fDeviation->GetDevZdRad(),
852 pos0.Phi() -fDeviation->GetDevAzRad());
853 }
854
855 // Calculate source position in camera, and convert to mm:
856 TVector2 v = MAstro::GetDistOnPlain(pos0, pos, -fGeom->GetCameraDist()*1000);
857
858 pos.RotateZ(-pos0.Phi());
859 pos.RotateY(-pos0.Theta());
860 pos.RotateZ(-TMath::Pi()/2); // exchange x and y
861 pos *= -fGeom->GetCameraDist()*1000/pos.Z();
862
863 TVector2 v = pos.XYvector();
864
865 if (fDeviation)
866 v -= fDeviation->GetDevXY()/fGeom->GetConvMm2Deg();
867
868 SetSrcPos(v);
869
870 // ================================================
871
872 if (fMode==kWobble)
873 {
874 // The trick here is that the anti-source position in case
875 // of the off-source regions is always the on-source positon
876 // thus a proper anti-source cut is possible.
877 fSrcPosAnti->SetXY(v);
878 if (fCallback)
879 v = Rotate(v, fCallback->GetNumPass(), fCallback->GetNumPasses());
880 else
881 v *= -1;
882 fSrcPosCam->SetXY(v);
883 }
884 else
885 {
886 // Because we don't process this three times like in the
887 // wobble case we have to find another way to determine which
888 // off-source region is the right anti-source position
889 // We do it randomly.
890 fSrcPosCam->SetXY(v);
891 if (fNumRandomOffPositions>1)
892 v = Rotate(v, gRandom->Integer(fNumRandomOffPositions), fNumRandomOffPositions);
893 else
894 v *= -1;
895 fSrcPosAnti->SetXY(v);
896 }
897
898
899
900 */
901
902
903
904
905/*
906 PointingData CalcPointingPos(const PointingSetup &setup, double _mjd, const Weather &weather, uint16_t timeout, bool tpoint=false)
907 {
908 PointingData out;
909 out.mjd = _mjd;
910
911 const double elong = Nova::kORM.lng * M_PI/180;
912 const double lat = Nova::kORM.lat * M_PI/180;
913 const double height = 2200;
914
915 const bool valid = weather.time+boost::posix_time::seconds(timeout) > Time();
916
917 const double temp = valid ? weather.temp : 10;
918 const double hum = valid ? weather.hum : 0.25;
919 const double press = valid ? weather.press : 780;
920
921 const double dtt = palDtt(_mjd); // 32.184 + 35
922
923 const double tdb = _mjd + dtt/3600/24;
924 const double dut = 0;
925
926 // prepare calculation: Mean Place to geocentric apperent
927 // (UTC would also do, except for the moon?)
928 double fAmprms[21];
929 palMappa(2000.0, tdb, fAmprms); // Epoche, TDB
930
931 // prepare: Apperent to observed place
932 double fAoprms[14];
933 palAoppa(_mjd, dut, // mjd, Delta UT=UT1-UTC
934 elong, lat, height, // long, lat, height
935 0, 0, // polar motion x, y-coordinate (radians)
936 273.155+temp, press, hum, // temp, pressure, humidity
937 0.40, 0.0065, // wavelength, tropo lapse rate
938 fAoprms);
939
940 out.source.ra = setup.source.ra * M_PI/ 12;
941 out.source.dec = setup.source.dec * M_PI/180;
942
943 if (setup.planet!=kENone)
944 {
945 // coordinates of planet: topocentric, equatorial, J2000
946 // One can use TT instead of TDB for all planets (except the moon?)
947 double ra, dec, diam;
948 palRdplan(tdb, setup.planet, elong, lat, &ra, &dec, &diam);
949
950 // ---- apparent to mean ----
951 palAmpqk(ra, dec, fAmprms, &out.source.ra, &out.source.dec);
952 }
953
954 if (setup.wobble_offset<=0 || tpoint)
955 {
956 out.pointing.dec = out.source.dec;
957 out.pointing.ra = out.source.ra;
958 }
959 else
960 {
961 const double dphi =
962 setup.orbit_period==0 ? 0 : 2*M_PI*(_mjd-setup.start)/setup.orbit_period;
963
964 const double phi = setup.wobble_angle + dphi;
965
966 const double cosdir = cos(phi);
967 const double sindir = sin(phi);
968 const double cosoff = cos(setup.wobble_offset);
969 const double sinoff = sin(setup.wobble_offset);
970 const double cosdec = cos(out.source.dec);
971 const double sindec = sin(out.source.dec);
972
973 const double sintheta = sindec*cosoff + cosdec*sinoff*cosdir;
974
975 const double costheta = sintheta>1 ? 0 : sqrt(1 - sintheta*sintheta);
976
977 const double cosdeltara = (cosoff - sindec*sintheta)/(cosdec*costheta);
978 const double sindeltara = sindir*sinoff/costheta;
979
980 out.pointing.dec = asin(sintheta);
981 out.pointing.ra = atan2(sindeltara, cosdeltara) + out.source.ra;
982 }
983
984 // ---- Mean to apparent ----
985 double r=0, d=0;
986 palMapqkz(out.pointing.ra, out.pointing.dec, fAmprms, &r, &d);
987
988 //
989 // Doesn't work - don't know why
990 //
991 // slaMapqk (radec.Ra(), radec.Dec(), rdpm.Ra(), rdpm.Dec(),
992 // 0, 0, (double*)fAmprms, &r, &d);
993 //
994
995 // -- apparent to observed --
996 palAopqk(r, d, fAoprms,
997 &out.sky.az, // observed azimuth (radians: N=0,E=90) [-pi, pi]
998 &out.sky.zd, // observed zenith distance (radians) [-pi/2, pi/2]
999 &out.apparent.ha, // observed hour angle (radians)
1000 &out.apparent.dec, // observed declination (radians)
1001 &out.apparent.ra); // observed right ascension (radians)
1002
1003 // ----- fix ambiguity -----
1004 if (out.sky.zd<0)
1005 {
1006 out.sky.zd = -out.sky.zd;
1007 out.sky.az += out.sky.az<0 ? M_PI : -M_PI;
1008 }
1009
1010 // Star culminating behind zenith and Az between ~90 and ~180deg
1011 if (out.source.dec<lat && out.sky.az>0)
1012 out.sky.az -= 2*M_PI;
1013
1014 out.mount = SkyToMount(out.sky);
1015
1016 return out;
1017 }
1018};
1019
1020*/
1021
Note: See TracBrowser for help on using the repository browser.