%NOTES:
%Run this program to genereate the mesh for the subduction interface.
%The outputs are interface_tri and interface_XYZ
%interface_XYZ contains the coordinates of all the nodes for the mesh
%interface_tri contains the triplets of nodes for all the triangles

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

%upper edge of interface segments
%enter (lat,long) coordinates of segment endpoints
%NOTE: The order is important! Enter segments from north to south.
%      Each segment is assigned to a different cell interface_edge{k}.
%      First row of interface_edge{k} is lat,long for the southern endpoint
%      and the second row is lat,long for the northern endpoint
interface_edge{1}=[23  119.4;25 120.25];
interface_edge{2}=[22.25 119.7;23  119.4];
interface_edge{3}=[21.7 120.5;22.25  119.7];
interface_edge{4}=[20 120.5;21.7 120.5];
%enter dip of each segment
dip=[10 10 12 15]*pi/180;
%maximum depth (km) of patches on interface
maxdepth=35;

%origin for cartesian coordinate system
origin = [22 120];
density=1;  %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;

%maximum extent of interface to the south and north in local xy coordinates
southern_end=-300; 
northern_end=400;

%intersection matrix
%Enter "1" if segments intersect. For example, if segment 1 intersects
%segment 2 then enter doesintersect(1,2)=1 and doesintersect(2,1)=1
%You may need to experiment with this to get it right.
doesintersect=zeros(length(interface_edge));
doesintersect(1,2)=1; doesintersect(2,1)=1;
doesintersect(1,3)=1; doesintersect(3,1)=1;
doesintersect(2,3)=1; doesintersect(3,2)=1;
doesintersect(3,4)=1; doesintersect(4,3)=1;

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

load TaiwanCoastline
taiwan = llh2localxy(fliplr(taiwan)', origin);

%change to Cartesian coordinate system (x,y,z)
%make strike, dip, and normal vectors for each segment
for k=1:length(interface_edge)
interface_edge{k} = llh2localxy(interface_edge{k}', origin);
strikevec=[interface_edge{k}(2,1)-interface_edge{k}(1,1) interface_edge{k}(2,2)-interface_edge{k}(1,2) 0];
strikevec=strikevec/norm(strikevec);
dipvec=[strikevec(2) -strikevec(1) -tan(dip(k))];
dipvec=dipvec/norm(dipvec);
normal{k}=cross(dipvec,strikevec);
strike(k)=pi/2-atan2(interface_edge{k}(2,2)-interface_edge{k}(1,2),interface_edge{k}(2,1)-interface_edge{k}(1,1));
end
%Note: first column of interface_edge{k} is x-coordinate, second column
%y-coordinate   

%generate a set of points along planes defined by the fault segments --
%many of these points do not actually lie on the interface
for k=1:length(interface_edge)
    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)) 0 sin(dip(k));0 1 0;-sin(dip(k)) 0 cos(dip(k))];
    xyrot=ROT*[X';Y';0*X'];
    ROT=[cos(strike(k)) sin(strike(k)) 0;-sin(strike(k)) cos(strike(k)) 0;0 0 1];
    xyz{k}=ROT*xyrot+repmat([interface_edge{k}(1,:)';0],1,size(xyrot,2));

    %these are needed to find intersection of edges
    slope(k)=(interface_edge{k}(2,2)-interface_edge{k}(1,2))/(interface_edge{k}(2,1)-interface_edge{k}(1,1));
    ycrd(k)=interface_edge{k}(1,2); 
    xcrd(k)=interface_edge{k}(1,1);
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 interface

%find intersection of edges
for k=1:length(interface_edge)-1
    for j=k+1:length(interface_edge)
    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:length(interface_edge)-1
    for j=k+1:length(interface_edge)
        %vector in direction of intersection
        intersect{k,j}=cross(normal{k},normal{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 
%NOTE: This probably only works for interface edges dipping to the east.
%Need to generalize this algorithm...
%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

%discard nodes from intersecting segments
for k=1:length(interface_edge)-1
    azimuth=(atan2(intersect{k,k+1}(2),intersect{k,k+1}(1))*180/pi);
    slopeint=intersect{k,k+1}(2)/intersect{k,k+1}(1);
    y=xyz{k}(2,:);
    x=xyz{k}(1,:);
    y1=interface_edge{k}(1,2);
    x1=interface_edge{k}(1,1);
    if azimuth < 90
        index= y < (y1 + slopeint*(x-x1));
    else
        index= y > (y1 + slopeint*(x-x1));
    end
    xyz{k}(:,index)=[];
    
    y=xyz{k+1}(2,:);
    x=xyz{k+1}(1,:);
    if azimuth < 90
        index = y > (y1 + slopeint*(x-x1));
    else
        index = y < (y1 + slopeint*(x-x1));
    end
    xyz{k+1}(:,index)=[];
end

%discard nodes from non-neighboring segments
for k=1:length(interface_edge)-1
    for j=k+1:length(interface_edge)
        if doesintersect(k,j)==1

    azimuth=(atan2(intersect{k,j}(2),intersect{k,j}(1))*180/pi);
    slopeint=intersect{k,j}(2)/intersect{k,j}(1);
    y=xyz{k}(2,:);
    x=xyz{k}(1,:);
    y1=edgeintersect{k,j}(1,2);
    x1=edgeintersect{k,j}(1,1);

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

%store all nodes in matrix XYZ
XYZ=[];
for k=1:length(interface_edge)
    XYZ=[XYZ xyz{k}];
end

%discard patches below maxdepth
index=abs(XYZ(3,:))>maxdepth;
XYZ(:,index)=[];
%discard patches outside southern and northern borders
index=XYZ(2,:)<southern_end;
XYZ(:,index)=[];
index=XYZ(2,:)>northern_end;
XYZ(:,index)=[];

%plot the mesh
figure; hold on
%tri=delaunay(XYZ(1,:),XYZ(2,:));
tri = delaunay(XYZ(1,:),XYZ(2,:),{'Qt','Qbb','Qc','Qz'});


%this disposes of unwanted triangles generated at edges of faces
sidelengths(:,1)=sqrt((XYZ(1,tri(:,1))-XYZ(1,tri(:,2))).^2+(XYZ(2,tri(:,1))-XYZ(2,tri(:,2))).^2+(XYZ(3,tri(:,1))-XYZ(3,tri(:,2))).^2);
sidelengths(:,2)=sqrt((XYZ(1,tri(:,2))-XYZ(1,tri(:,3))).^2+(XYZ(2,tri(:,2))-XYZ(2,tri(:,3))).^2+(XYZ(3,tri(:,2))-XYZ(3,tri(:,3))).^2);
sidelengths=max(sidelengths')';
avglength=mean(sidelengths);
index=sidelengths>avglength*2;
tri(index,:)=[];

trisurf(tri,XYZ(1,:),XYZ(2,:),XYZ(3,:),'FaceColor',rand(3,1))
axis equal
xlabel('easting')
ylabel('northing')
plot(taiwan(:,1),taiwan(:,2))
title('Interface Mesh')


%plot the nodes only
figure;hold on
for k=1:length(xyz)
    plot3(xyz{k}(1,:),xyz{k}(2,:),xyz{k}(3,:),'.','Color',rand(3,1))
end
axis equal
xlabel('easting')
ylabel('northing')
plot(taiwan(:,1),taiwan(:,2))
title('Interface Nodes')


%x=linspace(-400,600,10);
%y=linspace(-1000,600,10);
%[X,Y]=meshgrid(x,y);
%xyz=[X(:)';Y(:)';0*X(:)'];
%xp{1}=[XYZ(1,:)' XYZ(2,:)' XYZ(3,:)'];
%triangles=tri;
%for k=1:size(tri,1)
%    tris{1}=triangles(k,:);
%    u = tridisloc(xyz, xp, tris, [1 0 0], .25);
%    G(:,k)=u(:);
%end

%%%%%OUTPUTS%%%%%%
interface_tri=tri;
interface_XYZ=XYZ;

