module fraktals
  implicit none
  contains

!-----------------------------------------------------
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

    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 (modulus2(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) = iter
          endif

      enddo
    enddo 

end subroutine simple_julia
!-----------------------------------------------------
!
!
!
subroutine pickover_0(pic, count)
    implicit none
    integer, intent(inout), dimension (:,:)  ::  pic
    integer, intent(in)                      ::  count

    double precision      :: xa, ya, za, xb, yb, zb
    double precision      :: ka, kb, kc, kd
    integer               :: i, px, py

    ka = 2.24 ; kb = 0.43 ; kc = -0.65 ; kd = -2.43

    xa = 0.00 ; ya = 0.00 ; za = 0.0


    do i=1, count

      xb = sin(ka*ya) - za*cos(kb*xa)
      yb = za*sin(kc*xa) - cos(kd*ya)
      zb = sin(xa)

      print *, i, xb, yb, zb

      xa = xb ; ya = yb ; za = zb 

    enddo

end subroutine pickover_0
!-----------------------------------------------------
!           -- some support functions --
!-----------------------------------------------------

function dist0 (x, y)
  implicit none
  real, intent(in)    :: x, y
  real                :: dist0
  dist0 = ( x*x + y*y )
end function

!-----------------------------------------------------
function modulus2(pt)
  implicit none
  complex, intent(in)    ::   pt
  real                   ::   modulus2
  modulus2 = real(pt)*real(pt) + imag(pt)*imag(pt)
end
!-----------------------------------------------------

end module fraktals