!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Calcolo numerico del volume di un' ipersfera           !
!d-dimensionale tramite il metodo "hit and miss"        !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Program hypersphere

  Implicit none
  real,parameter :: pi = acos(-1.0)
  real, dimension(:), allocatable :: rnd
  real :: truevol, vol
  Integer :: i, Nmax, accept, dim,n
  character(len = 8) :: file_name
!  INTEGER :: sizer
!  integer, dimension(:), allocatable :: seed


!  call random_seed(sizer)
!  allocate(seed(sizer))
!  print *,'Qui il seme ha ',sizer,' componenti; inserirle (o scrivere "/") >'  
!  read(*,*)seed
!  call random_seed(put=seed)

  print*,'Inserire la dimensionalità della sfera in numero (<10):'
  read*, dim
  

  allocate(rnd(dim))
  print*,'Numero massimo di punti da generare:'
  read*, Nmax
  !volume analitico dell'ipersfera
  truevol = pi**(0.5*dim)/gamma((0.5*dim) + 1)

  write(unit=file_name,fmt="(i1,a)") dim,"dim.dat" 
   print*," scrivo sul file ", file_name," i risultati ogni 100 punti"
   open (unit=7,file=file_name,action="write",Status='Replace') !NON USARE unita' 5 o 6


   write(7,*)'# n punti, vol. numerico, vol. esatto, abs(error), log(n), log(abs(error))'
  do n=1,Nmax
  accept=0
  if (mod(n,100) == 0) then
    !metodo hit and miss
    Do i=1, n
        call random_number(rnd)
        rnd = (rnd - 0.5)*2.0   ! genero numero random tra -1 e 1
        If (dot_product(rnd,rnd) <= 1) accept = accept + 1
    end do

    !calcolo del volume numerico
    vol = (2.0**dim)*float(accept)/float(n)
    
    !salvataggio dati
    write(7,*) n,vol,truevol, abs(truevol - vol),log(float(n)),log(abs(truevol-vol))
 end if
  end do
  Close(7)
  print*,' volume sfera numerico e esatto:',vol, truevol

End program hypersphere