光参量振荡器Matlab仿真（全）

% OPO close all; format shorteng tic % Parameters mu=4*pi*1e-7; % magnetic permeability c=3*1e8; % speed of light in vacuum (m/s) j=1i; Ngrid=64; % Grid Size L=4e-3; % Physical Grid Size, Window size in meters [dx,dy,x,y,fx,fy]=Grid(Ngrid,L); %Initialize Grid NxN with width/length of L [X,Y]=meshgrid(x,y); % Meshgrid of spatial coord. [FX,FY]=meshgrid(fx,fy); % Meshgrid of freq. coord. % Wavelength of pump signal and idler lp=2.054e-6; %m ls=2.68e-6; li=1/(1/lp-1/ls); % lambda_p=lp; % lambda_s=ls; % lambda_i=li; % % Cavity Modes CrystalLength=15e-3; %Crystal Length M1=500e-3; % Radius of Curvature of input mirror M2=500e-3; % Radius of Curvature of output mirror Mir2CryLengthM1=5e-3; % Free space from M1 to crystal Mir2CryLengthM2=5e-3; % Free space from crystal to M2 T=305; %K - Temperature np=GaAsSell(lp,T); % pump index of refraction ns=GaAsSell(ls,T); % signal index of refraction ni=GaAsSell(li,T); % idler index of refraction % Prop Distance Z=CrystalLength+Mir2CryLengthM1+Mir2CryLengthM2; % Total Single Pass Cavity Prop distance zz=linspace(0,CrystalLength,30);% Prop. distance array in xtal dz=zz(2)-zz(1); % Calculates Cavity Parameters - location of beam waist, beam waist, sizes % at M1 and M2, and provides array z, Wz, Rz, to plot beam size and beam % curvature [zop,wop,wmp1,wmp2,zp,Wzp,Rzp]=CavityBeamSize(CrystalLength,Mir2CryLengthM1,Mir2CryLengthM2,M1,M2,lp,np); [zos,wos,wms1,wms2,zs,Wzs,Rzs]=CavityBeamSize(CrystalLength,Mir2CryLengthM1,Mir2CryLengthM2,M1,M2,ls,ns); [zoi,woi,wmi1,wmi2,zi,Wzi,Rzi]=CavityBeamSize(CrystalLength,Mir2CryLengthM1,Mir2CryLengthM2,M1,M2,li,ni); % Initialize Unit Gaussian Beam with width from cavity calculation % Width taken at Input Mirror pgausbeam=Field(wmp1,wmp1,X,Y,lp,0,M1); % Pump sgausbeam=Field(wms1,wms1,X,Y,ls,0,M1); % Signal igausbeam=Field(wmi1,wmi1,X,Y,li,0,M1); % Idler % Input Power/Energy for pump signal and idler % In Joules for energy (pulsed), In Watts for Power (cw) Ppin=300e-6; % Pump Psin=001e-80; % Signal Piin=00e-6; % Idler Energy_in=Ppin+Psin+Piin; % Input Pulse m=0; %0 for pulse, 1 for cw pulsewidth=55e-9; %FWHM of pump pulse or nt duration for CW beam tau_pulse=pulsewidth*sqrt(2*log(2)); % 1/e2 temporal pulsewidth trt=2*CrystalLength*ni/c+2*(Mir2CryLengthM1+Mir2CryLengthM2)/c; % roundtrip nt for idler lm=1.5; nt=-lm*tau_pulse:trt:lm*tau_pulse; % time array N=length(nt); % Pulse Shape if m==0 % Pulsed B=exp(-(nt./(tau_pulse/2)).^2); % Unit Gaussian pulse elseif m==1 % CW B=ones(1,length(nt)); end % Normalizing Time Slices [p,~]=NormInpBeam(m,pgausbeam,dx,dy,np,mu,c,Ppin,nt,B); % Pump [s,~]=NormInpBeam(m,sgausbeam,dx,dy,ns,mu,c,Psin,nt,B); % Signal [id,dt]=NormInpBeam(m,igausbeam,dx,dy,ni,mu,c,Piin,nt,B); % Idler % Mirror Curvature, % f=R/2 -> R=2f f1=M1/2; f2=M2/2; % Lens Lp1=Lens(lp,f1,X,Y); % M1 Pump Ls1=Lens(ls,f1,X,Y); % M1 Signal Li1=Lens(li,f1,X,Y); % M1 Idler Lp2=Lens(lp,f2,X,Y); % M2 Pump Ls2=Lens(ls,f2,X,Y); % M2 Signal Li2=Lens(li,f2,X,Y); % M2 Idler % Mirror Reflectivities % Input Mirror - M1 Rpi=0; % Pump Rsi=0.99; % Signal Rii=0.99; % Idler % Output Mirror - M2 Rpo=0.99; % Pump Rso=0.95; % Signal Rio=0.05; % Idler % Crystal Loss 1/m loss_p=0.00/1e-3; % Pump loss_s=0.00/1e-3; % Signal loss_i=0.00/1e-3; % Idler % Crystal Reflectivity pf=0.0; % Pump - Front side sf=0.0; % Signal - Front side idf=0.0; % Idler - Front side pb=0.0; % Pump - Back side sb=0.0; % Signal - Back side idb=0.0; % Idler - Back side deff=94*2/pi*1e-12; %pm/V (factor of 2/pi for QPM structures) % Initialize Temporal Outputs PumpOutputOC=zeros(1,length(nt)); SignalOutputOC=zeros(1,length(nt)); IdlerOutputOC=zeros(1,length(nt)); PumpOutputIC=zeros(1,length(nt)); SignalOutputIC=zeros(1,length(nt)); IdlerOutputIC=zeros(1,length(nt)); counter=0; % for t=1:length(nt) % Loops for the number of pulse iterations N if t==1 % First pass, time slice * spatially norm gaus beam * input mirror % reflectivity Ap=p(t)*pgausbeam*sqrt(1-Rpi); As=s(t)*sgausbeam*sqrt(1-Rsi); Ai=id(t)*igausbeam*sqrt(1-Rii); else % Subsequent Passes, remaining cavity energy + new input energy Ap=p(t)*pgausbeam*sqrt(1-Rpi)+Apcirc; As=s(t)*sgausbeam*sqrt(1-Rsi)+Ascirc; Ai=id(t)*igausbeam*sqrt(1-Rii)+Aicirc; end % %%%%%%%%%%%%%%%%% % %Forward Pass % %%%%%%%%%%%%%%%% % Free Space Prop from M1 to Crystal Ap=FSProp(Ap,lp,Mir2CryLengthM1,FX,FY); As=FSProp(As,ls,Mir2CryLengthM1,FX,FY); Ai=FSProp(Ai,li,Mir2CryLengthM1,FX,FY); % Crystal Reflection - Front Ap=Ap*sqrt(1-pf); As=As*sqrt(1-sf); Ai=Ai*sqrt(1-idf); % Nonlinear Crystal [Ap,As,Ai,cavp1,cavs1,cavi1]=NLO4mix(lp,ls,li,np,ns,ni,dz,zz,Ap,As,Ai,FX,FY,dx,dy,mu,c,loss_p,loss_s,loss_i,deff); % Crystal Reflection - Back Ap=Ap*sqrt(1-pb); As=As*sqrt(1-sb); Ai=Ai*sqrt(1-idb); % Free Space Prop from Crystal to M2 Ap=FSProp(Ap,lp,Mir2CryLengthM2,FX,FY); As=FSProp(As,ls,Mir2CryLengthM2,FX,FY); Ai=FSProp(Ai,li,Mir2CryLengthM2,FX,FY); % Displays Beam at Outcoupler mesh(X*1e3,Y*1e3,abs(Ap)); xlabel('X (mm)') ylabel('Y (mm)') F(t)=getframe;%202307221311 % Output Energy (Output Coupler - M2) PumpOutputOC(t)=(np/(2*mu*c))*sum(sum(abs((Ap.*sqrt(1-Rpo))).^2))*dx*dy; SignalOutputOC(t)=(ns/(2*mu*c))*sum(sum(abs((As.*sqrt(1-Rso))).^2))*dx*dy; IdlerOutputOC(t)=(ni/(2*mu*c))*sum(sum(abs((Ai.*sqrt(1-Rio))).^2))*dx*dy; %%%%%%%%%%%%% % Reverse Pass %%%%%%%%%%%%% % Refelcted Beam from Outcoupler Ap=(Ap.*Lp2*sqrt(Rpo)); As=(As.*Ls2*sqrt(Rso)); Ai=(Ai.*Li2*sqrt(Rio)); % Free Space Prop from M2 to Crystal Ap=FSProp(Ap,lp,Mir2CryLengthM2,FX,FY); As=FSProp(As,ls,Mir2CryLengthM2,FX,FY); Ai=FSProp(Ai,li,Mir2CryLengthM2,FX,FY); % Crystal Reflection - Back Ap=Ap*sqrt(1-pb); As=As*sqrt(1-sb); Ai=Ai*sqrt(1-idb); % Nonlinear Crystal [Ap,As,Ai,cavp2,cavs2,cavi2]=NLO4mix(lp,ls,li,np,ns,ni,dz,zz,Ap,As,Ai,FX,FY,dx,dy,mu,c,loss_p,loss_s,loss_i,deff); % Crystal Reflection - Front Ap=Ap*sqrt(1-pf); As=As*sqrt(1-sf); Ai=Ai*sqrt(1-idf); % Free Space Prop from Crystal to M1 Ap=FSProp(Ap,lp,Mir2CryLengthM2,FX,FY); As=FSProp(As,ls,Mir2CryLengthM2,FX,FY); Ai=FSProp(Ai,li,Mir2CryLengthM2,FX,FY); % Circulating Fields Apcirc=Ap.*Lp1.*sqrt(Rpi); Ascirc=As.*Ls1.*sqrt(Rsi); Aicirc=Ai.*Li1.*sqrt(Rii); % Output Energy (Input Coupler - M1) PumpOutputIC(t)=(np/(2*mu*c))*sum(sum(abs(Ap.*Lp1*sqrt(1-Rpi)-p(t)*pgausbeam*sqrt(Rpi)).^2))*dx*dy; SignalOutputIC(t)=(ns/(2*mu*c))*sum(sum(abs(As.*Ls1*sqrt(1-Rsi)-s(t)*sgausbeam*sqrt(Rsi)).^2))*dx*dy; IdlerOutputIC(t)=(ni/(2*mu*c))*sum(sum(abs(Ai.*Li1*sqrt(1-Rii)-id(t)*igausbeam*sqrt(Rii)).^2))*dx*dy; % Internal Fields if t==1 PCav=[cavp1,cavp2]; SCav=[cavs1,cavs2]; ICav=[cavi1,cavi2]; else Ptemp=[cavp1,cavp2]; Stemp=[cavs1,cavs2]; Itemp=[cavi1,cavi2]; ICav=[ICav,Itemp]; SCav=[SCav,Stemp]; PCav=[PCav,Ptemp]; end counter=counter+1 ; end % toc if m == 0 % Energy Outputs for Pulse operation EpOC=sum(PumpOutputOC)*dt; % Pump - Out Coupler EsOC=sum(SignalOutputOC)*dt;% Signal - Out Coupler EiOC=sum(IdlerOutputOC)*dt;% Idler - Out Coupler EpHR=sum(PumpOutputIC)*dt ;% Pump - In Coupler EsHR=sum(SignalOutputIC)*dt ;% Signal - In Coupler EiHR=sum(IdlerOutputIC)*dt; % Idler - In Coupler total=EpOC+EsOC+EiOC+EpHR+EsHR+EiHR ;% Total Energy elseif m == 1 % Power Outputs for or CW operation EpOC=PumpOutputOC(N); % Pump - Out Coupler EsOC=SignalOutputOC(N); % Signal - Out Coupler EiOC=IdlerOutputOC(N); % Idler - Out Coupler EpHR=PumpOutputIC(N); % Pump - In Coupler EsHR=SignalOutputIC(N); % Signal - In Coupler EiHR=IdlerOutputIC(N) ;% Idler - In Coupler total=EpOC+EsOC+EiOC+EpHR+EsHR+EiHR; % Total Power end %% Output from Output Mirror (M2) - Compared to undepleted input pump close all; for t=1:N A=p(t)*pgausbeam; un(t)=(np/(2*mu*c))*sum(sum(abs((A).^2)))*dx*dy; end % Normalized to max of undepleted beam plot(nt*1e9,un/max(un),'k-