%% base on GO ray tracing method,calculate cassegrain ray-path and reflect
%%ray,writing on 2014,10,9
%unit m,
%---------------------start------------------------------------------------
clear
clc
% define parameter
tic
%base parameter
global ideg dm mc fm
ideg=pi/180; %弧度与角度的比例系数
dm=13.7; %main reflector diameter,口径
mc=11; %magnification ratio,放大率，副面定点到卡焦点与到主焦点的距离比
%distance from subreflector vertex to
%secondary-focus divide distance from subreflector vertex to primary-focus
i_fm=0.3708; %focal diameter ratio，焦径比
fm=dm*i_fm; %main focal length，焦距
ds=1; %subreflector diameter 副面口径
h_svpf=0.3852; %distance from subreflector vertex to primary-focus 副面顶点到主焦点距离
h_svcf=mc*h_svpf; %副面顶点到卡焦点的距离
h_pvcf=fm-h_svpf-h_svcf;%主面顶点到卡焦点距离
%variable parameter
p0=[0 0 h_pvcf]; %馈源相位中心点坐标
gamma0=10; %发出射线与Z轴夹角
phi0=[90 90-gamma0 gamma0]*ideg; %方向余弦角
I0=cos(phi0); %入射射线方向余弦
l_0=I0(1);
m_0=I0(2);
n_0=I0(3);
% ray tracing-----------------------
% first reflect ------
x0=p0(1);
y0=p0(2);
z0=p0(3);
syms y1
[z1 nl_1 nm_1 nn_1]=fsubref((y1-y0)*l_0/m_0+x0,y1); %返回z1的符号变量形式
f11=z1-z0-(y1-y0)*n_0/m_0; %求解y1的方程表达式，以字符变量的形式
f1=strcat(char(f11),'=0'); %转化成字符串等式形式
y1=eval(solve(f1)); %求解y1值,,,,,,??????
x1=(y1-y0)*l_0/m_0+x0;
z1=eval(z1);
% nl_1=eval(nl_1);
% nm_1=eval(nm_1);
% nn_1=eval(nn_1);
nl_1=-1*eval(nl_1);
nm_1=-1*eval(nm_1);
nn_1=-1*eval(nn_1);
s1=2*(l_0*nl_1+m_0*nm_1+n_0*nn_1);
I1=[l_0-s1*nl_1 m_0-s1*nm_1 n_0-s1*nn_1]; %射线L1的方向余弦I1
l_1=I1(1);
m_1=I1(2);
n_1=I1(3);
% second reflect ------
syms y2
[z2 nl_2 nm_2 nn_2]=fprimref((y2-y1)*l_1/m_1+x1,y2); %返回z2的符号变量形式
f21=z2-z1-(y2-y1)*n_1/m_1; %求解y2的方程表达式，以字符变量的形式
f2=strcat(char(f21),'=0'); %转化成字符串等式形式
y2=eval(solve(f2)); %求解y2值
y2=y2(1);
x2=(y2-y1)*l_1/m_1+x1;
z2=eval(z2);
nl_2=eval(nl_2);
nm_2=eval(nm_2);
nn_2=eval(nn_2);
s2=2*(l_1*nl_2+m_1*nm_2+n_1*nn_2);
I2=[l_1-s2*nl_2 m_1-s2*nm_2 n_1-s2*nn_2]; %射线L2的方向余弦I2
l_2=I2(1);
m_2=I2(2);
n_2=I2(3);
%apeture plate point
syms y3
z3=fap((y3-y2)*l_2/m_2+x2,y3);
f31=z2-z2-(y3-y2)*n_2/m_2; %求解y3的方程表达式，以字符变量的形式
f3=strcat(char(f31),'=0'); %转化成字符串等式形式
y3=eval(solve(f3)); %求解y3值
y3=y3(1);
x3=(y3-y2)*l_2/m_2+x2;
z3=eval(z3);
%plot ray tracing 绘制
% plot3([x0 x1 x2 x3],[y0 y1 y2 y3],[z0 z1 z2 z3])
% plot3([x0 x1],[y0 y1],[z0 z1],'r')
% hold on
% plot3([x1 x2],[y1 y2],[z1 z2],'g')
% hold on
% plot3([x2 x3],[y2 y3],[z2 z3],'b')
plot([z0 z1],[y0 y1],'r')
hold on
plot([z1 z2],[y1 y2],'g')
hold on
plot([z2 z3],[y2 y3],'b')
plot([0 nn_1],[0 nm_1]),hold on
% time(S_time)
toc
t=toc
评论0