一段matlab绘制透镜的代码，欢迎探讨

clear;close all;  %消除工作变量/关闭窗口
R=100;d=10;%thick of lens   设置透镜初始参数
%plot the lens
lensPlot(R,d);   %调用函数绘制平凸透镜轮廓
%plot the optical rays
ND=5;             %进行调制
%Red light        红光追迹
DL=-15;DH=-5;n=1.515;
[XLmin,XRmax1]=raysPlot(DL,DH,ND,R,d,n,'r');
%Green light      绿光追迹
DL=-5;DH=5;n=1.519;
[XLmin,XRmax2]=raysPlot(DL,DH,ND,R,d,n,'g');
%Blue light        蓝光追迹
DL=5;DH=15;n=1.523;
[XLmin,XRmax2]=raysPlot(DL,DH,ND,R,d,n,'b');
XRmaxA=sort([XRmax1,XRmax2,XRmax3]);        %设置合适的坐标参数
Ymax=sqrt(R^2-(R-d)^2)+10;
axis([XLmin XRmaxA(3)-Ymax Ymax])
grid on;                                    %设置图形坐标及标题
title('simulate the optical ray through the Lens');
function lensPlot(R,d)                     %绘制透镜函数声明
%plot the lens
Hmax=sqrt(R^2-(R-d)^2);%Lens apeture       %绘制平凹透镜的“平”面
Nlens=3;
ylens=linspace(-Hmax,Hmax,Nlens);
xlens=(R-d)*ones(size(ylens));
plot(xlens,ylens);
hold on;
str1=sprintf('sqrt(%f^2-x^2)',R);        %绘制平凹透镜的“凹”面
fplot(str1,[R-d R]);
str2=sprintf('-sqrt(%f^2-x^2)',R);
fplot(str2,[R-d R]);
hold off;
function[XLmin,XRmax]=raysPlot(DL,DH,ND,R,d,n,str)  %光线追迹函数申明
%plot the optical rays through the plane-conves lens
%parameters introduce;
%DL/DH;
%ND:the number of rays
%str:set the colour of rays
if(nargin<6)
    disp('please confirm the number of input parameters!');
    return;
else if nargin==6
        str='b';%set the default colour:blue
else 
         disp('OK!Now Start Plotting....');
end        
Dmax=sqrt(R^2-(R-d+1)^2);%__________1mm
if(abs(DL)>=Dmax|abs(DH)>=Dmax)
    disp('The DL/DH is so big!\n');
    sprintf('The aperture of Lens is %%f\n',Dmax)
    return;
end
hold on;
NP=5;
h=linspace(DL,DH,ND);    %根据输入的DL/DH生成入射光线的高度数组
temp=find(h~=0);         %确定光线追迹的左右边界，激光线的起点和中点横坐标
theta1=asin(h(temp(1))/R);%左右边界数值作为函数输出
theta2=asin(n*h(temp(1))/R);
theta=theta2-theta1;
f=h(temp(1)/tan(theta));
XRmax=R+3/4*f;
XLmin=R-d-30;
for hindex=1:ND         %绘制平行光轴的光线部分
    xinterPoint=sqrt(R^2-H(hindex)^2);
    xleft=linspace(XLmin,xinterPoint,NP);
    yleft=h(hindex)*ones(size(xleft));
    plot(xleft,yleft,str);%plot the left rays
    
%caculate the changed angle for the optical ray
theta1=asin(h(hindex)/R);      %确定会聚光线的斜率
theta2=asin(n*(hindex)/R);
theta=theta2-theta1;
k=tan(theta);
xright=linspace(XinterPoint,XRmax,NP);   %绘制经透镜会聚的光线部分
Yright=K*xright+h(hindex)-k*XinterPoint;
plot(xright,yright,str);%plot the left rays
end
hold off;