MODULE MODEL
    use READ_DATA
    IMPLICIT NONE
    INTEGER,DIMENSION(:),ALLOCATABLE:: visited_node
    CONTAINS
    FUNCTION  DISTANCE(i,j) !UPDATED
        INTEGER i,j
        INTEGER a,b
        REAL*8 distance,pi
        pi=dacos(-1d0)
        distance = 0d0
        IF(distance_definition.eq.1)THEN
                ! Arc distance with PBC
                 a=abs(j-i) ; b=n_spins-abs((j-i))
                 distance = dfloat(min(a,b))*distance_factor
        ELSE IF (distance_definition.eq.2)THEN
                ! Chord distance with PBC
                 distance =dfloat(n_spins)*sin(pi*dfloat(abs(i-j))/dfloat(n_spins))/pi
        END IF
    END FUNCTION
    SUBROUTINE COUPLING_MATRIX() !UPDATED
        IMPLICIT NONE
        INTEGER :: i,j,l,k,number_of_interactions,visit,n_iter
        REAL*8 ::probability_number,A_min,A_max,A_mid,Z_min,Z_max,Z_mid,norm_ct
        REAL :: random1,random2,pi,r1279
        LOGICAL :: repeated,accepted_j,connected
        pi=acos(-1.)
        ALLOCATE(visited_node(n_spins))
        n_iter=0
        connected=.FALSE.
        IF(is_initial_seed.eqv..TRUE.)THEN
                OPEN(41,FILE='interact_metrix.dat')
        END IF
        IF (matrix_definition.eq.1)THEN
                DO WHILE ((connected.eqv..FALSE.).and.(n_iter.lt.1000))
                n_iter=n_iter+1
                !+++++++++++++++++++++++++++++++++++++++++++++++
                !  NORMALIZATION CONSTANT
                !+++++++++++++++++++++++++++++++++++++++++++++++
                norm_ct=0d0
                DO j=2,n_spins
                        norm_ct=norm_ct+distance(1,j)**(-power_interactions)
                END DO
 
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                !  COUPLING MATRIX
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++

                number_neighbours=0
                eff_neigh=0
                spin_interaction=0.

                number_of_interactions=int(real(coordination_number*n_spins)/real(2))
                DO l=1,number_of_interactions
                        ! Choose random i
                        random1=r1279()
                        i=mod(int(real(n_spins)*random1),n_spins)+1
                        accepted_j=.FALSE.
                        DO WHILE (accepted_j.eqv..FALSE.)
                                !Choose random j
                                random2=r1279()            
                                j=mod(int(real(n_spins)*random2),n_spins)+1
                                !Chech that is different form i and the
                                !interaction in not assigned
                                repeated=.FALSE.
                                DO k=1,number_neighbours(i)!If this interaction was assinged before, go to next pair
                                         IF(eff_neigh(i,k).eq.j)repeated=.TRUE.
                                END DO
                                IF(i.ne.j)THEN
                                        probability_number=(distance(i,j)**(-power_interactions))/norm_ct
                                END IF
                                random1=r1279()
                                IF((real(random1,8).lt.probability_number).and.(i.ne.j).and.(repeated.eqv..FALSE.))THEN
                                        !Valid interaction
                                        accepted_j=.TRUE.
                                        number_neighbours(i)=number_neighbours(i)+1
                                        number_neighbours(j)=number_neighbours(j)+1
                                        eff_neigh(i,number_neighbours(i))=j
                                        eff_neigh(j,number_neighbours(j))=i
                                        IF (interaction_definition.eq.1) THEN
                                                ! +-1 interaction
                                                ! distribution
                                                random1=r1279()
                                                IF (random1.lt.0.5) THEN
                                                        spin_interaction(i,number_neighbours(i))=1.
                                                        spin_interaction(j,number_neighbours(j))=1.
                                                ELSE
                                                        spin_interaction(i,number_neighbours(i))=-1.
                                                        spin_interaction(j,number_neighbours(j))=-1.
                                                END IF
                                        ELSE IF(interaction_definition.eq.2)THEN
                                                ! Gausian interaction distribution
                                                random1=r1279()
                                                random2=r1279()
                                                spin_interaction(i,number_neighbours(i))=sqrt(-2.*log(random1))*cos(2.*pi*random2)
                                                spin_interaction(j,number_neighbours(j))=spin_interaction(i,number_neighbours(i)) 
                                        ENDIF
                                        IF(is_initial_seed.eqv..TRUE.)THEN
                                                WRITE(41,*)i,j,spin_interaction(i,number_neighbours(i))
                                        END IF
                                END IF
                        END DO
                END DO
                visit=1
                visited_node=0
                CALL CHECK_CONNECTED_GRAPH(visit)
                !print*,'visited nodes',sum(visited_node)
                IF(sum(visited_node).eq.n_spins)THEN
                        connected=.TRUE.
                        print*,'connected graph found with',n_iter,'iterations'
                END IF
                END DO
                OPEN(40, FILE='seeds_used.dat', action='write',status='unknown',position='append')
                WRITE(40,*)read_seed,real(sum(number_neighbours))/real(n_spins),power_interactions,&
                &n_spins
                close(40)
        ELSE IF (matrix_definition.eq.2)THEN
                DO WHILE ((connected.eqv..FALSE.).and.(n_iter.lt.1000))
                n_iter=n_iter+1
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                !   COMPUTING LOWER BOUND FOR A
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                probability_number=0d0
                DO j=2,n_spins
                        probability_number=probability_number+(1d0/(distance(1,j)**power_interactions))
                END DO
                A_min=dfloat(coordination_number)/probability_number
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                !   COMPUTING Z
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                probability_number=0d0
                DO j=2,n_spins
                        probability_number=probability_number+(1d0-exp(-A_min/(distance(1,j)**power_interactions)))
                END DO
                Z_min=probability_number
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                !   COMPUTING A FOR MAXIMUM Z
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                A_max=A_min
                DO i=1,500
                        probability_number=0d0
                        DO j=2,n_spins
                                probability_number=probability_number+(1d0-exp(-A_max/(distance(1,j)**power_interactions)))
                        END DO
                        IF (probability_number.gt.dfloat(coordination_number)) THEN
                                Z_max=probability_number
                                exit
                        END IF
                        A_max=A_max+1d-1*A_max
                END DO
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                !   FINDING THE BEST A
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                DO i=1,500
                        A_mid=(A_max+A_min)/2d0

                        probability_number=0d0
                        DO j=2,n_spins
                                probability_number=probability_number+(1d0-exp(-A_mid/(distance(1,j)**power_interactions)))
                        END DO
                        Z_mid=probability_number

                        IF (Z_mid.gt.dfloat(coordination_number))THEN
                                A_max=A_mid
                                Z_max=Z_mid

                        ELSE IF (Z_mid.lt.dfloat(coordination_number))THEN
                                A_min=A_mid
                                Z_min=Z_mid
                        END IF
                        IF(abs(Z_max-Z_min).le.1d-5)THEN
                                exit
                        END IF
                END DO
                probability_number=0d0
                DO j=2,n_spins
                        probability_number=probability_number+(1d0-exp(-A_mid/(distance(1,j)**power_interactions)))
                END DO

                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                !  COUPLING MATRIX
                !+++++++++++++++++++++++++++++++++++++++++++++++++++++++
                number_neighbours=0
                eff_neigh=0
                spin_interaction=0.

                DO i=1,n_spins
                        DO j=i+1,n_spins
                                repeated=.FALSE.
                                DO k=1,number_neighbours(i)!If this interaction was assinged before, go to next pair
                                        IF(eff_neigh(i,k).eq.j) repeated=.TRUE.
                                END DO
                                IF(i.ne.j)THEN
                                        probability_number=1d0-exp(-A_mid/(distance(i,j)**power_interactions))
                                END IF
                                random1=r1279()
                                IF((real(random1,8).lt.probability_number).and.(i.ne.j).and.(repeated.eqv..FALSE.))THEN
                                        !Valid interaction
                                        number_neighbours(i)=number_neighbours(i)+1
                                        number_neighbours(j)=number_neighbours(j)+1
                                        eff_neigh(i,number_neighbours(i))=j
                                        eff_neigh(j,number_neighbours(j))=i
                                        IF (interaction_definition.eq.1) THEN
                                                ! +-1 interaction
                                                ! distribution
                                                random1=r1279()
                                                IF (random1.lt.0.5) THEN
                                                        spin_interaction(i,number_neighbours(i))=1.
                                                        spin_interaction(j,number_neighbours(j))=1.
                                                ELSE
                                                        spin_interaction(i,number_neighbours(i))=-1.
                                                        spin_interaction(j,number_neighbours(j))=-1.
                                                END IF
                                        ELSE IF(interaction_definition.eq.2)THEN
                                                ! Gausian interaction
                                                ! distribution
                                                random1=r1279()
                                                random2=r1279()
                                                spin_interaction(i,number_neighbours(i))=sqrt(-2.*log(random1))*cos(2.*pi*random2)
                                                spin_interaction(j,number_neighbours(j))=spin_interaction(i,number_neighbours(i))
 
                                        ENDIF
                                        IF(is_initial_seed.eqv..TRUE.)THEN
                                                WRITE(41,*)i,j,spin_interaction(i,number_neighbours(i))
                                        END IF
                                END IF
                        END DO
                END DO
                visit=1
                visited_node=0
                CALL CHECK_CONNECTED_GRAPH(visit)
                !stop
                !print*,'visited nodes',sum(visited_node)
                IF(sum(visited_node).eq.n_spins)THEN
                        connected=.TRUE.
                        print*,'connected graph with',n_iter,'iterations'
                END IF
                END DO
                OPEN(40, FILE='seeds_used.dat', action='write',status='unknown',position='append')
                WRITE(40,*)read_seed,A_mid,real(sum(number_neighbours))/real(n_spins),power_interactions,&
                &n_spins
                close(40)
        ENDIF
        IF(is_initial_seed.eqv..TRUE.)THEN
                WRITE(41,*)'#####################'
                close(41)
        END IF
    END SUBROUTINE
    RECURSIVE SUBROUTINE CHECK_CONNECTED_GRAPH(to_visit)
        IMPLICIT NONE
        INTEGER :: to_visit,i,next_visit
        visited_node(to_visit)=1
        DO i=1,number_neighbours(to_visit)
                !print*,to_visit,eff_neigh(to_visit,i),visited_node(eff_neigh(to_visit,i))
                IF (visited_node(eff_neigh(to_visit,i)).eq.0)then
                        next_visit=eff_neigh(to_visit,i)
                        CALL CHECK_CONNECTED_GRAPH(next_visit)
                END IF
        END DO
     END SUBROUTINE
END MODULE MODEL
