program newton
      implicit none
      complex c,z_ini,z,z2,z3_1
      integer pgopen,i,j,k,i_again,isymbol,n_gen,n_fifteenth
      integer idummy,pgband,indx,i_color
      real lmd_x_min,lmd_x_max,lmd_y_min,lmd_y_max,lmd_x,lmd_y,abs_z3_1
      real x_pos,y_pos,x1,x2,y1,y2,temp
      character*(1) ch

      isymbol = -1

      write (*,*) 'Plots Newton initial z such that z3-1 divserges.'
      write (*,*)
c     z_ini = (0.5,0.0)
c     z_ini = (0.0,0.0)

      if ( pgopen('/xwin') .le. 0 ) stop
      call pgpap(5.0,1.0)
c pgask(0) means not to ask confirmation when erase a page
      call pgask(0)


 100  write(*,*) 'Color (1) or balck-and-white (0) ?'
      read (*,*) i_color
 
      write (*,*) 'What is the numer of generation for the mapping ?'
      write (*,*) 'For best effect, the numer must not be too big.'
      read (*,*) n_gen

      call pgeras
      call pgsci(1)

      lmd_x_min = -1.0
      lmd_x_max =  1.0
      lmd_y_min = -1.0
      lmd_y_max =  1.0
 
c     lmd_x_min = -2.0
c     lmd_x_max =  0.5
c     lmd_y_min = -1.25
c     lmd_y_max =  1.25

 200  call pgenv(lmd_x_min, lmd_x_max, lmd_y_min, lmd_y_max, 1, 0)

      n_fifteenth = n_gen/15


      IF (i_color.eq.1) THEN

      do j=1,500
        lmd_x = lmd_x_min + j*(lmd_x_max-lmd_x_min)/500
        call pgbbuf
        do k=1,500
          lmd_y = lmd_y_min + k*(lmd_y_max-lmd_y_min)/500
          z_ini = (1.0,0.0)*lmd_x + (0.0,1.0)*lmd_y
          z = z_ini
          do i=1,n_gen
            z2 = z*z
            z = z - (z2*z-1)/(3*z2)
            z3_1 = z*z*z - 1
            abs_z3_1 = sqrt( conjg(z3_1)*z3_1 )
            if (abs_z3_1.lt.10e-5) then
              call pgsci(i/n_fifteenth)
              call pgpt(1,lmd_x,lmd_y,isymbol)
              exit
            endif
          end do
        end do
        call pgebuf
      end do
      call pgsci(1)
      call pgbox('BCNST',0.0,0,'BCNST',0.0,0)

      ELSE

      do j=1,1000
        lmd_x = lmd_x_min + j*(lmd_x_max-lmd_x_min)/1000
        do k=1,1000
          lmd_y = lmd_y_min + k*(lmd_y_max-lmd_y_min)/1000
          z_ini = (1.0,0.0)*lmd_x + (0.0,1.0)*lmd_y
          z = z_ini
c         call pgbbuf
          do i=1,n_gen
            z2 = z*z
            z = z - (z2*z-1)/(3*z2)
          end do
          z3_1 = z*z*z - 1
          abs_z3_1 = sqrt( conjg(z3_1)*z3_1 )
          if (abs_z3_1.gt.10e-8) then
            call pgpt(1,lmd_x,lmd_y,isymbol)
          endif
c         call pgebuf
        end do
      end do
      call pgsci(1)
      call pgbox('BCNST',0.0,0,'BCNST',0.0,0)

      ENDIF


      write (*,*) 'Enlarge some parts ? (1=YES;0=NO)'
      read (*,*) i_again
      if (i_again.eq.0) goto 102

      write (*,*) 'Now choose a range to enlarge using mouse pointer.'
      write (*,*)
      idummy = pgband (0, 1, 0.5, 0.5, x_pos, y_pos, CH)
      x1 = x_pos
      y1 = y_pos 
      idummy = pgband (2, 1, x1, y1, x_pos, y_pos, CH)
      x2 = x_pos
      y2 = y_pos

      if (x1.gt.x2) then
        temp = x1
        x1 = x2
        x2 = temp
      endif
      if (y1.gt.y2) then
        temp = y1
        y1 = y2
        y2 = temp
      endif

      write (*,*) 'Your choice is: ', ' r = (',x1,'~',x2,' y = (',y1,'~'
     &,y2,')'

      lmd_x_min = x1
      lmd_x_max = x2
      lmd_y_min = y1
      lmd_y_max = y2

      write (*,*) 'What is the numer of generation for the mapping ?'
      read (*,*) n_gen
      goto 200

 102  write (*,*) 'Again ? (1=YES;0=NO)'
      read (*,*) i_again
      if (i_again.eq.1) goto 100
      if (i_again.eq.0) stop
      goto 102

      end