%先输入为被旋转曲线的参数方程
%再输入转轴的标准方程
close all;clear all;clc;
disp('请输入首先被旋转曲线'); %输入提示信息
prompt1='请输入x=f(t)，例如x=2t+3，则输入数组[2 3]nx=';
f=input(prompt1);
prompt2='请输入y=g(t)，例如y=4t+5，则输入数组[4 5]ny=';
g=input(prompt2);
prompt3='请输入z=h(t)，例如y=6t+7，则输入数组[6 7]nz=';
h=input(prompt3);
disp(['您的被旋转曲线方程是：']);
disp(['x=(' num2str(f(1)) ')×t+(' num2str(f(2)) ')']);
disp(['y=(' num2str(g(1)) ')×t+(' num2str(g(2)) ')']);
disp(['z=(' num2str(h(1)) ')×t+(' num2str(h(2)) ')']);
disp('请输入首先旋转轴的直线标准方程：(x-x0)/m=(y-y0)/n=(z-z0)/p'); %输入提示信息
prompt4='请输入数组[m n p]n[m n p]=';
mnp=input(prompt4);
prompt5='请输入数组[x0 y0 z0]n[x0 y0 z0]=';
xyz=input(prompt5);
disp(['您的旋转轴直线方程是：']);
disp(['(x-(' num2str(xyz(1)) '))/(' num2str(mnp(1)) ')=']);
fprintf("%c",8) %删掉换行符号
disp(['(y-(' num2str(xyz(2)) '))/(' num2str(mnp(2)) ')=']);
fprintf("%c",8) %删掉换行符号
disp(['(z-(' num2str(xyz(3)) '))/(' num2str(mnp(3)) ')']);
F=[mnp(1)*f(1)+mnp(2)*g(1)+mnp(3)*h(1),mnp(2)*f(2)+mnp(2)*g(2)+mnp(3)*h(2)];
disp(['F=(' num2str(F(1)) ')t+(' num2str(F(2)) ')']) %调试用
if(F(1)==0)disp(['两直线垂直，方程退化，不予输出']);return;
end
F_1=[1/F(1) -F(2)/F(1)]; %F的逆为(1/a)t+(-b/a)
disp(['F^(-1)=(' num2str(F_1(1)) ')t+(' num2str(F_1(2)) ')']) %调试用
syms x y z;
F_1u=F_1(1)*(mnp(1)*x+mnp(2)*y+mnp(3)*z); %F^(-1)(mx+ny+pz)
Gamma=(x-xyz(1))^2+(y-xyz(2))^2+(z-xyz(3))^2 ...
-(f(1)*F_1u+f(2)-xyz(1))^2-(g(1)*F_1u+g(2)-xyz(2))^2-(h(1)*F_1u+h(2)-xyz(3))^2;
Gamma=expand(Gamma);
disp(['旋转曲面方程为：' char(Gamma) '=0']);
prompt6='请输入绘图区域[x0 x1 y0 y1 z0 z1]，如输入0，则为默认区间n';
plotrange=input(prompt6);
if length(plotrange)==6fimplicit3(Gamma,plotrange); %绘图
elsefimplicit3(Gamma); %绘图
end
xlabel('x');
ylabel('y');

%测试样例1：圆锥
[1 0]
[0 0]
[1 0]
[0 0 1]
[0 0 0]
0

%测试样例2：奇怪的图形
[2 3]
[4 5]
[6 7]
[1 2 3]
[4 5 6]
[-15 15 -15 15 0 20]