clear all;
N=76;
Qmax=140;
Tmax=160;

Coordinates=xlsread('Coordi07.xlsx');
d=xlsread('Demand07.xlsx');
P=xlsread('Profit07.xlsx');
T=zeros(76,76);

for i=1:N
    for j=1:N
        if i~=j
            x=(Coordinates(i,1)-Coordinates(j,1))^2;
            y=(Coordinates(i,2)-Coordinates(j,2))^2;
            T(i,j)=sqrt(x+y);
            T(j,i)=T(i,j);
        else
            T(j,i)=inf;
        end
    end
end
PP=zeros(N,N);
QQ=zeros(N,N);
Apostasis=zeros(N,N);
diadromi=zeros(N,N);
kostos=zeros(N,1);
Zitisi=zeros(N,1);
count=1;
L=1;
while count~=N && L<=2
K=0;
Q=0;
diadromi(L,1) =1;
A=0;
k=1;es
B=1;
P_diadromwn=0;
while A~=1
    
    [~,n]=max(P(:,1));
    
    if Q+d(n)<=Qmax && K+T(B,n)+ T(n,1)<=Tmax
   
    k=k+1;
    diadromi(L,k)=n;
    K=K+T(B,n);
    Q=Q+d(n);
    Apostasis(L,k)=T(B,n);
    QQ(L,k)=d(n);
    PP(L,k)=P(n);
    B=n;
    count=count+1;
    P_diadromwn=P_diadromwn+P(n);
    P(B,:)=0;
  
    T(:,B)=inf;
    
             
    else
       [~,n]=min(T(B,:));
          if Q+ d(n)<=Qmax && K + T(B,n)+ T(n,1)<=Tmax
              k=k+1;
              diadromi(L,k)=n;
              K=K+T(B,n);
              Q=Q+d(n);
              Apostasis(L,k)=T(B,n);
              QQ(L,k)=d(n);
              PP(L,k)=P(n);
              B=n;
              count=count+1;
             T(B,diadromi(L,k-1))=inf;
              T(:,B)=inf;
              P_diadromwn=P_diadromwn+P(n);
              P(B,:)=0;
               
          else
          
            k=k+1;
            Apostasis(L,k)=T(B,1);
            QQ(L,k)=0;
            PP(L,k)=0;
            diadromi(L,k)=1;
            kostos(L)=K+T(B,1);
            zitisi(L)=Q;
            A=1;
            P_diad(L)=P_diadromwn+P(B)
            L=L+1;
          end
      
  end
end

end
Synoliko_Profit=sum(P_diad);
for i=1:N
    for j=1:N
        if i~=j
            x=(Coordinates(i,1)-Coordinates(j,1))^2;
            y=(Coordinates(i,2)-Coordinates(j,2))^2;
            T(i,j)=sqrt(x+y);
            T(j,i)=T(i,j);
        else
            T(j,i)=inf;
        end
    end
end