program PLANETS

  implicit none

! Arguments

	integer, allocatable, dimension(:) :: seed
	integer                            :: taille_seed

	integer :: i
	integer :: n_part
	real(8) :: v_ini, alpha, v_med, ecart_v, a_med, ecart_a
	real(8) :: g, masse, time
	real(8) :: r_x, r_z, r, theta
	real(8) :: Pi = 3.141592654_8
	real(8) :: alea1, alea2
	real(8) :: r_z_med =0.0_8, r_x_med = 0.0_8, r_x_max = 0.0_8, r_z_max = 0.0_8
  
	real(8), allocatable, dimension(:, :) :: resu


! Initialisation
    
	time  = 10_8
	g     = 1.113_8
	masse = 1.0e-2_8
	
	v_med   = 20.0_8
	ecart_v = 1.0_8
	a_med   = Pi/2.0_8
	ecart_a = Pi/50.0_8
	
	
!	print*, "Donnez le nombre de particules et le temps de chute, v, alpha, g"
!	read*,  n_part
!	read*,  time
!	read*,  v_med
!	read*,  a_med
!	read*,  g

n_part = 100000
time = 1
v_med = 0.1_8
a_med = 90.0_8
g = 0.0_8

	
	call lancement_hasard
	
	allocate(resu(n_part, 4))  
	
	open(2, file = "resultat.dat")

!----------------------------------------------------------------------------

	do i = 1, n_part

		! Détermination de l'angle et da la vitesse de départ avec de l'aléatoire
		call normale(alea1, alea2)
		v_ini = v_med + ecart_v*alea1
		alpha = a_med + ecart_a*alea2
		alpha = alpha*(Pi/180.0_8)
	!	print*, v_ini, alpha
		
		! Calcul de r_x et r_z en tenant compte du temps time
		if (time .lt. (2_8/g)*v_ini*sin(alpha)) then
			r_x = v_ini*cos(alpha)*time
			r_z = -(g/2)*(time**2) + v_ini*sin(alpha)*time
		else
			r_x = (2_8/g)*(v_ini**2) * cos(alpha)*sin(alpha)
			r_z = 0.0_8
		endif
		
		! Calcul de r et theta, coordonnées sphérique au temps time
		r     = sqrt(r_x**2 + r_z**2)
		theta = atan2(r_z,r_x)
		
		! Sauvegarde des valeurs
		resu(i, 1) = r
		resu(i, 2) = theta
		resu(i, 3) = r_z
		resu(i, 4) = r_x
		
	!	print*, "v ini =", v_ini, "alpha =", alpha*180/Pi
	!	print*, "r_z =", resu(i, 3), "r_x =", resu(i, 4),                  &
	!	        "r =", resu(i, 1), "theta =", resu(i, 2)*180/Pi
		        
		write(2,*) r_x, r_z
	enddo
	
	do i = 1, n_part
		r_z_max = max(r_z_max, resu(i,3))
		r_x_max = max(r_x_max, resu(i,4))
		r_z_med = r_z_med +resu(i,3)
		r_x_med = r_x_med +resu(i,4)
	end do
	
	print*, "La hauteur moyenne est de : ", r_z_med/n_part
	print*, "La largeur moyenne est de : ", r_x_med/n_part
	print*, "La hauteur maximal est de : ", r_z_max
	print*, "La largeur maximal est de : ", r_x_max
!----------------------------------------------------------------------------

! Fin du programme on stoppe le générateur aléatoire, on déloue resu et on
! ferme resultat.dat

	call arret_hasard
  
	deallocate(resu)
	
	close(2)

!////////////////////////////////////////////////////////////////////////////

  
contains


!////////////////////////////////////////////////////////////////////////////



subroutine lancement_hasard

  implicit none
  
! Arguments
  integer :: k, ios


! Ouverture de fichier seed.dat
  open(unit = 1, file = "seed.dat", status = "old", iostat = ios)


! On appelle random et on alloue à seed
  call random_seed (size = taille_seed)
  allocate (seed(taille_seed))


! Remplissage de seed
  if (ios == 0) then
     read(1, *) (seed(k) , k=1, taille_seed)
  else
     close(1)
     open(unit = 1, file = "seed.dat")
     seed = (/ ( 2*k, k=1,taille_seed) /)
  endif

! On place seed
  call random_seed (put = seed)

  
end subroutine lancement_hasard


!////////////////////////////////////////////////////////////////////////////


subroutine arret_hasard

  implicit none
  
! Arguments
  integer :: k


! On récupère seed, on le place dans seed.dat, on déloue et on ferme les
! fichiers
  call random_seed( get = seed )

  close(1) ! Pour écrire sur la première ligne
  open(unit = 1, file = "seed.dat")
  write(1, *) (seed(k), k=1, taille_seed)

  deallocate(seed)

  close(1)

 
end subroutine arret_hasard


!////////////////////////////////////////////////////////////////////////////


real(kind = 8) function alea_continue(mini , maxi )

  implicit none
  
! Arguments
  real(kind = 8), intent(in) :: mini, maxi !intent(in) : trucs qui sont rentrés mais pas changeables
  										  				 !intent(out) : trucs qui doivent être modifiés
  

! Calcul du nombre aléatoire  
  call random_number(alea_continue)
  alea_continue = mini + (maxi-mini) * alea_continue


end function alea_continue


!////////////////////////////////////////////////////////////////////////////


subroutine normale(x1, x2)

	implicit none
  
! Argument
  real(8), intent(out) :: x1, x2
  real(8)              :: u1, u2, rho
  

! Attribution à x1, x2
  u1 = alea_continue(0.0_8, 1.0_8)
  u2 = alea_continue(0.0_8, 1.0_8)

  rho = sqrt(-2*log(u1))
  x1  = rho*cos(2*Pi*u2)
  x2 =  rho*sin(2*Pi*u2)



end subroutine


!////////////////////////////////////////////////////////////////////////////

	
end program PLANETS