verilog实现二倍频与三倍频KDP晶体二倍频与三倍频；基于耦合波方程组；已于实验结果校核.zip

%function [xiaolv_SHG,xiaolv_THG1,xiaolv_THG]=BBTHG_II_t(ENERGY,dtheta2,dtheta31,dtheta32) %function [xiaolv_SHG,xiaolv_THG1,xiaolv_THG]=BBTHG_II_t(dtheta31,dtheta32) %function [xiaolv_THG]=BBTHG_II_t(delta_miu) %%%%%%%%%%%%%% Input 1w pulse parameters %%%%%%%%%%%%%% clc,clear ENERGY=5046/1.5; dtheta2=150; %二倍频晶体失谐角，单位urad dtheta31=250; %三倍频晶体失谐角 ，单位urad dtheta32=-850; ENERGY=ENERGY/0.29; %输入基频光能量，单位J Nt=2^7; %时域切片数 Nx=2^6; crystal_type='KDP'; bandwidth=0.3e-9; %%%%%%%%%%%%%%%% Basic parameters %%%%%%%%%%%%%%%%% Lightspeed=2.9979e+8; %真空中的光速 m/s Epsilon0=8.8542e-12; %真空介电常数ε0 dp=0.2569e-12; %二阶有效非线性系数 m/v Lam1=1.053e-6; %输入基频光中心波长 m w1=2*pi*Lightspeed./Lam1; w2=2*w1; w3=3*w1; %%%%%%%%%%%%%%%% Input crystal parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%% Nz=20; %晶体切片数 crystal_L2=12e-3; %倍频晶体长度 crystal_L31=9e-3; %第一块混频晶体长度 crystal_L32=9e-3; %第一块混频晶体长度 theta20=41.17907*pi/180; %倍频相位匹配角 switch crystal_type case {'KDP'} theta30=59.0188*pi/180; %混频晶体匹配角 dp3=0.3442e-12; %三阶有效非线性系数 m/v [n1_T,n2_T,n3_T,vg1_T,vg2_T,vg3_T]=KDP_BBTHG(w1,theta30+dtheta31*1e-6); deltak_T1=(n1_T*w1+n2_T*w2-n3_T*w3)/Lightspeed; [n1_T_II,n2_T_II,n3_T_II,vg1_T_II,vg2_T_II,vg3_T_II]=KDP_BBTHG(w1,theta30+dtheta32*1e-6); deltak_T2=(n1_T_II*w1+n2_T_II*w2-n3_T_II*w3)/Lightspeed; case {'DKDP'} theta30=60.3569*pi/180; dp3=0.316e-12; %三阶有效非线性系数 m/v [n1_T,n2_T,n3_T,vg1_T,vg2_T,vg3_T]=DKDP_BBTHG(w1,theta30+dtheta31*1e-6); deltak_T1=(n1_T*w1+n2_T*w2-n3_T*w3)/Lightspeed; [n1_T_II,n2_T_II,n3_T_II,vg1_T_II,vg2_T_II,vg3_T_II]=DKDP_BBTHG(w1,theta30+dtheta32*1e-6); deltak_T2=(n1_T_II*w1+n2_T_II*w2-n3_T_II*w3)/Lightspeed; otherwise disp('Unknown crystal!') end Lt=4e-9; %时域取样范围3ns pulsewidth=2e-9; Lx=0.40; beamsize=0.29; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 初始化数据 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% tt=linspace(-Lt/2,Lt/2,Nt); dt=Lt/(Nt-1); xx=linspace(-Lx/2,Lx/2,Nx); dx=Lx/(Nx-1); [X,T]=meshgrid(xx,tt); ft=Nt*T/(Lt^2); fx=Nx*X/(Lx^2); ww=2*pi*ft; % 角频率 mm=2*pi*Lightspeed./(ww+w1); mm3=2*pi*Lightspeed./(ww+w3); %%%%%%%%%%%%%% phase modulation %%%%%%%%%%%% delta_miu=0e9*bandwidth/0.3699e-9; %50GHz %bandwidth=Lam1^2/Lightspeed*delta_miu*1e9; omega=5e9; sigma=delta_miu/omega/2; fi_t=sigma*sin(2*pi*omega*T); At=exp(-0.5*log(2)*(2*T/pulsewidth).^40).*exp(1i*fi_t).*exp(-0.5*log(2)*(2*X/beamsize).^16); % %%%%%%%%%% chirp stacked %%%%%%%%%%%% % T_single=125e-12; % chirp=pi*Lightspeed*bandwidth*T_single/2/log(2)/Lam1^2; % 啁啾系数 % Nsg=2; %超高斯阶数 % At_single=exp(-0.5*log(2)*(2*tt/T_single).^Nsg).*exp(-0.5*log(2)*i*chirp*(2*tt/T_single).^2); % R=90e-12; % At=At_single; % for n=1:1:10 %21 15个堆积成1ns，43个堆积成3ns % At=At+exp(-0.5*log(2)*(2*(tt-n*R)/T_single).^Nsg).*exp(-0.5*log(2)*i*chirp*(2*(tt-n*R)/T_single).^2); % At=At+exp(-0.5*log(2)*(2*(tt+n*R)/T_single).^Nsg).*exp(-0.5*log(2)*i*chirp*(2*(tt+n*R)/T_single).^2); % end %%%%%%%%%%%%% end %%%%%%%%%%%%%%% I10=Lightspeed*Epsilon0.*(abs(At)).^2./2; I1=I10*ENERGY/(sum(I10(:))*dt*dx); A1=sqrt(2*I1/(Lightspeed*Epsilon0)).*exp(1i*angle(At)); % v/m A2=zeros(size(A1)); A3=zeros(size(A1)); % far1w=(abs(fftshift(fft(A1)))).^2; %%%%%%%%%%%%%%%%%%%%%%%%%%%% 位相失配、空间走离、群速度色散 %%%%%%%%%%%%% [n1,n2,p2,vg1,vg2]=KDP_BBSHG(w1,theta20+dtheta2*1e-6); deltak=(2*n1*w1-n2*w2)/Lightspeed; GVD1=2*pi*(1/vg1-1/vg2)*ft; GVD2=2*pi*(1/vg2-1/vg2)*ft; %%%%%%%%%%%%%%%%%%%%%%%%%%% SHG %%%%%%%%%%%%%%%%%%%%%%%% dz=crystal_L2/Nz; A1_temp=A1/sqrt(n1);A2_temp=A2/sqrt(n2); for z=0:dz:crystal_L2-dz A1_temp=ifft2(ifftshift(fftshift(fft2(A1_temp)).*exp(1i*GVD1*dz/2))); A2_temp=ifft2(ifftshift(fftshift(fft2(A2_temp)).*exp(1i*GVD2*dz/2))); % ------4阶R-K法解常微分方程组，非线性频率转换 xishu1=1i*w1*dp./(n1*Lightspeed); xishu2=1i*w2*dp./(2*n2*Lightspeed); Gamaxishu1=1i*w1*n1*Epsilon0; Gamaxishu2=1i*w2*n2*Epsilon0; [K1,L1]=RK_SHG(xishu1,xishu2,Gamaxishu1,Gamaxishu2,deltak*z,A1_temp,A2_temp); [K2,L2]=RK_SHG(xishu1,xishu2,Gamaxishu1,Gamaxishu2,deltak*(z+dz/2),A1_temp+dz*K1/2,A2_temp+dz*L1/2); [K3,L3]=RK_SHG(xishu1,xishu2,Gamaxishu1,Gamaxishu2,deltak*(z+dz/2),A1_temp+dz*K2/2,A2_temp+dz*L2/2); [K4,L4]=RK_SHG(xishu1,xishu2,Gamaxishu1,Gamaxishu2,deltak*(z+dz),A1_temp+dz*K3,A2_temp+dz*L3); A1_temp=A1_temp+dz*(K1+2*K2+2*K3+K4)/6; A2_temp=A2_temp+dz*(L1+2*L2+2*L3+L4)/6; A1_temp=ifft2(ifftshift(fftshift(fft2(A1_temp)).*exp(1i*GVD1*dz/2))); A2_temp=ifft2(ifftshift(fftshift(fft2(A2_temp)).*exp(1i*GVD2*dz/2))); I21=n1.*Lightspeed*Epsilon0.*(abs(A1_temp)).^2./2; I22=n2.*Lightspeed*Epsilon0.*(abs(A2_temp)).^2./2; end A1=A1_temp*sqrt(n1);A2=A2_temp*sqrt(n2); E21_all=sum(I21)*dt; E22_all=sum(I22)*dt; xiaolv_SHG=E22_all/(E21_all+E22_all); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% THG_I %%%%%%%%%%%%%%%%%%%%%%%%$$$$$$$$$$$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%% 位相失配、空间走离、群速度色散 %%%%%%%%%%%%% deltak_T=deltak_T1; GVD1_T=2*pi*(1/vg1_T-1/vg3_T)*ft; GVD2_T=2*pi*(1/vg2_T-1/vg3_T)*ft; GVD3_T=2*pi*(1/vg3_T-1/vg3_T)*ft; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dz=crystal_L31/Nz; A1=A1/sqrt(n1_T);A2=A2/sqrt(n2_T);A3=A3/sqrt(n3_T); A1_temp=A1;A2_temp=A2;A3_temp=A3; qq=1; for z=0:dz:crystal_L31-dz A1_temp=ifft2(ifftshift(fftshift(fft2(A1_temp)).*exp(1i*GVD1_T*dz/2))); A2_temp=ifft2(ifftshift(fftshift(fft2(A2_temp)).*exp(1i*GVD2_T*dz/2))); A3_temp=ifft2(ifftshift(fftshift(fft2(A3_temp)).*exp(1i*GVD3_T*dz/2))); % ------4阶R-K法解常微分方程组，非线性频率转换 xishu1=1i*w1*dp3./(n1_T*Lightspeed); xishu2=1i*w2*dp3./(n2_T*Lightspeed); xishu3=1i*w3*dp3./(n3_T*Lightspeed); Gamaxishu1=1i*w1*n1_T*Epsilon0; Gamaxishu2=1i*w2*n2_T*Epsilon0; Gamaxishu3=1i*w3*n3_T*Epsilon0; Dkz=0; [K1,L1,M1]=RK_THG(xishu1,xishu2,xishu3,Gamaxishu1,Gamaxishu2,Gamaxishu3,Dkz+deltak_T*z,A1_temp,A2_temp,A3_temp); [K2,L2,M2]=RK_THG(xishu1,xishu2,xishu3,Gamaxishu1,Gamaxishu2,Gamaxishu3,Dkz+deltak_T*(z+dz/2),A1_temp+dz*K1/2,A2_temp+dz*L1/2,A3_temp+dz*M1/2); [K3,L3,M3]=RK_THG(xishu1,xishu2,xishu3,Gamaxishu1,Gamaxishu2,Gamaxishu3,Dkz+deltak_T*(z+dz/2),A1_temp+dz*K2/2,A2_temp+dz*L2/2,A3_temp+dz*M2/2); [K4,L4,M4]=RK_THG(xishu1,xishu2,xishu3,Gamaxishu1,Gamaxishu2,Gamaxishu3,Dkz+deltak_T*(z+dz),A1_temp+dz*K3,A2_temp+dz*L3,A3_temp+dz*M3); A1_temp=A1_temp+dz*(K1+2*K2+2*K3+K4)/6; A2_temp=A2_temp+dz*(L1+2*L2+2*L3+L4)/6; A3_temp=A3_temp+dz*(M1+2*M2+2*M3+M4)/6; A1_temp=ifft2(ifftshift(fftshift(fft2(A1_temp)).*exp(1i*GVD1_T*dz/2))); A2_temp=ifft2(ifftshift(fftshift(fft2(A2_temp)).*exp(1i*GVD2_T*dz/2))); A3_temp=ifft2(ifftshift(fftshift(fft2(A3_temp)).*exp(1i*GVD3_T*dz/2))); % I31=n1_T.*Lightspeed*Epsilon0.*(abs(A1_temp)).^2./2; % I32=n2_T.*Lightspeed*