! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2010
!
! ejemplo 2e
!
! numeros aleatorios - reproduce una distribucion arbitraria
! que se ingresa como un histograma
!
!
program eje2e
implicit none
  character (30) :: arch
  integer :: i,j,tot
  integer, dimension (2,30) :: aaa
  real, dimension (30) :: rrr
  real :: valr
  
! lee archivo de entrada
!
  write(*,*)'Archivo de entrada: '
  read(*,*)arch
  
  open(7,file=arch)
  tot=0
  aaa=0
  rrr=0.
  do i=1,30
  
! en el array aaa guarda f(a)
!
    read(7,*)j,aaa(1,i)
	
! en el array rrr guarda F(y)
!
	tot=tot+aaa(1,i)
	rrr(i)=float(tot)
  enddo
  close(7)
  
! renormaliza la funcion acumulativa
!
  do i=1,30
    rrr(i)=rrr(i)/float(tot)
  enddo
  
! inicializo la subrutina que genera numeros aleatorios
!
  write(*,*)'Ingrese semilla (int): '
  read(*,*)j
  valr=-float(j)
  
! itera 5e6 veces
!
  do j=1,5000000
  
! obtengo un numero aleatorio
!
    call random1(valr)
	i=1
	
! busco el bin del histograma que totaliza un
! valor de F(y) que sea mayor que el numero
! aleatorio encontrado
!
	do while(i < 30 .and. rrr(i) < valr)
	  i=i+1
    enddo	
	
! en el bin correspondiente de la nueva distribucion
! incremento el totalizador en 1
!
	aaa(2,i)=aaa(2,i)+1
  enddo
  
! renormaliza la distribucion para compararla
! con la original
!
  do i=1,30
    aaa(2,i)=int(float(aaa(2,i))*float(tot)/5000000.)
  enddo
  
! guarda resultados en archivo
!
  open(7,file='eje2e.out')
  do i=1,30
    write(7,*)i,aaa(1,i),aaa(2,i)
  enddo 
  close(7)
  
end program eje2e


! *******************************************************
! subrutina mejorada para generar numeros
! aleatorios con distribucion uniforme en [0,1)
! y portable (no usa random_number) con un generador
! lineal multiplicativo congruente.
!
! 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
! *******************************************************
