[已解决]Sullivan书中的傅里叶变化程序总运行不对，版主帮忙看看！

Sullivan书中2.2.c程序中傅里叶变化中入射波的变换计算部分，代码如下： QM@zy  
    /* Fourier Transform of the input pulse */ ry0YS\W  
    if (T<100) x.Tulo0/  
    { z['2  
     for (m=0; m<=2; m++) VE?Aa  
     { j%=X ps  
       real_in[m] = real_in[m] + cos (arg[m]*T)*ex[10]; in`|.#  
       imag_in[m] = imag_in[m] - sin (arg[m]*T)*ex[10]; ew<_2Xy"<  
     } aATNeAR  
    }入射波的位置是100，这里ex[10]是不是应该改成100？ 个人对书中程序觉得不正确的地方。 y+nX(@~f
]  
sq?js#C5  
>35W{d  
另：我是用MATLAB的，所以将傅里叶变换部分改成MATLAB程序运算，得不到正确的结果，各位帮忙看看。 NiD_v  
我的代码是在Sullivan书中一维有耗介质传播程序的中加入傅里叶变化部分，如下： 0G`@^`  
clc,clear kO4'|<  
% 空间时间步长 @W-0ybv  
delta_x =0.01; 8&3KVd`  
delta_t = delta_x/(2*3e8); i@mS8%|l  
% 空间时间网格数 2])e}&i  
max_time = 400; t8M\  
max_space = 100; *auT_*  
% 相关参数设置 c&C*'c-r  
sigma = 0.04; z0@BBXQ`  
epszero = 8.85419e-12; LZ RP}|  
muz=4.0*pi*1.0e-7;
`RXlqj#u
 
frequency=2.45e9; ~bsdy2&/q  
omega = 2*pi*frequency; gF
l@A}  
dielectric_constant = 4; izs=5
 
% 初始化电磁场变量 K #}t
\  
E = zeros(max_space,1); 8J:=@X^}  
H = E; U7HfDDh  
%传播波的傅里叶变化变量 b
>Ea_3T/  
real_pt = zeros(max_space,1); %实部变量 h.q9p!  
imag_pt = zeros(max_space,1); %虚部变量 w@pJ49  
ampn = zeros(max_space,1);   %振幅标量 aNP\Q23D  
phasen = zeros(max_space,1); %相位变量 J vq)%t8q>  
%入射波的傅里叶变化变量 7[I +1  
real_in = 0; D7=Irz!O\7  
imag_in = 0; X0]{8v%  
amp_in =0;
opTH6a  
phase_in =0; L8(2or  
% 介质位置设置 v2k@yxt(  
dielectric_begin =50; wsB  
dielectric_end = max_space; 8*Ty`G&v  
dielectic_space_a = ones(max_space,1); Up0kTL  
dielectic_space_b = ones(max_space,1); N5zx#g  
% 方程参数设置 wHh6y?g\  
eta = 0.5; c#TV2@  
for i = 1:max_space }{)Rnb@ >  
    if (i>=dielectric_begin & i <= dielectric_end)       5= T$h;O  
      loss_term = (delta_t*sigma) / (2*epszero*dielectric_constant); qiH)J- ~GZ  
      dielectric_space_a(i)=(1-loss_term)/(1+loss_term); HvmE'O8  
      dielectric_space_b(i)=0.5/(dielectric_constant *(1+loss_term));  'I$kDM mwh  
    else g'cVsO)S  
        dielectric_space_a(i)=1; NS-0-o|4#  
        dielectric_space_b(i)=eta; KW$.Yy  
    end oU se~  
end ~RLWr.pK  
%ABS needed for LHS of problem space do not work with Lossy medium *Qyu QF  
E_low_m2 = 0; M4(57b[`  
E_low_m1 = 0; (I/iD.A  
E_high_m2 = 0; /B)2L]6p  
E_high_m1 = 0; Mfnfp{.)  
%Initialize Pulse Variables（初始化脉冲变量） W8F@nY  
place_pulse = 2; sR/y|  
% spread = 12; 'x5p ?m  
% center_problem_space = max_space/2; *W;;L_V"   
% t0 = 40.0; &j,#5f(  
%Set up dielectric medium within the problem space q;*'V9#  
for n = 1:max_time `@ny!S|1/  
    for k = 2:max_space (4?^X  
        E(k) = dielectric_space_a(k)*E(k) + dielectric_space_b(k)*(H(k-1)-H(k)); <KfR)7I$0a  
    end S|af?IW  
% 傅里叶计算 +\$c_9|C+  
    for k = 1:max_space hu.p;A3p;  
        real_pt(k) = real_pt(k) + cos(2*pi*frequency*delta_t*n)*E(k); %对空间每一点计算实部值 A1r%cs  
        imag_pt(k) = imag_pt(k) - sin(2*pi*frequency*delta_t*n)*E(k); %对空间每一点计算虚部值 [w!
T  
_]-8gr-T  
    end  @Ne&%F?^Z  
% 入射波傅里叶计算 wY ??#pS  
    if n<100 Kj)sL0  
        real_in = real_in+cos(2*pi*frequency*delta_t*n)*E(place_pulse); 41P0)o  
        imag_in = imag_in-sin(2*pi*frequency*delta_t*n)*E(place_pulse); s\<UDW  
    end 2qojU%fiH  
  6lT< lzT  
    pulse = sin(2*pi*frequency*delta_t*n); w 62m}5eA  
    E(place_pulse) = E(place_pulse) + pulse; :{4C2qK>  
    \;KSx3o  
    E(1) = E_low_m2; ug.|ag'R  
    E_low_m2 = E_low_m1; |P`b"x  
    E_low_m1 = E(2); v|KIVBkbT  
    E(max_space) = E_high_m2; :W6'G@ p  
    E_high_m2 = E_high_m1; |}UA=? Xl  
    E_high_m1 = E(max_space-1); KDP"z  
    for h = 1:max_space-1 <B>hvuCoH  
        H(h) = H(h) +eta*(E(h)-E(h+1));
~`="tzr:  
    end ;K~=? k  
    plot(E) }zxf~41  
  &n .. u80C>sQ  
&*Xrh7K2e  

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

for k = 1:max_space  ?? RrrlfFms  
for (m=0; m<=2; m++) ?? *c/|/  
B/K=\qmm  
看看比两个的循环

版主可能没有看过Sullivan书中的C程序，Sullivan书中的C程序激励源是脉冲信号，然后对100MHZ，200MHZ，500MHZ，三个频率进行分析，m变量代表的是这个频率。我的程序只有一个频率，所以没有用m变量。Sullivan书中的C程序放在附件中，辛苦版主看看。我对书中的理论也看了很多次了，还是没找出问题。 tjm@+xs  

刚才把论坛的帖子翻了一遍，看到也有网友出现过这样的问题，帖子地址： ;dYpdy  
http://bbs.rfeda.cn/read-htm-tid-27893.html iwfv t^  
http://bbs.rfeda.cn/read-htm-tid-19840-page-1.html {N$G|bm]u<  
后面是怎么解决的没说啊。 q-}J0vu\K  
大家帮帮忙吧！ 0Q5ua`U  

自己认为ampn变量即幅度值最后应该除以总时间步长数即max_time变量。不除的话时间越长，幅度值越大，这样应该不合理吧！