24小时热门版块排行榜

返回列表

xpp8888

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 14
帖子: 8
在线: 10.8小时
虫号: 458562
注册: 2007-11-13
专业: 非线性光学与量子光学

[求助] 帮忙看一下程序解方程部分老说越界，不能正常运行已有1人参与

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 ]

回复此楼