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

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

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.7x=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 1if eval(dfr)~=0Dt=eval((-fr+sqrt(fr^2+2*dfr*de))/dfr);elseDt=eval(de/f);endtb=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==Nbreak;endj=j+1;
end