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

是仿照Sullivan那本书上的C代码写的，但是运行出来不正确，麻烦哪位仁兄帮忙找找错误？ /uR/,R++  
%%%%%%%%%%%%%%%%%%%单根微带线的仿真MATLAB编写%%%%%%%%%%%%%%%%%%% Eld[z{n"  
%%%%%%%%%%%%%%%%%%一、参数的设定%%%%%%%%%%%%%%%%%%%%% o6~JAvw  
IE=26; [N9yWuc  
JE=120; PP!-*~F0Jr  
KE=18; gE^pOn  
ktop=4;%%介质底板的厚度 :qB|~"9O  
%%%%%%%%%%%%%%%%%参数的初始化%%%%%%%%%%%%%%%%%%%%% }><[6Uz%  
gax=ones(IE,JE,KE); oqbz!dM(Z  
gay=ones(IE,JE,KE); L{r4hL [  
gaz=ones(IE,JE,KE); 7hPwa3D^  
gbx=zeros(IE,JE,KE);%%设为零 Dyo^O=0c  
gby=zeros(IE,JE,KE); ~G=E Q]a
  
gbz=zeros(IE,JE,KE); |`o1B;lc  
  0T(+
z)Ki  
fi1=zeros(1,IE); +zLw%WD[l  
fi2=ones(1,IE); B@dCCKc%/  
fi3=ones(1,IE); xb0,dZb  
gi1=zeros(1,IE); /|}yf/^9X  
 gi2=ones(1,IE); (jyufHm  
 gi3=ones(1,IE); Q}<QE:-&E  
(=c,b9cb  
    fj1=zeros(1,JE); 'PFjZGaKR  
    fj2=ones(1,JE); nsVLgTbx  
    fj3=ones(1,JE); FAM:; F30  
    o_k)x3I?  
    gj1=zeros(1,JE); #QcRN?s  
    gj2=ones(1,JE); +Q);t,  
    gj3=ones(1,JE); @+p(%  
p"jze3mF  
    fk1=zeros(1,KE); I 2OQ
  
    fk2=ones(1,KE); Pn.DeoHme  
    fk3=ones(1,KE); $YY{|8@kjv  
    4<E <sD  
    gk1=zeros(1,KE); -gt?5H h  
    gk2=ones(1,KE); sSGXd=":  
    gk3=ones(1,KE); L%\Wt1\[  
 #$2
/<  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BaIuOZ@,  
%微带线的设置% 8qw{e`c  
i_strip_st=10; QZ;DZMP  
i_strip_end=17; &hL2xx=  
j_strip_st=11; > cWE@P  
j_strip_end=110; 6P>}7R}  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uCuB>x&  
%%%%%介质和计算空间的参数%%%%%%%%% 6dz^%Ub  
eps0=8.85e-12; cqs.[0 z#B  
muo=4*pi*1.0e-7; I"3C/ pU2  
c=3*10^8; #Y   
ddx=0.0004; P3=#<Q.  
ddy=0.0004; k%O3\q  
ddz=0.0004; !plu;w  
dt=ddx/(2*c); ybFxz  
epsr_sub=4.4; mDbTOtD  
sigma=0.0;
,1Z([R*  
%%%%%%%%介质内计算时的参数设置%%%%%%%%%%% }KIS_krs  
for i=1:IE `\;Z&jlpT  
    for j=1:JE z^&$6c_  
        for k=1:ktop +'olC^?5 }  
            gax(i,j,k)=1.0/(epsr_sub+dt*sigma/eps0); (/]#G8  
            gay(i,j,k)=1.0/(epsr_sub+dt*sigma/eps0); aTeW#:m  
            gaz(i,j,k)=1.0/(epsr_sub+dt*sigma/eps0); cVxO\M  
        end 7D:rq 8$\  
    end .
%.7~Nu,  
end "cBqZzkk9j  
%%%%%%%%电压源和电阻的位置%%%%%%%%%% k_1@?&3  
   is=13; kb/BEJ  
   js=11; vbtZ5Gm   
   ks=3; e`7>QS;.  
   ir=13; gg%)#0Zi  
   jr=110; ^_P?EJ,)`  
   kr=3; xJ. kd Tr  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% z;<~j=lP  
%%%%%%%%%%%%%%介质与空气分界面上的参数设置%%%% qS!N\p~>  
for i=1:IE E7@Gpu,o  
    for j=1:JE 6\K\d_x  
        epsr_sub_surf=(epsr_sub+1)/2; u1#(~[.  
        gax(i,j,ktop)=1.0/(epsr_sub_surf+dt*sigma/eps0); }}~a4p>%  
        gay(i,j,ktop)=1.0/(epsr_sub_surf+dt*sigma/eps0);%%% m$v >r\*X  
    end uG6.(A1LM  
end i+~BVb  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2?Jw0Wq5D  
%%%在k=1的平面上加上金属面%%%%%%%%%%%%%%%%%%%%%%% GQA\JYw|oY  
for i=1:IE 0}`-vOLd-  
    for j=1:JE 9=T;Dxn  
        2@2d |  
        gax(i,j,1)=0.0; ,>kVVp
u  
        gay(i,j,1)=0.0; Y(kf<Wo  
    end W  
end w <"mS*Q  
%%%%%%%%%%%%%在衬底表面上加金属微带线%%%%%%%%%%%%%%%%%%%%%%% MOCcp s*  
for i=i_strip_st:i_strip_end 4Nt4(3Kf  
    for j=j_strip_st:j_strip_end uSQ#Y^V_  
        gax(i,j,ktop+1)=0.0; <)(W7#Ks  
        gay(i,j,ktop+1)=0.0; wik<#ke  
    end z,SI  
end +YW;63"
o  
%%%%%%%%%%%在电压源和电阻的位置，设置金属地面和导线的连接线%%%%%%%%%%%%%%%%%%% "Z,T%]  
for k=1:ks-1 .n YlYY'   
    .7b%7dQ<\  
    gax(13,11,k)=0; to&,d`k=-  
    gaz(13,11,k)=0; `W~  
    %%%%%%%%%%%%%%%%% s>L
.V2!$0  
   x`@`y7(  
    gax(13,110,k)=0; X_@|+d  
    gaz(13,110,k)=0; 3rMJC\h  
end ?#Z4Dg 9|  
for k=ks+1:ktop N_iy4W(NU  
gax(13,11,k)=0; I{[Z  
    gaz(13,11,k)=0; Q25VG5G  
    %%%%%%%%%%%%%%%%% +ls`;f  
   q jc4IW t~  
    gax(13,110,k)=0; 2/s42 FoG  
    gaz(13,110,k)=0; w"dKOdY  
end   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 修改了 vBF9!6X.  
    GCO: !,1  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% pXN'vP  
%%%%%%%%%%%%初始化PML的参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7[qL~BT+  
npml=5; Q{950$)L
 
for i=1:npml j_<!y(W  
    xnn=(npml-i)/npml; C:5d/9k  
    xn=0.33*xnn^3; {R@V  
    fi1(i)=xn; L2EQ 9i'[  
    fi1(IE-i)=xn; R*lq7n9  
    gi2(i)=1/(1+xn); +>!nqp  
    gi2(IE-i)=1/(1+xn); YMK ![ q-  
    gi3(i)=(1-xn)/(1+xn); Z/?{{}H+  
    gi3(IE-i)=(1-xn)/(1+xn); Fih pp<  
    %%%%%%%%%%%%%%%%%%%%%%%%%%% 7/QK"0  
    xnn=(npml-i-0.5)/npml; 'xkl|P>=],  
    xn=0.33*xnn^3; we^'R}d  
    gi1(i)=xn; S-gO  
    gi1(IE-i)=xn; T9Juq6|  
    fi2(i)=1/(1+xn); FibZT1-k  
    fi2(IE-i)=1/(1+xn); <a
nKw|  
    fi3(i)=(1+xn)/(1-xn); {en'8kS  
    fi3(IE-i)=(1+xn)/(1-xn); a4 N f\7  
end -6NoEmb)\'  
for j=1:npml pl1CPxSdO  
    xnn=(npml-j)/npml; &b5(Su  
    xn=0.33*xnn^3; a oU"  
    fj1(j)=xn; jED.0,+K!  
    fj1(JE-j)=xn; =|I
lORf<  
    gj2(j)=1/(1+xn); [{u3g4`}  
    gj2(JE-j)=1/(1+xn); -/{FGbpR;  
    gj3(j)=(1-xn)/(1+xn); {b4`\I@<  
    gj3(JE-j)=(1-xn)/(1+xn); 1Pw1TO"Z  
    %%%%%%%%%因为电场和磁场交错开0.5，那么xn的设定时也是错开0.5%%%%%%%%%%%%%%%%%% VlA]A,P}i  
    xnn=(npml-j-0.5)/npml; ;zD4#7=  
    xn=0.33*xnn^3; g:JSy  
    gj1(j)=xn; e]88 4FP  
    gj1(JE-j)=xn; ~).D\Q\  
    fj2(j)=1/(1+xn); \#dacQ2E@  
    fj2(JE-j)=1/(1+xn); `{Q'iydU  
    fj3(j)=(1+xn)/(1-xn); _r\M}lDh*  
    fj3(JE-j)=(1+xn)/(1-xn); OQ?N_zs,  
end !^su=c  
   for k=1:npml &U|c=$!\  
       xnn=(npml-k)/npml; !vRZh('R  
       xn=0.33*xnn^3; H#;*kc a4  
       fk1(KE-k)=xn; `
}=R  
       gk2(KE-k)=1/(1+xn); x@rQ7K>  
       gk3(KE-k)=(1-xn)/(1+xn); _Wg}#r  
       xnn=(npml-k-0.5)/npml; `DgK$QM  
       xn=0.33*xnn^3; OmBz's
p:  
       gk1(KE-k)=xn; X]@
"ZV[  
       fk2(KE-k)=1/(1+xn); c6 mS  
       fk3(KE-k)=(1-xn)/(1+xn); &Q?@VNi  
   end !aF~5P7%  
   %%%%%%%%%%%%%%%高斯脉冲响应的参数设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |w{Qwf!2  
   v0=2; $B%KkD  
   spread=150; U[?_|=~7  
   t0=3*T; Wmcd{MOS  
   t=0; Zc1x"j  
   nstep=1; j`RG Moq  
   5{V"!M+<  
   %pulse=v0*exp(-4*(((t-t0)/T)^2)); b1 w@toc  
   R=50; 1r$-Uh  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mWaij]1>  
   %%%%%%%%%%%%初始化电场和磁场%%%%%%%%%%%%%%%% TQ4L~8  
   jD9u(qAlH  
               dx=zeros(IE,JE,KE); }5oI` 9VT  
               dy=zeros(IE,JE,KE); X=]utn  
               dz=zeros(IE,JE,KE); 4{zy)GE|W  
               hx=zeros(IE,JE,KE); {P~rf&Ee  
               hy=zeros(IE,JE,KE); fy"}# 2  
               hz=zeros(IE,JE,KE); H@xS<=:lM  
               ex=zeros(IE,JE,KE); T*C25l;w  
               ey=zeros(IE,JE,KE); 9_:"`)]3B  
               ez=zeros(IE,JE,KE); 9c)#j&2?H  
               idx=zeros(IE,JE,KE); \'j(@b,  
               idy=zeros(IE,JE,KE); \N0vA~N.  
               idz=zeros(IE,JE,KE); ZuGd{p$  
               ix=zeros(IE,JE,KE); A[;deHg=  
               iy=zeros(IE,JE,KE); 65~E<)UJ  
               iz=zeros(IE,JE,KE); s~ 8g  
               ihx=zeros(IE,JE,KE); /C
sP@f_Gw  
               ihy=zeros(IE,JE,KE); LPt9+sauf1  
               ihz=zeros(IE,JE,KE); ;i6~iLY  
      oxc;DfJ_  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% H"AL@=  
   %%%%%%%%%%电磁场的迭代计算%%%%%%%%%%%%%%%%%%%% <L qJg
  
   9!Mh(KtQ  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% o6O-\d7
^M  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% o(jLirnk  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 循环开始的地方%%%%%%%%%%%%%%%%%%%%%%%%%%%% BM /FOY;  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 2n@`Og_0  
   T=1; pK3A/ry<  
for t=1:100 PtW2S 1?j  
    T=T+1; Io3-\Ff  
    pulse=sin(2*pi*1500*1e8*dt*T); t7F0[E'=5\  
       %pulse=exp((-2)*( (t0-T)/spread ).^2); Ju47}t%HB  
 %%%%%%%%%%%%计算Dx的值%%%%%%%%%%%%%%%%%= 2'S&%UyP  
 for i=1:IE-1 96V8R
<   
     for j=2:JE d%'#-w'  
         for k=2:KE P5_Ajb(@'  
             curl_h=hz(i,j,k)-hz(i,j-1,k)-(hy(i,j,k)-hy(i,j,k-1)); `Fr ,,Q81\  
             idx(i,j,k)=idx(i,j,k)+curl_h; uM[|>t  
             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)); 2\1+M)  
         end lF}@@e)N  
     end .i4aM;Qy  
 end ;Y*K!iFWH  
 %%%%%%%%%计算Dy的值%%%%%%%%%%%%%%%%%%%% 8~C}0
H  
 for i=2:IE BVb^xL  
     for j=1:JE-1 LsERcjwwK  
         for k=2:KE Yt
(FSb31H  
             curl_h=hx(i,j,k)-hx(i,j,k-1)-(hz(i,j,k)-hz(i-1,j,k)); ekyCZ8iai  
             idy(i,j,k)=idy(i,j,k)+curl_h; ,jg #^47I  
             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)); o1(;"5MM  
         end hTn"/|_SW  
     end d|NW&PG  
 end t(*n[7e  
  %%%%%%%%%计算Dz的值%%%%%%%%%%%%%%%%%%%% K& ^qn&  
   gaz0=(1-(dt*ddz/(2*eps0*epsr_sub*R*ddx*ddy)))/(1+(dt*ddz/(2*eps0*epsr_sub*R*ddx*ddy))); bOr11?
 
           gbz0=dt/(1+(dt*ddz/(2*eps0*epsr_sub*R*ddx*ddy)))/sqrt(muo*eps0); :zKW[sF  
           gap0=dt/(R*ddx*ddy)/(1+(dt*ddz/(2*R*eps0*epsr_sub*ddx*ddy)))/sqrt(muo*eps0); IVSC7SBiT  
  for i=2:IE 
=figat  
     for j=2:JE ^ul1{  
         for k=1:KE-1 :Pdh##k  
             curl_h=hy(i,j,k)-hy(i-1,j,k)-(hx(i,j,k)-hx(i,j-1,k)); "{D/a7]lC  
             idz(i,j,k)=idz(i,j,k)+curl_h; M+ %O-B  
             iiq `:G  
        if i==is 3/n?g7B  
           if j==js `Uz.9_6  
               if k==ks h76j|1gI  
          V[/9?5pM  
                  dz(i,j,k)=gaz0*dz(i,j,k)+gbz0*curl_h+gap0*pulse; 6L8nw+mEK  
             \T_ZcV  
               end }?Y -I> w  
           end La1:WYt  
        jJiuq#;T3  
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LvG$J*  
       % is,js,ks是源的位置，ir，jr，kr是电阻的位置，源是阻性电压源 7w)8s  
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% _r3Y$^!U  
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cDz@3So.b  
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 9/0H,qZc  
        elseif i==ir 3?FY?Q[  
               if j==jr a^J(TW/  
                  if k==kr K _VIk'RB  
          %;QK5L 
 
                  dz(i,j,k)=gaz0*dz(i,j,k)+gbz0*curl_h; 2Cp4aTGv#  
             j "<?9/r  
                  end yg}O9!MJ  
             end 49*f=gpGj2  
       c2g
[w;0"  
      %%%%%%%%%%%%%%%%  此处不对，因为源和负载的地方计算了两次         %   dz(i,j,k)=gaz0*dz(i,j,k)+gbz0*curl_h; R|qrK  
        else [m:cO6DM,  
         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)); _1gNU]"  
        end ek]JzD~w$  
         end !I?C8)  
     end #1J,!seJ  
  end Lbz/M_G  
ZzE(S  
i`F5  
  %%%%%%%%%%%计算电场值%%%%%%%%%%%%%%%%%%%     Qt4mg?X/  
  for i=1:IE-1 o4FHR+u<M  
      for j=2:JE-1 45.ks.  
          for k=2:KE-1 =H;n$ -P  
              ex(i,j,k)=gax(i,j,k)*(dx(i,j,k)-ix(i,j,k)); np^&cY]  
              ix(i,j,k)=ix(i,j,k)+gbx(i,j,k)*ex(i,j,k); <T[LugI  
          end 3'.3RKV  
      end fU$Jh/#":  
  end n}Z%D-b$  
  for i=2:IE-1 Lf%3-P  
      for j=1:JE-1 xFp$JN  
          for k=2:KE-1 3\m!  
              ey(i,j,k)=gay(i,j,k)*(dy(i,j,k)-iy(i,j,k)); Fc`IRPW<  
              iy(i,j,k)=iy(i,j,k)+gby(i,j,k)*ey(i,j,k); n`Pl:L*kG  
          end Y[7p
rjd  
      end 9)tb=  
  end 6t;;Fz  
  for i=2:IE-1 %W D^0U|  
      for j=2:JE-1 }cMkh  
          for k=1:KE-1 KU$,{Sn6@  
              ez(i,j,k)=gaz(i,j,k)*(dz(i,j,k)-iz(i,j,k)); :} =lE"2  
              iz(i,j,k)=iz(i,j,k)+gbz(i,j,k)*ez(i,j,k); ^oPFLez56  
          end ]}/Rl}_  
      end h623)C;  
  end L-?ty@-i  
  %%%%%%最外层表面的设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% apm%\dN  
  for i=1:IE-1 D@uVb4uK  
      for j=1:JE 72~L ?  
          ex(i,j,1)=0.0; :& Dv!z  
          ex(i,j,KE)=0.0; }TMO>eB'  
      end
OlyW/hd  
  end dnD@BQ  
  for i=1:IE HQ"T>xb  
      for j=JE-1 ZJs~,Q  
          ey(i,j,1)=0.0; +g?uvXC&  
          ey(i,j,KE)=0.0; &l2xh~L  
      end Y }VJ4!%U
 
  end [G8EX3  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ri4z^1\  
  for i=1:IE-1 A-4;$ QSm  
      for k=2:KE-1 c9EtUv~  
          ex(i,1,k)=0.0; AAa7)^R  
          ex(i,JE,k)=0.0; "W+>?u)  
      end xT&
~{,9  
  end ~

`BkCTT  
  for i=1:IE yrEh5v:  
      for k=1:KE-1 /A0_#g:2*#  
          ez(i,1,k)=0.0; $rW(*#C  
          ez(i,JE,k)=0.0; mF@7;dpr  
      end Ti)Me-g  
  end 5OWyxO3{  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% X9?)P5h=  
  for i=1:IE-1 >+7{PF+sB  
      for k=2:KE-1 9cB+x`+Lu  
          ey(1,j,k)=0.0; =!SV;^-q  
          ey(IE,j,k)=0.0; 0^>,  
      end n32"cFPpT  
  end #:BkDidt2v  
  for j=2:JE-1 \l(J6T
u  
      for k=1:KE-1 <Mvniz  
          ez(1,j,k)=0.0; u4FD}nV  
          ez(IE,j,k)=0.0; m BvO<?ec  
      end +:^l|6%}  
  end JqO1 a?H  
  %%%%%%%%%%磁场的计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mTu>S  
  for i=1:IE 3\}u#/Vb  
      for j=1:JE-1 rVP\F{Q4Tr  
          for k=1:KE-1 |qe;+)0>K  
              curl_e=ey(i,j,k+1)-ey(i,j,k)-(ez(i,j+1,k)-ez(i,j,k)); *CXc{{  
              ihx(i,j,k)=ihx(i,j,k)+curl_e; &X:;B'   
   .. O'98OH+u  
,(q] $eO
Z  

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

Sullivan书上讲的三维吸收边界有些问题，吸收效果不好，你试试用别的吸收边界吧!

但是好像不是吸收边界条件的问题，波形不在微带线上传播

难道用这本书的边界条件运行不出书上的程序？还是说，不怎么准备，但是图基本相似额》？？我也调了好久了，一直调不出来，还希望给点指点额··谢谢了··

全部检查对了，就能够运行出书上的图，还是希望楼主仔细检查下。

  太难了···

我想问下 楼主的pml这段代码 还需要修改么？？