CODE: function ctr
clc
clear
clf,
%syms t w;
T1=clock;
%DFB激光器参数
p.c=3*10^8;p.h=6.63*10^-34;
p.V=1.5*10^-16; %active region volume
p.q=1.6*10^-19; %电荷所带电量
p.Gn=2.125*10^-12; %增益斜率
p.af=5.5; %线宽增强因子
p.gm=1*10^-5; %自发辐射因子
p.eps=3*10^-23; %增益饱和参量
p.yt=0.4; %光场限制因子
p.tn=2*10^-9; %载流子寿命
p.tl=9*10^-12; %谐振腔的单次循环时间;
p.n0=4*10^23; %透明载流子密度
p.nt=9.9*10^23; %阈值载流子密度
p.th_e=0;
p.shixie=10*10^9; %失谐频率;
p.vg=8.57*10^7; %光在环形激光器中的传输速率;
p.miucw=p.vg/(1.55006*10^-6); %顺时针绕行的光信号的角频率;顺时针和逆时针的波长怎么给出?
p.w=p.shixie+p.miucw; %角频率;
p.tp=2*10^-12; %光子寿命
p.nps=2.7; %inversion factor;
%SRL的激光器的参数
p.yts=0.226; %光场限制因子;
p.vg=8.57*10^7; %光在环形激光器中的传输速率;
p.as=4.0; %线宽增强因子 4;
p.q=1.6*10^-19; %电荷所带电量;
p.n0s=7.6*10^23; %透明载流子密度;
p.tps=5.17*10^-12; %光子寿命;
p.tns=3.15*10^-9; %载流子寿命;
p.Vs=2.58*10^-16; %volume of the ring resonator cavity of SRL;
p.Gs=1.17*10^-19; %差分增益系数;
p.epss=3.16*10^-23; %self-gain自增益饱和参量;
p.epc=2*p.epss; %cross-gain交叉增益饱和参量;
p.tin=7.33*10^-12; %光在semiconductor ring laser的绕行时间=2*pi*R/p.vg,半径为200微米的laser出射波长1550nm;
p.miucw=p.vg/(1.55006*10^-6); %顺时针绕行的光信号的角频率;顺时针和逆时针的波长怎么给出?
p.miuccw=p.vg/(1.55008*10^-6); %逆时针绕行的光信号的角频率;
p.vs=p.vg/(1.550*10^-6); %光在环形腔的谐振角频率;
%*****控制项--电流--********************
%******DFB的控制电流******************
ith=0.012;%DFB的阈值电流;
rb=2.4; %几倍的阈值电流;
p.it=rb*ith; %DFB的偏置电流;
%****************SRL的偏置电流控制项************
etai=0.8;
ith=21.5*10^-3; %阈值电流;
I=etai*5.12*ith ;%semiconductor ring激光器的注入(偏置)电流;
%*****控制项--注入--********************
p.kinj=0.07;%注入系数;
p.tj=1*10^-9;%注入延迟时间;
p.shixie=10*10^9;%失谐频率;
%初始值**************************************
%t<=0时的值
fs=4*10^10;
ts=0;
tfs=20*10^-9;
%*****dde23************************
xt0=[9.9*10^23;4.02*10^20;0]; %t<=0时的值
xr0=[9.9*10^23;4.02*10^20;0;4.02*10^20;0];
%初始值**************************************
solt=dde23(@sctrt,xt0,[ts,tfs],[],p);
solr=dde23(@sctrr,xr0,[ts,tfs],[],p,solt);
%_____________________________________________________________
t=1*tf/5:1/fs:tf;
%***********************************
[yt,ypt]=deval(solt,t-p.tj);
Pt=p.V*p.h*(p.wt/(2*pi)).*yt(2,:)/(2*p.yt*p.tpt); %光功率
[ft,pwt]=spectrum(fs,Pt);
%***********************************
[yr,ypr]=deval(solr,t);
Pr=p.V*p.h*(p.wr/(2*pi)).*yr(2,:)/(2*p.yt*p.tpr); %光功率
[fr,pwr]=spectrum(fs,Pr);
sinsig=Pt-Pr;
figure(1)
subplot(211),plot(t*10^9,Pt,'black','LineWidth',1.3)
subplot(212),plot(t*10^9,Pr,'r','LineWidth',1.3)
xlabel('time t/ns');ylabel('output power p/W');
figure(2)
plot(Pt,Pr,'*black','MarkerSize',2)
xlabel('output power Pt/W');ylabel('output power Pr/W');
figure(3)
subplot(211),plot(yt(1,:)/p.nt,Pt,'black')
subplot(212),plot(yr(1,:)/p.nt,Pr,'black')
xlabel('n/nth');ylabel('output power P/W');
figure(4)
subplot(211),plot(ft,pwt,'black')
subplot(212),plot(fr,pwr,'r')
xlabel('frequency (Hz)');ylabel('power dBm');
title('spectrum')
%eyediagram(P,6,1*10^-10)
figure(5)
plot(t*10^9,yt(3,:),'black','LineWidth',1.5)
hold on
plot(t*10^9,yr(3,:),'r','LineWidth',1.5)
xlabel('time t/ns');ylabel('phase');
figure(6)
plot(yt(3,:),yr(3,:),'*black','MarkerSize',2)
xlabel('phase of transmitter');ylabel('phase of receiver');
%eyediagram(P,6,1*10^-10)
%*********
std(Pr-Pt)
dp=mean(abs(Pt/max(Pt)-Pr/(max(Pr))))
pht=abs(yt(3,:));phr=abs(yr(3,:));
dph=mean(abs(pht/max(pht)-phr/(max(phr))))
mean((Pr-Pt-mean(Pr-Pt)).^2)
%*********
%t=2*tf/5:1/fs:tf*0.8;
%y1=deval(solr,t);
%y1(2,:)=y1(2,:)-mean(y1(2,:));
%y2=deval(solr,t-1*tf/5);
%y2(2,:)=y2(2,:)-mean(y2(2,:));
%figure(7)
%[c,lag]=xcorr(y1(2,:),y2(2,:),'coeff');
%plot(10^9*lag/fs,c)
figure(7)
plot(t,sinsig)
fprintf('\n');
fprintf('\t');
fprintf('CPUTIME IS');
lll=etime(clock,T1);
disp(lll);
这是主程序,想把sctrt子程序的结果注入sctrr中求出信号功率等画图,但是出现问题,??? Attempted to access tspan(1); index out of bounds because numel(tspan)=0.
Error in ==> dde23 at 100
t0 = tspan(1);
Error in ==> ctr at 82
solt=dde23(@sctrt,xt0,[ts,tfs],[],p);
下面附上子程序:function xdot=sctrt(x,p)
i=p.it;%DFB的偏置电流;
%********transmitter laser with all-optical feedback*********
dnt=i/(p.V*p.q)-1/p.tn*x(1)-x(2)*p.Gn*(x(1)-p.n0)*(1+p.eps*x(2))^-1;
dst=1/p.tn*p.gm*p.yt*x(1)-1/p.tpt*x(2)+x(2)*p.Gn*(x(1)-p.n0)*((1+p.eps*x(2))^-1)*p.yt;
dft=1/2*p.af*(p.Gn*(x(1)-p.n0)*((1+p.eps*x(2))^-1)*p.yt-1/p.tpt);
xdot=[dnt;
dst;
dft];
另一个子程序:function xdot=sctrr(x,p,st)
kcj=p.kinj/p.tin;
if t-p.tj>=0
[xt,xpt]=deval(st,t-p.tj);
else
xt=[0,0,0];
end
%*receiver laser with all-optical feedback and transmitter light injection*
g2=(1-p.epss*x(2)^2-p.epc*x(4)^2);
g4=(1-p.epss*x(4)^2-p.epc*x(2)^2);
dnr=I/(p.q*p.Vs)-x(1)/p.tns-p.vg*p.Gs*(x(1)-p.n0s)*(g2*x(2)^2+g4*x(4)^2);
dsr=1/2*(p.yts*p.vg*p.Gs*(x(1)-p.n0s)*g2-1/p.tps)*x(2)+kcj*xt(2)*cos((x(3)-xt(3))-p.shixie*t);
dfr=1/2*p.as*(p.yts*p.vg*p.Gs*(x(1)-p.n0s)*g2-1/p.tps)-2*pi*(p.miucw-p.vs)-kcj*xt(2)/x(2)*sin((x(3)-xt(3))-p.shixie*t);
dsr2=1/2*(p.yts*p.vg*p.Gs*(x(1)-p.n0s)*g4-1/p.tps)*x(4);
dfr2=1/2*p.as*(p.yts*p.vg*p.Gs*(x(1)-p.n0s)*g4-1/p.tps)-2*pi*(p.miuccw-p.vs);
xdot=[dnr;
dsr;
dfr;
dsr2;
dfr2];
请大家帮我看看如何正确得出结果。如果能出来正确的结果愿意全部的金币奉上
[ Last edited by xiegangmai on 2013-5-9 at 15:25 ]