Gui

ReadGadget/IOGadget.go

@@ -0,0 +1,81 @@
+package ReadGadget
+
+import (
+ "binary"
+ "fmt"
+ "log"
+ "os"
+)
+
+type Gadget struct {
+ name string
+ format int
+ Part []Particule
+ header Header
+ log *log.Logger
+}
+
+func New(name string, format int) *Gadget {
+ res := new(Gadget)
+ res.name = name
+ res.format = format
+ return res
+}
+
+func (c *Gadget) SetLogger(val *log.Logger) {
+ c.log = val
+}
+
+func (c *Gadget) GetLogger() *log.Logger {
+ return c.log
+}
+
+type UnsupportedFormat int
+
+func (c UnsupportedFormat) Error() string {
+ return fmt.Sprint("Unsupported format: ", c)
+}
+
+type NotImplemented string
+
+func (c NotImplemented) Error() string {
+ return fmt.Sprint("You are trying to use a not implemented method: ", c)
+}
+
+func (c *Gadget) Read() error {
+ switch c.format {
+ case 1:
+ return c.read_format1()
+ case 2:
+ return c.read_format2()
+ case 3:
+ return c.read_format3()
+ default:
+ return UnsupportedFormat(c.format)
+ }
+
+}
+
+func (c *Gadget) read_format1() error {
+ var size, nb_files, pc, pc_new int = 0, 1, 0, 0
+ var err error = nil
+ var file *os.File
+
+ for i := 0; i < nb_files; i++ {
+ if file, err = os.Open(c.name); err != nil {
+ return err
+ }
+
+ pc = pc_new
+ }
+
+ return NotImplemented("Read()->read_format1()")
+}
+
+func (c *Gadget) read_format2() error {
+ return NotImplemented("Read()->read_format2()")
+}
+
+func (c *Gadget) read_format3() error {
+ return NotImplemented("Read()->read_format3()")
+}

octree/octree.go

@@ -1,14 +1,17 @@
 //-tags: main
 package octree
 
-import rg "github.com/ElricleNecro/WisiGo/ReadGadget"
-import "fmt"
-import m "math"
+import (
+ "fmt"
+ rg "github.com/ElricleNecro/WisiGo/ReadGadget"
+ m "math"
+)
 
 var (
  UseGoRoutine = false //Use of Go Routine for the tree creation?
  // If you set this to true, don't forget to set runtime.GOMAXPROCS
  // to the correct value.
+ Ggrav = 6.67e-11 //Gravitational constant express in International system of units (L^3 M^{-1}T^{-2}).
 )
 
 // This function is used to create a tree using a list of the particle in the node, the center of the cube and his size. Her essential case of use will be for create the root of the tree.
@@ -68,69 +71,31 @@ func (e *Node) Create(NbMin int32) {
  //We are calculating the particle contain in the Node :
  NbUtil := t1.setPart()
 
- println((*t1).level, NbUtil, NbUse)
- print((*t1).level, " ")
- for z := 0; z < 3; z++ {
-  print("] ", (*t1).Center[z]-(*t1).Size/2., "; ", (*t1).Center[z]+(*t1).Size/2., "] ")
- }
- println("")
+ //println((*t1).level, NbUtil, NbUse)
+ //print((*t1).level, " ")
+ //for z := 0; z < 3; z++ {
+ //print("] ", (*t1).Center[z]-(*t1).Size/2., "; ", (*t1).Center[z]+(*t1).Size/2., "] ")
+ //}
+ //println("")
 
  //If there is particle in the node, we are actualizing the Particle slice, launching calculation of the son and preparing the next brother.
  if NbUtil != 0 {
  (*t1).Part = e.Part[NbUse:(NbUse + NbUtil)]
+ (*t1).Mass = 0.0
+ for _, v := range (*t1).Part {
+ (*t1).Mass += float64(v.Mass)
+ }
 
- fmt.Println("Going down!")
+ //fmt.Println("Going down!")
  (*t1).Create(NbMin)
- fmt.Println("Going up!")
+ //fmt.Println("Going up!")
 
  NbUse += NbUtil
  t1 = &(*t1).Frere
  }
  }
 }
 
-type Search struct {
- Radius float64
- Part rg.Particule
-}
-
-type ByDist []Search
-
-func (e ByDist) Len() int {
- return len(e)
-}
-
-func (e ByDist) Swap(i, j int) {
- e[i], e[j] = e[j], e[i]
-}
-
-func (e ByDist) Less(i, j int) bool {
- return e[i].Radius < e[j].Radius
-}
-
-func Insert(e []Search, to Search) {
- var tmp, bak Search
- var i int
-
- // We search where to place the particule :
- for i, _ = range e {
- if e[i].Radius > to.Radius {
- tmp = e[i]
- e[i] = to
- i++
- break
- }
- }
-
- // We move all particule from where we place 'to' to the end of the array
- // (we don't keep those who get out from the array) :
- for ; i < len(e); i++ {
- bak = e[i]
- e[i] = tmp
- tmp = bak
- }
-}
-
 func (e *Node) Nearest(part rg.Particule, NbVois int) []Search {
  var res []Search = make([]Search, NbVois)
 
@@ -163,12 +128,57 @@ func (e *Node) SearchNeighbor(part rg.Particule, searchy []Search) {
  return
  }
  if e.Fils != nil {
- e.Fils.SearchNeighbor(part, searchy[:])
+ //e.Fils.SearchNeighbor(part, searchy[:])
+ for t1 := e.Fils; t1 != nil; t1 = t1.Frere {
+ t1.SearchNeighbor(part, searchy[:])
+ }
  } else {
  e.fill_neighboorhood(part, searchy[:])
  }
 }
 
+func (e *Node) GetBrother() *Node {
+ return e.Frere
+}
+
+func (e *Node) GetSon() *Node {
+ return e.Fils
+}
+
+func (e *Node) TotalPotential(opening float64) float64 {
+ var pot float64 = 0.
+ for _, v := range e.Part {
+ pot += e.Potential(v, opening)
+ }
+ return pot * 0.5
+}
+
+func (e *Node) Potential(part rg.Particule, opening float64) float64 {
+ //We can process the Node and it has sons:
+ if (float64)(e.Size)/e.Dist(part) >= opening && e.Fils != nil {
+ var pot float64 = 0.0
+ for t1 := e.Fils; t1 != nil; t1 = t1.Frere {
+ pot += t1.Potential(part, opening)
+ }
+ return pot
+ }
+
+ //The previous condition is wrong, so we have two possibilities:
+ if e.Fils == nil {
+ //we are in a leaf:
+ var pot float64 = 0.
+ for _, v := range e.Part {
+ if v.Dist(part) != 0.0 {
+ pot += Ggrav * (float64)(v.Mass/v.Dist(part))
+ }
+ }
+ return pot
+ } else {
+ //the node is too far or too small to be processed, we do some approximation:
+ return Ggrav * e.Mass / m.Sqrt((float64)((e.Center[0]-part.Pos[0])*(e.Center[0]-part.Pos[0])+(e.Center[1]-part.Pos[1])*(e.Center[1]-part.Pos[1])+(e.Center[2]-part.Pos[2])*(e.Center[2]-part.Pos[2]))) //e.Dist(part)
+ }
+}
+
 func (e *Node) setPart() (NbUse int32) {
  NbUse = 0
  for i, _ := range e.Part {

octree/octree_test.go

@@ -1,7 +1,6 @@
 package octree
 
 import (
- //"fmt"
  rg "github.com/ElricleNecro/WisiGo/ReadGadget"
  "math/rand"
  "os"
@@ -245,3 +244,80 @@ func TestSetPart(t *testing.T) {
  root.Create(1)
  root.setPart()
 }
+
+func TestPotential(t *testing.T) {
+ rdm := rand.New(rand.NewSource(int64(33)))
+ nb_part := 10
+ part := make([]rg.Particule, nb_part)
+ var testy rg.Particule
+
+ for i := 0; i < nb_part; i++ {
+ for j, _ := range part[i].Pos {
+ part[i].Pos[j] = 2.0*rdm.Float32() - 1.0
+ part[i].Id = int64(i + 1)
+ part[i].Mass = 1.0
+ }
+ }
+
+ center := [...]float32{0., 0., 0.}
+ root := New(part, center, 2.0)
+ Ggrav = 1.0
+
+ testy.Mass = 1.0
+ testy.Id = 0
+ for i, _ := range testy.Pos {
+ testy.Pos[i] = 0.
+ }
+
+ pot := root.Potential(testy, 0.0)
+ pot2 := bruteforcePotential(part, testy)
+
+ if pot != pot2 {
+ t.Fatal("Tree Code is different from brute force:", pot, " != ", pot2)
+ }
+}
+
+func TestTotalPotential(t *testing.T) {
+ rdm := rand.New(rand.NewSource(int64(33)))
+ nb_part := 10
+ part := make([]rg.Particule, nb_part)
+ var testy rg.Particule
+
+ for i := 0; i < nb_part; i++ {
+ for j, _ := range part[i].Pos {
+ part[i].Pos[j] = 2.0*rdm.Float32() - 1.0
+ part[i].Id = int64(i + 1)
+ part[i].Mass = 1.0
+ }
+ }
+
+ center := [...]float32{0., 0., 0.}
+ root := New(part, center, 2.0)
+ Ggrav = 1.0
+
+ testy.Mass = 1.0
+ testy.Id = 0
+ for i, _ := range testy.Pos {
+ testy.Pos[i] = 0.
+ }
+
+ pot := root.TotalPotential(0.0)
+ var pot2 float64 = 0.0
+ for _, v := range part {
+ pot2 += bruteforcePotential(part, v)
+ }
+
+ if pot != pot2 {
+ t.Fatal("Tree Code is different from brute force:", pot, " != ", pot2)
+ }
+}
+
+func bruteforcePotential(part []rg.Particule, testy rg.Particule) float64 {
+ var pot float64 = 0.0
+ for _, v := range part {
+ if v.Dist(testy) != 0.0 {
+ pot += Ggrav * (float64)(v.Mass/v.Dist(testy))
+ }
+ }
+ return pot
+}

octree/types.go

@@ -7,10 +7,18 @@ type Node struct {
  Parent, Fils, Frere *Node //These is the Main Node to Now : the parent, the first son and the next brother.
  Center [3]float32 //Center of the cube.
  Size float32 //Side size of the cube.
+ Mass float64 //Total Mass contain into the Cube.
  Part []rg.Particule //List of the particle contain in the node.
  level int64 //Level of the tree.
 }
 
 type WriteStringer interface {
  WriteString(s string) (ret int, err error)
 }
+
+type Search struct {
+ Radius float64
+ Part rg.Particule
+}
+
+type ByDist []Search

octree/types_func.go

@@ -0,0 +1,36 @@
+package octree
+
+func (e ByDist) Len() int {
+ return len(e)
+}
+
+func (e ByDist) Swap(i, j int) {
+ e[i], e[j] = e[j], e[i]
+}
+
+func (e ByDist) Less(i, j int) bool {
+ return e[i].Radius < e[j].Radius
+}
+
+func Insert(e []Search, to Search) {
+ var tmp, bak Search
+ var i int
+
+ // We search where to place the particule :
+ for i, _ = range e {
+ if e[i].Radius > to.Radius {
+ tmp = e[i]
+ e[i] = to
+ i++
+ break
+ }
+ }
+
+ // We move all particule from where we place 'to' to the end of the array
+ // (we don't keep those who get out from the array) :
+ for ; i < len(e); i++ {
+ bak = e[i]
+ e[i] = tmp
+ tmp = bak
+ }
+}