1130阅读
6回复

[求助]请教牛人帮忙看看分析一下程序

UID ：36059

注册：2009-06-25
登录：2011-06-20
发帖：46
等级：仿真新人

0楼发表于: 2009-09-20 09:59:57

是仿照Sullivan那本书上的C代码写的，但是运行出来不正确，麻烦哪位仁兄帮忙找找错误？
%%%%%%%%%%%%%%%%%%%单根微带线的仿真MATLAB编写%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%一、参数的设定%%%%%%%%%%%%%%%%%%%%%
IE=26;
JE=120;
KE=18;
ktop=4;%%介质底板的厚度
%%%%%%%%%%%%%%%%%参数的初始化%%%%%%%%%%%%%%%%%%%%%
gax=ones(IE,JE,KE);
gay=ones(IE,JE,KE);
gaz=ones(IE,JE,KE);
gbx=zeros(IE,JE,KE);%%设为零
gby=zeros(IE,JE,KE);
gbz=zeros(IE,JE,KE);

fi1=zeros(1,IE);
fi2=ones(1,IE);
fi3=ones(1,IE);
gi1=zeros(1,IE);
gi2=ones(1,IE);
gi3=ones(1,IE);

    fj1=zeros(1,JE);
    fj2=ones(1,JE);
    fj3=ones(1,JE);

    gj1=zeros(1,JE);
    gj2=ones(1,JE);
    gj3=ones(1,JE);

    fk1=zeros(1,KE);
    fk2=ones(1,KE);
    fk3=ones(1,KE);

    gk1=zeros(1,KE);
    gk2=ones(1,KE);
    gk3=ones(1,KE);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%微带线的设置%
i_strip_st=10;
i_strip_end=17;
j_strip_st=11;
j_strip_end=110;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%介质和计算空间的参数%%%%%%%%%
eps0=8.85e-12;
muo=4*pi*1.0e-7;
c=3*10^8;
ddx=0.0004;
ddy=0.0004;
ddz=0.0004;
dt=ddx/(2*c);
epsr_sub=4.4;
sigma=0.0;
%%%%%%%%介质内计算时的参数设置%%%%%%%%%%%
for i=1:IE
    for j=1:JE
        for k=1:ktop
            gax(i,j,k)=1.0/(epsr_sub+dt*sigma/eps0);
            gay(i,j,k)=1.0/(epsr_sub+dt*sigma/eps0);
            gaz(i,j,k)=1.0/(epsr_sub+dt*sigma/eps0);
        end
    end
end
%%%%%%%%电压源和电阻的位置%%%%%%%%%%
   is=13;
   js=11;
   ks=3;
   ir=13;
   jr=110;
   kr=3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%介质与空气分界面上的参数设置%%%%
for i=1:IE
    for j=1:JE
        epsr_sub_surf=(epsr_sub+1)/2;
        gax(i,j,ktop)=1.0/(epsr_sub_surf+dt*sigma/eps0);
        gay(i,j,ktop)=1.0/(epsr_sub_surf+dt*sigma/eps0);%%%
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%在k=1的平面上加上金属面%%%%%%%%%%%%%%%%%%%%%%%
for i=1:IE
    for j=1:JE

        gax(i,j,1)=0.0;
        gay(i,j,1)=0.0;
    end
end
%%%%%%%%%%%%%在衬底表面上加金属微带线%%%%%%%%%%%%%%%%%%%%%%%
for i=i_strip_st:i_strip_end
    for j=j_strip_st:j_strip_end
        gax(i,j,ktop+1)=0.0;
        gay(i,j,ktop+1)=0.0;
    end
end
%%%%%%%%%%%在电压源和电阻的位置，设置金属地面和导线的连接线%%%%%%%%%%%%%%%%%%%
for k=1:ks-1

    gax(13,11,k)=0;
    gaz(13,11,k)=0;
    %%%%%%%%%%%%%%%%%

    gax(13,110,k)=0;
    gaz(13,110,k)=0;
end
for k=ks+1:ktop
gax(13,11,k)=0;
    gaz(13,11,k)=0;
    %%%%%%%%%%%%%%%%%

    gax(13,110,k)=0;
    gaz(13,110,k)=0;
end   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 修改了

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%初始化PML的参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
npml=5;
for i=1:npml
    xnn=(npml-i)/npml;
    xn=0.33*xnn^3;
    fi1(i)=xn;
    fi1(IE-i)=xn;
    gi2(i)=1/(1+xn);
    gi2(IE-i)=1/(1+xn);
    gi3(i)=(1-xn)/(1+xn);
    gi3(IE-i)=(1-xn)/(1+xn);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    xnn=(npml-i-0.5)/npml;
    xn=0.33*xnn^3;
    gi1(i)=xn;
    gi1(IE-i)=xn;
    fi2(i)=1/(1+xn);
    fi2(IE-i)=1/(1+xn);
    fi3(i)=(1+xn)/(1-xn);
    fi3(IE-i)=(1+xn)/(1-xn);
end
for j=1:npml
    xnn=(npml-j)/npml;
    xn=0.33*xnn^3;
    fj1(j)=xn;
    fj1(JE-j)=xn;
    gj2(j)=1/(1+xn);
    gj2(JE-j)=1/(1+xn);
    gj3(j)=(1-xn)/(1+xn);
    gj3(JE-j)=(1-xn)/(1+xn);
    %%%%%%%%%因为电场和磁场交错开0.5，那么xn的设定时也是错开0.5%%%%%%%%%%%%%%%%%%
    xnn=(npml-j-0.5)/npml;
    xn=0.33*xnn^3;
    gj1(j)=xn;
    gj1(JE-j)=xn;
    fj2(j)=1/(1+xn);
    fj2(JE-j)=1/(1+xn);
    fj3(j)=(1+xn)/(1-xn);
    fj3(JE-j)=(1+xn)/(1-xn);
end
   for k=1:npml
       xnn=(npml-k)/npml;
       xn=0.33*xnn^3;
       fk1(KE-k)=xn;
       gk2(KE-k)=1/(1+xn);
       gk3(KE-k)=(1-xn)/(1+xn);
       xnn=(npml-k-0.5)/npml;
       xn=0.33*xnn^3;
       gk1(KE-k)=xn;
       fk2(KE-k)=1/(1+xn);
       fk3(KE-k)=(1-xn)/(1+xn);
   end
   %%%%%%%%%%%%%%%高斯脉冲响应的参数设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   v0=2;
   spread=150;
   t0=3*T;
   t=0;
   nstep=1;

   %pulse=v0*exp(-4*(((t-t0)/T)^2));
   R=50;
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%初始化电场和磁场%%%%%%%%%%%%%%%%

               dx=zeros(IE,JE,KE);
               dy=zeros(IE,JE,KE);
               dz=zeros(IE,JE,KE);
               hx=zeros(IE,JE,KE);
               hy=zeros(IE,JE,KE);
               hz=zeros(IE,JE,KE);
               ex=zeros(IE,JE,KE);
               ey=zeros(IE,JE,KE);
               ez=zeros(IE,JE,KE);
               idx=zeros(IE,JE,KE);
               idy=zeros(IE,JE,KE);
               idz=zeros(IE,JE,KE);
               ix=zeros(IE,JE,KE);
               iy=zeros(IE,JE,KE);
               iz=zeros(IE,JE,KE);
               ihx=zeros(IE,JE,KE);
               ihy=zeros(IE,JE,KE);
               ihz=zeros(IE,JE,KE);

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%电磁场的迭代计算%%%%%%%%%%%%%%%%%%%%

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 循环开始的地方%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
   T=1;
for t=1:100
    T=T+1;
    pulse=sin(2*pi*1500*1e8*dt*T);
       %pulse=exp((-2)*( (t0-T)/spread ).^2);
%%%%%%%%%%%%计算Dx的值%%%%%%%%%%%%%%%%%=
for i=1:IE-1
     for j=2:JE
         for k=2:KE
             curl_h=hz(i,j,k)-hz(i,j-1,k)-(hy(i,j,k)-hy(i,j,k-1));
             idx(i,j,k)=idx(i,j,k)+curl_h;
             dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h+gi1(i)*idx(i,j,k));
         end
     end
end
%%%%%%%%%计算Dy的值%%%%%%%%%%%%%%%%%%%%
for i=2:IE
     for j=1:JE-1
         for k=2:KE
             curl_h=hx(i,j,k)-hx(i,j,k-1)-(hz(i,j,k)-hz(i-1,j,k));
             idy(i,j,k)=idy(i,j,k)+curl_h;
             dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*(curl_h+gj1(i)*idy(i,j,k));
         end
     end
end
%%%%%%%%%计算Dz的值%%%%%%%%%%%%%%%%%%%%
   gaz0=(1-(dt*ddz/(2*eps0*epsr_sub*R*ddx*ddy)))/(1+(dt*ddz/(2*eps0*epsr_sub*R*ddx*ddy)));
           gbz0=dt/(1+(dt*ddz/(2*eps0*epsr_sub*R*ddx*ddy)))/sqrt(muo*eps0);
           gap0=dt/(R*ddx*ddy)/(1+(dt*ddz/(2*R*eps0*epsr_sub*ddx*ddy)))/sqrt(muo*eps0);
for i=2:IE
     for j=2:JE
         for k=1:KE-1
             curl_h=hy(i,j,k)-hy(i-1,j,k)-(hx(i,j,k)-hx(i,j-1,k));
             idz(i,j,k)=idz(i,j,k)+curl_h;

        if i==is
           if j==js
               if k==ks

                  dz(i,j,k)=gaz0*dz(i,j,k)+gbz0*curl_h+gap0*pulse;

               end
           end

       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % is,js,ks是源的位置，ir，jr，kr是电阻的位置，源是阻性电压源
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif i==ir
               if j==jr
                  if k==kr

                  dz(i,j,k)=gaz0*dz(i,j,k)+gbz0*curl_h;

                  end
             end

      %%%%%%%%%%%%%%%% 此处不对，因为源和负载的地方计算了两次         %   dz(i,j,k)=gaz0*dz(i,j,k)+gbz0*curl_h;
        else
         dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*(curl_h+gk1(k)*idz(i,j,k));
        end
         end
     end
end

%%%%%%%%%%%计算电场值%%%%%%%%%%%%%%%%%%%
for i=1:IE-1
      for j=2:JE-1
          for k=2:KE-1
              ex(i,j,k)=gax(i,j,k)*(dx(i,j,k)-ix(i,j,k));
              ix(i,j,k)=ix(i,j,k)+gbx(i,j,k)*ex(i,j,k);
          end
      end
end
for i=2:IE-1
      for j=1:JE-1
          for k=2:KE-1
              ey(i,j,k)=gay(i,j,k)*(dy(i,j,k)-iy(i,j,k));
              iy(i,j,k)=iy(i,j,k)+gby(i,j,k)*ey(i,j,k);
          end
      end
end
for i=2:IE-1
      for j=2:JE-1
          for k=1:KE-1
              ez(i,j,k)=gaz(i,j,k)*(dz(i,j,k)-iz(i,j,k));
              iz(i,j,k)=iz(i,j,k)+gbz(i,j,k)*ez(i,j,k);
          end
      end
end
%%%%%%最外层表面的设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:IE-1
      for j=1:JE
          ex(i,j,1)=0.0;
          ex(i,j,KE)=0.0;
      end
end
for i=1:IE
      for j=JE-1
          ey(i,j,1)=0.0;
          ey(i,j,KE)=0.0;
      end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:IE-1
      for k=2:KE-1
          ex(i,1,k)=0.0;
          ex(i,JE,k)=0.0;
      end
end
for i=1:IE
      for k=1:KE-1
          ez(i,1,k)=0.0;
          ez(i,JE,k)=0.0;
      end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:IE-1
      for k=2:KE-1
          ey(1,j,k)=0.0;
          ey(IE,j,k)=0.0;
      end
end
for j=2:JE-1
      for k=1:KE-1
          ez(1,j,k)=0.0;
          ez(IE,j,k)=0.0;
      end
end
%%%%%%%%%%磁场的计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:IE
      for j=1:JE-1
          for k=1:KE-1
              curl_e=ey(i,j,k+1)-ey(i,j,k)-(ez(i,j+1,k)-ez(i,j,k));
              ihx(i,j,k)=ihx(i,j,k)+curl_e;
   ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册