该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
syms p q r s t k R
syms u i j a
syms x y w
syms cx cy
syms p1
h=35;
l=11;
x=-30:1:30;
y=-30:1:30;
[X,Y]=meshgrid(x,y);
m=size(X);
Z=ones(m);
P1=ones(m);
P2=ones(m);
j=1;
for x1=-30:1:30
i=1;
for y1=-30:1:30
d=x-h/2+l;
g=(1/2).*(d+((y.^2+h^2/4)./d));
if x1<6.5
w=h/2-l+g+(g.^2-h^2/4).^(1/2);
elseif x1==6.5
w=x;
else
w=h/2-l+g-(g.^2-h^2/4).^(1/2);
end
u=w;
%f=1/84+(1/63-1/84).*((u+11)/35);
%m n取值后
f=1/84+(1/63-1/84).*(70*((u+l)./h).^4-224*((u+l)./h).^5+280*((u+l)./h).^6-160*((u+l)./h).^7+35*((u+l)./h).^8);
%J=(u.*(u + 232))./17640;
%m n取值后
J=(u.*(u.^8 - 81*u.^7 + 1116*u.^6 + 60564*u.^5 - 1318464*u.^4 - 27143424*u.^3 + 562028544*u.^2 + 15897378816*u + 1865020561164))./145921525312500;
theta=asin(J); %theta
p=cos(theta);
%p=simple(p); %cos
n=tan(theta);
%n=simple(n); %tan
%任一点球心坐标
%任一点球心坐标
fx=u-(1./f).*J;
fy=0;
%fz=(1./f).*p+(0.00002951*u.^3+0.006562*u.^2-0.0009396*u+0.0008536);
%Cn=1,拟合曲线取6阶
%fz=(1./f).*p+(0.00000000007519*u.^6+0.000000002433*u.^5+0.0000002793*u.^4+0.00001898*u.^3+0.006576*u.^2-0.000002726*u-0.0000009864);
%m n取值后
fz=(1./f).*p+(2.277e-07*u.^4 + 3.704e-05*u.^3 +0.006494*u.^2 -0.0002022.*u-0.00154);
Z=fz-((1./f).^2-(x-fx).^2-y.^2).^(1/2);
%%%%柱面球面表达
P1=diff(Z,'x');
%q=diff(Z,'y');
%s=diff(p1,'y');
P2=subs(P1,{x,y},{cx,cy});
%q1=subs(q,{x,y},{x1,y1});
%s1=subs(s,{x,y},{x1,y1});
i=i+1;
end
j=j+1;
end