function line_inter_3d(xe, ye, ze, fc, fs, fe, at)
%3d line interpolate
%input x end,y end,z end,feed velocity,speed rising time
%draw the interpolate points of this line
%INT16 xe, int16 ye;
%unit16 fc, uint16 fs, uint16 fe;
%uint16 a_t;

%變量聲明
int16 xe;
int16 ye;
int16 ze;
int16 fc;
int16 fs;
int16 fe;
int16 at;
int16 fi;
fi = int16(0);
double xi;
double yi;
double zi;
double dx;
double dy;
double dz;
char x_dir;
char y_dir;
char z_dir;
double L;
double La;
double Le;
double Lsum;
double Lxy;
double dL;
double dLxy;     %XY平面周期進給量
double cos_z;    %L與XY平面夾角余弦值
double sin_z;    %L與XY平面夾角正弦值
double tan_xy;   %Y與X夾角正切值
double sec_xy;   %Y與X夾角余弦值
int16 cParam;

%變量初始化
xi = double(0);
yi = double(0);
zi = double(0);
acc_flag = char(0);
ac_step = char(0);
end_step = char(0);
x_dir = char(0);
y_dir = char(0);
z_dir = char(0);
blnRun = char(1);
dx = double(0);
dy = double(0);
dz = double(0);
dL = double(0);
cos_z = double(0);
sin_z = double(0);
tan_xy = double(0);
sec_xy = double(0);
L = double(0);
La = double(0);
Le = double(0);
Lsum = double(0);
Lxy = double(0);
temp = int32(0);
acc_step = int16(0);
cParam = int16(125);

%設置繪圖模式為“添加”
hold on
grid on
box on

%判定x進給方向
if (xe < 0)
    xe = -xe;
    x_dir = 1;
end

%判定y進給方向
if (ye < 0)
    ye = -ye;
    y_dir = 1;
end

%判定y進給方向
if (ze < 0)
    ze = -ze;
    z_dir = 1;
end

%若本段無進給
if(xe==0 && ye==0 && Ze==0)
    blnRun = 0;
    return;
end

%求Lxy長度
Lxy = sqrt(double(xe*xe + ye*ye));
L   = sqrt(double(xe*xe + ye*ye + ze*ze));
if(ze ~= 0)
    %求L與XY平面夾角cos值
    cos_z = Lxy/L;

    %求L與XY平面夾角sin值
    sin_z = sqrt(double(1 - cos_z*cos_z));

    if(xe ~= 0)
        %計算tan值
        tan_xy  = ye/xe;
        %計算sec值
        sec_xy  = sqrt(double(tan_xy*tan_xy + 1));
    end
end

%速度預處理段,采用梯形加減速控制
if(fs < fc)
    %加速過程
    acc_flag = 0;
    ac = int16((fc - fs)/at);
elseif (fs == fc)
    %勻速過程
    acc_flag = 1;
    if(fe < fc)
        ac = int16((fc - fe)/at);
    end
elseif(fs > fc)
    %減速過程
    acc_flag = 2;
    ac = int16((fs - fc)/at);
end

%為了與DSP上插補時間一致，此處乘以時間系數
ac = ac*125;

while (blnRun)
    switch(acc_flag)
        case 0
            %加速段
            %梯形加速控制，加速度值不變
            %速度呈線性變化
            fi = fi + ac;
            Le = L - Lsum;

            if(fi >= fc)
                %升速階段完成，置恒速運行
                acc_flag = 1;
                La = Lsum + dL;
                fi = fc;
            end

            if(Le <= La)
                %越過加速區，進入減速區
                acc_flag = 2;
            end

            dL = fi/cParam;
            acc_step = acc_step + 1;
        case 1
            %勻速段
            Le = L - Lsum;

            if(Le <= La)
                acc_flag = 2;
            end
        case 2
            %減速段
            fi = fi - ac;
            dL = fi/cParam;
            acc_step = acc_step - 1;

            if(ac_step <= 1)
                acc_flag = 3;
                %if(dL == 0)
                 %   dL = 1;
                %end
                end_step = int16((La - Le)/dL);
                if(end_step <= 0)
                    blnRun = 0;
                end
            end
        case 3
            %低速段
            end_step = end_step - 1;
            if(end_step <= 0)
                blnRun = 0;
            end
    end
    
    %求qdLxy
    dLxy = dL*cos_z;
    
    if(ze ~= 0)
        dz = dL*sin_z;
    else
        dz = 0;
    end
    if (xe ~= 0)
        dx = int16(double(dLxy)/sec_xy);
        dy = int16(double(dx)*tan_xy);
    else
        dx = 0;
        dy = dLxy;
    end

    %重新計算dL值，減少舍入誤差
    dL = sqrt(double(dx*dx + dy*dy + dz*dz));

    xi = xi + dx;
    yi = yi + dy;
    zi = zi + dz;

    Lsum = Lsum + dL;
    if(Lsum >= L)
        dL = L - (Lsum - dL);
        blnRun  = 0;
    end
    
    if(x_dir == 1)
        dx = -dx;
    end
    if(y_dir == 1)
        dy = -dy;
    end
    if(z_dir == 1)
        dz = -dz;
    end

    double x;
    double y;
    double z;
    x = double(xi);
    y = double(yi);
    z = double(zi);
    for i = 1:1:dx
        x = x + 1;
        y = y + tan_xy;
        z = z + (sec_xy*sin_z/cos_z);
        %line(x, y);
        plot3(x, y, z, 'r*');
    end
end
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
axis([0, xe, 0, ye, 0, ze]);