curved1_s_cone.m

for fiber 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;
w3=0;w4=0;w=zeros(16,1);
xin=input('入射点的x坐标为:xi0=');
yin=input('入射点的y坐标为(<0):yi0=');
zin=input('入射点的z坐标为:zi0=');
phi=input('入射光与子午平面的夹角:phi=');

r=linspace(-r0,r0,202);
for tt=1:1:30  %角度变化;
    jd=0;
   for q=2:1:201  %位置变化;
   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,(xi0-xin)/l),max(0,(xi0-xin)/l),400);
    x=kl.*l+xin;
    y=kl.*m+yin;
    z=kl.*n+zin;
%     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 jdx<0&jdy>0  %判断光线是否溢出边界 1/4 circle
% 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];lin1=lin;
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 jdx<0&jdy>0 %判断光线是否溢出所定义的边界1/4
%      if jdx<0&jdy<0        %半圆边界;        
        disp('out of fiber1!');
        disp(jd);
        jdx=xi0;
        jdy=yi0;
        jdz=zi0;
        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];
    pin=[xi0,yi0,zi0];
    lsin=[l,m,n];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5对于s段的另一段分析%%%%%%%%%%%%%%%%%%%%%%%%%55
%input:[l,m,n];[xi0,yi0,zi0]
%c2='x^2+(y-2*R0)^2+z^2+R0^2-r0^2-2*R0*sqrt(x^2+(y-2*R0)^2)';
A=l*xi0+m*(yi0-2*R0)+n*zi0;
 B=xi0^2+(yi0-2*R0)^2+zi0^2+R0^2-r0^2;
 C=4*R0^2*(l^2+m^2);  
 D=2*(xi0*l+yi0*m-2*R0*m)*4*R0^2;
 E=(xi0^2+yi0^2+4*R0^2-4*R0*yi0)*4*R0^2;
 c=[1,4*A,2*B+4*A^2-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&xj(j)<0&yj(j)>R0-r0
            num(s+1)=j;     %记录实根的位置序号;
            dj1(s+1)=dj(j); %计算所有实交点到入射点的距离;
            s=s+1;
        end
    end
            
end      
if s==0
        disp('there is no real solutions');
   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];             %确定交点的位置;
% plot3(jdx,jdy,jdz,'r*');
disp(jd);


% if atan(jdy/jdx)>pi/2|atan(jdy/jdx)<0  %判断光线是否溢出所定义的边界 1/4
if jdx<0&jdy>2*R0
    disp('out of fiber2!');
    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*R0)^2+z^2+R0^2-r0^2-2*R0*sqrt(x^2+(y-2*R0)^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);
% plot3(xi0,yi0,zi0,'r*');  %进入s型光纤下半段的第一个交点;
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-2*R0)+n*zi0;
 B=xi0^2+(yi0-2*R0)^2+zi0^2+R0^2-r0^2;
 C=4*R0^2*(l^2+m^2);  
 D=2*(xi0*l+yi0*m-2*R0*m)*4*R0^2;
 E=(xi0^2+yi0^2+4*R0^2-4*R0*yi0)*4*R0^2;
 c=[1,4*A,2*B+4*A^2-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>2*R0         
        disp('out of fiber s !');
        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*R0)^2+z^2+R0^2-r0^2-2*R0*sqrt(x^2+(y-2*R0)^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)
%边界平面的坐标表示:yout=2*R0;
yout=2*R0;
xout=(yout-yi0)/m*l+xi0;
zout=(yout-yi0)/m*n+zi0;

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=n1*cos(thetat)-(n0*cos(thetain));

thetat=asin(n1*sin(thetain)/n0);
lt=(lin+p*lfin)/n1; %折射光线的方向向量;
%接受出射光的位置;
yout2=2*R0+300*r0;fout2=4;D1=(yout2-2*R0)/fout2;
kl=linspace(min(0,(yout2-yout)/m),max(0,(yout2-yout)/m),400);
x=kl.*lt(1)+xout;
y=kl.*lt(2)+yout;
z=kl.*lt(3)+zout;
plot3(x,y,z,'r');
xout2=(yout2-yout)/lt(2)*lt(1)+xout;
if abs(xout2)>(R0-D1/2)&abs(xout2)<(R0+D1/2)
    w(3,1)=w(3,1)+1;
end
zout2=(yout2-yout)/lt(2)*lt(3)+zout;
if abs(zout2)<(D1/2)
    w(4,1)=w(4,1)+1;
end

fout2=3.5;D1=(yout2-2*R0)/fout2;
if abs(xout2)>(R0-D1/2)&abs(xout2)<(R0+D1/2)
    w(5,1)=w(5,1)+1;
end
if abs(zout2)<(D1/2)
    w(6,1)=w(6,1)+1;
end

fout2=3;D1=(yout2-2*R0)/fout2;
if abs(xout2)>(R0-D1/2)&abs(xout2)<(R0+D1/2)
    w(7,1)=w(7,1)+1;
end
if abs(zout2)<(D1/2)
    w(8,1)=w(8,1)+1;
end

fout2=2.5;D1=(yout2-2*R0)/fout2;
if abs(xout2)>(R0-D1/2)&abs(xout2)<(R0+D1/2)
    w(9,1)=w(9,1)+1;
end
if abs(zout2)<(D1/2)
    w(10,1)=w(10,1)+1;
end

fout2=2;D1=(yout2-2*R0)/fout2;
if abs(xout2)>(R0-D1/2)&abs(xout2)<(R0+D1/2)
    w(11,1)=w(11,1)+1;
end
if abs(zout2)<(D1/2)
    w(12,1)=w(12,1)+1;
end
fout2=1.5;D1=(yout2-2*R0)/fout2;
if abs(xout2)>(R0-D1/2)&abs(xout2)<(R0+D1/2)
    w(13,1)=w(13,1)+1;
end
if abs(zout2)<(D1/2)
    w(14,1)=w(14,1)+1;
end

fout2=1;D1=(yout2-2*R0)/fout2;
if abs(xout2)>(R0-D1/2)&abs(xout2)<(R0+D1/2)
    w(15,1)=w(15,1)+1;
end
if abs(zout2)<(D1/2)
    w(16,1)=w(16,1)+1;
end

lxout2(q-1,tt)=xout2;
lzout2(q-1,tt)=zout2;
%某一个平面光锥中的光线在x,y,z方向的夹角;
jjy2(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;
jjx2(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;
jjz2(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+6;
deltax2(tt)=max(jjx2(:,tt))-min(jjx2(:,tt));
deltay2(tt)=max(jjy2(:,tt))-min(jjy2(:,tt));
deltaz2(tt)=max(jjz2(:,tt))-min(jjz2(:,tt));
end
disp (w3,w4);