lc=298; %comprimento da aresta do cubo metalico
l=lc+2; %comprimento da aresta do cubo coberto pelo polimero
t=1; %espessura das camadas do modelo polimerico
k=20; %número de malhas de metal
n=lc/t; %número de camadas
aux=(lc-1)/(k-1); %espaçamento entre cada malha
espa=475;

%coordenadas de onde vem o laser
cl(1,1)=(l+1)/2;
cl(1,2)=(l+1)/2;
cl(1,3)=1250;

for q=1:4
    for p=1:6
        cl(p+1+(q-1)*20,1)=1-espa;
        cl(p+1+(q-1)*20,2)=1-espa+(p-1)*250*1249/1250;
        cl(p+1+(q-1)*20,3)=1250-250*q;
    end

    for p=1:6
        cl(p+1+6+(q-1)*20,1)=l+espa;
        cl(p+1+6+(q-1)*20,2)=1-espa+(p-1)*250*1249/1250;
        cl(p+1+6+(q-1)*20,3)=1250-250*q;
    end

    for p=2:5
        cl(p+12+(q-1)*20,1)=1-espa+(p-1)*250*1249/1250;
        cl(p+12+(q-1)*20,2)=1-espa;
        cl(p+12+(q-1)*20,3)=1250-250*q;
    end

    for p=2:5
        cl(p+16+(q-1)*20,1)=1-espa+(p-1)*250*1249/1250;
        cl(p+16+(q-1)*20,2)=l+espa;
        cl(p+16+(q-1)*20,3)=1250-250*q;
    end
end

for p=1:20*4+1
    aux2=lc-1/2; %altura do cubo em relação à camada que esta a ser feita
    for m= 1:n-1/t %calculo dos pontos que o laser atinge para cada camada
        cont=0;
        %Todos os segmentos da rede metalica
        if (aux2<cl(p,3)-1)
            for i= 1:k
                a(i,1)=2;
                a(i,2)=2+(i-1)*aux;
                a(i,3)=aux2;
                a(i,4)=lc+1;
                a(i,5)=2+(i-1)*aux;
                a(i,6)=aux2;
            end

            for i= 1:k
                a(i+k,1)=2+(i-1)*aux;
                a(i+k,2)=2;
                a(i+k,3)=aux2;
                a(i+k,4)=2+(i-1)*aux;
                a(i+k,5)=lc+1;
                a(i+k,6)=aux2;
            end
            cont=2*k;
        end

         for i= 1:k
            a(i+2*cont,1)=2;
            a(i+2*cont,2)=2+(i-1)*aux;
            a(i+2*cont,3)=min(aux2,cl(p,3)-1);
            a(i+2*cont,4)=2;
            a(i+2*cont,5)=2+(i-1)*aux;
            a(i+2*cont,6)=0;
        end

        for i= 1:k
            a(i+k+2*cont,1)=lc+1;
            a(i+k+2*cont,2)=2+(i-1)*aux;
            a(i+k+2*cont,3)=min(aux2,cl(p,3)-1);
            a(i+k+2*cont,4)=lc+1;
            a(i+k+2*cont,5)=2+(i-1)*aux;
            a(i+k+2*cont,6)=0;
        end

        for i= 1:k
            a(i+2*k+2*cont,1)=2+(i-1)*aux;
            a(i+2*k+2*cont,2)=2;
            a(i+2*k+2*cont,3)=min(aux2,cl(p,3)-1);
            a(i+2*k+2*cont,4)=2+(i-1)*aux;
            a(i+2*k+2*cont,5)=2;
            a(i+2*k+2*cont,6)=0;
        end

        for i= 1:k
            a(i+3*k+2*cont,1)=2+(i-1)*aux;
            a(i+3*k+2*cont,2)=lc+1;
            a(i+3*k+2*cont,3)=min(aux2,cl(p,3)-1);
            a(i+3*k+2*cont,4)=2+(i-1)*aux;
            a(i+3*k+2*cont,5)=lc+1;
            a(i+3*k+2*cont,6)=0;
        end

        cont=cont+4*k; %cont guarda o numero de segmentos que a rede possui

        aux3=aux2-aux; %serve para ver quais são os segmentos laterais que estão a cima da camada que esta a ser considerada

        while (aux3>=1)
            if(aux3<cl(p,3)-1)
                a(cont+1,1)=2;
                a(cont+1,2)=2;
                a(cont+1,3)=aux3;
                a(cont+1,4)=2;
                a(cont+1,5)=lc+1;
                a(cont+1,6)=aux3; 

                a(cont+2,1)=lc+1;
                a(cont+2,2)=2;
                a(cont+2,3)=aux3;
                a(cont+2,4)=lc+1;
                a(cont+2,5)=lc+1;
                a(cont+2,6)=aux3;

                a(cont+3,1)=2;
                a(cont+3,2)=2;
                a(cont+3,3)=aux3;
                a(cont+3,4)=lc+1;
                a(cont+3,5)=2;
                a(cont+3,6)=aux3; 

                a(cont+4,1)=2;
                a(cont+4,2)=lc+1;
                a(cont+4,3)=aux3;
                a(cont+4,4)=lc+1;
                a(cont+4,5)=lc+1;
                a(cont+4,6)=aux3;

                cont=cont+4;
            end
            aux3=aux3-aux;
        end

        for i=1:l
            for j=1:l
                b(i,j,m)=1;%Começa por supor que todos os pontos da camada são atingidos pelo lazer
            end
        end
        for i=1:cont %para cada segmento determinado
            %calcula a projecão dos pontos extermos do segmento na camada
            aux4=cl(p,3)/(cl(p,3)-a(i,3));

            x1=cl(p,1)+aux4*(a(i,1)-cl(p,1));
            y1=cl(p,2)+aux4*(a(i,2)-cl(p,2));

            aux5=cl(p,3)/(cl(p,3)-a(i,6));

            x2=cl(p,1)+aux5*(a(i,4)-cl(p,1));
            y2=cl(p,2)+aux5*(a(i,5)-cl(p,2));
            %calcula quais os pontos que se encontram entre os pontos
            %calculados antes
            if(abs(x1-x2)>abs(y1-y2))
                if (x1-x2~=0)%Verica se a projecção nao é um unico ponto
                    if(x1>x2)%ordena os pontos
                        aux6=x1;
                        x1=x2;
                        x2=aux6;
                        aux6=y1;
                        y1=y2;
                        y2=aux6;
                    end
                    %determina o declive
                    d=(y2-y1)/(x2-x1);

                    %calculo dos pontos que não são atingidos pelo lazer
                    y1=y1+d*(round(x1)-x1);
                    y2=y2+d*(round(x2)-x2);
                    x1=round(x1);
                    x2=round(x2);
                    y=y1;
                    for j=x1:x2
                        %verifica se o ponto se encontra na zona pretendida
                        if(1<=j&&j<=l)
                            if(1<=round(y)&&round(y)<=l)
                                b(j,round(y),m)=0;%coloca o seu valor a 0
                            end
                        end
                        y=y+d;
                    end
                else
                    %verifica se o ponto se encontra na zona pretendida
                    if(1<=round(x1)&&round(x1)<=l)
                        if(1<=round(y1)&&round(y1)<=l)
                                b(round(x1),round(y1),m)=0;%coloca o seu valor a 0
                        end
                    end
                end
            else
                if (y1-y2~=0)%Verica se a projecção nao é um unico ponto
                    if(y1>y2)%ordena os pontos
                        aux6=x1;
                        x1=x2;
                        x2=aux6;
                        aux6=y1;
                        y1=y2;
                        y2=aux6;
                    end
                    %calculo do declive
                    d=(x2-x1)/(y2-y1);
                    %calculo dos pontos que não são atingidos pelo lazer
                    x1=x1+d*(round(y1)-y1);
                    x2=x2+d*(round(y2)-y2);
                    y1=round(y1);
                    y2=round(y2);
                    x=x1;

                    for j=y1:y2
                        %verifica se o ponto se encontra na zona pretendida
                        if(1<=round(x)&&round(x)<=l)
                            if(1<=j&&j<=l)
                                b(round(x),j,m)=0;%coloca o seu valor a 0
                            end
                        end
                        x=x+d;
                    end
                else
                        %verifica se o ponto se encontra na zona pretendida
                    if(1<=round(x1)&&round(x1)<=l)
                        if(1<=round(y1)&&round(y1)<=l)
                                b(round(x1),round(y1),m)=0;%coloca o seu valor a 0
                        end
                    end
                end
            end


        end
        aux2=aux2-t;%actualiza a altura da camada seguinte
    end


    for m=1:n-2/t
        for i=1:l
            final(i+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(i,1,m);
        end

        for i=1:l
            final(i+l+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(i,l,m);
        end

        for j=2:l-1
            final(j-1+2*l+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(1,j,m);
        end

        for j=2:l-1
            final(j-1+2*l+l-2+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(l,j,m);
        end
        
        j=0;
        for ii=1:k-1
            for i=3+round((ii-1)*aux):1+round(ii*aux);
                j = j + 1;
                final(j+4*(l-1)+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(i,2,m);
            end
        end
        
        j=0;
        for ii=1:k-1
            for i=3+round((ii-1)*aux):1+round(ii*aux);
                j = j + 1;
                final(j+l-2-k+4*(l-1)+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(i,l-1,m);
            end
        end
        
        j=0;
        for ii=1:k-1
            for i=3+round((ii-1)*aux):1+round(ii*aux);
                j = j + 1;
                final(j+2*(l-2-k)+4*(l-1)+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(2,i,m);
            end
        end

        j=0;
        for ii=1:k-1
            for i=3+round((ii-1)*aux):1+round(ii*aux);
                j = j + 1;
                final(j+3*(l-2-k)+4*(l-1)+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(l-1,i,m);
            end
        end
      
        
        for i=3:l-2
            final(i-2+4*(l-2-k)+4*(l-1)+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(i,3,m);
        end

        for i=3:l-2
            final(i-2+4*(l-2-k)+4*(l-1)+l-4+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(i,l-2,m);
        end

        for j=4:l-3
            final(j-3+4*(l-2-k)+6*(l-2)+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(3,j,m);
        end

        for j=4:l-3
            final(j-3+4*(l-2-k)+6*(l-2)+l-6+(m-1)*(8*(l-3)+4*(l-2-k)),p)=b(l-2,j,m);
        end
    end

    
    
    m=n-1/t;

    for i=1:l
        final(i+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(i,1,m);
    end

    for i=1:l
        final(i+l+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(i,l,m);
    end

    for j=2:l-1
        final(j-1+2*l+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(1,j,m);
    end

    for j=2:l-1
        final(j-1+2*l+l-2+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(l,j,m);
    end


            j=0;
    for ii=1:k-1
        for i=3+round((ii-1)*aux):1+round(ii*aux);
            j = j + 1;
            final(j+4*(l-1)+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(i,2,m);
        end
    end

    j=0;
    for ii=1:k-1
        for i=3+round((ii-1)*aux):1+round(ii*aux);
            j = j + 1;
            final(j+l-2-k+4*(l-1)+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(i,l-1,m);
        end
    end

    j=0;
    for ii=1:k-1
        for i=3+round((ii-1)*aux):1+round(ii*aux);
            j = j + 1;
            final(j+2*(l-2-k)+4*(l-1)+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(2,i,m);
        end
    end

    j=0;
    for ii=1:k-1
        for i=3+round((ii-1)*aux):1+round(ii*aux);
            j = j + 1;
            final(j+3*(l-2-k)+4*(l-1)+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(l-1,i,m);
        end
    end



    for i=3:l-2
        for j=3:l-2
            final((i-3)*(l-4)+(j-2)+4*(l-2-k)+4*(l-1)+(n-2/t)*(8*(l-3)+4*(l-2-k)),p)=b(i,j,m);
        end
    end
end



fileID = fopen('prob.lp','w');
fprintf(fileID,'min\n');
fprintf(fileID,'obj: x2');
for i=3:20*4+1
    fprintf(fileID,' + x%d', i);
end
fprintf(fileID,'\n\n');
fprintf(fileID,'subject to\n\n');
fprintf(fileID,'c1: x1 = 1\n\n');
for i=1:(n-2/t)*(8*(l-3)+4*(l-2-k))+(1/t)*(4*(l-1)+(l-4)*(l-4)+4*(l-2-k))
    fprintf(fileID,'c%d: %dx1', i+1, final(i,1));
    for j=2:20*4+1
        if (final(i,j)==1)
            fprintf(fileID,' + x%d', j);
        end 
    end
    fprintf(fileID,' >= 1\n\n');
end
fprintf(fileID,'bounds\n');
fprintf(fileID,'binary\n');
for i=2:20*4+1
    fprintf(fileID,'x%d\n', i);
end
fclose(fileID);