M=dlmread('Costa1.xyz',' ');
M=flip(M);
%M=sortrows(M);

N=dlmread('Roca.xyz',' ');
N=flip(N);
%N=sortrows(N);

%M([2,1],:)=M([1,2],:);

ID1=zeros(711,1);
ID2=zeros(16,1);


M1=[M(1,1) 3119571.25 0];
M2=[M(709,1) 3119571.25 0];

M=[M1;M;M2];

[C,ia,ib]=unique(M,'rows','stable');
i=true(size(M,1),1);
i(ia)=false;
M(i,:)=0;

%PERM=[1 3 2];

%M=M(:,PERM);
%M(:,3)=0;
%M(1,2)=800;
%M(166,2)=800;

L1=zeros(711,2);

L2=zeros(16,2);

for i=1:711
    L1(i,1)=i;
    L1(i,2)=i+1;
end

for i=1:16
    L2(i,1)=711+i;
    L2(i,2)=711+i+1;
end

L1(711,:)=[711 1];
L2(16,:)=[727 712];

S=1;
for i=1:711
    ID1(i,:)=i;
end

for i=1:16
    ID2(i,:)=711+i;
end

M=[ID1 M];
L1=[ID1 L1];

N=[ID2 N];
L2=[ID2 L2];
% GMsh geometry file is stored as .txt file
fileID = fopen('dominio.geo','w');

fprintf(fileID,'SetFactory("OpenCASCADE");\n');
%fprintf(fileID,'Mesh.Algorithm=7;\n');
%fprintf(fileID,'Mesh.Algorithm3D=7;\n');
fprintf(fileID,'//+\n');

for i=1:711
fprintf(fileID,'Point(%d)={%d,%d,%d};\n',M(i,:));
fprintf(fileID,'//+\n');
end

for i=1:16
fprintf(fileID,'Point(%d)={%d,%d,%d};\n',N(i,:));
fprintf(fileID,'//+\n');
end

fprintf(fileID,'///////////////////////////////////////////\n');

for i=1:711
fprintf(fileID,'Line(%d)={%d,%d};\n',L1(i,:));
fprintf(fileID,'//+\n');
end

for i=1:16
fprintf(fileID,'Line(%d)={%d,%d};\n',L2(i,:));
fprintf(fileID,'//+\n');
end

allOneString1 = sprintf('%.0f,' , ID1);
allOneString2 = sprintf('%.0f,' , ID2);

allOneString1 = allOneString1(1:end-1);% strip final comma
allOneString2 = allOneString2(1:end-1);% strip final comma
allOneString11= extractBefore(allOneString1,2 );% oeste
allOneString12= extractAfter(allOneString1,length(allOneString1)-3);% norte
allOneString13= extractBetween(allOneString1,length(allOneString1)-6,length(allOneString1)-4);% norte
allOneString13=char(allOneString13);
allOneString14= extractBetween(allOneString1,3,length(allOneString1)-7);
allOneString14=char(allOneString14);
allOneString14=strcat(allOneString14,allOneString2);%costa


fprintf(fileID,'Line Loop(1) = {%s};\n',allOneString1);
fprintf(fileID,'Line Loop(2) = {%s};\n',allOneString2);
fprintf(fileID,'//+\n');
fprintf(fileID,'Plane Surface(1) = {1,2};\n');
fprintf(fileID,'//+\n');

fprintf(fileID,'Characteristic Length{%s}=15;\n',allOneString14);
fprintf(fileID,'//+\n');
fprintf(fileID,'Physical Line("oeste",1)={%s};\n',allOneString11);
fprintf(fileID,'Physical Line("norte",2)={%s};\n',allOneString12);
fprintf(fileID,'Physical Line("este",3)={%s};\n',allOneString13);
fprintf(fileID,'Physical Line("costa",4)={%s};\n',allOneString14);
fprintf(fileID,'Physical Surface("mar",5)={1};\n');

