Swarm-NG  1.1
query.cpp
Go to the documentation of this file.
1 /*************************************************************************
2  * Copyright (C) 2011 by Saleh Dindar and the Swarm-NG Development Team *
3  * *
4  * This program is free software; you can redistribute it and/or modify *
5  * it under the terms of the GNU General Public License as published by *
6  * the Free Software Foundation; either version 3 of the License. *
7  * *
8  * This program is distributed in the hope that it will be useful, *
9  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
11  * GNU General Public License for more details. *
12  * *
13  * You should have received a copy of the GNU General Public License *
14  * along with this program; if not, write to the *
15  * Free Software Foundation, Inc., *
16  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
17  ************************************************************************/
18 
31 #include "query.hpp"
32 #include "kepler.h"
33 #include "bdb_query.hpp"
34 
35 
36 namespace swarm { namespace query {
37 
38 using log::body;
39 
40 
41 void get_Tsys(gpulog::logrecord &lr, double &T, int &sys)
42 {
43  //std::cerr << "msgid=" << lr.msgid() << "\n";
44  if(lr.msgid() < 0)
45  {
46  // system-defined events that have no (T,sys) heading
47  T = -1; sys = -1;
48  }
49  else
50  {
51  lr >> T >> sys;
52  }
53 }
54 
55 // Default output, if no handler is registered
56  std::ostream& record_output_default(std::ostream &out, gpulog::logrecord &lr)
57 {
58  double T; int sys;
59  get_Tsys(lr, T, sys);
60  out << lr.msgid() << "\t" << T << "\t" << sys;
61  return out;
62 }
63 
64 bool keplerian_output = false;
65 planets_coordinate_system_t planets_coordinate_system = jacobi;
66 
67 void set_keplerian_output(const planets_coordinate_system_t& coordinate_system)
68 { keplerian_output = true; planets_coordinate_system = coordinate_system; }
69 
70 void set_cartesian_output(const planets_coordinate_system_t& coordinate_system)
71 { keplerian_output = false; planets_coordinate_system = coordinate_system; }
72 
73 void set_coordinate_system(const planets_coordinate_system_t& coordinate_system)
74 { planets_coordinate_system = coordinate_system; }
75 
76 
77 struct keplerian_t {
78  double a, e, i, O, w , M;
79 };
80 
81 keplerian_t keplerian_for_cartesian ( const body& b, const body& center )
82 {
83  keplerian_t k;
84  double x = b.x - center.x;
85  double y = b.y - center.y;
86  double z = b.z - center.z;
87  double vx = b.vx - center.vx;
88  double vy = b.vy - center.vy;
89  double vz = b.vz - center.vz;
90  double mass = b.mass + center.mass;
91  calc_keplerian_for_cartesian(k.a, k.e, k.i, k.O, k.w, k.M, x, y,z, vx, vy, vz, mass );
92  return k;
93 }
94 
95 body center_of_mass(const body* bodies, const int nbod ){
96  body center;
97  center.x = center.y = center.z = center.vx = center.vy = center.vz = 0.;
98  center.mass = 0.;
99  int i=0;
100  while(i<nbod)
101  {
102  center.x += bodies[i].x*bodies[i].mass;
103  center.y += bodies[i].y*bodies[i].mass;
104  center.z += bodies[i].z*bodies[i].mass;
105  center.vx += bodies[i].vx*bodies[i].mass;
106  center.vy += bodies[i].vy*bodies[i].mass;
107  center.vz += bodies[i].vz*bodies[i].mass;
108  center.mass += bodies[i].mass;
109  i++;
110  }
111  center.x /= center.mass;
112  center.y /= center.mass;
113  center.z /= center.mass;
114  center.vx /= center.mass;
115  center.vy /= center.mass;
116  center.vz /= center.mass;
117 
118  return center;
119 }
120 
121 // EVT_SNAPSHOT
122  std::ostream& record_output_1(std::ostream &out, gpulog::logrecord &lr, body_range_t &body_range)
123 {
124  double time;
125  int nbod, sys, flags;
126  const body *bodies;
127  lr >> time >> sys >> flags >> nbod >> bodies;
128 
129  body center;
130  const body &star = bodies[0];
131 
132  switch(planets_coordinate_system)
133  {
134  case astrocentric:
135  center = star;
136  break;
137  case barycentric:
138  center = center_of_mass( bodies, nbod );
139  break;
140  case jacobi:
141  center = star;
142  break;
143  case origin:
144  center.x = center.y = center.z = center.vx = center.vy = center.vz = 0.; center.mass = star.mass;
145  break;
146  };
147 
148  if((time<=0.) && false) // was used for debugging at some point
149  {
150  if (keplerian_output)
151  { std::cerr << "# Output in Keplerian coordinates "; }
152  else { std::cerr << "# Output in Cartesian coordinates "; }
153  switch(planets_coordinate_system)
154  {
155  case astrocentric:
156  std::cerr << "(astrocentric) " << center.x << ' ' << center.y << ' '<< center.z << " " << center.vx << ' ' << center.vy << ' ' << center.vz;
157  break;
158  case barycentric:
159  std::cerr << "(barycentric) "<< center.x << ' ' << center.y << ' '<< center.z << " " << center.vx << ' ' << center.vy << ' ' << center.vz;
160  break;
161  case jacobi:
162  std::cerr << "(jacobi) " << center.x << ' ' << center.y << ' '<< center.z << " " << center.vx << ' ' << center.vy << ' ' << center.vz;
163  break;
164  case origin:
165  std::cerr << "(origin) " << center.x << ' ' << center.y << ' '<< center.z << " " << center.vx << ' ' << center.vy << ' ' << center.vz;
166  break;
167  }
168  std::cerr << "\n";
169  }
170 
171 
172  size_t bufsize = 1000;
173  char buf[bufsize];
174  for(int bod = 0; bod < nbod; bod++)
175  {
176  if(!body_range.in(bod)) { continue; }
177  const body &b = bodies[bod];
178  if( keplerian_output && (bod==0) ) { continue; }
179  if( (planets_coordinate_system==jacobi) && (bod==0) ) { continue; }
180  if( ( !keplerian_output && (planets_coordinate_system!=jacobi) && (bod> 0) ) ||
181  ( !keplerian_output && (planets_coordinate_system==jacobi) && (bod> 1) ) ||
182  ( keplerian_output && (bod> 1 ) ) ){ out << "\n"; }
183 
184  if( planets_coordinate_system == jacobi )
185  { center = center_of_mass ( bodies, bod); }
186 
187  if( keplerian_output && bod > 0)
188  {
189  if(planets_coordinate_system==barycentric) center.mass -= b.mass;
190  keplerian_t orbit = keplerian_for_cartesian( b, center );
191  if(planets_coordinate_system==barycentric) center.mass += b.mass;
192  const double rad2deg = 180./M_PI;
193  snprintf(buf, bufsize, "%10d %lg %6d %6d %lg % 9.5lg % 9.5lg % 9.5lg % 9.5lg % 9.5lg % 9.5lg %d", lr.msgid(), time, sys, bod, b.mass, orbit.a, orbit.e , orbit.i*rad2deg, orbit.O*rad2deg, orbit.w *rad2deg, orbit.M*rad2deg, flags);
194  }
195  if(!keplerian_output)
196  {
197  double x = b.x - center.x;
198  double y = b.y - center.y;
199  double z = b.z - center.z;
200  double vx= b.vx- center.vx;
201  double vy= b.vy- center.vy;
202  double vz= b.vz- center.vz;
203  snprintf(buf, bufsize, "%10d %lg %6d %6d %lg %9.5lg %9.5lg %9.5lg %9.5lg %9.5lg %9.5lg %d", lr.msgid(), time, sys, bod, b.mass, x, y, z, vx, vy, vz, flags);
204  }
205  out << buf; // << "\n";
206  }
207  return out;
208 }
209 
210 
211 // EVT_EJECTION
212 // WARNING: Logging for ejections not yet fully implemented
213  std::ostream& record_output_2(std::ostream &out, gpulog::logrecord &lr, body_range_t &body_range)
214 {
215  double T;
216  int sys, body_id = 0;
217  body b;
218  lr >> T >> sys >> b;
219 
220  if(!body_range.in(b.body_id)) return out;
221 
222  size_t bufsize = 1000;
223  char buf[bufsize];
224 
225  if( keplerian_output)
226  {
227  double mass = 1. + b.mass; // + center.mass;
228  keplerian_t orbit;
229  calc_keplerian_for_cartesian(orbit.a, orbit.e, orbit.i, orbit.O, orbit.w, orbit.M, b.x, b.y,b.z,b.vx,b.vy,b.vz,mass );
230  const double rad2deg = 180./M_PI;
231  snprintf(buf, bufsize, "%10d %lg %6d %6d %lg % 9.5lg % 9.5lg % 9.5lg % 9.5lg % 9.5lg % 9.5lg ***WARNING: NOT ACCURATE***", lr.msgid(), T, sys, body_id, b.mass, orbit.a, orbit.e , orbit.i*rad2deg, orbit.O*rad2deg, orbit.w *rad2deg, orbit.M*rad2deg);
232  }
233  else
234  {
235  snprintf(buf, bufsize, "%10d %lg %6d %6d %lg % 9.5lf % 9.5lf % 9.5lf % 9.5lf % 9.5lf % 9.5lf", lr.msgid(), T, sys, body_id, b.mass, b.x, b.y, b.z, b.vx, b.vy, b.vz);
236  }
237  out << buf;
238 
239  return out;
240 }
241 
242 
243 // EVT_TRANSIT
244 std::ostream& record_output_15(std::ostream &out, gpulog::logrecord &lr, body_range_t &body_range)
245 {
246  double T;
247  int sys, body_id;
248  double minb, vproj;
249  lr >> T >> sys >> body_id >> minb >> vproj;
250  if(!body_range.in(body_id)) return out;
251 
252  size_t bufsize = 1000;
253  char buf[bufsize];
254  snprintf(buf, bufsize, "%10d %lg %6d %6d %lg %lg Transit Experimental", lr.msgid(), T, sys, body_id, minb, vproj);
255  out << buf;
256 
257  return out;
258 }
259 
260 
261 // EVT_RV_OBS
262 std::ostream& record_output_11(std::ostream &out, gpulog::logrecord &lr, body_range_t &body_range)
263 {
264  double T;
265  int sys, body_id;
266  double vz;
267  lr >> T >> sys >> body_id >> vz;
268  if(!body_range.in(body_id)) return out;
269 
270  size_t bufsize = 1000;
271  char buf[bufsize];
272  snprintf(buf, bufsize, "%10d %lg %6d %6d %lg Rv Experimental", lr.msgid(), T, sys, body_id, vz);
273  out << buf;
274 
275  return out;
276 }
277 
278 
279 // EVT_OCCULTATION
280 std::ostream& record_output_16(std::ostream &out, gpulog::logrecord &lr, body_range_t &body_range)
281 {
282  double T;
283  int sys, body_id;
284  double minb, vproj;
285  lr >> T >> sys >> body_id >> minb >> vproj;
286  if(!body_range.in(body_id)) return out;
287 
288  size_t bufsize = 1000;
289  char buf[bufsize];
290  snprintf(buf, bufsize, "%10d %lg %6d %6d %lg %lg Occultation Experimental", lr.msgid(), T, sys, body_id, minb, vproj);
291  out << buf;
292 
293  return out;
294 }
295 
296 
297  std::ostream &output_record(std::ostream &out, gpulog::logrecord &lr, const body_range_t &_bod)
298 {
299  body_range_t bod = _bod;
300  int evtid = lr.msgid();
301 
302  // Make sure these stay synced with src/swarm/log/log.hpp
303  switch(evtid){
304  case 1: // standard system snapshot
305  return record_output_1(out,lr,bod);
306  case 2: // data one body upon ejection
307  return record_output_2(out,lr,bod);
308  case 3: // reserved for data for one pair of bodies upon close encounter/collision
309  return record_output_default(out,lr); // feature still missing
310  case 11: // star v_z at observation time
311  return record_output_11(out,lr,bod);
312  case 15: // near a transit of planet in front of star
313  return record_output_15(out,lr,bod);
314  case 16: // near an occultation of star in front of planet
315  return record_output_16(out,lr,bod);
316  default:
317  return record_output_default(out,lr);
318  }
319 }
320 
321 
322  // void execute(const std::string &datafile, time_range_t T, sys_range_t sys)
323 void execute_binary_query(const std::string &datafile, time_range_t T, sys_range_t sys, body_range_t bod) {
324  swarmdb db(datafile);
325  swarmdb::result r = db.query(sys, T);
326  // swarmdb::result r = db.query(sys, bod, T);
327  gpulog::logrecord lr;
328  while(lr = r.next())
329  {
330  output_record(std::cout, lr,bod );
331  std::cout << "\n";
332  }
333 }
334 
335 void execute(const std::string &datafile, time_range_t T, sys_range_t sys, body_range_t bod)
336 {
337  if(datafile.substr(datafile.length()-3,datafile.length()) == ".db"){
338  execute_bdb_query(datafile,T,sys,bod);
339  } else {
340  execute_binary_query(datafile,T,sys,bod);
341  }
342 }
343 
344  } } // end namespace swarm::query