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

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