% This part will generate the source of light beams.
% The light beams are parallel to each other and the central one pass through P point(view point),Pv(its position) .Their distances from each other are all a .
%..........................................................................
% The parameters that can bee adjusted are n1,n2 (the refractive index of 1st(source zone) & 2nd environment (sphere)),
% Pv , o & r (center of the sphere & it's radius), a (distance between light beams),
% aq (the width of the beams group)
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
clear all
n1=1 ; n2=4/3;             % for air & glass
Pv=[-4,1,1];
o=[0,0,0] ; r=1 ;          % the center of sphere & it's radius                 
nv=o-Pv ;          %the vector of light beams                                                       
a=0.2 ;
aq=5 ;             %indecating the number of prallel light beams placing in the plate
Lv1=[0.01,0.01,((dot(Pv,nv)-(Pv(1)+.01)*nv(1)-(Pv(2)+.01)*nv(2))/nv(3))-Pv(3)] ;    
t1=a/norm(Lv1) ;
Lv2=cross(nv,Lv1) ;            
t2=a/norm(Lv2) ;
for s1=1:(2*aq+1)           %Generating mesh of light beams
    PL1(:,:,s1)=[Pv(1)+(s1-aq-1)*t1*Lv1(1),Pv(2)+(s1-aq-1)*t1*Lv1(2),Pv(3)+(s1-aq-1)*t1*Lv1(3)] ;
    for s2=1:(2*aq+1)
        PL2(s2,:)=[PL1(1,1,s1)+(s2-aq-1)*t2*Lv2(1),PL1(1,2,s1)+(s2-aq-1)*t2*Lv2(2),PL1(1,3,s1)+(s2-aq-1)*t2*Lv2(3)] ;
    end
    PLT(:,:,s1)=[PL2] ;
end

%Primary step Of visualising ---- glass sphere

sphere(50) ;
shading interp ; colormap cool ;
grid off ;
hold on ; light('Color','w','Position',[-3 -1 0],'Style','infinite') ;
lighting phong ; alpha(.2) ;
set(findobj(gca,'type','surface'),'FaceLighting','phong','SpecularStrength',1,'DiffuseStrength',0.6,...
    'AmbientStrength',0.9,'SpecularExponent',200,...
    'SpecularColorReflectance',0.4,'BackFaceLighting','lit') ;



%-----------------------------------------------------------------------------------------------------------------

       %finding the contact points
   %Pc1 is the first contact point of each beam of light to sphere

L1=nv./norm(nv) ;
for i=1:2*aq+1
    for j=1:2*aq+1
        if (norm(PLT(j,:,i)-Pv) <= r)
          t01=[ 1/2/(L1(1)^2+L1(2)^2+L1(3)^2)*(2*L1(3)*o(3)-2*PLT(j,3,i)*L1(3)-2*PLT(j,2,i)*L1(2)+2*L1(1)*o(1)+2*L1(2)*o(2)-2*PLT(j,1,i)*L1(1)+2*(2*L1(3)*o(3)*L1(1)*o(1)+2*L1(3)*o(3)*L1(2)*o(2)-2*L1(3)*o(3)*PLT(j,1,i)*L1(1)+2*PLT(j,3,i)*L1(3)*PLT(j,2,i)*L1(2)-2*PLT(j,3,i)*L1(3)*L1(1)*o(1)-2*PLT(j,3,i)*L1(3)*L1(2)*o(2)+2*PLT(j,3,i)*L1(3)*PLT(j,1,i)*L1(1)-2*PLT(j,2,i)*L1(2)*L1(1)*o(1)+2*PLT(j,2,i)*L1(2)*PLT(j,1,i)*L1(1)+2*L1(1)*o(1)*L1(2)*o(2)-2*L1(2)*o(2)*PLT(j,1,i)*L1(1)-L1(1)^2*o(2)^2-L1(1)^2*o(3)^2-L1(1)^2*PLT(j,2,i)^2+L1(1)^2*r^2-L1(1)^2*PLT(j,3,i)^2-L1(2)^2*PLT(j,1,i)^2-L1(2)^2*o(3)^2-L1(2)^2*o(1)^2+L1(2)^2*r^2-L1(2)^2*PLT(j,3,i)^2-L1(3)^2*PLT(j,1,i)^2-L1(3)^2*o(2)^2-L1(3)^2*o(1)^2-L1(3)^2*PLT(j,2,i)^2+L1(3)^2*r^2-2*L1(3)*o(3)*PLT(j,2,i)*L1(2)+2*L1(1)^2*PLT(j,3,i)*o(3)+2*L1(1)^2*PLT(j,2,i)*o(2)+2*L1(2)^2*PLT(j,3,i)*o(3)+2*L1(2)^2*PLT(j,1,i)*o(1)+2*L1(3)^2*PLT(j,1,i)*o(1)+2*L1(3)^2*PLT(j,2,i)*o(2))^(1/2))] ;
          t02=[ 1/2/(L1(1)^2+L1(2)^2+L1(3)^2)*(2*L1(3)*o(3)-2*PLT(j,3,i)*L1(3)-2*PLT(j,2,i)*L1(2)+2*L1(1)*o(1)+2*L1(2)*o(2)-2*PLT(j,1,i)*L1(1)-2*(2*L1(3)*o(3)*L1(1)*o(1)+2*L1(3)*o(3)*L1(2)*o(2)-2*L1(3)*o(3)*PLT(j,1,i)*L1(1)+2*PLT(j,3,i)*L1(3)*PLT(j,2,i)*L1(2)-2*PLT(j,3,i)*L1(3)*L1(1)*o(1)-2*PLT(j,3,i)*L1(3)*L1(2)*o(2)+2*PLT(j,3,i)*L1(3)*PLT(j,1,i)*L1(1)-2*PLT(j,2,i)*L1(2)*L1(1)*o(1)+2*PLT(j,2,i)*L1(2)*PLT(j,1,i)*L1(1)+2*L1(1)*o(1)*L1(2)*o(2)-2*L1(2)*o(2)*PLT(j,1,i)*L1(1)-L1(1)^2*o(2)^2-L1(1)^2*o(3)^2-L1(1)^2*PLT(j,2,i)^2+L1(1)^2*r^2-L1(1)^2*PLT(j,3,i)^2-L1(2)^2*PLT(j,1,i)^2-L1(2)^2*o(3)^2-L1(2)^2*o(1)^2+L1(2)^2*r^2-L1(2)^2*PLT(j,3,i)^2-L1(3)^2*PLT(j,1,i)^2-L1(3)^2*o(2)^2-L1(3)^2*o(1)^2-L1(3)^2*PLT(j,2,i)^2+L1(3)^2*r^2-2*L1(3)*o(3)*PLT(j,2,i)*L1(2)+2*L1(1)^2*PLT(j,3,i)*o(3)+2*L1(1)^2*PLT(j,2,i)*o(2)+2*L1(2)^2*PLT(j,3,i)*o(3)+2*L1(2)^2*PLT(j,1,i)*o(1)+2*L1(3)^2*PLT(j,1,i)*o(1)+2*L1(3)^2*PLT(j,2,i)*o(2))^(1/2))] ;
          tc1=min(t01,t02) ;
          Pc1=real([PLT(j,1,i)+tc1*L1(1),PLT(j,2,i)+tc1*L1(2),PLT(j,3,i)+tc1*L1(3)]) ;       %Pc1 : Point of first contact     
        
        %..................................................................
        %formulating the equation of deflected light beam
        %First deflection :
        
          delC=2*Pc1-2*o ;
          beta1=asin((n1/n2)*(norm(cross(delC,-L1)))/(norm(delC)*norm(-L1))) ;
          ev1=o+0.46*r ;                          %evaluation amount (2)
          Pw1=fsolve(@solver1,ev1,optimset('fsolve'),Pc1,-delC,L1,beta1) ;
          L2=(Pw1-Pc1)./norm(Pw1-Pc1) ;
           L22=(2/(norm(-delC)^2)*dot(L2,-delC)*-delC)-L2 ;
          if dot(L22,L1) >= dot(L2,L1) 
             L2=L22./norm(L22) ;
         end 
          t11=[ 1/2/(L2(1)^2+L2(2)^2+L2(3)^2)*(2*L2(3)*o(3)-2*Pc1(3)*L2(3)-2*Pc1(2)*L2(2)+2*L2(1)*o(1)+2*L2(2)*o(2)-2*Pc1(1)*L2(1)+2*(2*L2(3)*o(3)*L2(1)*o(1)+2*L2(3)*o(3)*L2(2)*o(2)-2*L2(3)*o(3)*Pc1(1)*L2(1)+2*Pc1(3)*L2(3)*Pc1(2)*L2(2)-2*Pc1(3)*L2(3)*L2(1)*o(1)-2*Pc1(3)*L2(3)*L2(2)*o(2)+2*Pc1(3)*L2(3)*Pc1(1)*L2(1)-2*Pc1(2)*L2(2)*L2(1)*o(1)+2*Pc1(2)*L2(2)*Pc1(1)*L2(1)+2*L2(1)*o(1)*L2(2)*o(2)-2*L2(2)*o(2)*Pc1(1)*L2(1)-L2(1)^2*o(2)^2-L2(1)^2*o(3)^2-L2(1)^2*Pc1(2)^2+L2(1)^2*r^2-L2(1)^2*Pc1(3)^2-L2(2)^2*Pc1(1)^2-L2(2)^2*o(3)^2-L2(2)^2*o(1)^2+L2(2)^2*r^2-L2(2)^2*Pc1(3)^2-L2(3)^2*Pc1(1)^2-L2(3)^2*o(2)^2-L2(3)^2*o(1)^2-L2(3)^2*Pc1(2)^2+L2(3)^2*r^2-2*L2(3)*o(3)*Pc1(2)*L2(2)+2*L2(1)^2*Pc1(3)*o(3)+2*L2(1)^2*Pc1(2)*o(2)+2*L2(2)^2*Pc1(3)*o(3)+2*L2(2)^2*Pc1(1)*o(1)+2*L2(3)^2*Pc1(1)*o(1)+2*L2(3)^2*Pc1(2)*o(2))^(1/2))] ;
          t12=[ 1/2/(L2(1)^2+L2(2)^2+L2(3)^2)*(2*L2(3)*o(3)-2*Pc1(3)*L2(3)-2*Pc1(2)*L2(2)+2*L2(1)*o(1)+2*L2(2)*o(2)-2*Pc1(1)*L2(1)-2*(2*L2(3)*o(3)*L2(1)*o(1)+2*L2(3)*o(3)*L2(2)*o(2)-2*L2(3)*o(3)*Pc1(1)*L2(1)+2*Pc1(3)*L2(3)*Pc1(2)*L2(2)-2*Pc1(3)*L2(3)*L2(1)*o(1)-2*Pc1(3)*L2(3)*L2(2)*o(2)+2*Pc1(3)*L2(3)*Pc1(1)*L2(1)-2*Pc1(2)*L2(2)*L2(1)*o(1)+2*Pc1(2)*L2(2)*Pc1(1)*L2(1)+2*L2(1)*o(1)*L2(2)*o(2)-2*L2(2)*o(2)*Pc1(1)*L2(1)-L2(1)^2*o(2)^2-L2(1)^2*o(3)^2-L2(1)^2*Pc1(2)^2+L2(1)^2*r^2-L2(1)^2*Pc1(3)^2-L2(2)^2*Pc1(1)^2-L2(2)^2*o(3)^2-L2(2)^2*o(1)^2+L2(2)^2*r^2-L2(2)^2*Pc1(3)^2-L2(3)^2*Pc1(1)^2-L2(3)^2*o(2)^2-L2(3)^2*o(1)^2-L2(3)^2*Pc1(2)^2+L2(3)^2*r^2-2*L2(3)*o(3)*Pc1(2)*L2(2)+2*L2(1)^2*Pc1(3)*o(3)+2*L2(1)^2*Pc1(2)*o(2)+2*L2(2)^2*Pc1(3)*o(3)+2*L2(2)^2*Pc1(1)*o(1)+2*L2(3)^2*Pc1(1)*o(1)+2*L2(3)^2*Pc1(2)*o(2))^(1/2))] ;
          tc2=max(t11,t12) ;
          Pc2=Pc1+tc2*L2 ;
         
        %Second deflection :
        
          delC2=2*Pc2-2*o ;
          beta2=asin((n2/n1)*(norm(cross(-delC2,-L2)))/(norm(-delC2)*norm(-L2))) ;
          beta2=real(beta2) ;
          ev2=o+1.34*r ;          % Evaluation amount (2)
          Pw2=fsolve(@solver1,ev2,optimset('fsolve'),Pc2,delC2,L2,beta2) ;
          L3=(Pw2-Pc2)./norm(Pw2-Pc2) ;
          
          L32=(2/(norm(delC2)^2)*dot(L3,delC2)*delC2)-L3 ;
          if dot(L32,L2) >= dot(L3,L2)
             L3=L32./norm(L32) ;
         end 
        
        %Showing Results in the Figure
        
          PT1=[PLT(j,:,i)',Pc1'] ;
          PT2=[Pc1',Pc2'] ;
          PT3=[Pc2',(Pc2+3*L3)'] ;
          plot3(PT1(1,:),PT1(2,:),PT1(3,:),'b',PT2(1,:),PT2(2,:),PT2(3,:),'g',PT3(1,:),PT3(2,:),PT3(3,:),'r') ;
          hold on ;
          
      else
          PT2=[PLT(j,:,i)',(PLT(j,:,i)+8*L1)'] ;
          plot3(PT2(1,:),PT2(2,:),PT2(3,:),'b') ;
          hold on ;
      end 
    end 
end 
plot3(o(1),o(2),o(3),'.k','MarkerSize',5) ;
axis equal ;
hold off