You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('') and can be up to 35 characters long.
83 lines
3.0 KiB
83 lines
3.0 KiB
/*


* Copyright (c) 20192020 Thomas Kramer.


*


* This file is part of LibrEDA


* (see https://codeberg.org/libreda).


*


* This program is free software: you can redistribute it and/or modify


* it under the terms of the GNU Affero General Public License as


* published by the Free Software Foundation, either version 3 of the


* License, or (at your option) any later version.


*


* This program is distributed in the hope that it will be useful,


* but WITHOUT ANY WARRANTY; without even the implied warranty of


* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


* GNU Affero General Public License for more details.


*


* You should have received a copy of the GNU Affero General Public License


* along with this program. If not, see <http://www.gnu.org/licenses/>.


*/




//! Compute electrostatic repulsion between many particles using a quadtree for approximation.




use super::types::*;


use super::quadtree::QuadTree;




/// Compute electrostatic repulsion between many particles using a quadtree aproximation


/// for efficient computation.


///


/// * `particles`: All the particles.


/// * `threshold`: Accurracy setting for the quad tree. `0` is the most inaccurate `1` is very accurate but slow.


/// Decent values are around `0.5`.


pub fn calculate_electric_forces(particles: &Vec<Particle>, threshold: f64) > Vec<Cost> {


if particles.is_empty() {


Vec::new()


} else {


// Find upper right and lower left corner by finding maxima and minima of x and y coordinates.




let (lower_left, upper_right) =


{


// Use first particle as the starting value to find minimum and maximum.


let loc0 = (particles[0].location.x, particles[0].location.y);


// Find lowerleft and upperright corner of the bounding box around the particles.


particles


.iter()


.fold((loc0, loc0),


((xmin, ymin), (xmax, ymax)), p


{


let l = p.location;


((xmin.min(l.x), ymin.min(l.y)),


(xmax.max(l.x), ymax.max(l.y)))


},


)


};




// Find minimum side length for the squareshaped quadtree such that all particles fit in.


let side_length = {


let (x1, y1) = lower_left;


let (x2, y2) = upper_right;


(x2  x1).max(y2  y1)


} * 1.001;




// Create an empty quadtree.


let mut qt = QuadTree::new(lower_left, side_length);




// Insert particles into the tree.


for &p in particles.iter() {


qt.insert(p);


}




// Precalculate centers of charge.


qt.update_center_of_charge();




// Calculate the force on each particle.


let f_df = particles.iter()


.map(p {


qt.force_onto_particle(&p, threshold)


})


.collect();




// Return (force, force_derivative)


f_df


}


}
