RCM.zip_物理

%% 环境设置 clear all close all clc; %% 参数赋值 lamd=632.8*10^-3;% 入射光的波长,单位um nP=1.9649;% 棱镜折射率;金红石材料的折射率是1.9649;钽酸锂材料的折射率为2.1677 zP=0.5;% 棱镜厚度 zA=250*10^-3;% 空气层厚度 % nZ=2.41;% 衬底折射率 %%%%%%%%%%========KTP衬底折射率======= % % %%%%%-------nx切衬底折射率-------- nZ=1.76244;% % %%%%%-------ny切衬底折射率-------- % nZ=1.7718;%1.7715 % % %%%%-------nz切衬底折射率-------- % nZ=1.8642;%1.8644 % %%%%%%%%%%%%=========载入KTP损伤率======= load HeKTP300keV10e1016RBSCx.dat; ws=HeKTP300keV10e1016RBSCx; nweileia=ws(:,2);%%%%ws(:,2)=nx;ws(:,3)=ny;ws(:,4)=nz; zweileia=ws(:,1);%%%%%% longa=length(nweileia); long=longa+2; nweilei=zeros(long,1); zweilei=zeros(long,1); nweilei(1:long-2,1)=nweileia; zweilei(1:long-2,1)=zweileia; nweilei(long-1:long,1)=nZ; zweilei(long-1:long,1)=[zweileia(end)+0.1 zweileia(end)+0.2]'; % % % % % ttweilei=100*10^-3;%单位是微米 % % long=fix(zweileia(end)/ttweilei); % % zweilei=linspace(0,zweileia(end),long); % % nweilei=interp1(zweileia,nweileia,zweilei); % %%%%%%%%%%%============采用RCM给定假设折射率分布========== nG=1.7380+0.01;% 波导层折射率 P1=3;% P1：光学位垒的位置； P2=nG;% P2：波导表面的折射率； P3=0.01;%nZ-1.75;% P3：折射率的最大改变量； P4=P1/10;% P4：上升沿高斯曲线的标准偏差；P1/as，as值越小位垒折射率分布越平缓 P5=0;% P5：峰值的顶端宽度； P6=(10-P1)/28;% P6：下降沿高斯曲线的标准偏差。 P7=P1+P5; zweileiz=8;%单位是微米 ttweilei=100*10^-3;%单位是微米 long=fix(zweileiz/ttweilei); zweilei=linspace(0,zweileiz,long); nweilei=zeros(1,long); for i=1:long if zweilei(i)<P1 nweilei(i)=P2-P3*exp((-(zweilei(i)-P1)^2)/P4^2); else if zweilei(i)>P7 nweilei(i)=nZ-(nZ-P2+P3)*exp((-(zweilei(i)-P7)^2)/P6^2); else if zweilei(i)>=P1 && zweilei(i)<=P7 nweilei(i)=P2-P3; end end end end % % % %%%%%%%%%%%============人为假设的折射率变化========== % nG=1.7380+0.01;% 波导层折射率 % nGmax=1.73;%最小折射率 % z11=0.4;% P1：光学位垒的前位置； % t12=0.6; % z12=z11+t12;% P1：光学位垒的位置； % t13=0.3;% 光学位垒厚度 % % zweileiz=8;%单位是微米 % ttweilei=100*10^-3;%单位是微米 % long=fix(zweileiz/ttweilei); % zweilei=linspace(0,zweileiz,long); % nweilei=zeros(1,long); % zhishu=6; % for i=1:long % if zweilei(i)<z11 % nweilei(i)=nG; % else if zweilei(i)>=z11 && zweilei(i)<z12 % nweilei(i)=((nG-nGmax)/(z11^zhishu-z12^zhishu))*zweilei(i)^zhishu+nG-((nG-nGmax)/(z11^zhishu-z12^zhishu))*z11^zhishu; % % bb=log(nG/nGmax)/(z11-z12); % % nweilei(i)=nG*exp(-bb)*exp(bb*zweilei(i)); % % else if zweilei(i)>=z12 && zweilei(i)<=z12+t13 % nweilei(i)=nGmax; % else if zweilei(i)>z12+t13 % nweilei(i)=nZ; % end % end % end % end % end % % % %%%%==============薄膜材料折射率========= % P1=7; % zweileiz=8;%单位是微米, % ttweilei=50*10^-3;%单位是微米 % long=fix(zweileiz/ttweilei); % zweilei=linspace(0,zweileiz,long); % nweilei=zeros(1,long); % for i=1:long % if zweilei(i)<P1 % nweilei(i)=nG; % else if zweilei(i)>=P1 % nweilei(i)=nZ; % % end % end % end % %%%%%%=============光直接进入衬底============== % zweilei=[0 1 2 3 4]; % nweilei=[nZ nZ nZ nZ nZ]; % long=5; % % figure(1) plot(zweilei,nweilei) legend('位垒折射率分布') xlabel('depth/um') ylabel('refractive index') [nzd,jddd]=min(nweilei); ddd=zweilei(jddd); % MM=((nG^2-nzd^2)^(1/2))*(4*ddd/lamd); % disp(['波导厚度d=' num2str(ddd)]) % disp(['模式数M=' num2str(MM)]) % %%%%%%%%%%%%%===========存储结果======= % nzheshelv=[zweilei' nweilei']; % save('RCM.dat','nzheshelv','-ascii') %%%%%%%%符号说明%%%%%%% % iE1zz=exp(i*k1z*z1); % iE1zf=exp(-i*k1z*z1); % % E1zz=exp(k1z*z1); % E1zf=exp(-k1z*z1); % % iK12z=1+i*k1z/k2z; % iK12f=1-i*k1z/k2z; % % K12z=1+k1z/k2z; % K12f=1-k1z/k2z; %%%%%%%%%%j到j+1%%%%%%%%%%%%% n=zeros(1,long+1); z=zeros(1,long+1); n(1)=nP;%棱镜折射率 n(2)=1.0;%空气折射率 z(1)=zP;%棱镜厚度 z(2)=zA;%空气层厚度 n(3:long+1)=nweilei(2:end); zweileia=zeros(1,long-1); for i=1:long-1 zweileia(i)=zweilei(i+1)-zweilei(i); end z(3:long+1)=zweileia; ii=(-1)^(1/2); longnz=length(n); xxs=zeros(2,2*longnz-2); bbs=zeros(2,2*longnz-2); k0=2*pi/lamd; longst=500; sita=linspace(53*pi/180,70*pi/180,longst); % % sita=linspace(-80*pi/180,80*pi/180,longst); beitaa=k0.*nP.*sin(sita); fanshexishu=zeros(1,longst); for i=1:longst beita=beitaa(i); for j=1:longnz-1 n1k0=(n(j)*k0)^2; n2k0=(n(j+1)*k0)^2; k1z=(n1k0-beita^2)^(1/2); k2z=(n2k0-beita^2)^(1/2); iE1zz=1;%exp(ii*k1z*z(j)); iE1zf=1;%exp(-ii*k1z*z(j)); E1zz=1;%exp(k1z*z(j)); E1zf=1;%exp(-k1z*z(j)); iE2zz=exp(ii*k2z*z(j+1)); iE2zf=exp(-ii*k2z*z(j+1)); E2zz=exp(k2z*z(j+1)); E2zf=exp(-k2z*z(j+1)); iK21z=1+ii*k2z/k1z; iK21f=1-ii*k2z/k1z; K21z=1+k2z/k1z; K21f=1-k2z/k1z; aa(j)=0.5*iE1zf*K21z*iE2zz; bb(j)=0.5*iE1zf*K21f*iE2zf; cc(j)=0.5*iE1zz*K21f*iE2zz; dd(j)=0.5*iE1zz*K21z*iE2zf; % Ej=exp(ii*k1z*z(j)); % % aa(j)=0.5*(1+k2z/k1z)/Ej; % bb(j)=0.5*(1-k2z/k1z)/Ej; % cc(j)=0.5*(1-k2z/k1z)*Ej; % dd(j)=0.5*(1+k2z/k1z)*Ej; xxs(1,2*j-1)=aa(j); xxs(1,2*j)=bb(j); xxs(2,2*j-1)=cc(j); xxs(2,2*j)=dd(j); end bbs(1:2,1:2)=xxs(1:2,1:2); for j=1:longnz-2 bbs(1:2,2*(j+1)-1:2*(j+1))=bbs(1:2,2*j-1:2*j)*xxs(1:2,2*(j+1)-1:2*(j+1)); end MM=bbs(1:2,end-1:end); fanshexishu(i)=(abs(MM(2,1)/MM(1,1))).^2; end for i=1:length(fanshexishu) if abs(fanshexishu(i)-1)<=10^(-4) fanshexishu(i)=1; end end % % dfs1=zeros(1,length(fanshexishu)); % for i=1:length(fanshexishu)-1 % dfs1(i)=fanshexishu(i+1)-fanshexishu(i); % end % % dfsmax=max(dfs1); % dfs=zeros(1,length(fanshexishu)); % for i=1:length(fanshexishu) % if dfs1(i)<dfsmax*0.001 % dfs(i)=0; % else if dfs1(i)>=dfsmax*0.001 % dfs(i)=dfs1(i); % end % end % end % % %%%%%%%%%%寻找极小值点，即寻找波导模式角 IndMin=find(diff(sign(diff(fanshexishu)))>0)+1; % % IndMax=find(diff(sign(diff(data)))<0)+1; % figure; hold on; box on; % plot(1:length(data),data); % figure(2) % jiaodu=sita.*180/pi; % plot(jiaodu,fanshexishu,jiaodu(IndMin),fanshexishu(IndMin),'r^'); figure(3) AA=(45+35/60)*pi/180; neffa=beitaa./k0;%nP.*sin(sita); % % staa=asin(sin(sita)./nP)+AA; % % neffa=nP.*sin(staa); % % neffa=sin(sita).*cos(AA)+((nP^2-(sin(sita)).^2).^(1/2)).*sin(AA); plot(neffa,fanshexishu,neffa(IndMin),fanshexishu(IndMin),'r^'); % %%%%%%%%%%===========存储结果======= % nfs=[neffa' fanshexishu']; % save('RCMRCMP4fs8.dat','nfs','-ascii') % % %%%%%%%%%%%%%===========存储结果======= % nzheshelv=[zweilei' nweilei']; % save('RCRCMP4zheshelv8.dat','nzheshelv','-ascii') % % % % % % plot(IndMin,fanshexishu(IndMin),'r^') % plot(IndMax,data(IndMax),'k*') % legend('曲线','波谷点','波峰点') % title('计算离散节点的波峰波谷信息', 'FontWeight', 'Bold'); % % % xcv=neffa(IndMin); % kx=((k0*nG)^2-(k0*xcv(end-1))^2)^(1/2) % % a2=((k0*xcv(end-1))^2-(k0*nZ)^2)^(1/2) % % a0=((k0*xcv(end-1))^2-(k0)^2)^(1/2) % % kx*P1 % % % atan(a2/kx)+atan(a0/kx)+pi