Swarm-NG  1.1
ensemble.hpp
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 
25 #pragma once
26 
27 #include "allocators.hpp"
28 #include "coalescedstructarray.hpp"
29 #include <config.h>
30 
31 namespace swarm {
32 
33 template<class N>
34 
36 GENERIC N square(const N& x) { return x*x; }
37 
38 
43 template<typename _value_type, int _N, int _CHUNK_SIZE>
44 struct CoalescedMemberArray {
45  typedef _value_type value_type;
46  const static int N = _N;
47  const static int CHUNK_SIZE = _CHUNK_SIZE;
48 
49 
50  GENERIC value_type getitem(const size_t& i)const{
51  return _array[i][0];
52  }
53  GENERIC void setitem(const size_t& i, const value_type& v){
54  _array[i][0] = v ;
55  }
56 
57  GENERIC const value_type& operator[](const size_t& i)const{
58  return _array[i][0];
59  }
60  GENERIC value_type& operator[](const size_t& i){
61  return _array[i][0];
62  }
63 
64 private:
65  value_type _array[N][CHUNK_SIZE];
66 
67 };
68 
69 
109 template< int _CHUNK_SIZE, int _NUM_BODY_ATTRIBUTES = NUM_PLANET_ATTRIBUTES, int _NUM_SYS_ATTRIBUTES = NUM_SYSTEM_ATTRIBUTES >
110 class EnsembleBase {
111  public:
112 
114  static const int CHUNK_SIZE = _CHUNK_SIZE;
115 
117  static const int NUM_BODY_ATTRIBUTES = _NUM_BODY_ATTRIBUTES;
118 
120  static const int NUM_SYS_ATTRIBUTES = _NUM_SYS_ATTRIBUTES;
121 
131  struct Body {
132  struct Component {
134  double _pos[CHUNK_SIZE];
136  double _vel[CHUNK_SIZE];
138  GENERIC double& pos() { return _pos[0]; }
140  GENERIC double& vel() { return _vel[0]; }
142  GENERIC const double& pos() const { return _pos[0]; }
144  GENERIC const double& vel() const { return _vel[0]; }
145  } component[3];
146 
147  double _mass[CHUNK_SIZE];
148  typedef CoalescedMemberArray<double,NUM_BODY_ATTRIBUTES,CHUNK_SIZE>
149  attributes_t;
150  attributes_t _attributes;
151 
152  attributes_t& attributes(){ return _attributes; }
153 
155  GENERIC double& mass() { return _mass[0]; }
156 
158  GENERIC const double& mass() const { return _mass[0]; }
159 
163  GENERIC Component& operator[] (const int & i) { return component[i]; };
164 
166  GENERIC const Component& operator[] (const int & i) const { return component[i]; };
167 
168  GENERIC Component& getitem(const int& i){
169  return operator[](i);
170  }
171 
172  GENERIC const double& x() const { return component[0].pos(); }
173  GENERIC double& x() { return component[0].pos(); }
174  GENERIC const double& y() const { return component[1].pos(); }
175  GENERIC double& y() { return component[1].pos(); }
176  GENERIC const double& z() const { return component[2].pos(); }
177  GENERIC double& z() { return component[2].pos(); }
178 
179  GENERIC const double& vx() const { return component[0].vel(); }
180  GENERIC double& vx() { return component[0].vel(); }
181  GENERIC const double& vy() const { return component[1].vel(); }
182  GENERIC double& vy() { return component[1].vel(); }
183  GENERIC const double& vz() const { return component[2].vel(); }
184  GENERIC double& vz() { return component[2].vel(); }
185 
187  GENERIC double& attribute(const int& i) { return _attributes[i]; }
188  GENERIC const double& attribute(const int& i) const { return _attributes[i]; }
189  GENERIC int num_attributes() const { return NUM_BODY_ATTRIBUTES; };
190 
192  GENERIC double distance_to_origin_squared() {
193  return square(operator[](0).pos())
194  + square(operator[](1).pos())
195  + square(operator[](2).pos());
196  }
197 
199  GENERIC double speed_squared() {
200  return square(operator[](0).vel())
201  + square(operator[](1).vel())
202  + square(operator[](2).vel());
203  }
204 
205  GENERIC double distance_to_origin() { return sqrt(distance_to_origin_squared()); }
206  GENERIC double speed() { return sqrt(speed_squared()); }
207 
209  GENERIC void get(double& x,double & y, double & z
210  , double & vx, double & vy, double & vz) {
211  x = operator[](0).pos();
212  y = operator[](1).pos();
213  z = operator[](2).pos();
214  vx = operator[](0).vel();
215  vy = operator[](1).vel();
216  vz = operator[](2).vel();
217  }
218 
219  };
220 
227  struct Sys {
228  double _time[CHUNK_SIZE];
229 
230  struct {
231  int state; int id;
232  } _bunch[CHUNK_SIZE];
233 
234  typedef CoalescedMemberArray<double,NUM_SYS_ATTRIBUTES,CHUNK_SIZE>
235  attributes_t;
236  attributes_t _attributes;
237 
238  attributes_t& attributes(){ return _attributes; }
239 
241  GENERIC double& time() { return _time[0]; }
242  GENERIC const double& time() const { return _time[0]; }
243 
245  GENERIC int& state() { return _bunch[0].state; }
246  GENERIC const int& state()const { return _bunch[0].state; }
247 
249  GENERIC int& id() { return _bunch[0].id; }
250  GENERIC const int& id() const { return _bunch[0].id; }
251 
253  GENERIC const double& attribute(const int& i) const { return _attributes[i]; }
254  GENERIC double& attribute(const int& i) { return _attributes[i]; }
255  GENERIC int num_attributes() const { return NUM_SYS_ATTRIBUTES; };
256 
269  enum ActivationStates {
270  SYSTEM_ACTIVE = 0,
271  SYSTEM_INACTIVE = 1,
272  SYSTEM_NEEDS_EXAM = 2,
273  SYSTEM_DISABLED = -1
274  };
275  };
276 
277 
278 
279 
287  class SystemRef {
289  const int _nbod;
291  const int _number;
293  Body* _body;
295  Sys* _sys;
296  public:
297  typedef SystemRef second_type;
298 
300  GENERIC SystemRef(const int& nbod,const int& number,Body* body,Sys* sys):_nbod(nbod),_number(number),_body(body),_sys(sys){}
301 
302 
304  GENERIC Body& operator[](const int & i ) const { return _body[i]; };
305 
306  GENERIC Body& getitem(const int& i){
307  return operator[](i);
308  }
309 
311  GENERIC typename Sys::attributes_t& attributes() const { return _sys[0].attributes(); }
312 
314  GENERIC double& time() const { return _sys[0].time(); }
315 
317  GENERIC bool is_active() const { return _sys[0].state() == Sys::SYSTEM_ACTIVE; }
318 
320  GENERIC bool is_inactive() const { return _sys[0].state() == Sys::SYSTEM_INACTIVE; }
321 
323  GENERIC bool is_disabled() const { return _sys[0].state() == Sys::SYSTEM_DISABLED; }
324 
326  GENERIC bool is_enabled() const { return !is_disabled(); }
327 
329  GENERIC int& state() const { return _sys[0].state(); }
330 
332  GENERIC void set_active() const { _sys[0].state() = Sys::SYSTEM_ACTIVE; }
333 
335  GENERIC void set_inactive() const { _sys[0].state() = Sys::SYSTEM_INACTIVE; }
336 
338  GENERIC void set_disabled() const { _sys[0].state() = Sys::SYSTEM_DISABLED; }
339 
347  GENERIC int& id() const { return _sys[0].id(); }
349  GENERIC const int& nbod() const{ return _nbod; }
351  GENERIC const int& number() const{ return _number; }
352 
354  GENERIC double& attribute(const int& i) const { return _sys[0].attribute(i); }
355  GENERIC int num_attributes() const { return NUM_SYS_ATTRIBUTES; };
356  GENERIC int num_body_attributes() const { return NUM_BODY_ATTRIBUTES; };
357 
360  GENERIC double distance_between(const int& i , const int & j ) const {
361  return sqrt(distance_squared_between(i,j));
362  }
364  GENERIC double distance_squared_between(const int& i , const int & j )const {
365  const Body& b1 = _body[i], & b2 = _body[j];
366  return square(b1[0].pos()-b2[0].pos())
367  + square(b1[1].pos()-b2[1].pos())
368  + square(b1[2].pos()-b2[2].pos());
369  }
370 
376  GENERIC void copyTo(const SystemRef& s ) {
377  for(int i = 0; i < _nbod; i++){
378  s[i][0].pos() = _body[i][0].pos();
379  s[i][1].pos() = _body[i][1].pos();
380  s[i][2].pos() = _body[i][2].pos();
381  s[i][0].vel() = _body[i][0].vel();
382  s[i][1].vel() = _body[i][1].vel();
383  s[i][2].vel() = _body[i][2].vel();
384  s[i].mass() = _body[i].mass();
385  for(int j = 0; j < NUM_BODY_ATTRIBUTES; j++){
386  s[i].attribute(j) = _body[i].attribute(j);
387  }
388  }
389  for(int j = 0; j < NUM_SYS_ATTRIBUTES; j++){
390  s.attribute(j) = attribute(j);
391  }
392  s.time() = time();
393  s.state() = state();
394  s.id() = id();
395  }
396 
397  GENERIC double total_energy()const{
398  return calc_total_energy();
399  }
400 
402  GENERIC double calc_total_energy() const {
403  double K = 0.0, U = 0.0;
404 
405  for(int i = 0; i < _nbod; i++)
406  K += 0.5 * _body[i].mass() * _body[i].speed_squared();
407 
408  for(int i = 0; i < _nbod; i++)
409  for(int j = 0; j < i; j++)
410  U += - _body[i].mass() * _body[j].mass() / distance_between(i,j);
411 
412  return K + U;
413  }
414  };
415 
419  class SystemRefConst {
420  SystemRef _ref;
421  public:
422 
423  // Constructor
424  GENERIC SystemRefConst(const SystemRef& ref):_ref(ref){}
425 
426  // Accessor
427  GENERIC const Body& operator[](const int & i ) const { return _ref[i]; }
428  GENERIC const double& time() { return _ref.time(); }
429  GENERIC const int& number() { return _ref.number(); }
430  GENERIC const double& time() const { return _ref.time(); }
431  GENERIC bool is_active() const { return _ref.state() == Sys::SYSTEM_ACTIVE; }
432  GENERIC bool is_disabled() const { return _ref.is_disabled(); }
433  GENERIC bool is_enabled() const { return !_ref.is_disabled(); }
434  GENERIC const int& state() const { return _ref.state(); }
435  GENERIC const int& id() const { return _ref.id(); }
436  GENERIC const double& attribute(const int& i) const { return _sys[0].attribute(i); }
437  GENERIC int num_attributes() const { return NUM_SYS_ATTRIBUTES; };
438  GENERIC int num_body_attributes() const { return NUM_BODY_ATTRIBUTES; };
439  GENERIC const int& nbod()const{ return _ref.nbod(); }
440  GENERIC double distance_squared_between(const int& i , const int & j ) { return _ref.distance_squared_between(i,j); }
441  GENERIC double distance_between(const int& i , const int & j ) { return _ref.distance_between(i,j); }
442  GENERIC void copyTo( const SystemRef& r ) { _ref.copyTo( r ) ; }
443  GENERIC double calc_total_energy() const { return _ref.calc_total_energy(); }
444  };
445 
446 
448  GENERIC static size_t body_element_count(const int& nbod,const int& nsys){
449  return (nsys + CHUNK_SIZE-1) / CHUNK_SIZE * nbod ;
450  }
451 
453  GENERIC static size_t sys_element_count(const int& nsys){
454  return (nsys + CHUNK_SIZE-1) / CHUNK_SIZE ;
455  }
456 
458  typedef CoalescedStructArray< Body, double, CHUNK_SIZE> BodyArray;
460  typedef CoalescedStructArray< Sys, double, CHUNK_SIZE> SysArray;
461  typedef typename BodyArray::PItem PBody;
462  typedef typename SysArray::PItem PSys;
463 
464  protected:
466  int _nbod;
468  int _nsys;
470  BodyArray _body;
472  SysArray _sys;
473 
477  GENERIC EnsembleBase():_nbod(0),_nsys(0),_body(0,0),_sys(0,0) {};
478 
480  GENERIC explicit EnsembleBase(const int& nbod, const int& nsys,PBody body_array, PSys sys_array):_nbod(nbod),_nsys(nsys),_body(body_array,body_element_count(nbod,nsys)),_sys(sys_array,sys_element_count(nsys)){}
481 
482  public:
483 
484  typedef struct SystemRef value_type;
485  typedef int key_type;
486  typedef int difference_type;
487  typedef int size_type;
488 
490  GENERIC const int& nbod()const{ return _nbod; }
492  GENERIC const int& nsys()const{ return _nsys; }
494  BodyArray& bodies() { return _body; }
496  SysArray& systems() { return _sys; }
497 
498 
512  GENERIC SystemRef operator[] (const int & i) {
513  const int sysinblock= i % CHUNK_SIZE;
514  const int blockid = i / CHUNK_SIZE;
515  const int idx = blockid * _nbod * CHUNK_SIZE + sysinblock;
516  return SystemRef(_nbod,i,&_body[idx], &_sys[i] ) ;
517  };
518 
519  GENERIC SystemRef getitem(const int& i){
520  return operator[](i);
521  }
522 
523  GENERIC SystemRefConst operator[] (const int & i) const {
524  return SystemRefConst( const_cast<EnsembleBase*>(this)->operator[](i) ) ;
525  };
526 
528 
529  GENERIC double& mass(const int& sys, const int & bod){
530  return operator[] ( sys )[bod].mass();
531  }
532 
533  GENERIC double& p(const int& sys, const int & bod, const int& c){
534  return operator[] ( sys )[bod][c].pos();
535  }
536 
537  GENERIC double& v(const int& sys, const int & bod, const int& c){
538  return operator[] ( sys )[bod][c].vel();
539  }
540 
541  GENERIC double& x(const int& sys, const int& bod ) { return p(sys,bod,0); }
542  GENERIC double& y(const int& sys, const int& bod ) { return p(sys,bod,1); }
543  GENERIC double& z(const int& sys, const int& bod ) { return p(sys,bod,2); }
544 
545  GENERIC double& vx(const int& sys, const int& bod ) { return v(sys,bod,0); }
546  GENERIC double& vy(const int& sys, const int& bod ) { return v(sys,bod,1); }
547  GENERIC double& vz(const int& sys, const int& bod ) { return v(sys,bod,2); }
548 
549  GENERIC double& time( const int & sys ) {
550  return operator[] ( sys ).time();
551  }
552 
553  GENERIC int& state(const int& sys){
554  return operator[] ( sys ).state();
555  }
556 
557  GENERIC int& flags(const int& sys){
558  return operator[] ( sys ).state();
559  }
560 
562 
563  GENERIC const double& mass(const int& sys, const int & bod)const{
564  return operator[] ( sys )[bod].mass();
565  }
566 
567  GENERIC const double& p(const int& sys, const int & bod, const int& c)const{
568  return operator[] ( sys )[bod][c].pos();
569  }
570 
571  GENERIC const double& v(const int& sys, const int & bod, const int& c)const{
572  return operator[] ( sys )[bod][c].vel();
573  }
574 
575  GENERIC const double& x(const int& sys, const int& bod )const { return p(sys,bod,0); }
576  GENERIC const double& y(const int& sys, const int& bod )const { return p(sys,bod,1); }
577  GENERIC const double& z(const int& sys, const int& bod )const { return p(sys,bod,2); }
578 
579  GENERIC const double& vx(const int& sys, const int& bod ) const { return v(sys,bod,0); }
580  GENERIC const double& vy(const int& sys, const int& bod ) const { return v(sys,bod,1); }
581  GENERIC const double& vz(const int& sys, const int& bod ) const { return v(sys,bod,2); }
582 
583  GENERIC const double& time( const int & sys ) const {
584  return operator[] ( sys ).time();
585  }
586 
587  GENERIC const int& state(const int& sys) const{
588  return operator[] ( sys ).state();
589  }
590 
591  GENERIC const int& flags(const int& sys) const{
592  return operator[] ( sys ).state();
593  }
594 
595  GENERIC void set_time( const int& sys, const double& time ) {
596  SystemRef s = operator[] ( sys );
597  s.time() = time;
598  }
599 
600  GENERIC bool is_active(const int& sys) const {
601  return operator[] ( sys ).is_active();
602  }
603  GENERIC bool is_inactive(const int& sys) const {
604  return ! operator[] ( sys ).is_inactive();
605  }
606  GENERIC bool is_disabled(const int& sys) const {
607  return ! operator[] ( sys ).is_disabled();
608  }
609  GENERIC bool is_enabled(const int& sys) const {
610  return ! operator[] ( sys ).is_disabled();
611  }
612  GENERIC void set_active(const int& sys) {
613  operator[] ( sys ).set_active();
614  }
615  GENERIC void set_inactive(const int& sys) {
616  operator[] ( sys ).set_inactive();
617  }
618  GENERIC void set_disabled(const int& sys) {
619  operator[] ( sys ).set_disabled();
620  }
621 
622 
623  GENERIC void set_body( const int& sys, const int& bod, const double& mass_planet
624  , const double& x, const double & y, const double& z
625  , const double& vx, const double& vy, const double& vz) {
626  SystemRef s = operator[] ( sys );
627  s[bod][0].pos() = x ;
628  s[bod][1].pos() = y ;
629  s[bod][2].pos() = z ;
630 
631  s[bod][0].vel() = vx ;
632  s[bod][1].vel() = vy ;
633  s[bod][2].vel() = vz ;
634 
635  s[bod].mass() = mass_planet;
636 
637  }
638 
639  GENERIC void get_body(const int & sys, const int & bod, double & m
640  , double& x, double & y, double& z
641  , double& vx, double& vy, double& vz) const {
642  SystemRefConst s = operator[] ( sys );
643 
644  x = s[bod][0].pos();
645  y = s[bod][1].pos();
646  z = s[bod][2].pos();
647 
648  vx = s[bod][0].vel();
649  vy = s[bod][1].vel();
650  vz = s[bod][2].vel();
651 
652  m = s[bod].mass();
653  }
654 
655 
656  // Utilities
657  //
658 
659 
665  GENERIC void get_barycenter(const int& sys, double& x, double& y, double& z, double& vx, double& vy, double& vz, const int& max_body_id = MAX_NBODIES) const
666  {
667 
668  x = 0.; y = 0.; z = 0.; vx = 0.; vy = 0.; vz = 0.;
669  double mass_sum = 0.;
670  const int stop_at_body = std::min(nbod()-1,max_body_id);
671  for(int bod=0;bod<=stop_at_body;++bod)
672  {
673  double m = mass(sys,bod);
674  x += m* this->x(sys,bod);
675  y += m* this->y(sys,bod);
676  z += m* this->z(sys,bod);
677  vx += m* this->vx(sys,bod);
678  vy += m* this->vy(sys,bod);
679  vz += m* this->vz(sys,bod);
680  mass_sum += m;
681  }
682  x /= mass_sum;
683  y /= mass_sum;
684  z /= mass_sum;
685  vx /= mass_sum;
686  vy /= mass_sum;
687  vz /= mass_sum;
688  };
689 
690 
692  GENERIC double calc_total_energy( int sys ) const {
693  SystemRefConst s = operator[] ( sys );
694  return s.calc_total_energy();
695  }
696 
699  GENERIC void calc_total_energy(double* E) const {
700  for (int sys = 0; sys != nsys(); sys++)
701  E[sys] = calc_total_energy(sys);
702  }
703 
705  struct range_t {
706  double average, min, max,median;
707  range_t(const double& a,const double& m, const double& M, const double& d)
708  :average(a),min(m),max(M),median(d){}
709 
712  template<class RandomAccessIterator>
713  static range_t calculate(RandomAccessIterator begin, RandomAccessIterator end){
714  size_t n = end - begin;
715  double min = *begin;
716  double max = *begin;
717  double sum = *begin;
718  for(RandomAccessIterator i = begin+1; i != end; i++) {
719  double v = *i;
720  if( v < min ) min = v;
721  if( v > max ) max = v;
722  sum += v;
723  }
724  std::nth_element(begin,end, begin+n/2);
725  return range_t(sum/n,min,max,*(begin+n/2));
726  }
727  };
728 
732  GENERIC range_t time_ranges() const {
733  std::vector<double> times(nsys());
734  for(int i= 0; i < nsys(); i++)
735  times[i] = operator[](i).time();
736 
737  return range_t::calculate(times.begin(),times.end());
738  }
739 
740 };
741 
743 typedef EnsembleBase< ENSEMBLE_CHUNK_SIZE > ensemble;
744 
745 
746 }