62 #include <siscone/siscone_error.h> 63 #include <siscone/circulator.h> 64 #include "protocones.h" 108 if (hc!=NULL)
delete hc;
121 if (protocones.size()!=0)
124 multiple_centre_done.clear();
127 set_particle_list(_particle_list);
160 for (p_idx=0;p_idx<n_part;p_idx++){
163 build(&plist[p_idx], 2.0*R);
168 if (vicinity_size==0){
169 protocones.push_back(*parent);
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;
187 }
while (!update_cone());
190 return proceed_with_stability();
208 int CSphstable_cones::init_cone(){
216 prepare_cocircular_lists();
224 centre = vicinity[first_cone];
226 centre_idx = first_cone;
231 compute_cone_contents();
243 int CSphstable_cones::test_cone(){
282 if (parent->E >= child->E){
285 cone_candidate = cone;
286 if (cone.ref.not_empty()){
287 hc->insert(&cone_candidate, parent, child,
false,
false);
291 cone_candidate += *parent;
292 hc->insert(&cone_candidate, parent, child,
true,
false);
294 cone_candidate = cone;
295 cone_candidate += *child;
296 hc->insert(&cone_candidate, parent, child,
false,
true);
299 cone_candidate += *parent;
300 hc->insert(&cone_candidate, parent, child,
true,
true);
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 ?
'+' :
'-') <<
")";
324 if (centre_idx==vicinity_size)
326 if (centre_idx==first_cone)
334 #ifdef DEBUG_STABLE_CONES 335 cout <<
" old_enter";
341 centre->is_inside->cone =
true;
344 dpt += fabs(child->px)+fabs(child->py)+fabs(child->pz);
348 centre = vicinity[centre_idx];
356 if (cocircular_check()){
357 #ifdef DEBUG_STABLE_CONES 358 cout <<
" Co-circular case detected" << endl;
360 return update_cone();
367 if ((centre->side) && (cone.ref.not_empty())){
368 #ifdef DEBUG_STABLE_CONES 376 centre->is_inside->cone =
false;
379 dpt += fabs(child->px)+fabs(child->py)+fabs(child->pz);
387 if ((dpt>
PT_TSHOLD*(fabs(cone.px)+fabs(cone.py)+fabs(cone.pz))) && (cone.ref.not_empty())){
388 recompute_cone_contents();
390 if (cone.ref.is_empty()){
395 #ifdef DEBUG_STABLE_CONES 408 int CSphstable_cones::proceed_with_stability(){
412 for (i=0;i<=hc->mask;i++){
414 elm = hc->hash_array[i];
421 #ifdef USE_QUADTREE_FOR_STABILITY_TEST 423 if (quadtree->circle_intersect(elm->eta, elm->phi, R2)==elm->ref)
439 #ifdef DEBUG_STABLE_CONES 440 nb_hash_cones = hc->n_cones;
441 nb_hash_occupied = hc->n_occupied_cells;
447 return protocones.size();
474 void CSphstable_cones::prepare_cocircular_lists() {
489 if ( siscone::abs_dphi((*search())->angle, here_pntr->
angle) <
491 && search() != here()) {
492 (*search())->cocircular.push_back(here_pntr);
502 if ( siscone::abs_dphi((*search())->angle, here_pntr->
angle) <
504 && search() != here()) {
505 (*search())->cocircular.push_back(here_pntr);
512 }
while (here() != vicinity.begin());
522 void CSphstable_cones::test_cone_cocircular(
CSphmomentum & borderless_cone,
523 list<CSphmomentum *> & border_list) {
529 angl_dir1/=angl_dir1.
_norm;
530 angl_dir2/=angl_dir2.
_norm;
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));
541 sort(border_vect.begin(), border_vect.end());
546 start(border_vect.begin(), border_vect.begin(),border_vect.end());
553 test_stability(candidate, border_vect);
559 mid()->is_in =
false;
560 }
while (++mid != start);
563 candidate = borderless_cone;
564 while (++mid != start) {
567 candidate += *(mid()->mom);
568 test_stability(candidate, border_vect);
571 }
while (++start != end);
576 candidate += *(mid()->mom);
577 test_stability(candidate, border_vect);
587 void CSphstable_cones::test_stability(
CSphmomentum & candidate,
const vector<CSphborder_store> & border_vect) {
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)) {
600 if (stable) hc->insert(&candidate);
610 bool CSphstable_cones::cocircular_check(){
618 if (centre->cocircular.empty())
return false;
621 if ((centre->side) && (cone.ref.not_empty())){
626 centre->is_inside->cone =
false;
629 dpt += fabs(child->px)+fabs(child->py)+fabs(child->pz);
636 list<siscone::Cvicinity_inclusion *> removed_from_cone;
637 list<siscone::Cvicinity_inclusion *> put_in_border;
638 list<CSphmomentum *> border_list;
642 border_list.push_back(parent);
645 centre->cocircular.push_back(centre);
649 for(list<CSphvicinity_elm *>::iterator it = centre->cocircular.begin();
650 it != centre->cocircular.end(); it++) {
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);
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);
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))
685 multiple_centre_done.push_back(pair<siscone::Creference,siscone::Creference>(borderless_cone.
ref,
689 double local_dpt = fabs(cone_removal.
px) + fabs(cone_removal.
py);
690 double total_dpt = dpt + local_dpt;
692 recompute_cone_contents_if_needed(borderless_cone, total_dpt);
693 if (total_dpt == 0) {
696 cone = borderless_cone + cone_removal;
700 test_cone_cocircular(borderless_cone, border_list);
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;
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;
737 void CSphstable_cones::compute_cone_contents() {
739 start(vicinity.begin()+first_cone, vicinity.begin(), vicinity.end());
753 if (!(*here())->side) ((*here())->is_inside->cone) = 1;
759 if ((*here())->side) ((*here())->is_inside->cone) = 0;
760 }
while (here != start);
765 recompute_cone_contents();
776 void CSphstable_cones::recompute_cone_contents(){
787 for (i=0;i<vicinity_size;i++){
789 if ((vicinity[i]->side) && (vicinity[i]->is_inside->cone))
790 cone += *vicinity[i]->v;
804 void CSphstable_cones::recompute_cone_contents_if_needed(
CSphmomentum & this_cone,
807 if (this_dpt >
PT_TSHOLD*(fabs(this_cone.
px)+fabs(this_cone.
py))) {
808 if (cone.ref.is_empty()) {
819 for (
unsigned int i=0;i<vicinity_size;i++){
821 if ((vicinity[i]->side) && (vicinity[i]->is_inside->cone))
822 this_cone += *vicinity[i]->v;
853 for (i=0;i<n_part;i++){
855 if (is_closer(&cone_centre, &(plist[i]), tan2R))
856 intersection+=plist[i].
ref;
CSph3vector centre
centre of the cone
double _norm
particle spatial norm (available ONLY after a call to build_norm)
sph_hash_element * next
pointer to the next element
void set_position(const circulator< T > &other)
set just the position without resetting the begin and end elements
class for storing a border momentum (in context of co-circularity checks).
list of element in the vicinity of a parent.
siscone::Creference ref
reference number for the vector
a circulator that is foreseen to take as template member either a pointer or an iterator; ...
double cocircular_range
amount by which the angle can be varied while maintaining this point within co-circularity margin ...
double angle
angle with parent
CSphstable_cones()
default ctor
void get_angular_directions(CSph3vector &angular_dir1, CSph3vector &angular_dir2)
for this direction, compute the two reference directions used to measure angles
unsigned int ref[3]
actual data for the reference
~CSphstable_cones()
default dtor
bool not_empty()
test non-emptyness
#define PT_TSHOLD
program name
bool is_stable
true if stable w.r.t. "border particles"
void init(std::vector< CSphmomentum > &_particle_list)
initialisation
base class for dynamic coordinates management
list of cones candidates.
base class for managing the spatial part of Cmomentum (defined after)
element in the vicinity of a parent.
information on store cones candidates.
references used for checksums.
int get_stable_cones(double _radius)
compute stable cones.