diff --git a/ReadGadget/IOGadget.go b/ReadGadget/IOGadget.go new file mode 100644 index 0000000..6b615ff --- /dev/null +++ b/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()") +} diff --git a/octree/octree.go b/octree/octree.go index d5fea4e..f0e8f12 100644 --- a/octree/octree.go +++ b/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,20 +71,24 @@ 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 @@ -89,48 +96,6 @@ func (e *Node) Create(NbMin int32) { } } -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 { diff --git a/octree/octree_test.go b/octree/octree_test.go index acad83b..10e3048 100644 --- a/octree/octree_test.go +++ b/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 +} diff --git a/octree/types.go b/octree/types.go index ca65217..03778a5 100644 --- a/octree/types.go +++ b/octree/types.go @@ -7,6 +7,7 @@ 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. } @@ -14,3 +15,10 @@ type Node struct { type WriteStringer interface { WriteString(s string) (ret int, err error) } + +type Search struct { + Radius float64 + Part rg.Particule +} + +type ByDist []Search diff --git a/octree/types_func.go b/octree/types_func.go new file mode 100644 index 0000000..ffe4748 --- /dev/null +++ b/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 + } +}