program mira

real(kind=8) :: a,b,Xmax,Ymax,eps,M
real(kind=8), allocatable, dimension(:) :: x,y,gridX,gridY
integer :: N,KX,KY,L,i,j,k,z


print *,'Donnez des valeurs aux coefficients a et b : '
read(*,*) a
read(*,*) b
print *,"Saisir le nombre d'iterations"
read(*,*) N

!print *,"Nombre de epsilons à considérer"
!read(*,*) L


ALLOCATE(x(N))
ALLOCATE(y(N))


print *,"Saisir les conditions initiales "
print *,"X0=? "
read(*,*) x(1)
print *,"YO=?"
read(*,*) y(1)

!print *,"Saisirilon"
!print *,"Epsilon=?"
!read(*,*) eps


open (unit = 1, file = "xy.txt")
write (1,*) '********************'
write (1,*) 'Conditions initiales'
write (1,*) 'a=',a
write (1,*) 'b=',b
write (1,*) 'N=',N
write (1,*) 'X0=',x(1)
write (1,*) 'Y0=',y(1)
write (1,*) '********************'
write (1,*) ' '
write (1,*) 'x                              Y'

do i=1,N

x(i+1) = b*y(i)+a*x(i)+((2*x(i)*x(i)*(1-a))/(1+(x(i)*x(i))))
y(i+1) = -x(i) + a*x(i+1) + ((2*x(i+1)*x(i+1)*(1-a))/(1+(x(i+1)*x(i+1))))
open (unit = 1, file = "xy.txt")
write (1,*) x(i),y(i)
!open (unit = 2, file = "x.txt")
!write (2,*) x(i)
!open (unit = 3, file = "y.txt")
!write (3,*) y(i)
end do



L=10
eps=0.1


do z=1,10!L
Xmax=int(2*MAXVAL(abs(x))) !max en X
Ymax=int(2*MAXVAL(abs(y))) !max en Y


KX=0
KY=0
eps=eps+0.01


print *, 'Epsilon=',eps

M=0

KX=int(2*Xmax/eps) ! NB de CARRES COTE X (taille eps)
KY=int(2*Ymax/eps) ! NB de CARRES COTE Y (taille eps)

ALLOCATE(gridX(KX))
ALLOCATE(gridY(KX))
gridX(1)=-Xmax
gridY(1)=-Ymax

do i=1,KX !compute gridX
gridX(i+1)=gridX(i)+eps
end do

do i=1,KY !compute gridY
gridY(i+1)=gridY(i)+eps
end do

!do i=1,KX !print gridY.txt
!open (unit = 4, file = "gridx.txt")
!write (4,*) gridX(i)
!end do

!do i=1,KY !print gridY.txt
!open (unit = 5, file = "gridy.txt")
!write (5,*) gridY(i)
!end do

!à ce stade, nous avons une figure fractale décrite par N couples (x,y), ainsi qu'un cadrillage autour de la fractale.
!maintenant, l'idée est de chercher à savoir combien de carrés de taille eps sont nécessaire pour recouvrir la fractale,
!pour ce faire, on fait tendre eps vers 0, et on compte le nombre de carrés qui contiennent au moins un des couples (x,y).
!Si aucun point ne passe dans le carré, ce dernier ne compte pas.
!On va donc parcourir tous le cadrillage, carré par carré, et tester chaque couple (x,y).

do i=1,KX
    do j=1,KY
        do k=1,N
            if (x(k)>gridX(i) .AND. x(k)<gridX(i+1) .AND. y(k)>gridY(j) .AND. y(k)<gridY(j+1)) then
                M=M+1


            exit
            end if

        end do
    end do
end do

print *, 'M=',M



DEALLOCATE(gridX)
DEALLOCATE(gridY)
end do
DEALLOCATE(x)
DEALLOCATE(y)
end program mira

!///////////////////////////////////////////////////////////////////////!