program fern
      real x,y,newx,newy,rand_num,ran2
      integer maxit,i,iseed,i_switch,pgopen,i_symbol
      external ran2

      if (pgopen('/xwin').le.0) stop
      call pgpap(6.0,1.0)
      call pgenv(-5.5,5.5,-0.5,10.5,1,-1)
      call pgsci(1)
      call pgscr(16,0.1,0.7,0.05)
      i_symbol = -1

 100  write(*,*) 'What is the total number of iteration ?'
      read (*,*) maxit
      write(*,*) 'Give a random seed integer :'
      read (*,*) iseed

      call pgsci(16)

      DO 1000, II=1,MAXIT

      x = ran2(iseed)
      y = ran2(iseed)

      call pgbbuf

c     do i=1,maxit
      do i=1,20
        rand_num = ran2(iseed)
        if (rand_num.lt.0.01) then
          newx = 0.0
          newy = 0.16*y
          x = newx
          y = newy
        else
          if (rand_num.lt.0.86) then
            newx = 0.85*x + 0.04*y
            newy = -0.04*x + 0.85*y + 1.6
            x = newx
            y = newy
          else
            if (rand_num.lt.0.93) then
              newx = 0.2*x - 0.26*y
              newy = 0.23*x + 0.22*y + 1.6
              x = newx
              y = newy
            else
              newx = -0.15*x + 0.28*y
              newy = 0.26*x + 0.24*y + 0.44
              x = newx
              y = newy
            endif
          endif
        endif
      enddo
      call pgpt(1,x,y,i_symbol)

      call pgebuf

 1000 CONTINUE


 104  write(*,*) 'More points to draw ? (1=YES;0=NO)'
      read (*,*) i_switch
      if (i_switch.eq.0) goto 102
      if (i_switch.eq.1) goto 100
      goto 104
 102  continue

      call pgclos

      end