% Case 3: multiple constant slope fans on a faceted topography
% Tzu-Yin Kasha Chen, Dec. 2020
clearvars; %close all; clc

Case = 'd';
% 'a': multiple fans with the same slope on a steep inclined plane
% 'b': multiple fans with various slopes on a steep inclined plane
% 'c': bajadas, fans with the same slope on a series of faceted valleys
% 'd': bajadas, fans with with various slopes on a series of faceted valleys

%% case 3(a): multiple fans with the same slope on a steep inclined plane
if Case=='a'
    % boundary surface
    xMin = -60; xMax = 60;
    yMin = 0; yMax = 30;
    zMin = 0; zMax = 20;
    dMesh = 0.1;
    [xMesh,yMesh] = meshgrid(xMin:dMesh:xMax,yMin:dMesh:yMax);
    tanBeta = 1; % slope of the inclined wall
    yToe = 15;
    zRock = zeros(size(xMesh))+zMin;
    zRock(yMesh>yToe & yMesh<yToe+zMax/tanBeta) = (yMesh(yMesh>yToe & yMesh<yToe+zMax/tanBeta)-yToe)*tanBeta;
    zRock(yMesh>=yToe+zMax/tanBeta) = zMax;
    
    % fan apex and slope
    xApex = -50+(20:20:100); yApex = yToe+[12,6,8,10,4];
    zApex = interp2(xMesh,yMesh,zRock,xApex,yApex);
    tanAlpha = 0.5;
    K = zeros(size(xApex));
    
    % computational solution with visibility method
    [zFan_visi,kFan,xyzkApex] = FanTopo_slope(xMesh,yMesh,zRock,xApex,yApex,zApex,tanAlpha*ones(size(xApex)),K,0);
    V_visi = sum(zFan_visi(:)-zRock(:))*dMesh^2;
    kParent = find(isnan(xyzkApex(:,4)));
    kParent = [kParent;length(xyzkApex)+length(kParent)+1];
    for i = 1:length(kParent)-1
            kFan(kFan>=kParent(i) & kFan<kParent(i+1))=i;
    end
    
    % analytical solution of the weld lines
    phi = asin(tanAlpha/tanBeta);
    B = tanBeta/tanAlpha;
    xWeldline = zeros(length(xApex)-1,50);
    yWeldline = zeros(length(xApex)-1,50);
    zWeldline = zeros(length(xApex)-1,50);
    xd = zeros(length(xApex)-1,1);
    yd = zeros(length(xApex)-1,1);
    xu = zeros(length(xApex)-1,1);
    yu = zeros(length(xApex)-1,1);
    zu = zeros(length(xApex)-1,1);
    for i = 1:length(xApex)-1
        x1 = xApex(i); x2 = xApex(i+1);
        y1 = yApex(i); y2 = yApex(i+1);
        z1 = zApex(i); z2 = zApex(i+1);
        xu(i) = (x1+x2+(y1-y2)*cot(phi))/2;
        yu(i) = (y1+y2+(x1-x2)*tan(phi))/2;
        zu(i) = (yu(i)-yToe)*tanBeta;
        syms x y
        eqToe1 = (y-yToe)^2-2*(y1-yToe)*(y-yToe)+(1-B^2)*(y1-yToe)^2+(x-x1)^2==0;
        eqToe2 = (y-yToe)^2-2*(y2-yToe)*(y-yToe)+(1-B^2)*(y2-yToe)^2+(x-x2)^2==0;
        [solxd,solyd] = vpasolve([eqToe1,eqToe2],[x,y],[xMin,xMax;-inf,yu(i)]);
        xd(i) = double(solxd);
        yd(i) = double(solyd);
        yWeldline(i,:) = linspace(yu(i),yd(i),length(yWeldline));
        for j = 1:length(yWeldline)
            yt = yWeldline(i,j);
            S = vpasolve((z1-z2)/tanAlpha-sqrt((x-x1).^2+(yt-y1).^2)+sqrt((x-x2).^2+(yt-y2).^2)==0,x,[x1,x2]);
            xWeldline(i,j) =  double(S);
        end
        zWeldline(i,:) = z1-tanAlpha*sqrt((xWeldline(i,:)-x1).^2+(yWeldline(i,:)-y1).^2);
    end
    %% comparison
    figure
    axis equal
    hold on
    axis off
    pcolor(xMesh,yMesh,kFan)
    c = colorcube(64);
    c = c([20;24;35;26;38;23],:);
    colormap(c)
    caxis([1 max(unique(kFan(~isnan(kFan))))])
    shading flat
    zFan = zFan_visi;
    zFan(isnan(zFan)) = zRock(isnan(zFan));
    axis([xMin+5 xMax yMin yMax zMin zMax])
    contour(xMesh,yMesh,zFan,0:2:zMax,'-','color',[0.5,0.5,0.5])
    axis equal
    hold on
    axis off
    plot(xWeldline',yWeldline','k-','linewidth',1)
    %% case 3(b): multiple fans with various slopes on a steep inclined plane
elseif Case=='b'
    % boundary surface
    xMin = -60; xMax = 60;
    yMin = 0; yMax = 30;
    zMin = 0; zMax = 20;
    dMesh = 0.1;
    [xMesh,yMesh] = meshgrid(xMin:dMesh:xMax,yMin:dMesh:yMax);
    tanBeta = 1; % slope of the inclined wall
    yToe = 15;
    zRock = zeros(size(xMesh))+zMin;
    zRock(yMesh>yToe & yMesh<yToe+zMax/tanBeta) = (yMesh(yMesh>yToe & yMesh<yToe+zMax/tanBeta)-yToe)*tanBeta;
    zRock(yMesh>=yToe+zMax/tanBeta) = zMax;
    
    % fan apex and slope
    xApex = -60+10+(20:20:80); yApex = yToe+[12,6,8,10];
    zApex = interp2(xMesh,yMesh,zRock,xApex,yApex);
    tanAlpha = 0.5*[1,0.6,0.8,0.7];
    K = tanAlpha*0;
    
    % computational solution with visibility method
    [zFan_visi,kFan,xyzkApex] = FanTopo_slope(xMesh,yMesh,zRock,xApex,yApex,zApex,tanAlpha,K,0);
    V_visi = sum(zFan_visi(:)-zRock(:))*dMesh^2;
    kParent = find(isnan(xyzkApex(:,4)));
    kParent = [kParent;length(xyzkApex)+length(kParent)+1];
    for i = 1:length(kParent)-1
            kFan(kFan>=kParent(i) & kFan<kParent(i+1))=i;
    end
    %%
    % analytical solution of the weld lines
    xWeldline = zeros(length(xApex)-1,50);
    yWeldline = zeros(length(xApex)-1,50);
    zWeldline = zeros(length(xApex)-1,50);
    xd = zeros(length(xApex)-1,1);
    yd = zeros(length(xApex)-1,1);
    xu = zeros(length(xApex)-1,1);
    yu = zeros(length(xApex)-1,1);
    zu = zeros(length(xApex)-1,1);
    for i = 1:length(xApex)-1
        x1 = xApex(i); x2 = xApex(i+1);
        y1 = yApex(i); y2 = yApex(i+1);
        z1 = zApex(i); z2 = zApex(i+1);
        S1 = tanAlpha(i); S2 = tanAlpha(i+1);
        B1 = tanBeta/S1; B2 = tanBeta/S2;
        phi1 = asin(S1/tanBeta);
        phi2 = asin(S2/tanBeta);
        xu(i) = (x1*tan(phi1)+x2*tan(phi2)+(y1-y2))/(tan(phi1)+tan(phi2));
        yu(i) = y1-(xu(i)-x1)*tan(phi1);
        zu(i) = (yu(i)-yToe)*tanBeta;
        syms x y
        eqToe1 = (y-yToe)^2-2*(y1-yToe)*(y-yToe)+(1-B1^2)*(y1-yToe)^2+(x-x1)^2==0;
        eqToe2 = (y-yToe)^2-2*(y2-yToe)*(y-yToe)+(1-B2^2)*(y2-yToe)^2+(x-x2)^2==0;
        [solxd,solyd] = vpasolve([eqToe1,eqToe2],[x,y],[xMin,xMax;-inf,yu(i)]);
        xd(i) = double(solxd);
        yd(i) = double(solyd);
        yWeldline(i,:) = linspace(yu(i),yd(i),length(yWeldline));
        for j = 1:length(yWeldline)
            yt = yWeldline(i,j);
            S = vpasolve((z1-z2)-S1*sqrt((x-x1).^2+(yt-y1).^2)+S2*sqrt((x-x2).^2+(yt-y2).^2)==0,x,[xMin,xMax]);
            xWeldline(i,j) =  double(S);
        end
        zWeldline(i,:) = z1-S1*sqrt((xWeldline(i,:)-x1).^2+(yWeldline(i,:)-y1).^2);
    end
    %% comparison
    figure
    axis equal
    hold on
    axis off
    pcolor(xMesh,yMesh,kFan)
    shading flat
    c = colorcube(64);
    c = c([20;24;35;26;38;23],:);
    colormap(c)
    caxis([1 max(kFan(:))])
    view(0,90)
    axis([xMin+5 xMax yMin yMax zMin zMax])
    axis off
    axis equal
    contour(xMesh,yMesh,zFan_visi,0:2:zMax,'-','color',[0.5,0.5,0.5])
    plot(xWeldline',yWeldline','k-','linewidth',1)
    %% case 3(c): bajadas, fans with the same slope on a series of faceted valleys
elseif Case =='c'
    % boundary surface
    xMin = -60; xMax = 60;
    yMin = 5; yMax = 35;
    zMin = 0; zMax = 30;
    dMesh = 0.1;
    [xMesh,yMesh] = meshgrid(xMin:dMesh:xMax,yMin:dMesh:yMax);
    tanBeta = 1; % slope of the inclined wall
    yToe = 15;
    psi = atan(3/4); % rotation angle of contours in the valley
    zRock = zeros(size(xMesh))+zMin;
    zRock(yMesh>yToe & yMesh<yToe+zMax/tanBeta) = (yMesh(yMesh>yToe & yMesh<yToe+zMax/tanBeta)-yToe)*tanBeta;
    %     zRock(yMesh>=yToe+zMax/tanBeta) = zMax;
    dValleyWallRight = (xMesh-0)*sin(psi)+(yMesh-yToe)*cos(psi);
    zValleyWallRight = dValleyWallRight*tanBeta;
    zValleyWallLeft = fliplr(zValleyWallRight);
    zRock(xMesh>=0 & yMesh>yToe & zRock>zValleyWallRight) = zValleyWallRight(xMesh>=0 & yMesh>yToe & zRock>zValleyWallRight);
    zRock(xMesh<=0 & yMesh>yToe & zRock>zValleyWallLeft) = zValleyWallLeft(xMesh<=0 & yMesh>yToe & zRock>zValleyWallLeft);
    zRock(xMesh>=-30 & xMesh<=-10) = zRock(xMesh>=-10 & xMesh<=10);
    zRock(xMesh>=-50 & xMesh<=-30) = zRock(xMesh>=-10 & xMesh<=10);
    zRock(xMesh>=10 & xMesh<=30) = zRock(xMesh>=-10 & xMesh<=10);
    zRock(xMesh>=30 & xMesh<=50) = zRock(xMesh>=-10 & xMesh<=10);
    
    % fan apex and slope
    xApex = [-40,-20,0,20,40]; yApex = yToe+[10,20,12,14,17];
    zApex = interp2(xMesh,yMesh,zRock,xApex,yApex);
    tanAlpha = 0.5;
    
    % computational solution with visibility method
    [zFan_visi,kFan,xyzkApex] = FanTopo_slope(xMesh,yMesh,zRock,xApex,yApex,zApex,tanAlpha*ones(size(xApex)),zeros(size(xApex)),0);
    V_visi = sum(zFan_visi(:)-zRock(:))*dMesh^2;
    kParent = find(isnan(xyzkApex(:,4)));
    kParent = [kParent;length(xyzkApex)+length(kParent)+1];
    for i = 1:length(kParent)-1
            kFan(kFan>=kParent(i) & kFan<kParent(i+1))=i;
    end
    
    % analytical solution
    zFan_analytical = zeros(size(xMesh))*nan;
    phi = asin(tanAlpha/tanBeta); % angle between the fan margin and the contour of the inclined wall
    for i = 1:length(xApex)
        psiRidge = (pi-psi)/2;
        yApexC(i) = (yApex(i)/tan(phi+psi)+yToe/tan(psiRidge))/(1/tan(psiRidge)+1/tan(phi+psi));
        xApexC(i) = xApex(i)+(yApexC(i)-yToe)/tan(psiRidge);
        zApexC(i) = zApex(i)-(yApex(i)-yApexC(i))/sin(phi+psi)*tanAlpha;
        yApexB(i) = yApexC(i);
        xApexB(i) = xApex(i)-(yApexC(i)-yToe)/tan(psiRidge);
        zApexB(i) = zApexC(i);
    end
    xWeldline = zeros(length(xApex)-1,50);
    yWeldline = zeros(length(xApex)-1,50);
    zWeldline = zeros(length(xApex)-1,50);
    xd = zeros(length(xApex)-1,1);
    yd = zeros(length(xApex)-1,1);
    xu = zeros(length(xApex)-1,1);
    yu = zeros(length(xApex)-1,1);
    zu = zeros(length(xApex)-1,1);
    for i = 1:length(xApex)-1
        x1 = xApexC(i); x2 = xApexB(i+1);
        y1 = yApexC(i); y2 = yApexB(i+1);
        z1 = zApexC(i); z2 = zApexB(i+1);
        x10 = xApex(i); x20 = xApex(i+1);
        y10 = yApex(i); y20 = yApex(i+1);
        z10 = zApex(i); z20 = zApex(i+1);
        S1 = tanAlpha; S2 = tanAlpha;
        B1 = tanBeta/S1; B2 = tanBeta/S2;
        phi1 = asin(S1/tanBeta);
        phi2 = asin(S2/tanBeta);
        xu(i) = (x1*tan(phi1)+x2*tan(phi2)+(y1-y2))/(tan(phi1)+tan(phi2));
        yu(i) = y1-(xu(i)-x1)*tan(phi1);
        zu(i) = (yu(i)-yToe)*tanBeta;
        syms x y
        eqToe1 = (y-yToe)^2-2*(y1-yToe)*(y-yToe)+(1-B1^2)*(y1-yToe)^2+(x-x1)^2==0;
        eqToe2 = (y-yToe)^2-2*(y2-yToe)*(y-yToe)+(1-B2^2)*(y2-yToe)^2+(x-x2)^2==0;
        [solxd,solyd] = vpasolve([eqToe1,eqToe2],[x,y],[xMin,xMax;-inf,yu(i)]);
        xd(i) = double(solxd);
        yd(i) = double(solyd);
        if z1<z2
            eqWeld1 = (z10-z2)-S1*sqrt((x-x10).^2+(y-y10).^2)+S2*sqrt((x-x2).^2+(y-y2).^2)==0;
            eqVisi1 = (x-x10)*(y1-y10)-(x1-x10)*(y-y10)==0;
            [solx,soly] = vpasolve([eqWeld1,eqVisi1],[x,y],[xMin,xMax;-inf,yu(i)]);
        else
            eqWeld2 = (z1-z20)-S1*sqrt((x-x1).^2+(y-y1).^2)+S2*sqrt((x-x20).^2+(y-y20).^2)==0;
            eqVisi2 = (x-x20)*(y2-y20)-(x2-x20)*(y-y20)==0;
            [solx,soly] = vpasolve([eqWeld2,eqVisi2],[x,y],[xMin,xMax;-inf,yu(i)]);
        end
        if double(soly)<yd(i)
            yWeldline(i,:) = linspace(yu(i),yd(i),length(yWeldline));
            n1 = length(yWeldline);
        elseif z1<z2
            eqToe1 = (y-yToe)^2-2*(y10-yToe)*(y-yToe)+(y10-yToe)^2-z10^2/S1^2+(x-x10)^2==0;
            [solxd,solyd] = vpasolve([eqToe1,eqToe2],[x,y],[xMin,xMax;-inf,yu(i)]);
            n1 = round(length(yWeldline)*(double(soly)-yu(i))/(double(solyd)-yu(i)));
            yWeldline(i,1:n1) = linspace(yu(i),double(soly),n1);
            yWeldline(i,n1+1:end) = linspace(double(soly)+yWeldline(i,2)-yWeldline(i,1),double(solyd),length(yWeldline)-n1);
            plot([x10,x1,double(solx)],[y10,y1,double(soly)],'r.')
            plot(double(solxd),double(solyd),'r.')
        else
            eqToe2 = (y-yToe)^2-2*(y20-yToe)*(y-yToe)+(y20-yToe)^2-z20^2/S2^2+(x-x20)^2==0;
            [solxd,solyd] = vpasolve([eqToe1,eqToe2],[x,y],[xMin,xMax;-inf,yu(i)]);
            n1 = round(length(yWeldline)*(double(soly)-yu(i))/(double(solyd)-yu(i)));
            yWeldline(i,1:n1) = linspace(yu(i),double(soly),n1);
            yWeldline(i,n1+1:end) = linspace(double(soly)+yWeldline(i,2)-yWeldline(i,1),double(solyd),length(yWeldline)-n1);
            plot([x20,x2,double(solx)],[y20,y2,double(soly)],'r.')
            plot(double(solxd),double(solyd),'r.')
        end
        for j = 1:length(yWeldline)
            yt = yWeldline(i,j);
            if j==n1+1 && z1<z2
                x1=x10; y1=y10; z1=z10;
            elseif j==n1+1
                x2=x20; y2=y20; z2=z20;
            end
            S = vpasolve((z1-z2)-S1*sqrt((x-x1).^2+(yt-y1).^2)+S2*sqrt((x-x2).^2+(yt-y2).^2)==0,x,[xMin,xMax]);
            xWeldline(i,j) =  double(S);
            zWeldline(i,j) = z1-S1*sqrt((xWeldline(i,j)-x1).^2+(yWeldline(i,j)-y1).^2);
        end
    end
    
    % comparison
    figure
    axis equal
    hold on
    axis off
    pcolor(xMesh,yMesh,kFan)
    c = colorcube(64);
    c = c([20;24;35;26;38;23],:);
    colormap(c)
    caxis([1 max(kFan(:))])
    shading flat
    zFan = zFan_visi;
    zFan(isnan(zFan)) = zRock(isnan(zFan));
    axis([xMin+5 xMax yMin yMax zMin zMax])
    contour(xMesh,yMesh,zFan,0:2:zMax,'-','color',[0.5,0.5,0.5])
    axis equal
    hold on
    axis off
    plot(xWeldline',yWeldline','k-','linewidth',1)
    
    %% case 3(d): bajadas, fans with various slopes on a series of faceted valleys
elseif Case =='d'
    %%
    % boundary surface
    xMin = -60; xMax = 60;
    yMin = 5; yMax = 35;
    zMin = 0; zMax = 30;
    dMesh = 0.1;
    [xMesh,yMesh] = meshgrid(xMin:dMesh:xMax,yMin:dMesh:yMax);
    tanBeta = 1; % slope of the inclined wall
    yToe = 15;
    psi = atan(3/4); % rotation angle of contours in the valley
    zRock = zeros(size(xMesh))+zMin;
    zRock(yMesh>yToe & yMesh<yToe+zMax/tanBeta) = (yMesh(yMesh>yToe & yMesh<yToe+zMax/tanBeta)-yToe)*tanBeta;
    dValleyWallRight = (xMesh-0)*sin(psi)+(yMesh-yToe)*cos(psi);
    zValleyWallRight = dValleyWallRight*tanBeta;
    zValleyWallLeft = fliplr(zValleyWallRight);
    zRock(xMesh>=0 & yMesh>yToe & zRock>zValleyWallRight) = zValleyWallRight(xMesh>=0 & yMesh>yToe & zRock>zValleyWallRight);
    zRock(xMesh<=0 & yMesh>yToe & zRock>zValleyWallLeft) = zValleyWallLeft(xMesh<=0 & yMesh>yToe & zRock>zValleyWallLeft);
    zRock(xMesh>=-30 & xMesh<=-10) = zRock(xMesh>=-10 & xMesh<=10);
    zRock(xMesh>=-50 & xMesh<=-30) = zRock(xMesh>=-10 & xMesh<=10);
    zRock(xMesh>=10 & xMesh<=30) = zRock(xMesh>=-10 & xMesh<=10);
    zRock(xMesh>=30 & xMesh<=50) = zRock(xMesh>=-10 & xMesh<=10);
    
    % fan apex and slope
    xApex = [-40,-20,0,20,40]; yApex = yToe+[10,15,12,14,17];
    zApex = interp2(xMesh,yMesh,zRock,xApex,yApex);
    tanAlpha = 0.5*[0.9,1.2,0.9,1.2,1];
    K = tanAlpha*0;
    
    % computational solution with visibility method
    [zFan_visi,kFan,xyzkApex] = FanTopo_slope(xMesh,yMesh,zRock,xApex,yApex,zApex,tanAlpha,K,0);
    V_visi = sum(zFan_visi(:)-zRock(:))*dMesh^2;
    kParent = find(isnan(xyzkApex(:,4)));
    kParent = [kParent;length(xyzkApex)+length(kParent)+1];
    for i = 1:length(kParent)-1
            kFan(kFan>=kParent(i) & kFan<kParent(i+1))=i;
    end
    
    % analytical solution
    for i = 1:length(xApex)
        phi = asin(tanAlpha(i)/tanBeta); % angle between the fan margin and the contour of the inclined wall
        psiRidge = (pi-psi)/2;
        yApexC(i) = (yApex(i)/tan(phi+psi)+yToe/tan(psiRidge))/(1/tan(psiRidge)+1/tan(phi+psi));
        xApexC(i) = xApex(i)+(yApexC(i)-yToe)/tan(psiRidge);
        zApexC(i) = zApex(i)-(yApex(i)-yApexC(i))/sin(phi+psi)*tanAlpha(i);
        yApexB(i) = yApexC(i);
        xApexB(i) = xApex(i)-(yApexC(i)-yToe)/tan(psiRidge);
        zApexB(i) = zApexC(i);
    end
    xWeldline = zeros(length(xApex)-1,50);
    yWeldline = zeros(length(xApex)-1,50);
    zWeldline = zeros(length(xApex)-1,50);
    xd = zeros(length(xApex)-1,1);
    yd = zeros(length(xApex)-1,1);
    xu = zeros(length(xApex)-1,1);
    yu = zeros(length(xApex)-1,1);
    zu = zeros(length(xApex)-1,1);
    for i = 1:length(xApex)-1
        x1 = xApexC(i); x2 = xApexB(i+1);
        y1 = yApexC(i); y2 = yApexB(i+1);
        z1 = zApexC(i); z2 = zApexB(i+1);
        x10 = xApex(i); x20 = xApex(i+1);
        y10 = yApex(i); y20 = yApex(i+1);
        z10 = zApex(i); z20 = zApex(i+1);
        S1 = tanAlpha(i); S2 = tanAlpha(i+1);
        B1 = tanBeta/S1; B2 = tanBeta/S2;
        phi1 = asin(S1/tanBeta);
        phi2 = asin(S2/tanBeta);
        xu(i) = (x1*tan(phi1)+x2*tan(phi2)+(y1-y2))/(tan(phi1)+tan(phi2));
        yu(i) = y1-(xu(i)-x1)*tan(phi1);
        zu(i) = (yu(i)-yToe)*tanBeta;
        syms x y
        eqToe1 = (y-yToe)^2-2*(y1-yToe)*(y-yToe)+(1-B1^2)*(y1-yToe)^2+(x-x1)^2==0;
        eqToe2 = (y-yToe)^2-2*(y2-yToe)*(y-yToe)+(1-B2^2)*(y2-yToe)^2+(x-x2)^2==0;
        [solxd,solyd] = vpasolve([eqToe1,eqToe2],[x,y],[xMin,xMax;-inf,yu(i)]);
        xd(i) = double(solxd);
        yd(i) = double(solyd);
        if z1<z2
            eqWeld1 = (z10-z2)-S1*sqrt((x-x10).^2+(y-y10).^2)+S2*sqrt((x-x2).^2+(y-y2).^2)==0;
            eqVisi1 = (x-x10)*(y1-y10)-(x1-x10)*(y-y10)==0;
            [solx,soly] = vpasolve([eqWeld1,eqVisi1],[x,y],[xMin,xMax;-inf,yu(i)]);
        else
            eqWeld2 = (z1-z20)-S1*sqrt((x-x1).^2+(y-y1).^2)+S2*sqrt((x-x20).^2+(y-y20).^2)==0;
            eqVisi2 = (x-x20)*(y2-y20)-(x2-x20)*(y-y20)==0;
            [solx,soly] = vpasolve([eqWeld2,eqVisi2],[x,y],[xMin,xMax;-inf,yu(i)]);
        end
        if isempty(soly)
            yWeldline(i,:) = linspace(yu(i),yd(i),length(yWeldline));
            n1 = length(yWeldline);
        elseif double(soly)<yd(i)
            yWeldline(i,:) = linspace(yu(i),yd(i),length(yWeldline));
            n1 = length(yWeldline);
        elseif z1<z2
            eqToe1 = (y-yToe)^2-2*(y10-yToe)*(y-yToe)+(y10-yToe)^2-z10^2/S1^2+(x-x10)^2==0;
            [solxd,solyd] = vpasolve([eqToe1,eqToe2],[x,y],[xMin,xMax;-inf,yu(i)]);
            n1 = round(length(yWeldline)*(double(soly)-yu(i))/(double(solyd)-yu(i)));
            yWeldline(i,1:n1) = linspace(yu(i),double(soly),n1);
            yWeldline(i,n1+1:end) = linspace(double(soly)+yWeldline(i,2)-yWeldline(i,1),double(solyd),length(yWeldline)-n1);
            plot([x10,x1,double(solx)],[y10,y1,double(soly)],'r.')
            plot(double(solxd),double(solyd),'r.')
        else
            eqToe2 = (y-yToe)^2-2*(y20-yToe)*(y-yToe)+(y20-yToe)^2-z20^2/S2^2+(x-x20)^2==0;
            [solxd,solyd] = vpasolve([eqToe1,eqToe2],[x,y],[xMin,xMax;-inf,yu(i)]);
            n1 = round(length(yWeldline)*(double(soly)-yu(i))/(double(solyd)-yu(i)));
            yWeldline(i,1:n1) = linspace(yu(i),double(soly),n1);
            yWeldline(i,n1+1:end) = linspace(double(soly)+yWeldline(i,2)-yWeldline(i,1),double(solyd),length(yWeldline)-n1);
            plot([x20,x2,double(solx)],[y20,y2,double(soly)],'r.')
            plot(double(solxd),double(solyd),'r.')
        end
        for j = 1:length(yWeldline)
            yt = yWeldline(i,j);
            if j==n1+1 && z1<z2
                x1=x10; y1=y10; z1=z10;
            elseif j==n1+1
                x2=x20; y2=y20; z2=z20;
            end
            S = vpasolve((z1-z2)-S1*sqrt((x-x1).^2+(yt-y1).^2)+S2*sqrt((x-x2).^2+(yt-y2).^2)==0,x,[xMin,xMax]);
            xWeldline(i,j) =  double(S);
            zWeldline(i,j) = z1-S1*sqrt((xWeldline(i,j)-x1).^2+(yWeldline(i,j)-y1).^2);
        end
    end
    
    % comparison
    figure
    axis equal
    hold on
    axis off
    zFan_visi(isnan(zFan_visi))=zRock(isnan(zFan_visi));
    pcolor(xMesh,yMesh,kFan)
    c = colorcube(64);
    c = c([20;24;35;26;38;23],:);
    colormap(c)
    caxis([1 6])
    shading flat
    axis([xMin+5 xMax yMin yMax zMin zMax])
    contour(xMesh,yMesh,zFan_visi,0:2:zMax,'-','color',[0.5,0.5,0.5])
    axis equal
    hold on
    axis off
    plot(xWeldline',yWeldline','k-','linewidth',1)
end