! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2010
!
! ejemplo 2c
!
! numeros aleatorios - subrutina portable utilizando
! un generador lineal multiplicativo congruente 
! (Park and Miller, 1988)
!
! distribucion uniforme en [0,1)
!
!
program eje2c
implicit none
  integer :: i,j,k,old
  integer, dimension(10) :: arr
  real :: valr
  
! pongo a cero el array
!
  arr=0
  
! inicializo la subrutina que genera numeros aleatorios
!
  write(*,*)'Ingrese semilla (int): '
  read(*,*)j
  valr=-float(j)
  
! itero j veces
!
  write(*,*)'Ingrese nro. iteraciones: '
  read(*,*)j
  
  do i=1,j
  
! obtengo un numero aleatorio
!
    call random1(valr)

! lo guardo en un histograma
!
	k=int(valr*10.)+1
	arr(k)=arr(k)+1
  enddo
  
! encuentro el bin con el mayor valor
!
  old=-12
  k=0
  do i=1,10
    if(arr(i) > old)then
	  old=arr(i)
	  k=i
	endif
  enddo
  
! muestro el resultado
!
  open(7,file='eje2c.out')
  do i=1,10
    write(*,*)i,arr(i),float(arr(i))/float(arr(k))
    write(7,*)i,arr(i),float(arr(i))/float(arr(k))
  enddo
  close(7)
end program eje2c


! *******************************************************
! subrutina mejorada para generar numeros
! aleatorios con distribucion uniforme en [0,1)
! y portable (no usa random_number) con un generador
! lineal multiplicativo congruente.
!
! se utiliza el algoritmo de schrage para evitar overflow
!
! la subrutina devuelve en val el numero
! aleatorio. Si en la entrada val es negativo
! la subrutina se inicializa
!
subroutine random1(val)
implicit none
  integer, parameter :: ma=2147483647, aa=16807, qq=127773, rr=2836
  real, parameter :: fac=1./ma
  integer :: i,hi,lo
  integer, save :: idum
  integer, dimension(97), save :: dd
  real :: val
  
! inicializa
!
  if(val < 0.)then
    idum=-int(val)
	do i=1,97
	  hi=idum/qq
	  lo=mod(idum,qq)
	  idum=aa*lo-rr*hi
	  if(idum < 0)idum=idum+ma
	  dd(i)=idum
	enddo
  endif
  
! punto de entrada general
!
  i=1+(97*idum)/ma
  if(i > 97 .or. i< 0)pause
  idum=dd(i)
  val=float(idum)*fac
  hi=idum/qq
  lo=mod(idum,qq)
  idum=aa*lo-rr*hi
  if(idum < 0)idum=idum+ma
  dd(i)=idum
  return
end subroutine random1
! *******************************************************