package main

import (
	"fmt"
	"github.com/veandco/go-sdl2/sdl"
	"math"
	"math/rand"
)

type Particule struct {
	// À l'instant présent
	m, // Masse
	x0, y0, z0, // Position
	vx0, vy0, vz0, // Vitesse
	ax0, ay0, az0, // Accélération
	px0, py0, pz0, // Pente

	// À l'instant dt suivant
	x1, y1, z1, // Position
	vx1, vy1, vz1, // Vitesse
	ax1, ay1, az1, // Accélération
	fx1, fy1, fz1, // Force s'excerçant sur la particule

	// Durant l'intervalle de temps dt précédent
	px_, py_, pz_ float64 // Pente

}

const N = 4           // Nombre de particules
var P [N]Particule    // Liste des N particules
const epsilon = 1e-15 // Erreur tolérée à chaque itération
var G = 100.0         // Constante des forces de gravité

var dt, t, // Intervalle de temps et temps cumulé
	r, r2, // Distance entre deux particules et son carré
	rx, ry, rz, // Vecteur de distance entre deux particules
	ux, uy, uz, // Vecteur unitaire entre deux particules
	fx, fy, fz, // Force
	Ec, Ep, // Energie cinétique et energie potentielle
	mp float64 // Maximum des variations de la pente

var dilatation float64 = 1
var Sym = false
var Has = false
var zoom = float64(1)
var rmin = float64(2)
var rmin0 = float64(2)

const fd = float64(1.000001)

func affichePresentation() {
	fmt.Println("Modèle newtonien d'attraction de N corps ")
	fmt.Println("-----------------------------------------")
	fmt.Println("Commandes :")
	fmt.Println("Touche c -> (Clear) Efface l'écran ")
	fmt.Println("Touche d -> (Dilate) Multiplie toutes les distances de 1.000001 ")
	fmt.Println("Touche s -> (Symétrise) Resymétrise les particules ")
	fmt.Println("Touche h -> (Hasard) Perturbation des vitesses au hasard ")
	fmt.Println("Touche F1 -> Flip-flop - symétrise en permanence ")
	fmt.Println("Touche F2 -> Flip-flop - Hasard en permanence ")
	fmt.Println("Touche + -> Zoom ")
	fmt.Println("Touche - -> Dézoom ")
	fmt.Println("Touche espace -> Affiche le temps et l'énergie totale ")
	fmt.Println("Touche echape -> Quitte le programme ")
	fmt.Println("-----------------------------------------")
}
func symetrisation() {
	a := &P[0]
	b := &P[1]
	a.x0 = (a.x0 + 1000 - b.x0) / 2
	a.y0 = (a.y0 + 1000 - b.y0) / 2
	a.z0 = (a.z0 + 1000 - b.z0) / 2
	a.x1 = (a.x1 + 1000 - b.x1) / 2
	a.y1 = (a.y1 + 1000 - b.y1) / 2
	a.z1 = (a.z1 + 1000 - b.z1) / 2
	a.z1 = (a.z1 + 1000 - b.z1) / 2
	a.vx0 = (a.vx0 - b.vx0) / 2
	a.vy0 = (a.vy0 - b.vy0) / 2
	a.vz0 = (a.vz0 - b.vz0) / 2
	a.vx1 = (a.vx1 - b.vx1) / 2
	a.vy1 = (a.vy1 - b.vy1) / 20
	a.vz1 = (a.vz1 - b.vz1) / 2
	b.x0 = 1000 - a.x0
	b.y0 = 1000 - a.y0
	b.z0 = 1000 - a.z0
	b.x1 = 1000 - a.x1
	b.y1 = 1000 - a.y1
	b.z1 = 1000 - a.z1
	b.vx0 = -a.vx0
	b.vy0 = -a.vy0
	b.vz0 = -a.vz0
	b.vx1 = -a.vx1
	b.vy1 = -a.vy1
	b.vz1 = -a.vz1

	a = &P[2]
	b = &P[3]
	a.x0 = (a.x0 + 1000 - b.x0) / 2
	a.y0 = (a.y0 + 1000 - b.y0) / 2
	a.z0 = (a.z0 + 1000 - b.z0) / 2
	a.x1 = (a.x1 + 1000 - b.x1) / 2
	a.y1 = (a.y1 + 1000 - b.y1) / 2
	a.z1 = (a.z1 + 1000 - b.z1) / 2
	a.z1 = (a.z1 + 1000 - b.z1) / 2
	a.vx0 = (a.vx0 - b.vx0) / 2
	a.vy0 = (a.vy0 - b.vy0) / 2
	a.vz0 = (a.vz0 - b.vz0) / 2
	a.vx1 = (a.vx1 - b.vx1) / 2
	a.vy1 = (a.vy1 - b.vy1) / 20
	a.vz1 = (a.vz1 - b.vz1) / 2
	b.x0 = 1000 - a.x0
	b.y0 = 1000 - a.y0
	b.z0 = 1000 - a.z0
	b.x1 = 1000 - a.x1
	b.y1 = 1000 - a.y1
	b.z1 = 1000 - a.z1
	b.vx0 = -a.vx0
	b.vy0 = -a.vy0
	b.vz0 = -a.vz0
	b.vx1 = -a.vx1
	b.vy1 = -a.vy1
	b.vz1 = -a.vz1
}
func hasard() {
	for i := 0; i < N; i++ {
		P[i].vx0 += (rand.Float64() - 0.5) / 1e7
		P[i].vy0 += (rand.Float64() - 0.5) / 1e7
		P[i].vz0 += (rand.Float64() - 0.5) / 1e7
	}
}

func main() {
	affichePresentation()

	// initialise une fenêtre de fond noire 1000×1000
	sdl.Init(sdl.INIT_EVERYTHING)
	defer sdl.Quit()
	w, _ := sdl.CreateWindow("modèle", 0, 0, 1000, 1000, sdl.WINDOW_SHOWN)
	defer w.Destroy()
	/*
		W, _ := sdl.CreateRenderer(w, -1, sdl.RENDERER_ACCELERATED)
		defer W.Destroy()
		W.Clear()
	*/
	ws, err := w.GetSurface()
	if err != nil {
		panic(err)
	}

	var np, // Nombre de pas
		k, // Nombre d'iterations par pas
		i, j int // Indice des particules

	var p = new(Particule) // Référence d'une particule
	var q = new(Particule) // Référence d'une particule

	// Etats des particules à l'instant 0
	P[0] = Particule{1, 400, 500, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
	P[1] = Particule{1, 600, 500, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
	P[2] = Particule{1, 500, 500, -100, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
	P[3] = Particule{1, 500, 500, 100, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

	/*	Changement d'échelle laissant inchangé la forme des trajectoires.
		fm := float64(1)
		for i = 0; i < N; i++ {
			P[i].x0 *= fm
			P[i].y0 *= fm
			P[i].z0 *= fm
			P[i].vx0 /= math.Sqrt(fm)
			P[i].vy0 /= math.Sqrt(fm)
			P[i].vz0 /= math.Sqrt(fm)
		}
	*/
	t = 0.0    // Temps=0
	dt = 1e-10 // intervalle initiale

	// Boucle sans fin, exécute une infinité de pas
	_B := true
	for np = 0; _B; np++ {

		for _e := sdl.PollEvent(); _e != nil; _e = sdl.PollEvent() {
			switch _t := _e.(type) {
			case *sdl.KeyboardEvent:
				if _t.Keysym.Sym == sdl.K_ESCAPE {
					_B = false
				}
				if _t.Keysym.Sym == sdl.K_SPACE && _t.State == 1 {
					// Calcule de l'énergie totale Ep+Ec
					Ec = 0
					for i = 0; i < N; i++ {
						p = &P[i]
						Ec += p.m * (p.vx1*p.vx1 + p.vy1*p.vy1 + p.vz1*p.vz1) / 2
					}
					Ep = 0
					for i = 1; i < N; i++ {
						p = &P[i]
						for j = 0; j < i; j++ {
							q = &P[j]

							rx = q.x1 - p.x1
							ry = q.y1 - p.y1
							rz = q.z1 - p.z1

							r = math.Sqrt(rx*rx + ry*ry + rz*rz)

							Ep += -q.m / r
						}
					}
					Ep *= G
					fmt.Printf("Temps: %.0f \tEnergie totale: ", t)
					fmt.Println(float64(Ec + Ep))
				}
				if _t.Keysym.Sym == sdl.K_d && _t.State == 1 {
					dilatation *= fd
					for i = 0; i < N; i++ {
						P[i].x0 *= fd
						P[i].y0 *= fd
						P[i].z0 *= fd
						P[i].x1 *= fd
						P[i].y1 *= fd
						P[i].z1 *= fd
					}
					fmt.Println("Dilatation:  %f : ", dilatation)
				}
				if _t.Keysym.Sym == sdl.K_c && _t.State == 1 {
					ws.FillRect(nil, 0x00000000)
				}
				if _t.Keysym.Sym == sdl.K_s && _t.State == 1 {
					symetrisation()
					fmt.Println("Symétrisation ")
				}
				if _t.Keysym.Sym == sdl.K_F1 && _t.State == 1 {
					Sym = !Sym
					if Sym {
						fmt.Println("Symétrisation permanente ")
					} else {
						fmt.Println("Arrêt de la symétrisation ")
					}
				}
				if _t.Keysym.Sym == sdl.K_h && _t.State == 1 {
					hasard()
					fmt.Println("Perturbation ")
				}
				if _t.Keysym.Sym == sdl.K_F2 && _t.State == 1 {
					Has = !Has
					if Has {
						fmt.Println("Perturbation permanente ")
					} else {
						fmt.Println("Arrêt de la perturbation ")
					}
				}
				if _t.Keysym.Sym == sdl.K_KP_PLUS && _t.State == 1 {
					zoom /= 0.75
					fmt.Printf("Echelle: %.1f \n", zoom)
				}
				if _t.Keysym.Sym == sdl.K_KP_MINUS && _t.State == 1 {
					zoom *= 0.75
					fmt.Printf("Echelle: %.1f \n", zoom)
				}

			}
		}

		if Sym {
			symetrisation()
		}
		if Has {
			hasard()
		}

		// Répéter 3 fois
		for k = 0; k < 3; k++ {

			// Pour chaque particule faire
			for i = 0; i < N; i++ {
				p = &P[i]

				// V1 = V0 + A0*dt + P0*dt^2/2
				// X1 = X0 + V0*dt + A0*dt^2/2 + P0*dt^3/6

				p.vx1 = p.vx0 + dt*(p.ax0+dt*p.px0/2)
				p.vy1 = p.vy0 + dt*(p.ay0+dt*p.py0/2)
				p.vz1 = p.vz0 + dt*(p.az0+dt*p.pz0/2)

				p.x1 = p.x0 + dt*(p.vx0+dt*(p.ax0+dt*p.px0/3)/2)
				p.y1 = p.y0 + dt*(p.vy0+dt*(p.ay0+dt*p.py0/3)/2)
				p.z1 = p.z0 + dt*(p.vz0+dt*(p.az0+dt*p.pz0/3)/2)

			}

			// Pour chaque particule
			for i = 0; i < N; i++ {
				p = &P[i]

				// f1 = 0

				p.fx1 = 0
				p.fy1 = 0
				p.fz1 = 0
			}

			// Pour chaque ensemble de deux particules {p,q}
			for i = 1; i < N; i++ {
				p = &P[i]
				for j = 0; j < i; j++ {
					q = &P[j]

					// R = q.X1-p.X1
					// r = norm(R)
					// r2 = r*r
					// U = R/r
					// F = U*G*q.m*p.m/r2
					// p.F = p.F + F
					// q.F = q.F - F

					rx = q.x1 - p.x1
					ry = q.y1 - p.y1
					rz = q.z1 - p.z1

					r = math.Sqrt(rx*rx + ry*ry + rz*rz)
					r2 = r * r

					if r < rmin {
						rmin = r
					}

					ux = rx / r
					uy = ry / r
					uz = rz / r

					fx = ux * G * q.m * p.m / r2
					fy = uy * G * q.m * p.m / r2
					fz = uz * G * q.m * p.m / r2

					p.fx1 += fx
					p.fy1 += fy
					p.fz1 += fz

					q.fx1 -= fx
					q.fy1 -= fy
					q.fz1 -= fz

				}
			}

			mp = 1e-10 // mp = infiniment petit

			// Pour chaque particule
			for i = 0; i < N; i++ {
				p = &P[i]

				// A1 = F1/m
				// P0 = (A1-A0)/dt

				p.ax1 = p.fx1 / p.m
				p.ay1 = p.fy1 / p.m
				p.az1 = p.fz1 / p.m

				p.px0 = (p.ax1 - p.ax0) / dt
				p.py0 = (p.ay1 - p.ay0) / dt
				p.pz0 = (p.az1 - p.az0) / dt

				// R = P0-Ppred
				// r = Norm(R)
				// if r>mp then mp=r

				rx = q.px0 - p.px_
				ry = q.py0 - p.py_
				rz = q.pz0 - p.pz_

				r = math.Sqrt(rx*rx + ry*ry + rz*rz)
				if r > mp {
					mp = r
				}
			}
			dt = math.Pow(24*epsilon/mp, 0.3333) // Recalcul de dt = (24*epsilon/mp)^(1/3)

			//---------------------------------------------------------------
		}
		t += dt // Temps t=t+dt

		// Pour chaque particule
		for i = 0; i < N; i++ {
			p = &P[i]

			// X0 = X1
			// V0 = V1
			// A0 = A1
			// Ppred = P0

			p.x0 = p.x1
			p.y0 = p.y1
			p.z0 = p.z1

			p.vx0 = p.vx1
			p.vy0 = p.vy1
			p.vz0 = p.vz1

			p.ax0 = p.ax1
			p.ay0 = p.ay1
			p.az0 = p.az1

			p.px_ = p.px0
			p.py_ = p.py0
			p.pz_ = p.pz0

		}

		// Pour chaque particule
		for i = 0; i < N; i++ {
			p = &P[i]

			// Afficher le point X
			/*
				W.SetDrawColor(255, 255, 255, 255)
				W.DrawPoint(int32(p.x1), int32(p.y1))
			*/

			ws.FillRect(&sdl.Rect{int32((p.x1-500)*zoom + 500), int32((p.y1-500)*zoom + 500), 1, 1}, 0xffffffff)
		}

		x := math.Mod(float64(np), 10000.0)
		if x == 0 {
			// W.Present() // Actualiser l'affichage de la fenêtre
			w.UpdateSurface()
			// fmt.Print("t=", t, " | ", Ep+Ec, "\n")
			if rmin0 < 2 && rmin >= 2 {
				fmt.Printf("Sortie collision, distance minimum : %.2f \n", rmin0)
				rmin0 = 2
			}
			if rmin < 2 && rmin0 >= 2 {
				fmt.Printf("Entré collision (distance < 2) à l'intant %.0f \n", t)
				rmin0 = rmin
			}
			if rmin < 2 {
				if rmin0 > rmin {
					rmin0 = rmin
				}
				rmin = 2
			}

		}

	}
	ws.Free()
	w.Destroy()
	sdl.Quit()
}