等精度分割曲线（问题篇）

clc;clear;
syms t t1 Eds

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

hold on
%%绘制分割曲线图案
for tx=0.1:0.01:0.7
   x=cos(tx)*sin(tx/2);
   y=sin(tx)*sin(tx/2);
   plot(x,y,'o'); 
end
hold off

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

x1=10*cos(t1)*sin(t1/2);
y1=10*sin(t1)*sin(t1/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));

dx1=simplify(diff(x1,t1));
dy1=simplify(diff(y1,t1));
ddx1=simplify(diff(dx1,t1));
ddy1=simplify(diff(dy1,t1));

dd=[ddx,ddy];
dd1=[ddx1,ddy1];

f=norm(dd,2)%%2范数
f1=norm(dd1,2)%%2范数
ff=-f;

ft=sqrt((8*Err)/f);%%给的一个初值，由于ddc是按初值给的，可能按这个初值去求dt不合适
%%会产生误差，会在该ti+1和ti之间的误差超差，即得在区间内重新求出最大的c的范数对应tmax，
%%这时tmax为该误差下分割最大的tmax,在重新按误差算出dt与ti相加小于tmax即可，不然只能dt0=dt1
ft1=sqrt((8*Err)/f1);
t0=0.1;
t=t0;
dt=eval(ft);
j=1;
while 1
    
ti=t0+dt;

ta=double(t0);
tb=double(ti);

goalf=inline(ff);
tmax=fminbnd(goalf,ta,tb);%%求函数的最小值，即求曲线的曲率半径最小对应的t值，和其R值
t=tmax;
goalC=eval(f);
goalt=sqrt((8*Err)/goalC);

    if (ta+goalt)>tmax
        sd(j)=tb;
        fd(j)=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),ta,tb);
        t0=tb;
        t=t0;
        ft=sqrt((8*Err)/f);
        dt=eval(ft);
    else
        dt=goalt;
        sd(j)=t0+dt;
        fd(j)=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),ta,tb);
        t0=t0+dt;
        ft=sqrt((8*Err)/f);
        dt=eval(ft);
    end
   
    if abs(sd(max(size(sd)))-0.7)>eps&&j==7
        ft1=sqrt((8*Err)/f1);
        eqn=0.7==t1+ft1;
        tr=solve(eqn,t1)
        Err2=sd(max(size(sd))-1)-tr;
        Err=(1-Err2)*Err;
        t0=0.1;
        t=t0;
        ft=sqrt((8*Err)/f);
        dt=eval(ft);
        if Err2<eps
            break;
        end
    end
    if sd(max(size(sd)))<0.7&&j==7
        break; 
    end
    
    if j==7
        j=0;
        sd=[];
        fd=[]; 
    end
     j=j+1;
end