!------------------------!
 program binomialthemovie
 implicit none

 integer :: i,n
 real(8) :: p1

 write(*,*)'P(1) [P(0)=1-P(1)]';read*,p1
 write(*,*)'Highest n';read*,n

 do i=1,n
   call initframe(i)
   call makeframe
   call histogram(i,p1)
   call closeframe
 end do

 end program binomialthemovie
!----------------------------!   

!--------------------------!
 subroutine histogram(n,p1)
!--------------------------!
 implicit none

 integer :: n,i,k
 real(8) :: p0,p1,fn,fk,fnk,fmax,f(0:n)

 p0=1.d0-p1
 fn=0.d0
 do i=2,n
   fn=fn+log(dble(i))
 end do
 fmax=0.d0
 do k=0,n
   fk=0.d0
   do i=2,k
     fk=fk+log(dble(i))
   end do
   fnk=0.d0
   do i=2,n-k
     fnk=fnk+log(dble(i))
   enddo
   f(k)=exp(k*log(p1)+(n-k)*log(p0)+fn-fk-fnk)
   if (f(k) > fmax) fmax=f(k)
 end do

 write(1,*)'gsave 0 0 0.3 setrgbcolor'
 write(1,*)'3 setlinewidth'
 do k=0,n
   if (n > 25) write(1,*)75./real(n),' setlinewidth'
   write(1,10)100.*real(k)/real(n),' 0 moveto 0 ',100.*real(f(k)/fmax),' rlineto stroke'
 enddo
 10 format(F8.3,a,F8.3,a)
 write(1,*)'grestore'
 write(1,*)'/Helvetica findfont 6 scalefont setfont'
 write(1,*)'4 102 moveto (N =',n,') show' 

 end subroutine histogram
!------------------------!

!-------------------!
 subroutine makeframe
!-------------------!
 implicit none

 write(1,*)'15 10 translate 0.72 0.72 scale'
 write(1,*)'0.4 setlinewidth'
 write(1,*)'-5 0 moveto 110 0 lineto stroke'
 write(1,*)'0 -1 moveto 0 110 lineto stroke'
 write(1,*)'100 -1 moveto 100 0 lineto stroke'
 write(1,*)'/Helvetica findfont 4 scalefont setfont'
 write(1,*)'-1 -5 moveto (0) show 99.5 -5 moveto (1) show'
 write(1,*)'46 -5 moveto (<x>) show'
 write(1,*)'-14 70 moveto (P(<x>)) show'

 end subroutine makeframe
!------------------------!

!-----------------------------!
 subroutine initframe(frameno)
!-----------------------------!
 implicit none

 integer :: i1,i2,i3,frameno
 character(16) :: fname

 i3=frameno/100
 i2=mod(frameno,100)/10
 i1=mod(frameno,10)
 fname='anim000.ps'
 fname(5:5)=achar(48+i3)
 fname(6:6)=achar(48+i2)
 fname(7:7)=achar(48+i1)
 open(1,file=fname,status='replace')
 write(1,*)'%!'
 write(1,*)'1 1 scale 0 0 translate 0.5 0.5 moveto'
 write(1,*)'612 0 rlineto 0 612 rlineto -612 0 rlineto 0 -612 rlineto'
 write(1,*)'closepath stroke'
 write(1,*)'6.12 6.12 scale'
   
 end subroutine initframe
!------------------------!

!---------------------!
 subroutine closeframe
!---------------------!
 implicit none

 write(1,*)'showpage'
 close(1)

 end subroutine closeframe
!-------------------------!