MATLAB---绘制Ferguson曲线

function FergusonCurve()
%本程序的功能是构造一个用的Ferguson曲线
N = 9;
SamplPs = CollectSPFergusonCir(N);%采集型值点
t0 = [0,1/3];
tn = [0,-1/3];%规定端点条件
SamplPTs = ThreeTangEq(SamplPs,t0,tn);%为每个型值点计算切矢量
RenderFergusonCurve(SamplPs,SamplPTs);%绘制Ferguson曲线

function SamplPs = CollectSPFergusonCir(N)
%N表示样点个数

step = pi/N;%为编程方便,把N作为间隔数
xita = 0:step:pi;%等间距取角度向量
%根据参数方程计算圆周上的点
SamplPs(1,:) = cos(xita);%x的坐标放在矩阵的第一行
SamplPs(2,:) = sin(xita);%y的坐标放在矩阵的第二行
%根据参数方程计算圆周上的点
% plot(SamplPs(1,:),SamplPs(2,:),'k*')%测试程序正确性用的
% axis equal;

function RenderFergusonCurve(SamplPs,SamplPTs)
%输入型值点SamplPs和型值点处的切矢量SamplPTs

%step1 汇总各个曲线段上的密化点
N = length(SamplPs(1,:));%得到密化点数目，参见Test3TEq
N = N-1;%每两个点之间有一个曲线段，这就是曲线段的数目
DensePs = [];%定义空矩阵
for i=1:N
    r0 = SamplPs(:,i)';r1 = SamplPs(:,i+1)';
    r0c = SamplPTs(i,:);r1c = SamplPTs(i+1,:);
    Ps = FergCurvSegDenPs(r0,r1,r0c,r1c);%在当前曲线段上构造密化点
    DensePs = [DensePs;Ps];
    %把当前曲线段上采集的密化点数组放入全局数组
end
plot(DensePs(:,1),DensePs(:,2),'linewidth',2)

hold on
plot(SamplPs(1,:),SamplPs(2,:),'k*')%添加型值点便于观察
N=N+1;
for i=1:N%逐个绘制切线段
    x = [SamplPs(1,i),SamplPs(1,i) + SamplPTs(i,1)];
    y = [SamplPs(2,i),SamplPs(2,i) + SamplPTs(i,2)];
    %每个切线段的起始点是[SamplPs(1,i),SamplPs(2,i)]
    %终了点是[SamplPs(1,i) + SamplPTs(i,1),SamplPs(2,i) + SamplPTs(i,2)]
    plot(x,y,'r','linewidth',1);
    %把切线绘制为红色,宽度为2
end
hold off
axis equal

function FergusonCurv2()

N = 12;
SamplPs = CollectSPFergusonCi2(N);%采集型值点
t0 = [0,1/3];
tn = [0,-1/3];%规定端点条件
SamplPTs = ThreeTangEq(SamplPs,t0,tn);%为每个型值点计算切矢量
RenderFergusonCurv2(SamplPs,SamplPTs);%绘制Ferguson曲线

function RenderFergusonCurv2(SamplPs,SamplPTs)
%输入型值点SamplPs和型值点处的切矢量SamplPTs

%step1 汇总各个曲线段上的密化点
N = length(SamplPs(1,:));%得到密化点数目，参见Test3TEq
N = N-1;%每两个点之间有一个曲线段，这就是曲线段的数目
DensePs = [];%定义空矩阵
for i=1:N
    r0 = SamplPs(:,i)';r1 = SamplPs(:,i+1)';
    r0c = SamplPTs(i,:);r1c = SamplPTs(i+1,:);
    Ps = FergCurvSegDenPs(r0,r1,r0c,r1c);%在当前曲线段上构造密化点
    DensePs = [DensePs;Ps];
    %把当前曲线段上采集的密化点数组放入全局数组
end
plot(DensePs(:,1),DensePs(:,2),'linewidth',2)

hold on
plot(SamplPs(1,:),SamplPs(2,:),'k*')%添加型值点便于观察
hold off
axis equal

function ts = ThreeTangEq(SamplPs,t0,tn)
%Sampls输入型值点;t0,tn是首末点的切矢量

N = length(SamplPs(1,:));
%SamPls(1,:)表示取矩阵Sampls第一行的所有元素
%注意矩阵一行或者一列的全部元素形成一个向量
%length()就是测一个向量的长度,也就是元素个数
%注意矩阵SamPs存储采集的型值点,所有型值点的x坐标放在第一行
%参见本文件夹中的函数SamplPs = CollectSPFergusonCir(N)
%因此N = length(SamPls(1,:))得到的就是本程序使用的型值点的个数
CoefA = eye(N,N);%建立一个N*N的单位矩阵
%对这个矩阵进行更改后就是系数矩阵
%更改第2行到第N-1行
for i=2:N-1
    CoefA(i,i-1)=1;CoefA(i,i)=4;CoefA(i,i+1)=1;
end
EqB2D = zeros(N,2);
%注意每一行表示一个向量,共有N行
EqB2D(1,:) = t0;
for i=2:N-1
    EqB2D(i,:) = 3*(SamplPs(:,i+1)'-SamplPs(:,i-1)');
    %V'表示取向量V的转置;SamPls(:,i+1)'就是把列向量表示的点变为行向量表示的点
end
EqB2D(N,:) = tn;
ts=CoefA\EqB2D;%求解线性方程组得到每个点的切矢量