module julias

  implicit none
  contains

!===============================================================
!-
!      some problems with color mapping, need more work
!-
subroutine simple_julia(pic, cx, cy, maxiter)
    implicit none
    integer, intent(inout), dimension (:,:)  ::  pic
    real, intent(in)                         ::  cx, cy
    integer, intent(in)                      ::  maxiter

    integer      ::  ix, iy, width, height
    real         ::  fx, fy
    complex      ::  Z, C
    integer      ::  iter
    logical      ::  over_iter

    width  = ubound(pic, 1)
    height = ubound(pic, 2)
    C = complex(cx, cy)
    print *, "Const = ", C
    ! ready ? ok, clear the picture
    pic = 0
    do ix = 1, width
      fx = (float(ix) / (float(width)/4.0) - 2.0)
      do iy = 1, height
        fy = (float(iy) / (float(height)/4.0) - 2.0)
          ! ------ traitement du pixel
          iter = 0  ; over_iter = .FALSE.
          Z = complex(fx, fy)
          do while (real(Z)*real(Z) + imag(Z)*imag(Z) .LT. 4.0) 
            Z = (Z * Z) + C
            iter = iter + 1
            if (iter .GE. maxiter) then
              over_iter = .TRUE.
              exit
            endif
          end do
          if (over_iter) then
            pic(ix, iy) = 0
          else
            pic(ix, iy) = mod(iter*13, 256)
          endif
      enddo      ! iy
    enddo        ! ix

end subroutine simple_julia
!===============================================================

subroutine julia_colormapped(pic, cx, cy, maxiter)
    use pixrgb
    type(t_pixrgb), intent(inout), dimension (:,:)  ::  pic
    real, intent(in)                                ::  cx, cy
    integer, intent(in)                             ::  maxiter

    integer      ::  ix, iy, width, height
    real         ::  fx, fy
    complex      ::  Z, C
    integer      ::  iter
    logical      ::  over_iter

    width  = ubound(pic, 1)
    height = ubound(pic, 2)
    C = complex(cx, cy)
    ! print *, "Color julia, const = ", C
    do ix = 1, width
      fx = (float(ix) / (float(width*2)/10.0) - 2.5)
      do iy = 1, height
        fy = (float(iy) / (float(height*2)/10.0) - 2.5)
          ! ------ traitement du pixel
          iter = 0  ; over_iter = .FALSE.
          Z = complex(fx, fy)
          do while ((real(Z)*real(Z) + imag(Z)*imag(Z)) .LT. 4.0) 
            Z = (Z * Z) + C
            iter = iter + 1
            if (iter .GE. maxiter) then
              over_iter = .TRUE.
              exit
            endif
          end do
          if (over_iter) then
            pic(ix, iy)%r = 0
            pic(ix, iy)%g = mod(abs(int(real(Z) *140)), 255)
            pic(ix, iy)%b = mod(abs(int(aimag(Z)*140)), 255)
          else
            pic(ix, iy)%r = mod(iter*22, 255)
            pic(ix, iy)%g = mod(iter*59, 255)
            pic(ix, iy)%b = mod(iter*21, 255)
          endif
      enddo      ! iy
    enddo        ! ix

end subroutine 
!===============================================================
end module