Dist_Nciclo

! Cálculo da distântica para a rede circular
implicit complex(a-h,o-z)

real, dimension(2500,2500)      :: u, plu, pru
complex, dimension(2500,2500)   :: rho, rhoI

real, dimension(2500)           :: ProbM !Vector for calculation of the average probability

integer, parameter              :: T = 100
integer, parameter              :: N = 20
integer, parameter              :: dim_R = 2*N !Dimension of the network
real, parameter                 :: c = 1.0/sqrt(2.0)

integer, parameter              :: Ndiv = 20
real, parameter                 :: p_min = 0.000001
real, parameter                 :: p_max = 1.0
real, parameter                 :: Pu = 1/N !1/N

real                            :: prob_medio
real                            :: p, dist

open(1,file="Dados_Dist_N-cycle.dat")

print *, "Simulation Started ..."
print *, "Parameters: "
print *,"   >>> T = " , T
print *, "  >>> N = ", N
print *, "  >>> Nº division = ",Ndiv
print *, "  >>> Minimum rate (p) = ", p_min
print *, "  >>> Maximum rate (p) = ", p_max
print *, " "
do i = 2,N-1
    k = 2*i-1
    u(k-2,k) = c
    u(k+3,k) = c
    u(k-2,k+1) = c
    u(k+3,k+1) = -c
end do

!=> x=1: x-1 -> n (no zero)
u(2*N-1,1) = c
u(4,1) = c
u(2*N-1,2) = c
u(4,2) = -c

!=> x=n: x+1 -> 1
u(2*N-3,2*N-1) = c
u(2,2*N-1) = c
u(2*N-3,2*N) = c
u(2,2*N) = -c

do i = 2,dim_R,2
    do j = 1,dim_R
        plu(i-1,j) = u(i-1,j)
        pru(i,j) = u(i,j)
    end do
end do

print *, " >>> Steps on time and calculation of the probabilities ..."
p = p_min
do l = 1,Ndiv

    do k = 1,dim_R
        do j = 1,dim_R
            rhoI(k,j) = (0.0,0.0)
            rho(k,j) = (0.0,0.0)
        end do
    end do

    do i = 1,dim_R
        ProbM(i) = 0.0
    end do

    rhoI(1,1) = (0.5,0.0)
    rhoI(1,2) = (0.0,-0.5)
    rhoI(2,1) = (0.0,0.5)
    rhoI(2,2) = (0.5,0.0)

    do m = 1,T
        call caminhada(rho,rhoI,u,PLU,PRU,p,N)

        do i = 1,N
            k = 2*i-1
            pnl = rho(k,k)
            pnr = rho(k+1,k+1)
            prob = pnl + pnr

            ProbM(i) = ProbM(i) + prob
        end do
    end do
    dist = 0.0

    do i = 1,N
        prob_medio = ProbM(i)/T
        dist = dist + abs(prob_medio - Pu)
    end do

    write(1,*) p, dist
    print *, l, p, dist

    p = p + (p_max - p_min)/Ndiv
end do
close(1)
print *, " ... End of the Simulation"
end                                     !--[[ Fim do Algoritmo ]]--!


subroutine caminhada(rho,rhoI,u,plu,pru,p,N)

    implicit complex(a-h,o-z)

    real                            :: p
    real, dimension(2500,2500)      :: u, PLU, PRU
    complex,dimension(2500,2500)    :: rhoI, rho, aux1, aux2

    !Cálculo feito sem os sítios n = 1 e n = N
    Nmin = 3
    Nmax = 2*N-2

    do i = Nmin, Nmax
        do j = i, Nmax
            aux1(i,j) = (0.0,0.0)
            aux2(i,j) = (0.0,0.0)

            do k=i-3,i+3
                do l = j-3,j+3
                    aux1(i,j) = aux1(i,j) + u(i,k)*rhoI(k,l)*u(j,l)
                    aux2(i,j) = aux2(i,j) + PLU(i,k)*rhoI(k,l)*PLU(j,l) + PRU(i,k)*rhoI(k,l)*PRU(j,l)
                end do
            end do

            rho(i,j) = (1.0-p)*aux1(i,j) + p*aux2(i,j)
            rho(j,i) = conjg(rho(i,j))

        end do
    end do

    ! sítio n = 1
    do i = 1,2
        do j = 1,2*N
            aux1(i,j) = (0.0,0.0)
            aux2(i,j) = (0.0,0.0)

            ! A soma é feita em duas etapas: [3,4] e [2N-1,2N]

            ! Para n := [3,4]
            do k = 3,4
                do l = 1, 2*N
                    aux1(i,j) = aux1(i,j) + u(i,k)*rhoI(k,l)*u(j,l)
                    aux2(i,j) = aux2(i,j) + PLU(i,k)*rhoI(k,l)*PLU(j,l) + PRU(i,k)*rhoI(k,l)*PRU(j,l)
                end do
            end do

            ! Para n := [2N-1,2N]
            do k = 2*N-1,2*N
                do l = 1,2*N
                    aux1(i,j) = aux1(i,j) + u(i,k)*rhoI(k,l)*u(j,l)
                    aux2(i,j) = aux2(i,j) + PLU(i,k)*rhoI(k,l)*PLU(j,l) + PRU(i,k)*rhoI(k,l)*PRU(j,l)
                end do
            end do


            rho(i,j) = (1.0-p)*aux1(i,j) + p*aux2(i,j)
            rho(j,i) = conjg(rho(i,j))
        end do
    end do

    ! sítio n = N

    Nmin = 2*N-1
    Nmax = 2*N

    do i = Nmin, Nmax
        do j = 1, Nmax
            aux1(i,j) = (0.0,0.0)
            aux2(i,j) = (0.0,0.0)

            ! A soma é feita em duas etapas: [1,2] e [2N-3,2N-2]

            !Para n := [1,2]
            do k = 1,2
                do l = 1,2*N
                    aux1(i,j) = aux1(i,j) + u(i,k)*rhoI(k,l)*u(j,l)
                    aux2(i,j) = aux2(i,j) + PLU(i,k)*rhoI(k,l)*PLU(j,l) + PRU(i,k)*rhoI(k,l)*PRU(j,l)
                end do
            end do

            ! Para n := [2N-3,2N-2]
            do k = 2*N-3, 2*N-2
                do l = 1,2*N
                    aux1(i,j) = aux1(i,j) + u(i,k)*rhoI(k,l)*u(j,l)
                    aux2(i,j) = aux2(i,j) + PLU(i,k)*rhoI(k,l)*PLU(j,l) + PRU(i,k)*rhoI(k,l)*PRU(j,l)
                end do
            end do


            rho(i,j) = (1.0-p)*aux1(i,j) + p*aux2(i,j)
            rho(j,i) = conjg(rho(i,j))
        end do
    end do


    do i = 1,2*N
        do j = i, 2*N
            rhoI(i,j) = rho(i,j)
            rhoI(j,i) = rho(j,i)
        end do
    end do
end subroutine