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 |
|
---|
12 | using namespace std;
|
---|
13 |
|
---|
14 | // ------------------------------------------------------------------------
|
---|
15 |
|
---|
16 | void 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 |
|
---|
59 | void 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 |
|
---|
122 | class MRotation : public TRotation
|
---|
123 | {
|
---|
124 | public:
|
---|
125 | MRotation() : TRotation(1, 0, 0, 0, -1, 0, 0, 0, 1)
|
---|
126 | {
|
---|
127 | }
|
---|
128 | };
|
---|
129 |
|
---|
130 |
|
---|
131 | int 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 |
|
---|
721 | TRotation & 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 |
|
---|
735 | TRotation & 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 |
|
---|
764 | void 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 |
|
---|