curved1_3d_cone2.m

calculating focus ratio degradation

%分析与子午平面不同夹角的平面光锥进行分析;
%会聚光入射的情况,与实际光相符;
clear all;
format long;
pi=3.1415926;
r0=160*10^-6;
ratio=500;
R0=r0*ratio;
n1=1.515;
n2=1.500;
n0=1.000;
xin=input('入射点的x坐标为:xi0=');
yin=input('入射点的y坐标为(convergence >0):yi0=');
zin=input('入射点的z坐标为:zi0=');
phi=input('入射光与子午平面的夹角:phi=');
w1=0;w2=0;

r=linspace(-r0,r0,202);
for tt=1:1:1  %角度变化
    jd=0;
for q=2:1:101    %位置变化 
   px=R0+r(q)*cos(phi/180*pi);
   pz=r(q)*sin(phi/180*pi);
   py=0;
   ab0=[(px-xin)/sqrt((px-xin)^2+(py-yin)^2+(pz-zin)^2),(py-yin)/sqrt((px-xin)^2+(py-yin)^2+(pz-zin)^2),(pz-zin)/sqrt((px-xin)^2+(py-yin)^2+(pz-zin)^2)];
    l=-ab0(1);
    m=-ab0(2);
    n=-ab0(3);
    xi0=-yin/m*l+xin;
    yi0=0;        
    zi0=-yin/m*n+zin;
    lin=[l,m,n];%incident- light directional vector\
    lfin=[0,1,0]; %法线的方向从入射介质到折射介质;
    cai=lin*lfin'; %入射角的余弦
    lt=0;
    thetain=acos(cai);
    cat=sqrt(1-n0^2/n1^2+(n0*cai/n1)^2);
    p=n1*cat-n0*cai;
    % p=n1*cos(thetat)-(n0*cos(thetain));

    thetat=asin(n0*sin(thetain)/n1);
    lt=(lin+p*lfin)/n1; %折射光线的方向向量;
    disp(acos(lt)/pi*180);
    kl=linspace(min(0,(-yi0-yin)/m),max(0,(-yi0-yin)/m),400);
    x=kl.*l+xi0;
    y=kl.*m+yi0;
    z=kl.*n+zi0;
    plot3(x,y,z,'b');
    hold on;
    grid on;

    lin1=lin;lin=lt;
    l=lt(1);m=lt(2);n=lt(3);

    A=l*xi0+m*yi0+n*zi0;
    B=xi0^2+yi0^2+zi0^2+R0^2-r0^2;
    C=4*R0^2*(l^2+m^2);
    D=2*(xi0*l+yi0*m)*4*R0^2;
    E=(xi0^2+yi0^2)*4*R0^2;
    c=[1,4*A,4*A^2+2*B-C,4*A*B-D,B^2-E];
    kj=roots(c);
    xj=kj*l+xi0;
    yj=kj*m+yi0;
    zj=kj*n+zi0;
    jd=[xj,yj,zj];
    dj=(xj-xi0).^2+(yj-yi0).^2+(zj-zi0).^2;
    dirj=[(xj-xi0)/l,(yj-yi0)/m,(zj-zi0)/n];
    disp(size(xj));
    s=0;num=0;dj1=0;
    for j=1:1:length(kj)        %筛选实根,并取沿光线传播方向的实解;
        if isreal(xj(j))==1&isreal(yj(j))==1&isreal(zj(j))==1%&dj(j)~=min(dj)
            if dirj(j)>=0
            num(s+1)=j;     %记录实根的位置序号;
            dj1(s+1)=dj(j); %计算所有实交点到入射点的距离;
            s=s+1;
        end
     end
            
    end      
    if s==0
        disp('there is no real solution');
    break;                     %无实根;退出对该光线的追踪;
    else
        disp(num) 
    for k=1:1:s
     if dj(num(k))==min(dj1)  %实根中光程最小者为真解;
         jdx=xj(num(k));
         jdy=yj(num(k));
         jdz=zj(num(k));
     end
 end
end
jd=[jdx,jdy,jdz];             %确定交点的位置;
disp(jdz);
% if atan(jdy/jdx)>pi/2|atan(jdy/jdx)<0  %判断光线是否溢出所定义的边界 1/4
if jdx<0&jdy<0
    disp('out of fiber1!');
    disp(jd);
    break;
end
 kl=linspace(min(0,(jdx-xi0)/l),max(0,(jdx-xi0)/l),400);
x=kl.*l+xi0;
y=kl.*m+yi0;
z=kl.*n+zi0;
% plot3(x,y,z,'g');


%求反射光线:
c1='x^2+y^2+z^2+R0^2-r0^2-2*R0*sqrt(x^2+y^2)';
dx=diff(c1,'x');
dy=diff(c1,'y');
dz=diff(c1,'z');
lx=subs(dx,{'R0','x','y'},{R0,jd(1),jd(2)});
ly=subs(dy,{'R0','x','y'},{R0,jd(1),jd(2)});
lz=subs(dz,{'z'},{jd(3)});
%法线的方向:
lfx=lx/sqrt(lx^2+ly^2+lz^2);
lfy=ly/sqrt(lx^2+ly^2+lz^2);
lfz=lz/sqrt(lx^2+ly^2+lz^2);
%反射定律求反射光线
lin=[l,m,n];
lf=[lfx,lfy,lfz];
ca=lin*lf';
if acos(ca)<asin(1/n1)  %不满足全反射定律，光线溢出纤芯;
    disp('light is out because of refraction1');
    break;
end
lr=lin-2*ca*lf;
ar=acos(lr)./pi*180;
%反射光线作为下次入射光线的参数：
xi0=jd(1);
yi0=jd(2);
zi0=jd(3);
l=lr(1);
m=lr(2);
n=lr(3);

%%%%%%%%%%%%%%%%%以下求得交点中都要排除入射点(xi0,yi0,zi0)伪解;%%%%%%%%%%%%%%%%%%%%%%%%%
time=1000;  %初始定义光线的反射次数;
for t=1:1:time
    num=0;
    dj1=0;
    dj=0;
    A=l*xi0+m*yi0+n*zi0;
    B=xi0^2+yi0^2+zi0^2+R0^2-r0^2;
    C=4*R0^2*(l^2+m^2);
    D=2*(xi0*l+yi0*m)*4*R0^2;
    E=(xi0^2+yi0^2)*4*R0^2;
    c=[1,4*A,4*A^2+2*B-C,4*A*B-D,B^2-E];
    kj=roots(c);
    xj=kj*l+xi0;
    yj=kj*m+yi0;
    zj=kj*n+zi0;
    jd=[xj,yj,zj];
    dj=(xj-xi0).^2+(yj-yi0).^2+(zj-zi0).^2;
    dirj=[(xj-xi0)/l,(yj-yi0)/m,(zj-zi0)/n];
    disp(size(xj));
    s=0;num=0;dj1=0;
    for j=1:1:length(kj)        %筛选实根,并取沿光线传播方向的实解;
        if isreal(xj(j))==1&isreal(yj(j))==1&isreal(zj(j))==1
            if dirj(j)>=0&isreal(dj(j))==1&dj(j)>10^-13
                num(s+1)=j;     %记录实根的位置序号;
                dj1(s+1)=dj(j); %计算所有实交点到入射点的距离;
                s=s+1;
            end
        end
            
    end      
    if s==0
        disp('there is no real solution2');
        break;                     %无实根;退出对该光线的追踪;
    else
        disp(num) 
        for k=1:1:s
            if dj(num(k))==min(dj1)  %实根中光程最小者为真解;
                jdx=xj(num(k));
                jdy=yj(num(k));
                jdz=zj(num(k));
            end
        end
    end
    jd=[jdx,jdy,jdz];             %确定交点的位置;
    disp(jdz);
%     if atan(jdy/jdx)>pi/2|atan(jdy/jdx)<0  %判断光线是否溢出所定义的边界1/4
     if jdx<0&jdy<0        %半圆边界;        
        disp('out of fiber!');
        break;
     end
kl=linspace(min(0,(jdx-xi0)/l),max(0,(jdx-xi0)/l),400);
x=kl.*l+xi0;
y=kl.*m+yi0;
z=kl.*n+zi0;
% plot3(x,y,z,'r');

%求反射光线:
    c1='x^2+y^2+z^2+R0^2-r0^2-2*R0*sqrt(x^2+y^2)';
    dx=diff(c1,'x');
    dy=diff(c1,'y');
    dz=diff(c1,'z');
    lx=subs(dx,{'R0','x','y'},{R0,jd(1),jd(2)});
    ly=subs(dy,{'R0','x','y'},{R0,jd(1),jd(2)});
    lz=subs(dz,{'z'},{jd(3)});
%法线的方向:
    lfx=lx/sqrt(lx^2+ly^2+lz^2);
    lfy=ly/sqrt(lx^2+ly^2+lz^2);
    lfz=lz/sqrt(lx^2+ly^2+lz^2);
%反射定律求反射光线
    lin=[l,m,n];
    lf=[lfx,lfy,lfz];
    ca=lin*lf';
    if acos(ca)<asin(1/n1)  %不满足全反射定律，光线溢出纤芯;
        disp('light is out because of refraction');
        break;
    end
    lr=lin-2*ca*lf;
    ar=acos(lr)./pi*180;
%反射光线作为下次入射光线的参数：
    xi0=jd(1);
    yi0=jd(2);
    zi0=jd(3);
    l=lr(1);
    m=lr(2);
    n=lr(3);
    lin=[l,m,n];
end

%(xi0,yi0,zi0),lin=(l,m,n),(jdx,jdy,jdy)
%边界平面的坐标表示:(x=0);
xout=-yi0/m*l+xi0;
zout=-yi0/m*n+zi0;
yout=0;

kl=linspace(min(0,(yout-yi0)/m),max(0,(yout-yi0)/m),400);
x=kl.*l+xi0;
y=kl.*m+yi0;
z=kl.*n+zi0;
% plot3(x,y,z,'g');

%求出射光线
lin=[l,m,n];   %incident- light directional vector\
lfin=[0,-1,0]; %法线的方向从入射介质到折射介质;
cai=lin*lfin'; %入射角的余弦
thetain=acos(cai);
cat=sqrt(1-n1^2/n0^2+(n1*cai/n0)^2);
p=n0*cat-n1*cai;
% p=n0*cos(thetat)-(n1*cos(thetain));

thetat=asin(n1*sin(thetain)/n0);
lt=(lin+p*lfin)/n0; %折射光线的方向向量;
lt=lt./sqrt(lt(1)^2+lt(2)^2+lt(3)^2);
%接受出射光的位置;
yout1=-300*r0;fout0=4;D=-yout1/fout0;
kl=linspace(min(0,(yout1-yout)/m),max(0,(yout1-yout)/m),400);
x=kl.*lt(1)+xout;
y=kl.*lt(2)+yout;
z=kl.*lt(3)+zout;
plot3(x,y,z,'r');
xout1=(yout1-yout)/lt(2)*lt(1)+xout;
if abs(xout1)>(R0-D/2)&abs(xout1)<(R0+D/2)
    w1=w1+1;
end
zout1=(yout1-yout)/lt(2)*lt(3)+zout;
if abs(zout1)<(D/2)
    w2=w2+1;
end
lxout(q-1,tt)=xout1;
lzout(q-1,tt)=zout1;
%某一个平面光锥中的光线在x,y,z方向的夹角;
jjy(q-1,tt)=acos((0*lt(1)+1*lt(2)+0*lt(3))/sqrt(0^2+1^2+0^2)/sqrt(lt(1)^2+lt(2)^2+lt(3)^2))/pi*180;
jjx(q-1,tt)=acos((1*lt(1)+0*lt(2)+0*lt(3))/sqrt(1^2+0^2+0^2)/sqrt(lt(1)^2+lt(2)^2+lt(3)^2))/pi*180;
jjz(q-1,tt)=acos((0*lt(1)+0*lt(2)+1*lt(3))/sqrt(0^2+0^2+1^2)/sqrt(lt(1)^2+lt(2)^2+lt(3)^2))/pi*180;
end
phi=phi+10;
deltax(tt)=max(jjx(:,tt))-min(jjx(:,tt));
deltay(tt)=max(jjy(:,tt))-min(jjy(:,tt));
deltaz(tt)=max(jjz(:,tt))-min(jjz(:,tt));
end
disp w1,w2;