implemented methods calculating the center of mass, the total mass of given trees

1 files changed, 128 insertions, 30 deletions

diff --git a/quadtree.go b/quadtree.go
index 4a84843..dbbf229 100644
--- a/quadtree.go
+++ b/quadtree.go
@@ -3,6 +3,8 @@ package structs
 import (
 	"fmt"
 	"io/ioutil"
+	"log"
+	"math"
 	"os/exec"
 )
 
@@ -205,53 +207,149 @@ func (n Node) GetAllStars() []Star2D {
 }
 
 // CalcCenterOfMass calculates the center of mass for every node in the treee
-func (n *Node) CalcCenterOfMass() {
-	meanPosX := 0.0
-	meanPosY := 0.0
-	totalMass := 0.0
+func (n *Node) CalcCenterOfMass() Vec2 {
+	fmt.Printf("[CalcCenterOfMass] %v\n", n.Boundry.Width)
 
-	// if the subtree is not empty
+	if n.Boundry.Width < 10 {
+		return Vec2{}
+	}
+
+	// if the subtrees are not empty
 	if n.Subtrees != ([4]*Node{}) {
+		for i := 0; i < len(n.Subtrees); i++ {
+			// get the center of mass for the subtree
+			subtreeCenterOfMassX := n.Subtrees[i].CalcCenterOfMass().X
+			subtreeCenterOfMassY := n.Subtrees[i].CalcCenterOfMass().Y
 
-		// iterate over all the subtrees
-		for index, _ := range n.Subtrees {
+			// get the totalmass of the subtree
+			subtreeTotalMass := n.Subtrees[i].TotalMass
 
-			// calculate the center of mass of the subtree
-			n.Subtrees[index].CalcCenterOfMass()
+			// calculate the center of mass of both
+			n.CenterOfMass.X = subtreeCenterOfMassX * subtreeTotalMass
+			n.CenterOfMass.Y = subtreeCenterOfMassY * subtreeTotalMass
+		}
+	}
 
-			if n.Subtrees[index].TotalMass == 0 {
-				n.Subtrees[index].CalcTotalMass()
-			}
+	if n.Star != (Star2D{}) {
+		n.CenterOfMass = n.Star.C
+	}
 
-			meanPosX += n.Subtrees[index].CenterOfMass.X * n.Subtrees[index].TotalMass
-			meanPosY += n.Subtrees[index].CenterOfMass.Y * n.Subtrees[index].TotalMass
-			totalMass += n.Subtrees[index].TotalMass
+	n.CenterOfMass.X = n.CenterOfMass.X / n.TotalMass
+	n.CenterOfMass.Y = n.CenterOfMass.Y / n.TotalMass
+
+	return n.CenterOfMass
+}
+
+// CalcTotalMass calculates the total mass for every node in the tree
+func (n *Node) CalcTotalMass() float64 {
+
+	// if the subtrees are not empty
+	if n.Subtrees != ([4]*Node{}) {
+		for i := 0; i < len(n.Subtrees); i++ {
+			n.TotalMass += n.Subtrees[i].CalcTotalMass()
 		}
 	}
-	meanPosX = meanPosX / totalMass
-	meanPosY = meanPosY / totalMass
 
-	n.TotalMass = totalMass
-	n.CenterOfMass = Vec2{meanPosX, meanPosY}
+	// if the star in the subtree is not empty
+	if n.Star != (Star2D{}) {
+		n.TotalMass += n.Star.M
+	}
+
+	return n.TotalMass
 }
 
-// CalcTotalMass calculates the total mass for every node in the tree
-func (n *Node) CalcTotalMass() {
-	totalMass := 0.0
+func (n Node) CalcAllForces(star Star2D, theta float64) Vec2 {
+	// initialize a variable storing the overall force
+	var localForce Vec2 = Vec2{}
+	log.Printf("[CalcAllforces] Boundary Width: %f", n.Boundry.Width)
+
+	// calculate the local theta
+	var tmpX float64 = math.Pow(star.C.X-n.Star.C.X, 2)
+	var tmpY float64 = math.Pow(star.C.Y-n.Star.C.Y, 2)
+	var distance float64 = math.Sqrt(tmpX + tmpY)
+
+	log.Printf("[CalcAllforces] n.Boundary.Width=%f", n.Boundry.Width)
+	log.Printf("[CalcAllforces] distance=%f", distance)
+	var localtheta float64 = n.Boundry.Width / distance
+	log.Printf("[CalcAllforces] localtheta=%f", localtheta)
 
 	// if the subtree is not empty
 	if n.Subtrees != ([4]*Node{}) {
+		if localtheta < theta {
+			log.Printf("[CalcAllforces] not recursing deeper ++++++++++++++++++++++++ ")
+			// don't recurse further into the tree
+			// calculate the forces inbetween the star and the node
+
+			nodeStar := Star2D{
+				C: Vec2{
+					X: n.CenterOfMass.X,
+					Y: n.CenterOfMass.Y,
+				},
+				V: Vec2{
+					X: 0,
+					Y: 0,
+				},
+				M: n.TotalMass,
+			}
 
-		// iterate over all the subtrees
-		for index, _ := range n.Subtrees {
-
-			totalMass += n.Subtrees[index].Star.M
+			if star != nodeStar {
+				log.Printf("[CalcAllforces] (%v) < +++++++ > (%v)\n", star, nodeStar)
+				force := CalcForce(star, nodeStar)
+				localForce.X += force.X
+				localForce.Y += force.Y
+			}
 
-			// calculate the mass of the subtree and add it to the totalmass
-			n.Subtrees[index].CalcTotalMass()
-			totalMass += n.Subtrees[index].TotalMass
+		} else {
+			log.Printf("[CalcAllforces] recursing deeper ----------------------------")
+			// recurse deeper
+			for i := 0; i < len(n.Subtrees); i++ {
+				force := n.Subtrees[i].CalcAllForces(star, theta)
+				localForce.X += force.X
+				localForce.Y += force.Y
+			}
+		}
+	} else {
+		if n.Star != (Star2D{}) {
+			if star != n.Star {
+				log.Printf("[CalcAllforces] not recursing deeper ====================== ")
+				log.Printf("[CalcAllforces] (%v) < ------- > (%v)\n", star, n.Star)
+				force := CalcForce(star, n.Star)
+				localForce.X += force.X
+				localForce.Y += force.Y
+			}
 		}
 	}
 
-	n.TotalMass = totalMass
+	log.Printf("[CalcAllforces] localforce: %v +-+-+-+-+-+-+-+-+-+- ", localForce)
+	return localForce
+
+}
+
+// CalcForce calculates the force exerted on s1 by s2 and returns a vector representing that force
+func CalcForce(s1 Star2D, s2 Star2D) Vec2 {
+	G := 6.6726 * math.Pow(10, -11)
+
+	// calculate the force acting
+	var combinedMass float64 = s1.M * s2.M
+	var distance float64 = math.Sqrt(math.Pow(math.Abs(s1.C.X-s2.C.X), 2) + math.Pow(math.Abs(s1.C.Y-s2.C.Y), 2))
+	fmt.Printf("\t- combinedMass: %f\n", combinedMass)
+	fmt.Printf("\t- distance: %f\n", distance)
+	fmt.Printf("\t- distance squared: %f\n", math.Pow(distance, 2))
+	fmt.Printf("\t- combinedMass / distance: %f\n", combinedMass/math.Pow(distance, 2))
+
+	var scalar float64 = G * ((combinedMass) / math.Pow(distance, 2))
+	fmt.Printf("\t- Scalar: %f\n", scalar)
+
+	// define a unit vector pointing from s1 to s2
+	var vector Vec2 = Vec2{s2.C.X - s1.C.X, s2.C.Y - s1.C.Y}
+	var UnitVector Vec2 = Vec2{vector.X / distance, vector.Y / distance}
+	fmt.Printf("\t- Vector: %v\n", vector)
+	fmt.Printf("\t- UnitVector: %v\n", UnitVector)
+
+	// multiply the vector with the force to get a vector representing the force acting
+	var force Vec2 = UnitVector.Multiply(scalar)
+	fmt.Printf("\t- force: %v\n", force)
+
+	// return the force exerted on s1 by s2
+	return force
 }