📄 dandiaohu.m
字号:
%function [P0x,P1x,x1,y1,r0,I,j,G,H,R]=barc(P0x,P1x,k0) %P0x=-1;P1x=2;k0=1.25;
close all
clear all
clc
P0x=0;
P1x=2;
k0=1.54;
h=0.01;
u1=0:h:2;
u=P0x:h:P1x; %
U=u.^2;
U=u1.^2;
plot(u1,U,'k');axis equal %原曲线
hold;
P0y=P0x^2;
P1y=P1x^2;
plot(P0x,P0y,'o',P1x,P1y,'o') %P0 P1
t=linspace(P0x,P1x);
T=(P0y+P1y)/2-(P1x-P0x)/(P1y-P0y)*(t-(P0x+P1x)/2);
%plot(t,T); axis equal %P0 P1中垂线
z=diff(U)/h;
v0x=cos(atan(z(1)));v0y=sin(atan(z(1)));
v1x=cos(atan(z(int32(abs((P1x-P0x)/h)))));v1y=sin(atan(z(int32(abs((P1x-P0x)/h)))));
V0x=P0x+v0x;V0y=P0y+v0y;
V1x=P1x+v1x;V1y=P1y+v1y;
plot(V0x,V0y,'*'); %单位向量v0端点
%plot(V1x,V1y,'*'); %单位向量v1端点
w=linspace(-1,2); %
W=(V0y+V1y)/2-(V1x-V0x)/(V1y-V0y)*(w-(V0x+V1x)/2);
%plot(w,W); %P0+v0 P1+v1中垂线
fun1=['(',num2str(V0y),'+',num2str(V1y),')/2-(',num2str(V1x),'-(',num2str(V0x),'))/(',num2str(V1y),'-',num2str(V0y),')*(e-(',num2str(V0x),'+',num2str(V1x),')/2)=f'];
fun0=['(',num2str(P0y),'+',num2str(P1y),')/2-(',num2str(P1x),'-(',num2str(P0x),'))/(',num2str(P1y),'-',num2str(P0y),')*(e-(',num2str(P0x),'+',num2str(P1x),')/2)=f'];
[e,f]=solve(fun1,fun0,'e,f');
%plot(double(e),double(f),'+'); %虚线圆圆心
r=sqrt((e-P0x)^2+(f-P0y)^2);%虚线圆半径
o=0:pi/100:2*pi;
x=e+r*cos(o);y=f+r*sin(o);
%plot(double(x),double(y),'--'); %画虚线圆
%k0=4; %小圆曲率
r0=1/k0;
x1=P0x-r0*v0y;y1=P0y+r0*v0x;
%plot(x1,y1,'+') %小圆圆心
fun2=['(g-(',num2str(double(e)),'))^2+(h-',num2str(double(f)),')^2-',num2str(double(r)),'^2=0'];
[g,h]=solve(fun2,['(g-(',num2str(double(x1)),'))^2+(h-',num2str(double(y1)),')^2-',num2str(double(r0)),'^2=0'],'g,h');
for i=1:length(g)
if (abs(double(g(i))-P0x)>0.01)&&(abs(double(h(i))-P0y)>0.01),
plot(double(g(i)),double(h(i)),'*');%小圆和虚线圆求交得两个交点剔掉一个
G=double(g(i));H=double(h(i));
end
end
O=3*pi/2-atan((x1-P0x)/(y1-P0y)):pi/100:3*pi/2+atan((G-x1)/(y1-H)); % ？
x2=x1+r0*cos(O);y2=y1+r0*sin(O);
plot(x2,y2); %弧a0
x4=linspace(-2,1);y4=(H-y1)/(G-x1)*(x4-G)+H;
%plot(x4,y4) %线(J c0)
x5=-0.5:0.1:2;y5=P1y-v1x/v1y*(x5-P1x);
%plot(x5,y5) %线（P1，P1+（v1垂直））
fun3=['j=(',num2str(H),'-',num2str(double(y1)),')/(',num2str(G),'-(',num2str(double(x1)),'))*(I-(',num2str(G),'))+',num2str(H),''];
[I,j]=solve(['j=',num2str(P1y),'-',num2str(v1x),'/',num2str(v1y),'*(I-',num2str(P1x),')'],fun3,'I,j');
%plot(double(I),double(j),'o') %弧a1的圆心（大圆)
O1=3*pi/2+atan((G-x1)/(y1-H)):pi/100:pi*2-atan((j-P1y)/(P1x-I));
R=sqrt((I-P1x)^2+(j-P1y)^2); %大圆半径
X=I+R*cos(O1);Y=j+R*sin(O1);
plot(double(X),double(Y)); %弧a1
fun4=['Jy=(',num2str(H),'-',num2str(double(y1)),')/(',num2str(G),'-(',num2str(double(x1)),'))*(Jx-',num2str(G),')+',num2str(H),''];
[Jx,Jy]=solve(fun4,['(Jx)^2-Jy=0'],'Jx,Jy');
for i=1:length(Jx)
if (double(Jx(i))>0)&&(double(Jx(i))<2)
plot(double(Jx(i)),double(Jy(i)),'d');%线（J C0）与原曲线求交得三个交点 ？ Jq
JX=double(Jx(i));
end
end
%求弧a0与原曲线的最大error
x7=P0x:JX/20:JX;
D0=0;
for i=1:length(x7)
y7(i)=x7(i).^2;
D0=max(D0,abs((sqrt((x7(i)-x1)^2+(y7(i)-y1)^2)-r0)));
end
%求弧a1与原曲线的最大error
x8=JX:(P1x-JX)/20:P1x;
D1=0;
for i=1:length(x8)
y8(i)=x8(i).^2;
D1=max(D1,abs((sqrt((x8(i)-double(I))^2+(y8(i)-double(j))^2)-double(R))));
end
if D0>0.8
P0x=P0x;
P1x=G;
P0y=P0y;
P1y=H;
h=0.01;
U=u.^2;
k0=1.85;
plot(P0x,P0y,'o',P1x,P1y,'o') %P0 P1
t=linspace(P0x,P1x);
T=(P0y+P1y)/2-(P1x-P0x)/(P1y-P0y)*(t-(P0x+P1x)/2);
%plot(t,T); axis equal %P0 P1中垂线
z=diff(U)/h;
v0x=cos(atan(z(1)));v0y=sin(atan(z(1)));
v1x=cos(atan((G-x1)/(y1-H)));v1y=sin(atan((G-x1)/(y1-H)));
V0x=P0x+v0x;V0y=P0y+v0y;
V1x=P1x+v1x;V1y=P1y+v1y;
plot(V0x,V0y,'*'); %单位向量v0端点
plot(V1x,V1y,'*'); %单位向量v1端点
w=linspace(-1,2); %
W=(V0y+V1y)/2-(V1x-V0x)/(V1y-V0y)*(w-(V0x+V1x)/2);
%plot(w,W); %P0+v0 P1+v1中垂线
fun1=['(',num2str(V0y),'+',num2str(V1y),')/2-(',num2str(V1x),'-(',num2str(V0x),'))/(',num2str(V1y),'-',num2str(V0y),')*(e-(',num2str(V0x),'+',num2str(V1x),')/2)=f'];
fun0=['(',num2str(P0y),'+',num2str(P1y),')/2-(',num2str(P1x),'-(',num2str(P0x),'))/(',num2str(P1y),'-',num2str(P0y),')*(e-(',num2str(P0x),'+',num2str(P1x),')/2)=f'];
[e,f]=solve(fun1,fun0,'e,f');
plot(double(e),double(f),'+'); %虚线圆圆心
r=sqrt((e-P0x)^2+(f-P0y)^2);%虚线圆半径
o=0:pi/100:2*pi;
x=e+r*cos(o);y=f+r*sin(o);
plot(double(x),double(y),'--'); %画虚线圆
%k0=4; %小圆曲率
r0=1/k0;
x1=P0x-r0*v0y;y1=P0y+r0*v0x;
plot(x1,y1,'+') %小圆圆心
fun2=['(g-(',num2str(double(e)),'))^2+(h-',num2str(double(f)),')^2-',num2str(double(r)),'^2=0'];
[g,h]=solve(fun2,['(g-(',num2str(double(x1)),'))^2+(h-',num2str(double(y1)),')^2-',num2str(double(r0)),'^2=0'],'g,h');
for i=1:length(g)
if (abs(double(g(i))-P0x)>0.01)&&(abs(double(h(i))-P0y)>0.01),
plot(double(g(i)),double(h(i)),'*');%小圆和虚线圆求交得两个交点剔掉一个
G=double(g(i));H=double(h(i));
end
end
O=3*pi/2-atan((x1-P0x)/(y1-P0y)):pi/100:3*pi/2+atan((G-x1)/(y1-H)); % ？
x2=x1+r0*cos(O);y2=y1+r0*sin(O);
plot(x2,y2); %弧a0
x4=linspace(-2,1);y4=(H-y1)/(G-x1)*(x4-G)+H;
%plot(x4,y4) %线(J c0)
x5=-0.5:0.1:2;y5=P1y-v1x/v1y*(x5-P1x);
%plot(x5,y5) %线（P1，P1+（v1垂直））
fun3=['j=(',num2str(H),'-',num2str(double(y1)),')/(',num2str(G),'-(',num2str(double(x1)),'))*(I-(',num2str(G),'))+',num2str(H),''];
[I,j]=solve(['j=',num2str(P1y),'-',num2str(v1x),'/',num2str(v1y),'*(I-',num2str(P1x),')'],fun3,'I,j');
%plot(double(I),double(j),'o') %弧a1的圆心（大圆)
O1=3*pi/2+atan((G-x1)/(y1-H)):pi/100:pi*2-atan((j-P1y)/(P1x-I));
R=sqrt((I-P1x)^2+(j-P1y)^2); %大圆半径
X=I+R*cos(O1);Y=j+R*sin(O1);
plot(double(X),double(Y)); %弧a1
end
fun4=['Jy=(',num2str(H),'-',num2str(double(y1)),')/(',num2str(G),'-(',num2str(double(x1)),'))*(Jx-',num2str(G),')+',num2str(H),''];
[Jx,Jy]=solve(fun4,['(Jx)^2-Jy=0'],'Jx,Jy');
for i=1:length(Jx)
if (double(Jx(i))>0)&&(double(Jx(i))<2)
plot(double(Jx(i)),double(Jy(i)),'d');%线（J C0）与原曲线求交得三个交点 ？ Jq
JX=double(Jx(i));
end
end
%求弧a0与原曲线的最大error
x7=P0x:JX/20:JX;
D0=0;
for i=1:length(x7)
y7(i)=x7(i).^2;
D0=max(D0,abs((sqrt((x7(i)-x1)^2+(y7(i)-y1)^2)-r0)));
end
%求弧a1与原曲线的最大error
x8=JX:(P1x-JX)/20:P1x;
D1=0;
for i=1:length(x8)
y8(i)=x8(i).^2;
D1=max(D1,abs((sqrt((x8(i)-double(I))^2+(y8(i)-double(j))^2)-double(R))));
end
if D1>0.05
P0x=G;
P1x=P1x;
P0y=H;
P1y=P1y;
h=0.01;
U=u.^2;
k0=0.4;
plot(P0x,P0y,'o',P1x,P1y,'o') %P0 P1
t=linspace(-1.5,1);
T=(P0y+P1y)/2-(P1x-P0x)/(P1y-P0y)*(t-(P0x+P1x)/2);
%plot(t,T); axis equal %P0 P1中垂线
z=diff(U)/h;
v0x=cos(atan(z(int32((JX)/h))));v0y=sin(atan(z(int32((JX)/h))));
v1x=cos(atan(z(int32((P1x)/h))));v1y=sin(atan(z(int32((P1x)/h))));
V0x=P0x+v0x;V0y=P0y+v0y;
V1x=P1x+v1x;V1y=P1y+v1y;
plot(V0x,V0y,'*'); %单位向量v0端点
plot(V1x,V1y,'*'); %单位向量v1端点
w=linspace(-1.5,2); %
W=(V0y+V1y)/2-(V1x-V0x)/(V1y-V0y)*(w-(V0x+V1x)/2);
%plot(w,W); %P0+v0 P1+v1中垂线
fun1=['(',num2str(V0y),'+',num2str(V1y),')/2-(',num2str(V1x),'-(',num2str(V0x),'))/(',num2str(V1y),'-',num2str(V0y),')*(e-(',num2str(V0x),'+',num2str(V1x),')/2)=f'];
fun0=['(',num2str(P0y),'+',num2str(P1y),')/2-(',num2str(P1x),'-(',num2str(P0x),'))/(',num2str(P1y),'-',num2str(P0y),')*(e-(',num2str(P0x),'+',num2str(P1x),')/2)=f'];
[e,f]=solve(fun1,fun0,'e,f');
%plot(double(e),double(f),'+'); %虚线圆圆心
r=sqrt((e-P0x)^2+(f-P0y)^2);%虚线圆半径
o=0:pi/100:2*pi;
x=e+r*cos(o);y=f+r*sin(o);
%plot(double(x),double(y),'--'); %画虚线圆
%k0=4; %小圆曲率
r0=1/k0;
x1=P0x-r0*v0y;y1=P0y+r0*v0x;
%plot(x1,y1,'+') %小圆圆心
fun2=['(g-(',num2str(double(e)),'))^2+(h-(',num2str(double(f)),'))^2-',num2str(double(r)),'^2=0'];
[g,h]=solve(fun2,['(g-(',num2str(double(x1)),'))^2+(h-',num2str(double(y1)),')^2-',num2str(double(r0)),'^2=0'],'g,h');
for i=1:length(g)
if (abs(double(g(i))-P0x)>0.01)&&(abs(double(h(i))-P0y)>0.01),
plot(double(g(i)),double(h(i)),'*');%小圆和虚线圆求交得两个交点剔掉一个
G=double(g(i));H=double(h(i));
end
end
O=3*pi/2-atan((x1-P0x)/(y1-P0y)):pi/100:3*pi/2+atan((G-x1)/(y1-H)); % ？
x2=x1+r0*cos(O);y2=y1+r0*sin(O);
plot(x2,y2); %弧a0
x4=linspace(-2,1);y4=(H-y1)/(G-x1)*(x4-G)+H;
%plot(x4,y4) %线(J c0)
x5=-0.5:0.1:2;y5=P1y-v1x/v1y*(x5-P1x);
%plot(x5,y5) %线（P1，P1+（v1垂直））
fun3=['j=(',num2str(H),'-',num2str(double(y1)),')/(',num2str(G),'-(',num2str(double(x1)),'))*(I-(',num2str(G),'))+',num2str(H),''];
[I,j]=solve(['j=',num2str(P1y),'-',num2str(v1x),'/',num2str(v1y),'*(I-',num2str(P1x),')'],fun3,'I,j');
%plot(double(I),double(j),'o') %弧a1的圆心（大圆)
O1=3*pi/2+atan((G-x1)/(y1-H)):pi/100:pi*2-atan((j-P1y)/(P1x-I));
R=sqrt((I-P1x)^2+(j-P1y)^2); %大圆半径
X=I+R*cos(O1);Y=j+R*sin(O1);
plot(double(X),double(Y)); %弧a1
end
fun4=['Jy=(',num2str(H),'-',num2str(double(y1)),')/(',num2str(G),'-(',num2str(double(x1)),'))*(Jx-',num2str(G),')+',num2str(H),''];
[Jx,Jy]=solve(fun4,['(Jx)^2-Jy=0'],'Jx,Jy');
for i=1:length(Jx)
if (double(Jx(i))>0)&&(double(Jx(i))<2)
plot(double(Jx(i)),double(Jy(i)),'d');%线（J C0）与原曲线求交得三个交点 ？ Jq
JX=double(Jx(i));
end
end
%求弧a0与原曲线的最大error
x7=P0x:JX/40:JX;
D0=0;
for i=1:length(x7)
y7(i)=x7(i).^2;
D0=max(D0,abs((sqrt((x7(i)-x1)^2+(y7(i)-y1)^2)-r0)));
end
%求弧a1与原曲线的最大error
x8=JX:(P1x-JX)/40:P1x;
D1=0;
for i=1:length(x8)
y8(i)=x8(i).^2;
D1=max(D1,abs((sqrt((x8(i)-double(I))^2+(y8(i)-double(j))^2)-double(R))));
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -