%NOTES:
%Before running this program, run make_interface_mesh.m
%Run this program to genereate the mesh for the faces of the blocks.
%The output is face_XYZ{block,face} and faces_tri{block,face}
% where "block" indicates the block number and "face" indicates the 
% the face number.
%face_XYZ contains the coordinates of the nodes of the faces.

clear dip xyz

%%%%%%INPUT SECTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%enter (lat,long) coordinates of segment endpoints
%NOTE: The order is important! Enter faces in counter-clockwise direction
%      starting with the northern-most segment.
%      Each block is assigned to a different cell block_face_edges{k}.
%      Each face is assigned to a different matrix block_face_edges{k}(:,:,j)
%      First row of block_face_edges{k}(:,:,j) is lat,long for the southern 
%      endpoint and the second row is lat,long for the northern endpoint
%NOTE: More blocks and faces can be added by adding more lines like below

%BLOCK 1
block_face_edges{1}(:,:,1)=[23.45 120.75;24.8 121.35];
block_face_edges{1}(:,:,2)=[21.9 120.2;23.45 120.75];
%enter dip of each face
dip{1}=[45 45]*pi/180;
isloop(1)=0; %if first and last face connect (faces make closed loop), 
             %then enter 1, otherwise enter 0 
%Enter faces that should be meshed together. For example if you have four
%faces and the first three are meshed together and the fourth by itself,
%enter faces{k,1}=[1 2 3] and faces{k,1}=[4] where k is the block number
faces{1,1}=[1 2]; 
%For neighoring faces, specify which face lies to the north of the 
%intersecting line projected on to the xy plane. For example, if face 4 lies 
%north of face 3, enter face_intersections{k}(4,3)=1 where k is the block
%number. If face 3 lies north of face 4, enter face_intersections{k}(3,4)=1.
%(This should be automated eventually, but for now, do it by hand)
face_intersections{1}=zeros(size(block_face_edges{1}));
face_intersections{1}(1,2)=1;
face_intersections{1}(2,1)=0;
%Specify faces that are connected. For example, if face 2 is connected to
%face 3, enter face_connections{1}(2,3)=1
face_connections{1}=zeros(size(block_face_edges{1}));
face_connections{1}(1,2)=1;
%Specify whether nodes on interface at bottom of block lie above or below
%the plane that defines the face. If bottom nodes lie above the face, set
%above_face{k}(j)=1 where k is block number and j is face number. For
%example, if bottom nodes lie above the second face on the third block, set
%above_face{3}(2)=1
above_face{1}=zeros(1,size(block_face_edges{1},3));
above_face{1}(1)=1;
above_face{1}(2)=1;
%Specify whether nodes on interface at bottom of block lie on the concave-in
%or concave-out side of an intersection. For example, if interface points on 
%bottom of block lie on the cancave-in side of the intersection of faces 1
%and 2, set concave_in{k}(1,2)=1 where k is block number
concave_in{1}=zeros(1,size(block_face_edges{1},3));
concave_in{1}(1,2)=0;

%BLOCK 2
block_face_edges{2}(:,:,1)=[23 120.4;23.45 120.75];
block_face_edges{2}(:,:,2)=[22.62 119.55;23 120.4];
block_face_edges{2}(:,:,3)=[21.9 120.2;23.45 120.75];
dip{2}=[55 45 45]*pi/180;
isloop(2)=1;
faces{2,1}=[1 2];
faces{2,2}=[3];
face_intersections{2}=zeros(size(block_face_edges{2}));
face_intersections{2}(1,2)=1;
face_intersections{2}(2,1)=0;
face_intersections{2}(3,1)=0;
face_intersections{2}(1,3)=0;
face_connections{2}=zeros(size(block_face_edges{2}));
face_connections{2}(1,2)=1;
face_connections{2}(3,1)=1;
above_face{2}=zeros(1,size(block_face_edges{2},3));
above_face{2}(1)=1;
above_face{2}(2)=1;
above_face{2}(3)=0;
concave_in{2}=zeros(1,size(block_face_edges{2},3));
concave_in{2}(1,2)=0;
concave_in{2}(3,1)=1;

density=2;  %scales density of triangles -- for example if you want 
            %to increase the number of triangles by 20%, set density=1.2. 
            %If you want to reduce to 50%, set density=0.5;


%%%%%%END OF INPUT SECTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%change to Cartesian coordinate system (x,y,z)
%make strike, dip, and normal vectors for each segment
for k=1:length(block_face_edges)
    for j=1:size(block_face_edges{k},3)
        block_face_edges{k}(:,:,j) = llh2localxy(block_face_edges{k}(:,:,j)', origin);
        strikevec=[block_face_edges{k}(2,1,j)-block_face_edges{k}(1,1,j) block_face_edges{k}(2,2,j)-block_face_edges{k}(1,2,j) 0];
        strikevec=strikevec/norm(strikevec);
        dipvec=[strikevec(2) -strikevec(1) -tan(dip{k}(j))];
        dipvec=dipvec/norm(dipvec);
        normal{k,j}=cross(dipvec,strikevec);
        strike(k,j)=pi/2-atan2(block_face_edges{k}(2,2,j)-block_face_edges{k}(1,2,j),block_face_edges{k}(2,1,j)-block_face_edges{k}(1,1,j));
        strikevecs_faces{k,j}=strikevec;
        dipvecs_faces{k,j}=dipvec;
        
%generate a set of points along planes defined by the fault segments --
%many of these points do not actually lie on the faces
        x=linspace(0,500,40*density);
        y=linspace(-500,500,80*density);
        [x,y]=meshgrid(x,y);
        X=x(:);Y=y(:);
        ROT=[cos(dip{k}(j)) 0 sin(dip{k}(j));0 1 0;-sin(dip{k}(j)) 0 cos(dip{k}(j))];
        xyrot=ROT*[X';Y';0*X'];
        ROT=[cos(strike(k,j)) sin(strike(k,j)) 0;-sin(strike(k,j)) cos(strike(k,j)) 0;0 0 1];
        xyz{k,j}=ROT*xyrot+repmat([block_face_edges{k}(1,:,j)';0],1,size(xyrot,2));
    
        %dispose of points lying below interface
        INDEX=[];
        for m=1:size(xyz{k,j},2)
            dist=((xyz{k,j}(1,m)-interface_XYZ(1,:)).^2+(xyz{k,j}(2,m)-interface_XYZ(2,:)).^2+(xyz{k,j}(3,m)-interface_XYZ(3,:)).^2);
            [dummy,index]=min(dist);
            if interface_XYZ(3,index)>=xyz{k,j}(3,m); INDEX=[INDEX m];end 
        end
        xyz{k,j}(:,INDEX)=[];
        %discard patches outside southern and northern borders
        index=xyz{k,j}(2,:)<southern_end;
        xyz{k,j}(:,index)=[];
        index=xyz{k,j}(2,:)>northern_end;
        xyz{k,j}(:,index)=[];
        
    end

end
%the rest of the code is concerned with selecting the points defined above
%that lie on the interface
    
%Now dispose of points not lying on faces
    
%loop over number of blocks

for loop=1:length(block_face_edges)

%get intersections of edges
for k=1:size(block_face_edges{loop},3)
    slope(k)=(block_face_edges{loop}(2,2,k)-block_face_edges{loop}(1,2,k))/(block_face_edges{loop}(2,1,k)-block_face_edges{loop}(1,1,k));
    ycrd(k)=block_face_edges{loop}(1,2,k);
    xcrd(k)=block_face_edges{loop}(1,1,k);
end

%find intersection of faces
for k=1:size(block_face_edges{loop},3)-1
    for j=k+1:size(block_face_edges{loop},3)
    tempx=(-ycrd(k)+slope(k)*xcrd(k)+ycrd(j)-slope(j)*xcrd(j))/(slope(k)-slope(j));
    tempy=ycrd(k)+slope(k)*(tempx-xcrd(k));
    %coordinates of intersection
    edgeintersect{k,j}=[tempx tempy];
    end
end

%get intersections of planes
for k=1:size(block_face_edges{loop},3)-1
    for j=k+1:size(block_face_edges{loop},3)
        %vector in direction of intersection
        intersect{k,j}=cross(normal{loop,k},normal{loop,j});
        if intersect{k,j}(3)>0
            intersect{k,j}=-intersect{k,j};
        end
    end
end

    
%dispose of points lying north or south of intersection 
%Here I find the equation for the line (projected into the xy plane) 
%passing through the intersection of the planes and identify points lying 
%above and below the line

for k=1:size(block_face_edges{loop},3)-1
    if face_connections{loop}(k,k+1)==1

    slopeint=intersect{k,k+1}(2)/intersect{k,k+1}(1);
    y=xyz{loop,k}(2,:);
    x=xyz{loop,k}(1,:);
    y1=block_face_edges{loop}(1,2,k);
    x1=block_face_edges{loop}(1,1,k);

    if face_intersections{loop}(k,k+1)==1
        index = y < (y1 + slopeint*(x-x1));
    else
        index = y > (y1 + slopeint*(x-x1));
    end
    xyz{loop,k}(:,index)=[];
    
    y=xyz{loop,k+1}(2,:);
    x=xyz{loop,k+1}(1,:);
    if face_intersections{loop}(k+1,k)==1
        index = y < (y1 + slopeint*(x-x1));
    else
        index = y > (y1 + slopeint*(x-x1));
    end
    xyz{loop,k+1}(:,index)=[];
    end %if isnan
end


%now close the loop
%get intersections of planes at end of loop
if isloop(loop)==1
k=size(block_face_edges{loop},3);
j=1;
intsect=cross(normal{loop,k},normal{loop,j});
if intsect(3)>0; intsect=-intsect; 
end
   if face_intersections{loop}(k,j)==1 
        index=xyz{loop,k}(2,:)<(block_face_edges{loop}(2,2,k)+...
        intsect(2)/intsect(1)*(xyz{loop,k}(1,:)-...
        block_face_edges{loop}(2,1,k)));
    else
        index=xyz{loop,k}(2,:)>(block_face_edges{loop}(2,2,k)+...
        intsect(2)/intsect(1)*(xyz{loop,k}(1,:)-...
        block_face_edges{loop}(2,1,k)));
    end
    xyz{loop,k}(:,index)=[];
    if face_intersections{loop}(j,k)==1
        index=xyz{loop,1}(2,:)<(block_face_edges{loop}(2,2,k)+...
        intsect(2)/intsect(1)*(xyz{loop,1}(1,:)-...
        block_face_edges{loop}(2,1,k)));
    else
        index=xyz{loop,1}(2,:)>(block_face_edges{loop}(2,2,k)+...
        intsect(2)/intsect(1)*(xyz{loop,1}(1,:)-...
        block_face_edges{loop}(2,1,k)));
    end
    xyz{loop,1}(:,index)=[];
end %isloop

end %loop over number of blocks

%store nodes for each face in each block in cell array face_XYZ
for loop1=1:size(faces,1)
    for loop2=1:size(faces,2)
    temp=[];
    for loop3=1:length(faces{loop1,loop2})
        temp=[temp xyz{loop1,faces{loop1,loop2}(loop3)}];
    end
        face_XYZ{loop1,loop2}=temp;
    end
end


for loop1=1:size(face_XYZ,1)
    for loop2=1:size(face_XYZ,2)
        if ~isempty(face_XYZ{loop1,loop2})
        tri = delaunay(face_XYZ{loop1,loop2}(1,:),face_XYZ{loop1,loop2}(2,:),{'Qt','Qbb','Qc','Qz'});
       
        %this disposes of unwanted triangles generated at edges of faces
        sidelengths=sqrt((face_XYZ{loop1,loop2}(1,tri(:,1))-face_XYZ{loop1,loop2}(1,tri(:,2))).^2+(face_XYZ{loop1,loop2}(2,tri(:,1))-face_XYZ{loop1,loop2}(2,tri(:,2))).^2+(face_XYZ{loop1,loop2}(3,tri(:,1))-face_XYZ{loop1,loop2}(3,tri(:,2))).^2)';
        sidelengths=[sidelengths sqrt((face_XYZ{loop1,loop2}(1,tri(:,2))-face_XYZ{loop1,loop2}(1,tri(:,3))).^2+(face_XYZ{loop1,loop2}(2,tri(:,2))-face_XYZ{loop1,loop2}(2,tri(:,3))).^2+(face_XYZ{loop1,loop2}(3,tri(:,2))-face_XYZ{loop1,loop2}(3,tri(:,3))).^2)'];
        sidelengths=max(sidelengths')';
        avglength=median(sidelengths);
        index=sidelengths>avglength*2;
        tri(index,:)=[];
        
        %cell array containing node triplets for each triangle
        faces_tri{loop1,loop2}=tri;
        end %if
    end %loop2
end %loop1

%Now, assign sections of interface to blocks
for k=1:size(face_XYZ,1)
    numfaces=length(dip{k});
        j=1;
        normal=cross(dipvecs_faces{k,j},strikevecs_faces{k,j});
        Z=(-(interface_XYZ(1,:)-block_face_edges{k}(1,1,j))*normal(1)-(interface_XYZ(2,:)-block_face_edges{k}(1,2,j))*normal(2))/normal(3);
        if above_face{k}(j)==1
            index=interface_XYZ(3,:)>Z; %keep these nodes
        else
            index=interface_XYZ(3,:)<Z; %keep these nodes
        end
        
        for j=1:numfaces-1
        normal=cross(dipvecs_faces{k,j+1},strikevecs_faces{k,j+1});
        Z=(-(interface_XYZ(1,:)-block_face_edges{k}(1,1,j+1))*normal(1)-(interface_XYZ(2,:)-block_face_edges{k}(1,2,j+1))*normal(2))/normal(3);
        if above_face{k}(j+1)==1
            index2=interface_XYZ(3,:)>Z; %keep these nodes
        else
            index2=interface_XYZ(3,:)<Z; %keep these nodes
        end
        
        if face_connections{k}(j,j+1)==1
            if concave_in{k}(j,j+1)==1
                index=index.*index2;
            else
                index=index+index2;
                index=index>0;
            end
        else
            index=index.*index2;
        end %if face_connections
  
    end%j
    index=logical(index);
    block_bottom{k}=interface_XYZ(:,index);
    
       tri = delaunay(block_bottom{k}(1,:),block_bottom{k}(2,:),{'Qt','Qbb','Qc','Qz'});
       
        %this disposes of unwanted triangles generated at edges of faces
        sidelengths=sqrt((block_bottom{k}(1,tri(:,1))-block_bottom{k}(1,tri(:,2))).^2+(block_bottom{k}(2,tri(:,1))-block_bottom{k}(2,tri(:,2))).^2+(block_bottom{k}(3,tri(:,1))-block_bottom{k}(3,tri(:,2))).^2)';
        sidelengths=[sidelengths sqrt((block_bottom{k}(1,tri(:,2))-block_bottom{k}(1,tri(:,3))).^2+(block_bottom{k}(2,tri(:,2))-block_bottom{k}(2,tri(:,3))).^2+(block_bottom{k}(3,tri(:,2))-block_bottom{k}(3,tri(:,3))).^2)'];
        sidelengths=max(sidelengths')';
        avglength=median(sidelengths);
        index=sidelengths>avglength*2;
        tri(index,:)=[];
             
        bottoms_tri{k}=tri;

end%k

%plot each block with separate color
figure; hold on
for loop1=1:size(face_XYZ,1)
    color=rand(3,1);
    for loop2=1:size(face_XYZ,2)
        if ~isempty(face_XYZ{loop1,loop2})
            trisurf(faces_tri{loop1,loop2},face_XYZ{loop1,loop2}(1,:),face_XYZ{loop1,loop2}(2,:),face_XYZ{loop1,loop2}(3,:),'FaceColor',color)
        end
     end
    trisurf(bottoms_tri{loop1},block_bottom{loop1}(1,:),block_bottom{loop1}(2,:),block_bottom{loop1}(3,:),'FaceColor',color)
end
axis equal
xlabel('easting')
ylabel('northing')
plot(taiwan(:,1),taiwan(:,2))



