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

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