mian_raytracing.rar_RayTracing_ray tracing_几何光学_射线光学_射线追踪图

%% 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