给定弓高误差，求离散曲线步长

clc;clear;
syms t Eds

%%给定弓高误差，求离散曲线步长
Err=0.01;


x=10*cos(t)*sin(t/2);
y=10*sin(t)*sin(t/2);

L=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),0.1,0.7)%%计算弧长

dx=simplify(diff(x,t));
dy=simplify(diff(y,t));

ddx=simplify(diff(dx,t));
ddy=simplify(diff(dy,t));

dc=[dx,dy];
ddc=[ddx,ddy];
tc=dc.*ddc
Ltc=simplify(sqrt(tc(1)^2+tc(2)^2));
Lc=simplify((sqrt(dc(1)^2+dc(2)^2))^3);
k=simplify(Ltc/Lc)%%以上步骤求曲率


fr=simplify(sqrt(dx^2+dy^2));

dfr=diff(fr,t);

R=1/k;

f=inline(R);
x=fminunc(f,0.1,0.7);%%求函数的最小值，即求曲线的曲率半径最小对应的t值，和其R值
t=x;
Rt=eval(R)%%已知t把符号公式计算出数值解出来其得出曲率R
%%根据弓高误差求弦长步长公式知，得求出最小的曲率半径R，即然后才能得出弦长步长。
hold on

%%绘制分割曲线图案
for t=0.1:0.01:0.7
   x=cos(t)*sin(t/2);
   y=sin(t)*sin(t/2);
   plot(x,y,'o'); 
end
hold off
%%弓高与曲率半径和弦长步长公式，得满足弓高误差前提下弦长步长
de=2*sqrt(2*Rt*Err-Err^2)
%%求分段数
Ni=L/de;
%%取正
N=6;

%%以下验证按该弦长分段出来，算出来的t值，tn到tn-1这段的弧长误差大了，不过弓高误差是满足的。可用前一章知识进行解决。
t0=0.1;
t=t0;
j=1;
ta=t0;
tb=t0;
while 1
    if eval(dfr)~=0
        Dt=eval((-fr+sqrt(fr^2+2*dfr*de))/dfr);
    else
        Dt=eval(de/f);
    end
    tb=tb+Dt;
    t=tb;
    st(j)=tb;
    ft(j)=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),ta,tb);
    ta=tb;
    if j==N
        break;
    end
     j=j+1;
end