!-----------------------------------------------------------------------------!
!This program iterates the logistic map for a specified number, npoint+nskip, !
!steps, for nmu values of the control parameter between mu1 and mu2 (nskip is !
!the number of skipped steps before graphing begins. A ps file 'lmap.ps' with !
!a graph showing the region mu=mu1,mu2, x=x1,x2 is produced.                  !
!-----------------------------------------------------------------------------!
!The y-axis (the values x_n for fixed mu) are discretized as lpix pixels that !
!are filled if a at least one point x_n falls within it.                      !
!-----------------------------------------------------------------------------!
!Input data is read from the file 'read.in', containing                       !
! mu1,mu2,x1,x2                                                               !
! nmu,npoint,nskip                                                            !
!-----------------------------------------------------------------------------!
 program logisticmap
 implicit none

 integer, parameter :: lpix=1000

 integer    :: i,j,nmu,pix
 integer(8) :: npoint,nskip
 logical    :: filled(0:lpix-1)
 real(8)    :: mu,mu1,mu2,x,x1,x2,dx 

 open(1,file='read.in',status='old')
 read(1,*)mu1,mu2,x1,x2
 read(1,*)nmu,npoint,nskip
 close(1)

 call initgraph(mu1,mu2,x1,x2,nmu,lpix)
 dx=dble(lpix)/(x2-x1)
 do j=0,nmu-1
   mu=mu1+dble(j)*(mu2-mu1)/dble(nmu)
   x=0.5d0
   do i=1,nskip
      x=mu*x*(1.d0-x)
   end do
   filled=.false.
   do i=1,npoint
      x=mu*x*(1.d0-x)
      if (x >= x1 .and. x < x2) then
        pix=int((x-x1)*dx)
        filled(pix)=.true.
      end if
   end do
   call plotmu(j,lpix,filled)
 end do
 call closegraph
 
 end program logisticmap
!-----------------------!

!-------------------------------------------------------!
!Plots column ix, using a 'filled' vector with ny pixels!
!-------------------------------------------------------!
 subroutine plotmu(ix,ny,filled)
 implicit none

 integer :: i,ix,ny
 logical :: filled(0:ny-1),last

 last=.false.
 do i=0,ny-1
    if (last.and.(.not.filled(i))) then
      write(1,*)ix,i,' l'
    else if ((.not.last).and.(filled(i))) then
      write(1,*)ix,i,' m'
    end if
    last=filled(i)
 end do
 if (filled(ny-1)) write(1,*)ix,ny,' l'
 write(1,*)'stroke'

 end subroutine plotmu
!---------------------!

!--------------------------------------------------------------------!
!Open the ps file, draws the coordinate system, and defines operators!
!--------------------------------------------------------------------!
 subroutine initgraph(x1,x2,y1,y2,nx,ny)
 implicit none

 integer :: nx,ny
 real(8) :: x1,x2,y1,y2

 open(1,file='lmap.ps',status='replace')
 1 format(a)
 write(1,1)'%!'
 write(1,1)'80 200 translate 4.5 4.5 scale 0.3 setlinewidth'
 write(1,1)'-2 0 moveto 102 0 rlineto stroke'
 write(1,1)'-2 100 moveto 102 0 rlineto stroke'
 write(1,1)'0 -2 moveto 0 102 rlineto stroke'
 write(1,1)'100 -2 moveto 0 102 rlineto stroke'
 write(1,1)'/TimesRoman findfont 5 scalefont setfont'
 write(1,2)'-5 -7 moveto (',x1,') show'
 write(1,2)'95 -7 moveto (',x2,') show'
 write(1,2)'-15 99 moveto (',y2,') show'
 write(1,2)'-15 -1 moveto (',y1,') show'
 2 format(a,f5.3,a)
 write(1,1)'36 102 moveto (Logistic Map) show'
 write(1,1)'-6 49 moveto (x) show'
 write(1,1)'/Symbol findfont 5 scalefont setfont'
 write(1,1)'48 -6 moveto (m) show'
 write(1,*)100./real(nx),100./real(ny),' scale'
 write(1,*)'1.02 setlinewidth 0.3 0 0 setrgbcolor'
 write(1,*)'0.5 0 translate'
 write(1,*)'/l {lineto} def'
 write(1,*)'/m {moveto} def'

 end subroutine initgraph
!------------------------!

!----------------------------------------------------!
!Adds the 'showpage' statement and closes the ps file!.
!----------------------------------------------------!
 subroutine closegraph

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

 end subroutine closegraph
!-------------------------!