SISCone  2.0.6
protocones.cpp
1 // File: protocones.cpp //
3 // Description: source file for stable cones determination (Cstable_cones) //
4 // This file is part of the SISCone project. //
5 // WARNING: this is not the main SISCone trunk but //
6 // an adaptation to spherical coordinates //
7 // For more details, see http://projects.hepforge.org/siscone //
8 // //
9 // Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez //
10 // //
11 // This program is free software; you can redistribute it and/or modify //
12 // it under the terms of the GNU General Public License as published by //
13 // the Free Software Foundation; either version 2 of the License, or //
14 // (at your option) any later version. //
15 // //
16 // This program is distributed in the hope that it will be useful, //
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
19 // GNU General Public License for more details. //
20 // //
21 // You should have received a copy of the GNU General Public License //
22 // along with this program; if not, write to the Free Software //
23 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
24 // //
25 // $Revision:: 311 $//
26 // $Date:: 2011-10-05 23:27:09 +0200 (Wed, 05 Oct 2011) $//
28 
29 /*******************************************************
30  * Introductory note: *
31  * Since this file has many member functions, we have *
32  * structured them in categories: *
33  * INITIALISATION FUNCTIONS *
34  * - ctor() *
35  * - ctor(particle_list) *
36  * - dtor() *
37  * - init(particle_list) *
38  * ALGORITHM MAIN ENTRY *
39  * - get_stable_cone(radius) *
40  * ALGORITHM MAIN STEPS *
41  * - init_cone() *
42  * - test_cone() *
43  * - update_cone() *
44  * - proceed_with_stability() *
45  * ALGORITHM MAIN STEPS FOR COCIRCULAR SITUATIONS *
46  * - cocircular_pt_less(v1, v2) *
47  * - prepare_cocircular_list() *
48  * - test_cone_cocircular() *
49  * - test_stability(candidate, border_list) *
50  * - updat_cone_cocircular() *
51  * RECOMPUTATION OF CONE CONTENTS *
52  * - compute_cone_contents() *
53  * - recompute_cone_contents() *
54  * - recompute_cone_contents_if_needed() *
55  * VARIOUS TOOLS *
56  * - circle_intersect() *
57  * - is_inside() *
58  * - abs_dangle() *
59  *******************************************************/
60 
61 #include <siscone/defines.h>
62 #include <siscone/siscone_error.h>
63 #include <siscone/circulator.h>
64 #include "protocones.h"
65 #include <math.h>
66 #include <iostream>
67 #include <algorithm>
68 
69 namespace siscone_spherical{
70 
71 using namespace std;
72 
73 /**********************************************************************
74  * CSphstable_cones implementation *
75  * Computes the list of stable comes from a particle list. *
76  * This class does the first fundamental task of te cone algorithm: *
77  * it is used to compute the list of stable cones given a list *
78  * of particles. *
79  **********************************************************************/
80 
82 // INITIALISATION FUNCTIONS //
83 // - ctor() //
84 // - ctor(particle_list) //
85 // - dtor() //
86 // - init(particle_list) //
88 
89 // default ctor
90 //--------------
92  nb_tot = 0;
93  hc = NULL;
94 }
95 
96 // ctor with initialisation
97 //--------------------------
98 CSphstable_cones::CSphstable_cones(vector<CSphmomentum> &_particle_list)
99  : CSphvicinity(_particle_list){
100 
101  nb_tot = 0;
102  hc = NULL;
103 }
104 
105 // default dtor
106 //--------------
108  if (hc!=NULL) delete hc;
109 }
110 
111 /*
112  * initialisation
113  * - _particle_list list of particles
114  * - _n number of particles
115  *********************************************************************/
116 void CSphstable_cones::init(vector<CSphmomentum> &_particle_list){
117  // check already allocated mem
118  if (hc!=NULL){
119  delete hc;
120  }
121  if (protocones.size()!=0)
122  protocones.clear();
123 
124  multiple_centre_done.clear();
125 
126  // initialisation
127  set_particle_list(_particle_list);
128 }
129 
130 
132 // ALGORITHM MAIN ENTRY //
133 // - get_stable_cone(radius) //
135 
136 /*
137  * compute stable cones.
138  * This function really does the job i.e. computes
139  * the list of stable cones (in a seedless way)
140  * - _radius: radius of the cones
141  * The number of stable cones found is returned
142  *********************************************************************/
144  int p_idx;
145 
146  // check if everything is correctly initialised
147  if (n_part==0){
148  return 0;
149  }
150 
151  R = _radius;
152  R2 = R*R;
153  tan2R = tan(R);
154  tan2R *= tan2R;
155 
156  // allow hash for cones candidates
157  hc = new sph_hash_cones(n_part, R);
158 
159  // browse all particles
160  for (p_idx=0;p_idx<n_part;p_idx++){
161  // step 0: compute the child list CL.
162  // Note that this automatically sets the parent P
163  build(&plist[p_idx], 2.0*R);
164 
165  // special case:
166  // if the vicinity is empty, the parent particle is a
167  // stable cone by itself. Add it to protocones list.
168  if (vicinity_size==0){
169  protocones.push_back(*parent);
170  continue;
171  }
172 
173 #ifdef DEBUG_STABLE_CONES
174  cout << endl << endl;
175  cout << "plot 'particles.dat' u 2:1 pt 1 ps 3" << endl;
176  cout << "set label 1 'x' at " << parent->_phi << ", " << parent->_theta << endl;
177 #endif
178 
179  // step 1: initialise with the first cone candidate
180  init_cone();
181 
182  do{
183  // step 2: test cone stability for that pair (P,C)
184  test_cone();
185 
186  // step 3: go to the next cone child candidate C
187  } while (!update_cone());
188  }
189 
190  return proceed_with_stability();
191 }
192 
193 
195 // ALGORITHM MAIN STEPS //
196 // - init_cone() //
197 // - test_cone() //
198 // - update_cone() //
199 // - proceed_with_stability() //
201 
202 /*
203  * initialise the cone.
204  * We take the first particle in the angular ordering to compute
205  * this one
206  * return 0 on success, 1 on error
207  *********************************************************************/
208 int CSphstable_cones::init_cone(){
209  // The previous version of the algorithm was starting the
210  // loop around vicinity elements with the "most isolated" child.
211  // given the nodist method to calculate the cone contents, we no
212  // longer need to worry about which cone comes first...
213  first_cone=0;
214 
215  // now make sure we have lists of the cocircular particles
216  prepare_cocircular_lists();
217 
218  //TODO? deal with a configuration with only degeneracies ?
219  // The only possibility seems a regular hexagon with a parent point
220  // in the centre. And this situation is by itself unclear.
221  // Hence, we do nothing here !
222 
223  // init set child C
224  centre = vicinity[first_cone];
225  child = centre->v;
226  centre_idx = first_cone;
227 
228  // build the initial cone (nodist: avoids calculating distances --
229  // just deduces contents by circulating around all in/out operations)
230  // this function also sets the list of included particles
231  compute_cone_contents();
232 
233  return 0;
234 }
235 
236 
237 /*
238  * test cones.
239  * We check if the cone(s) built with the present parent and child
240  * are stable
241  * return 0 on success 1 on error
242  *********************************************************************/
243 int CSphstable_cones::test_cone(){
244  siscone::Creference weighted_cone_ref;
245 
246  // depending on the side we are taking the child particle,
247  // we test different configuration.
248  // Each time, two configurations are tested in such a way that
249  // all 4 possible cases (parent or child in or out the cone)
250  // are tested when taking the pair of particle parent+child
251  // and child+parent.
252 
253  // here are the tests entering the first series:
254  // 1. check if the cone is already inserted
255  // 2. check cone stability for the parent and child particles
256 
257  //UPDATED(see below): if (centre->side){
258  //UPDATED(see below): // test when both particles are not in the cone
259  //UPDATED(see below): // or when both are in.
260  //UPDATED(see below): // Note: for the totally exclusive case, test emptyness before
261  //UPDATED(see below): cone_candidate = cone;
262  //UPDATED(see below): if (cone.ref.not_empty()){
263  //UPDATED(see below): hc->insert(&cone_candidate, parent, child, false, false);
264  //UPDATED(see below): }
265  //UPDATED(see below):
266  //UPDATED(see below): cone_candidate = cone;
267  //UPDATED(see below): cone_candidate+= *parent + *child;
268  //UPDATED(see below): hc->insert(&cone_candidate, parent, child, true, true);
269  //UPDATED(see below): } else {
270  //UPDATED(see below): // test when 1! of the particles is in the cone
271  //UPDATED(see below): cone_candidate = cone + *parent;
272  //UPDATED(see below): hc->insert(&cone_candidate, parent, child, true, false);
273  //UPDATED(see below):
274  //UPDATED(see below): cone_candidate = cone + *child;
275  //UPDATED(see below): hc->insert(&cone_candidate, parent, child, false, true);
276  //UPDATED(see below): }
277  //UPDATED(see below):
278  //UPDATED(see below): nb_tot+=2;
279 
280  // instead of testing 2 inclusion/exclusion states for every pair, we test the 4 of them
281  // when the parent has an energy bigger than the child
282  if (parent->E >= child->E){
283  // test when both particles are not in the cone
284  // Note: for the totally exclusive case, test emptiness before
285  cone_candidate = cone;
286  if (cone.ref.not_empty()){
287  hc->insert(&cone_candidate, parent, child, false, false);
288  }
289 
290  // test when 1! of the particles is in the cone
291  cone_candidate += *parent;
292  hc->insert(&cone_candidate, parent, child, true, false);
293 
294  cone_candidate = cone;
295  cone_candidate += *child;
296  hc->insert(&cone_candidate, parent, child, false, true);
297 
298  // test when both are in.
299  cone_candidate += *parent;
300  hc->insert(&cone_candidate, parent, child, true, true);
301 
302  nb_tot += 4;
303  }
304 
305 
306  return 0;
307 }
308 
309 
310 /*
311  * update the cone
312  * go to the next child for that parent and update 'cone' appropriately
313  * return 0 if update candidate found, 1 otherwise
314  ***********************************************************************/
315 int CSphstable_cones::update_cone(){
316 #ifdef DEBUG_STABLE_CONES
317  cout << "call 'circles_plot.gp' '" << centre->centre.px << "' '"
318  << centre->centre.py << "' '" << centre->centre.pz << "'" << endl
319  << "pause -1 '(" << centre->angle << " " << (centre->side ? '+' : '-') << ")";
320 #endif
321 
322  // get the next child and centre
323  centre_idx++;
324  if (centre_idx==vicinity_size)
325  centre_idx=0;
326  if (centre_idx==first_cone)
327  return 1;
328 
329  // update the cone w.r.t. the old child
330  // only required if the old child is entering inside in which
331  // case we need to add it. We also know that the child is
332  // inside iff its side is -.
333  if (!centre->side){
334 #ifdef DEBUG_STABLE_CONES
335  cout << " old_enter";
336 #endif
337  // update cone
338  cone += (*child);
339 
340  // update info on particles inside
341  centre->is_inside->cone = true;
342 
343  // update stability check quantities
344  dpt += fabs(child->px)+fabs(child->py)+fabs(child->pz);
345  }
346 
347  // update centre and child to correspond to the new position
348  centre = vicinity[centre_idx];
349  child = centre->v;
350 
351  // check cocircularity
352  // note that if cocirculaity is detected (i.e. if we receive 1
353  // in the next test), we need to recall 'update_cone' directly
354  // since tests and remaining part of te update has been performed
355  //if (cocircular_check())
356  if (cocircular_check()){
357 #ifdef DEBUG_STABLE_CONES
358  cout << " Co-circular case detected" << endl;
359 #endif
360  return update_cone();
361  }
362 
363  // update the cone w.r.t. the new child
364  // only required if the new child was already inside in which
365  // case we need to remove it. We also know that the child is
366  // inside iff its side is +.
367  if ((centre->side) && (cone.ref.not_empty())){
368 #ifdef DEBUG_STABLE_CONES
369  cout << " new exit";
370 #endif
371 
372  // update cone
373  cone -= (*child);
374 
375  // update info on particles inside
376  centre->is_inside->cone = false;
377 
378  // update stability check quantities
379  dpt += fabs(child->px)+fabs(child->py)+fabs(child->pz); //child->perp2();
380  }
381 
382  // check that the addition and subtraction of vectors does
383  // not lead to too much rounding error
384  // for that, we compute the sum of pt modifications and of |pt|
385  // since last recomputation and once the ratio overpasses a threshold
386  // we recompute vicinity.
387  if ((dpt>PT_TSHOLD*(fabs(cone.px)+fabs(cone.py)+fabs(cone.pz))) && (cone.ref.not_empty())){
388  recompute_cone_contents();
389  }
390  if (cone.ref.is_empty()){
391  cone = CSphmomentum();
392  dpt=0.0;
393  }
394 
395 #ifdef DEBUG_STABLE_CONES
396  cout << "'" << endl;
397 #endif
398 
399  return 0;
400 }
401 
402 
403 /*
404  * compute stability of all enumerated candidates.
405  * For all candidate cones which are stable w.r.t. their border particles,
406  * pass the last test: stability with quadtree intersection
407  ************************************************************************/
408 int CSphstable_cones::proceed_with_stability(){
409  int i;
410  sph_hash_element *elm;
411 
412  for (i=0;i<=hc->mask;i++){
413  // test ith cell of the hash array
414  elm = hc->hash_array[i];
415 
416  // browse elements therein
417  while (elm!=NULL){
418  // test stability
419  if (elm->is_stable){
420  // stability is not ensured by all pairs of "edges" already browsed
421 #ifdef USE_QUADTREE_FOR_STABILITY_TEST
422  // => testing stability with quadtree intersection
423  if (quadtree->circle_intersect(elm->eta, elm->phi, R2)==elm->ref)
424 #else
425  // => testing stability with the particle-list intersection
426  if (circle_intersect(elm->centre)==elm->centre.ref)
427 #endif
428  protocones.push_back(CSphmomentum(elm->centre,1.0));
429  }
430 
431  // jump to the next one
432  elm = elm->next;
433  }
434  }
435 
436  // free hash
437  // we do that at this level because hash eats rather a lot of memory
438  // we want to free it before running the split/merge algorithm
439 #ifdef DEBUG_STABLE_CONES
440  nb_hash_cones = hc->n_cones;
441  nb_hash_occupied = hc->n_occupied_cells;
442 #endif
443 
444  delete hc;
445  hc=NULL;
446 
447  return protocones.size();
448 }
449 
450 
452 // ALGORITHM MAIN STEPS FOR COCIRCULAR SITUATIONS //
453 // - cocircular_pt_less(v1, v2) //
454 // - prepare_cocircular_list() //
455 // - test_cone_cocircular() //
456 // - test_stability(candidate, border_vect) //
457 // - updat_cone_cocircular() //
459 
461 //NEVER USED
462 //bool cocircular_pt_less(CSphmomentum *v1, CSphmomentum *v2){
463 // return v1->perp2() < v2->perp2();
464 //}
465 
466 /*
467  * run through the vicinity of the current parent and for each child
468  * establish which other members are cocircular... Note that the list
469  * associated with each child contains references to vicinity
470  * elements: thus two vicinity elements each associated with one given
471  * particle may appear in a list -- this needs to be watched out for
472  * later on...
473  **********************************************************************/
474 void CSphstable_cones::prepare_cocircular_lists() {
476  vicinity.begin(),
477  vicinity.end());
478 
480 
481  do {
482  CSphvicinity_elm* here_pntr = *here();
483  search.set_position(here);
484 
485  // search forwards for things that should have "here" included in
486  // their cocircularity list
487  while (true) {
488  ++search;
489  if ( siscone::abs_dphi((*search())->angle, here_pntr->angle) <
490  here_pntr->cocircular_range
491  && search() != here()) {
492  (*search())->cocircular.push_back(here_pntr);
493  } else {
494  break;
495  }
496  }
497 
498  // search backwards
499  search.set_position(here);
500  while (true) {
501  --search;
502  if ( siscone::abs_dphi((*search())->angle, here_pntr->angle) <
503  here_pntr->cocircular_range
504  && search() != here()) {
505  (*search())->cocircular.push_back(here_pntr);
506  } else {
507  break;
508  }
509  }
510 
511  ++here;
512  } while (here() != vicinity.begin());
513 }
514 
515 /*
516  * Testing cocircular configurations in p^3 time,
517  * rather than 2^p time; we will test all contiguous subsets of points
518  * on the border --- note that this is till probably overkill, since
519  * in principle we only have to test situations where up to a
520  * half-circle is filled (but going to a full circle is simpler)
521  ******************************************************************/
522 void CSphstable_cones::test_cone_cocircular(CSphmomentum & borderless_cone,
523  list<CSphmomentum *> & border_list) {
524  // in spherical coordinates, we don't have a universal x-y axis system
525  // to measure the angles. So we first determine one minimising
526  // the uncertainties
527  CSph3vector angl_dir1, angl_dir2;
528  centre->centre.get_angular_directions(angl_dir1, angl_dir2);
529  angl_dir1/=angl_dir1._norm;
530  angl_dir2/=angl_dir2._norm;
531 
532  // now we have te reference axis, create the CSphborder_store structure
533  vector<CSphborder_store> border_vect;
534  border_vect.reserve(border_list.size());
535  for (list<CSphmomentum *>::iterator it = border_list.begin();
536  it != border_list.end(); it++) {
537  border_vect.push_back(CSphborder_store(*it, centre->centre, angl_dir1, angl_dir2));
538  }
539 
540  // get them into order of angle
541  sort(border_vect.begin(), border_vect.end());
542 
543  // set up some circulators, since these will help us go around the
544  // circle easily
546  start(border_vect.begin(), border_vect.begin(),border_vect.end());
548 
549  // test the borderless cone
550  CSphmomentum candidate = borderless_cone;
551  //candidate.build_etaphi();
552  if (candidate.ref.not_empty())
553  test_stability(candidate, border_vect);
554 
555  do {
556  // reset status wrt inclusion in the cone
557  mid = start;
558  do {
559  mid()->is_in = false;
560  } while (++mid != start);
561 
562  // now run over all inclusion possibilities with this starting point
563  candidate = borderless_cone;
564  while (++mid != start) {
565  // will begin with start+1 and go up to start-1
566  mid()->is_in = true;
567  candidate += *(mid()->mom);
568  test_stability(candidate, border_vect);
569  }
570 
571  } while (++start != end);
572 
573  // mid corresponds to momentum that we need to include to get the
574  // full cone
575  mid()->is_in = true;
576  candidate += *(mid()->mom);
577  test_stability(candidate, border_vect);
578 }
579 
580 
587 void CSphstable_cones::test_stability(CSphmomentum & candidate, const vector<CSphborder_store> & border_vect) {
588 
589  // this almost certainly has not been done...
590  //candidate.build_etaphi();
591 
592  bool stable = true;
593  for (unsigned i = 0; i < border_vect.size(); i++) {
594  if (is_closer(&candidate, border_vect[i].mom,tan2R) ^ (border_vect[i].is_in)) {
595  stable = false;
596  break; // it's unstable so there's no point continuing
597  }
598  }
599 
600  if (stable) hc->insert(&candidate);
601 }
602 
603 /*
604  * check if we are in a situation of cocircularity.
605  * if it is the case, update and test in the corresponding way
606  * return 'false' if no cocircularity detected, 'true' otherwise
607  * Note that if cocircularity is detected, we need to
608  * recall 'update' from 'update' !!!
609  ***************************************************************/
610 bool CSphstable_cones::cocircular_check(){
611  // check if many configurations have the same centre.
612  // if this is the case, branch on the algorithm for this
613  // special case.
614  // Note that those situation, being considered separately in
615  // test_cone_multiple, must only be considered here if all
616  // angles are on the same side (this avoid multiple counting)
617 
618  if (centre->cocircular.empty()) return false;
619 
620  // first get cone into status required at end...
621  if ((centre->side) && (cone.ref.not_empty())){
622  // update cone
623  cone -= (*child);
624 
625  // update info on particles inside
626  centre->is_inside->cone = false;
627 
628  // update stability check quantities
629  dpt += fabs(child->px)+fabs(child->py)+fabs(child->pz); //child->perp2();
630  }
631 
632 
633  // now establish the list of unique children in the list
634  // first make sure parent and child are in!
635 
636  list<siscone::Cvicinity_inclusion *> removed_from_cone;
637  list<siscone::Cvicinity_inclusion *> put_in_border;
638  list<CSphmomentum *> border_list;
639 
640  CSphmomentum cone_removal;
641  CSphmomentum border = *parent;
642  border_list.push_back(parent);
643 
644  // make sure child appears in the border region
645  centre->cocircular.push_back(centre);
646 
647  // now establish the full contents of the cone minus the cocircular
648  // region and of the cocircular region itself
649  for(list<CSphvicinity_elm *>::iterator it = centre->cocircular.begin();
650  it != centre->cocircular.end(); it++) {
651 
652  if ((*it)->is_inside->cone) {
653  cone_removal += *((*it)->v);
654  (*it)->is_inside->cone = false;
655  removed_from_cone.push_back((*it)->is_inside);
656  }
657 
658  // if a point appears twice (i.e. with + and - sign) in the list of
659  // points on the border, we take care not to include it twice.
660  // Note that this situation may appear when a point is at a distance
661  // close to 2R from the parent
662  if (!(*it)->is_inside->cocirc) {
663  border += *((*it)->v);
664  (*it)->is_inside->cocirc = true;
665  put_in_border.push_back((*it)->is_inside);
666  border_list.push_back((*it)->v);
667  //cout << " adding particle " << (*it)->v->_theta << ", " << (*it)->v->_phi << " to the border list" << endl;
668  }
669  }
670 
671 
672  // figure out whether this pairing has been observed before
673  CSphmomentum borderless_cone = cone;
674  borderless_cone -= cone_removal;
675  bool consider = true;
676  for (unsigned int i=0;i<multiple_centre_done.size();i++){
677  if ((multiple_centre_done[i].first ==borderless_cone.ref) &&
678  (multiple_centre_done[i].second==border.ref))
679  consider = false;
680  }
681 
682  // now prepare the hard work
683  if (consider) {
684  // record the fact that we've now seen this combination
685  multiple_centre_done.push_back(pair<siscone::Creference,siscone::Creference>(borderless_cone.ref,
686  border.ref));
687 
688  // first figure out whether our cone momentum is good
689  double local_dpt = fabs(cone_removal.px) + fabs(cone_removal.py);
690  double total_dpt = dpt + local_dpt;
691 
692  recompute_cone_contents_if_needed(borderless_cone, total_dpt);
693  if (total_dpt == 0) {
694  // a recomputation has taken place -- so take advantage of this
695  // and update the member cone momentum
696  cone = borderless_cone + cone_removal;
697  dpt = local_dpt;
698  }
699 
700  test_cone_cocircular(borderless_cone, border_list);
701  }
702 
703 
704  // relabel things that were in the cone but got removed
705  for(list<siscone::Cvicinity_inclusion *>::iterator is_in = removed_from_cone.begin();
706  is_in != removed_from_cone.end(); is_in++) {
707  (*is_in)->cone = true;
708  }
709 
710  // relabel things that got put into the border
711  for(list<siscone::Cvicinity_inclusion *>::iterator is_in = put_in_border.begin();
712  is_in != put_in_border.end(); is_in++) {
713  (*is_in)->cocirc = false;
714  }
715 
716  // we're done with everything -- return true to signal to user that we've
717  // been through the co-circularity rigmarole
718  return true;
719 }
720 
721 
723 // RECOMPUTATION OF CONE CONTENTS //
724 // - compute_cone_contents() //
725 // - recompute_cone_contents() //
726 // - recompute_cone_contents_if_needed() //
728 
737 void CSphstable_cones::compute_cone_contents() {
739  start(vicinity.begin()+first_cone, vicinity.begin(), vicinity.end());
740 
742 
743  // note that in the following algorithm, the cone contents never includes
744  // the child. Indeed, if it has positive sign, then it will be set as
745  // outside at the last step in the loop. If it has negative sign, then the
746  // loop will at some point go to the corresponding situation with positive
747  // sign and set the inclusion status to 0.
748 
749  do {
750  // as we leave this position a particle enters if its side is
751  // negative (i.e. the centre is the one at -ve angle wrt to the
752  // parent-child line
753  if (!(*here())->side) ((*here())->is_inside->cone) = 1;
754 
755  // move on to the next position
756  ++here;
757 
758  // as we arrive at this position a particle leaves if its side is positive
759  if ((*here())->side) ((*here())->is_inside->cone) = 0;
760  } while (here != start);
761 
762  // once we've reached the start the 'is_inside' information should be
763  // 100% complete, so we can use it to calculate the cone contents
764  // and then exit
765  recompute_cone_contents();
766  return;
767 
768 }
769 
770 
771 /*
772  * compute the cone momentum from particle list.
773  * in this version, we use the 'pincluded' information
774  * from the CSphvicinity class
775  */
776 void CSphstable_cones::recompute_cone_contents(){
777  unsigned int i;
778 
779  // set momentum to 0
780  cone = CSphmomentum();
781 
782  // Important note: we can browse only the particles
783  // in vicinity since all particles in the cone are
784  // withing a distance 2R w.r.t. parent hence in vicinity.
785  // Among those, we only add the particles for which 'is_inside' is true !
786  // This methos rather than a direct comparison avoids rounding errors
787  for (i=0;i<vicinity_size;i++){
788  // to avoid double-counting, only use particles with + angle
789  if ((vicinity[i]->side) && (vicinity[i]->is_inside->cone))
790  cone += *vicinity[i]->v;
791  }
792 
793  // set check variables back to 0
794  dpt = 0.0;
795 }
796 
797 
798 /*
799  * if we have gone beyond the acceptable threshold of change, compute
800  * the cone momentum from particle list. in this version, we use the
801  * 'pincluded' information from the CSphvicinity class, but we don't
802  * change the member cone, only the locally supplied one
803  */
804 void CSphstable_cones::recompute_cone_contents_if_needed(CSphmomentum & this_cone,
805  double & this_dpt){
806 
807  if (this_dpt > PT_TSHOLD*(fabs(this_cone.px)+fabs(this_cone.py))) {
808  if (cone.ref.is_empty()) {
809  this_cone = CSphmomentum();
810  } else {
811  // set momentum to 0
812  this_cone = CSphmomentum();
813 
814  // Important note: we can browse only the particles
815  // in vicinity since all particles in the this_cone are
816  // withing a distance 2R w.r.t. parent hence in vicinity.
817  // Among those, we only add the particles for which 'is_inside' is true !
818  // This methos rather than a direct comparison avoids rounding errors
819  for (unsigned int i=0;i<vicinity_size;i++){
820  // to avoid double-counting, only use particles with + angle
821  if ((vicinity[i]->side) && (vicinity[i]->is_inside->cone))
822  this_cone += *vicinity[i]->v;
823  }
824 
825  }
826  // set check variables back to 0
827  this_dpt = 0.0;
828  }
829 
830 }
831 
832 
834 // VARIOUS TOOLS //
835 // - circle_intersect() //
836 // - is_inside() //
837 // - abs_dangle() //
839 
840 
841 /*
842  * circle intersection.
843  * computes the intersection with a circle of given centre and radius.
844  * The output takes the form of a checkxor of the intersection's particles
845  * - cx circle centre x coordinate
846  * - cy circle centre y coordinate
847  * return the checkxor for the intersection
848  ******************************************************************/
849 siscone::Creference CSphstable_cones::circle_intersect(CSph3vector &cone_centre){
850  siscone::Creference intersection;
851  int i;
852 
853  for (i=0;i<n_part;i++){
854  // really check if the distance is less than R
855  if (is_closer(&cone_centre, &(plist[i]), tan2R))
856  intersection+=plist[i].ref;
857  }
858 
859  return intersection;
860 }
861 
862 }
CSph3vector centre
centre of the cone
Definition: hash.h:48
double _norm
particle spatial norm (available ONLY after a call to build_norm)
Definition: momentum.h:135
sph_hash_element * next
pointer to the next element
Definition: hash.h:51
void set_position(const circulator< T > &other)
set just the position without resetting the begin and end elements
Definition: circulator.h:51
class for storing a border momentum (in context of co-circularity checks).
Definition: protocones.h:55
list of element in the vicinity of a parent.
Definition: vicinity.h:83
siscone::Creference ref
reference number for the vector
Definition: momentum.h:142
a circulator that is foreseen to take as template member either a pointer or an iterator; ...
Definition: circulator.h:36
double cocircular_range
amount by which the angle can be varied while maintaining this point within co-circularity margin ...
Definition: vicinity.h:64
double angle
angle with parent
Definition: vicinity.h:62
double py
y-momentum
Definition: momentum.h:132
void get_angular_directions(CSph3vector &angular_dir1, CSph3vector &angular_dir2)
for this direction, compute the two reference directions used to measure angles
Definition: momentum.cpp:161
unsigned int ref[3]
actual data for the reference
Definition: reference.h:72
bool not_empty()
test non-emptyness
Definition: reference.cpp:81
#define PT_TSHOLD
program name
Definition: defines.h:61
bool is_stable
true if stable w.r.t. "border particles"
Definition: hash.h:49
void init(std::vector< CSphmomentum > &_particle_list)
initialisation
Definition: protocones.cpp:116
base class for dynamic coordinates management
Definition: momentum.h:158
list of cones candidates.
Definition: hash.h:61
base class for managing the spatial part of Cmomentum (defined after)
Definition: momentum.h:54
element in the vicinity of a parent.
Definition: vicinity.h:52
double px
x-momentum
Definition: momentum.h:131
information on store cones candidates.
Definition: hash.h:46
references used for checksums.
Definition: reference.h:43
int get_stable_cones(double _radius)
compute stable cones.
Definition: protocones.cpp:143
The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated for SISCone by  Doxygen 1.8.12