“ %linked with comsol, this program computes the opto-acoustic coupling
%coefficients and phononic energy spectra of an <elliptical> photonic
%waveguide
%note: the spatial characteristics of the excited phononic modes: for large
%cross-sectional dimension (micronscale), the difference from the eigenmode
%may be more pronounced (due to the densely-distributed eigmodes)
%note: the effect of shape may also be more pronounced at sub-wavelength
%scales!!!
%%
tic;
clear all;
% computing and geometric parameters
%Circle:
str1='LWRod_d2h_fsbs_si_Rev';
str2='EWRod_d2h_fsbs_si_Rev';
% str1='LWRod_Circle_Aniso_ang4';
% str2='EWRod_Circle_Aniso_fwd_4';
% str1='LWRod_Ellipse_Aniso_ang16_45';
% str2='EWRod_Ellipse_Aniso_fwd_ang16_45';
% str1='LW_ellipse_circle';
% str2='EW_ellipse_circle';
% wavlen=1.55e-6;
% dbNr=3.5;
% % scz1=0.6977;
% Note: in the case of circle, a, b should be taken a value slightly different
% from 1.0 to counteract the influence of degeneration (each time the FEM
% is run, different result may be obtained)
a=sqrt(1)*1.2;
b=sqrt(1)*1.2;
radius=1.55e-6/4;
% radius=0.19375e-06;
% scz2=-2*(2*pi*dbNr/(wavlen)*scz1)/(pi/(radius));
% % scz2=0;
dxy1=1.0*1.0e-8;
dxy2=0.5*1.0e-8;
% dxy1=0.1e-8;
% dxy2=0.2e-8;
freq_opt=1.9355e14;
Qf=1000; %mechanical Q-factor
num_eigen_ew=0;
Num_ew=100; %total number of elastic eigenmodes
% optical wave model initialization:
md_lw=mphopen(str1);
% md_lw=mphload(str1);
% md_lw.param.set('ellipse_a',num2str(a));
% md_lw.param.set('ellipse_b',num2str(b));
% md_lw.param.set('dbRadius',horzcat(num2str(radius), ' [m]'));
% md_lw.param.set('scz',num2str(scz1));
% md_lw.study('std1').run;
% mphsave(md_lw);
%normalization of electric field:
% X0=-3e-6:(1e-8):3e-6;
% Y0=-3e-6:(1e-8):3e-6;
X0=-a*radius:dxy1:a*radius;
Y0=-b*radius:dxy1:b*radius;
[xx,yy]=meshgrid(X0,Y0);
Coords=[xx(:),yy(:)]';
Solnum=1; % solnum of optical eigenmode
e1 = mphinterp(md_lw,'ewfd.Ex','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e1=reshape(e1,size(xx));
e2 = mphinterp(md_lw,'ewfd.Ey','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e2=reshape(e2,size(xx));
e3 = mphinterp(md_lw,'ewfd.Ez','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e3=reshape(e3,size(xx));
h1 = mphinterp(md_lw,'ewfd.Hx','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h1=reshape(h1,size(xx));
h2 = mphinterp(md_lw,'ewfd.Hy','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h2=reshape(h2,size(xx));
h3 = mphinterp(md_lw,'ewfd.Hz','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h3=reshape(h3,size(xx));
% 在基础代码的这部分替换原有的Pz计算:
% 方法一:坡印廷矢量方法
% P^{(i)} = 2∫∫ d²r ẑ·([e^{(i)}]^* × h^{(i)})
% 计算电场和磁场的叉乘 [e]^* × h
cross_x = conj(e2).*h3 - conj(e3).*h2;
cross_y = conj(e3).*h1 - conj(e1).*h3;
cross_z = conj(e1).*h2 - conj(e2).*h1;
% 取z分量 ẑ·([e]^* × h) = cross_z
integrand = cross_z;
% 计算功率积分
PzI = 2 * trapz(Y0, trapz(X0, integrand, 2));
% 验证功率的虚部(理论上应为实数)
if abs(imag(PzI)) > 1e-10 * abs(real(PzI))
warning('坡印廷方法:功率积分有显著虚部 %.2e + %.2ei', real(PzI), imag(PzI));
PzI = real(PzI); % 取实部
end
Nrf = sqrt(1/PzI); % normalization factor of electric field
% compact spatial domain for interpolation of electric fields
X0=-(a*radius):dxy1:(a*radius);
Y0=-(b*radius):dxy1:(b*radius);
% X0=-(a*radius*0.55):dxy:(a*radius*0.55);
% Y0=-(b*radius*0.55):dxy:(b*radius*0.55);
[xx,yy]=meshgrid(X0,Y0);
Coords=[xx(:),yy(:)]';
e1 = mphinterp(md_lw,'ewfd.Ex','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e1=reshape(e1,size(xx));
e2 = mphinterp(md_lw,'ewfd.Ey','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e2=reshape(e2,size(xx));
e3 = mphinterp(md_lw,'ewfd.Ez','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e3=reshape(e3,size(xx));
d1 = mphinterp(md_lw,'ewfd.Dx','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d1=reshape(d1,size(xx));
d2 = mphinterp(md_lw,'ewfd.Dy','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d2=reshape(d2,size(xx));
d3 = mphinterp(md_lw,'ewfd.Dz','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d3=reshape(d3,size(xx));
e1=Nrf*e1;
e2=Nrf*e2;
e3=Nrf*e3;
d1=Nrf*d1;
d2=Nrf*d2;
d3=Nrf*d3;
%
% save E_Data e1 e2 e3 d1 d2 d3 xx yy
% e1=conj(e1);
% e2=conj(e2);
% e3=conj(e3);
% d1=conj(d1);
% d2=conj(d2);
% d3=conj(d3);
save E1_data e1 xx yy
save E2_data e2 xx yy
save E3_data e3 xx yy
save D1_data d1 xx yy
save D2_data d2 xx yy
save D3_data d3 xx yy
toc
%%
tic
%load the elastic wave model, param setting and run:
md_ew=mphopen(str2);
% md_ew=mphload(str2);
% md_ew.param.set('ellipse_a',num2str(a));
% md_ew.param.set('ellipse_b',num2str(b));
% md_ew.param.set('dbRadius',horzcat(num2str(radius), ' [m]'));
% md_ew.param.set('scz',num2str(scz2));
% md_ew.study('std1').run;
% str2='EWRod_Circle_Aniso_ang0';
% md_ew=mphopen(str2);
%
% dth=2*pi/100;
% theta=-pi:dth:(pi);
%
% for i=1:length(theta)
% Coords1(:,i)=radius*[cos(theta(i)) sin(theta(i))]';
% end
%
% SolNum_ew=10;
%
% mi_kernel = mphinterp(md_ew,'mi_kernel','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',SolNum_ew,'coord',Coords1);
%
% egd_int=mpheval(md_ew,'egd_int','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
%
% Nrf_ew=sqrt(1/abs(egd_int.d1(1,1)/2));
%
% mi_kernel=Nrf_ew*mi_kernel;
%
% figure; hold on
% plot(theta, real(mi_kernel), 'r');
% plot(theta, imag(mi_kernel), 'g');
coupcoef_data=zeros(Num_ew,3); % array to save the opto-acoustic coupling coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f_count=0;
tic
% for SolNum_ew=1001:2:2000
% for SolNum_ew=1:Num_ew
% for SolNum_ew=1:1
for SolNum_ew=1:2:Num_ew
% for SolNum_ew=5
data=mpheval(md_ew,'freq','Complexfun','on','ComplexOut','on','Solnum',SolNum_ew);
freq=abs(real(data.d1(:,1)))
if abs(freq)>1e8
f_count=f_count+1
rp_coef=mpheval(md_ew,'rp_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
% rp_coef=0;
% es_coef=mpheval(md_ew,'es_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
egd_int=mpheval(md_ew,'egd_int','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
es_coef=mpheval(md_ew,'es_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
% es_coef_2=mpheval(md_ew,'es_coef_2','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
% es_coef_3=mpheval(md_ew,'es_coef_3','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
% es_coef_4=mpheval(md_ew,'es_coef_4','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
Nrf_ew=sqrt(1/abs(egd_int.d1(1,1)/2));
Imb=Nrf_ew*rp_coef.d1(1,1) %pay attention to the proportional factor
% Imb=0;
Ies=Nrf_ew*es_coef.d1(1,1)
% Ies_1=Nrf_ew*es_coef_1.d1(1,1)
% Ies_2=Nrf_ew*es_coef_2.d1(1,1)
% Ies_3=Nrf_ew*es_coef_3.d1(1,1)
% Ies_4=Nrf_ew*es_coef_4.d1(1,1)
% Ies=0;
% coupcoef_data(f_count,:)=[freq Imb Ies];
coupcoef_data(f_count,:)=[freq Imb Ies];
end
end
toc
%generating the phonon spectral data:
nf=0;
nQf=320;
ndf=800;
for i=1:f_count
freq0=abs(coupcoef_data(i,1));
eg_freq=0;
freq1=freq0-nQf*freq0/Qf;
freq2=freq0+nQf*freq0/Qf;
df=(freq2-freq1)/ndf;
for j=1:(ndf+1)
freq=freq1+(j-1)*df;
eg_freq=0;
for k=1:f_count
% eg_freq=eg_freq+abs(coupcoef_data(k,2)+coupcoef_data(k,3))^2*freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; %derivation???
% % eg_freq=eg_freq+abs(coupcoef_data(k,2))^2*freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2;
eg_freq=eg_freq+abs(coupcoef_data(k,2)+coupcoef_data(k,3))^2*coupcoef_data(k,1)^2/(coupcoef_data(k,1)^2+4*Qf^2*(freq-coupcoef_data(k,1))^2);
% coupcoef_data(k,1)^2/(coupcoef_data(k,1)^2+4*Qf^2*(freq-coupcoef_data(k,1))^2);
% freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; %derivation???
end
nf=nf+1;
eg_freq_data(nf,1)=abs(freq);
eg_freq_data(nf,2)=eg_freq/4;
end
end
% eg_freq_data(:,2)=2*2*pi*freq_opt*Qf*eg_freq_data(:,2);
eg_freq_data(:,2)=2*pi*freq_opt*Qf*eg_freq_data(:,2);
[a,b]=sort(eg_freq_data(:,1));
eg_freq_data(:,1)=eg_freq_data(b,1);
eg_freq_data(:,2)=eg_freq_data(b,2);
figure;hold on
filename='d2h_fsbs_si_opt1_Rev.xlsx';
xlswrite(filename,eg_freq_data);
plot(eg_freq_data(:,1),log10(eg_freq_data(:,2)),'r')
% [a,b]=sort(eg_freq_data(:,1));
% figure;hold on
% plot(eg_freq_data(b,1),log10(eg_freq_data(b,2)),'r')
% save data
toc”和“ %linked with comsol, this program computes the opto-acoustic coupling
%coefficients and phononic energy spectra of an <elliptical> photonic
%waveguide
%note: the spatial characteristics of the excited phononic modes: for large
%cross-sectional dimension (micronscale), the difference from the eigenmode
%may be more pronounced (due to the densely-distributed eigmodes)
%note: the effect of shape may also be more pronounced at sub-wavelength
%scales!!!
%%
tic;
clear all;
% computing and geometric parameters
%Circle:
str1='LWRod_d2h_fsbs_si_Rev';
str2='EWRod_d2h_fsbs_si_Rev';
% str1='LWRod_Circle_Aniso_ang4';
% str2='EWRod_Circle_Aniso_fwd_4';
% str1='LWRod_Ellipse_Aniso_ang16_45';
% str2='EWRod_Ellipse_Aniso_fwd_ang16_45';
% str1='LW_ellipse_circle';
% str2='EW_ellipse_circle';
% wavlen=1.55e-6;
% dbNr=3.5;
% % scz1=0.6977;
% Note: in the case of circle, a, b should be taken a value slightly different
% from 1.0 to counteract the influence of degeneration (each time the FEM
% is run, different result may be obtained)
a=sqrt(1)*1.2;
b=sqrt(1)*1.2;
radius=1.55e-6/4;
% radius=0.19375e-06;
% scz2=-2*(2*pi*dbNr/(wavlen)*scz1)/(pi/(radius));
% % scz2=0;
dxy1=1.0*1.0e-8;
dxy2=0.5*1.0e-8;
% dxy1=0.1e-8;
% dxy2=0.2e-8;
freq_opt=1.9355e14;
Qf=1000; %mechanical Q-factor
num_eigen_ew=0;
Num_ew=100; %total number of elastic eigenmodes
% optical wave model initialization:
md_lw=mphopen(str1);
% md_lw=mphload(str1);
% md_lw.param.set('ellipse_a',num2str(a));
% md_lw.param.set('ellipse_b',num2str(b));
% md_lw.param.set('dbRadius',horzcat(num2str(radius), ' [m]'));
% md_lw.param.set('scz',num2str(scz1));
% md_lw.study('std1').run;
% mphsave(md_lw);
%normalization of electric field:
% X0=-3e-6:(1e-8):3e-6;
% Y0=-3e-6:(1e-8):3e-6;
X0=-a*radius:dxy1:a*radius;
Y0=-b*radius:dxy1:b*radius;
[xx,yy]=meshgrid(X0,Y0);
Coords=[xx(:),yy(:)]';
Solnum=1; % solnum of optical eigenmode
e1 = mphinterp(md_lw,'ewfd.Ex','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e1=reshape(e1,size(xx));
e2 = mphinterp(md_lw,'ewfd.Ey','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e2=reshape(e2,size(xx));
e3 = mphinterp(md_lw,'ewfd.Ez','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e3=reshape(e3,size(xx));
h1 = mphinterp(md_lw,'ewfd.Hx','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h1=reshape(h1,size(xx));
h2 = mphinterp(md_lw,'ewfd.Hy','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h2=reshape(h2,size(xx));
h3 = mphinterp(md_lw,'ewfd.Hz','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); h3=reshape(h3,size(xx));
% 在基础代码的这部分替换原有的Pz计算:
% 方法二:电场旋度方法
% P^{(i)} = 1/(-iωμ₀) ∫∫ d²r { [e^*]·[ẑ×(∇×e)] + [e^*]·[∇×(ẑ×e)] }
% 常数定义
omega = 2*pi*freq_opt; % 角频率 [rad/s]
mu0 = 4*pi*1e-7; % 真空磁导率 [H/m]
% 计算电场梯度
dx = X0(2) - X0(1); % 正确的间距计算
dy = Y0(2) - Y0(1);
[dE1_dx, dE1_dy] = gradient(e1, dx, dy);
[dE2_dx, dE2_dy] = gradient(e2, dx, dy);
[dE3_dx, dE3_dy] = gradient(e3, dx, dy);
% 计算ẑ×(∇×e)
% ∇×e = (∂Ez/∂y, -∂Ez/∂x, ∂Ey/∂x - ∂Ex/∂y)
z_cross_curl_e_x = -(-dE3_dx); % = ∂Ez/∂x
z_cross_curl_e_y = -(-dE3_dy); % = ∂Ez/∂y (注意符号)
% 计算∇×(ẑ×e)
% ẑ×e = (-e2, e1, 0)
% ∇×(ẑ×e) = (0, 0, ∂e1/∂x + ∂e2/∂y)
curl_z_cross_e_z = dE1_dx + dE2_dy;
% 被积函数计算
term1 = conj(e1).*z_cross_curl_e_x + conj(e2).*z_cross_curl_e_y;
term2 = conj(e3).*curl_z_cross_e_z;
integrand = term1 + term2;
% 积分并计算功率
P_integral = trapz(Y0, trapz(X0, integrand, 1));
P_total = real(1/(-1i*omega*mu0) * P_integral);
Nrf = sqrt(1/P_total); % normalization factor of electric field
% compact spatial domain for interpolation of electric fields
X0=-(a*radius):dxy1:(a*radius);
Y0=-(b*radius):dxy1:(b*radius);
% X0=-(a*radius*0.55):dxy:(a*radius*0.55);
% Y0=-(b*radius*0.55):dxy:(b*radius*0.55);
[xx,yy]=meshgrid(X0,Y0);
Coords=[xx(:),yy(:)]';
e1 = mphinterp(md_lw,'ewfd.Ex','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e1=reshape(e1,size(xx));
e2 = mphinterp(md_lw,'ewfd.Ey','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e2=reshape(e2,size(xx));
e3 = mphinterp(md_lw,'ewfd.Ez','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); e3=reshape(e3,size(xx));
d1 = mphinterp(md_lw,'ewfd.Dx','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d1=reshape(d1,size(xx));
d2 = mphinterp(md_lw,'ewfd.Dy','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d2=reshape(d2,size(xx));
d3 = mphinterp(md_lw,'ewfd.Dz','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',Solnum,'coord',Coords); d3=reshape(d3,size(xx));
e1=Nrf*e1;
e2=Nrf*e2;
e3=Nrf*e3;
d1=Nrf*d1;
d2=Nrf*d2;
d3=Nrf*d3;
%
% save E_Data e1 e2 e3 d1 d2 d3 xx yy
% e1=conj(e1);
% e2=conj(e2);
% e3=conj(e3);
% d1=conj(d1);
% d2=conj(d2);
% d3=conj(d3);
save E1_data e1 xx yy
save E2_data e2 xx yy
save E3_data e3 xx yy
save D1_data d1 xx yy
save D2_data d2 xx yy
save D3_data d3 xx yy
toc
%%
tic
%load the elastic wave model, param setting and run:
md_ew=mphopen(str2);
% md_ew=mphload(str2);
% md_ew.param.set('ellipse_a',num2str(a));
% md_ew.param.set('ellipse_b',num2str(b));
% md_ew.param.set('dbRadius',horzcat(num2str(radius), ' [m]'));
% md_ew.param.set('scz',num2str(scz2));
% md_ew.study('std1').run;
% str2='EWRod_Circle_Aniso_ang0';
% md_ew=mphopen(str2);
%
% dth=2*pi/100;
% theta=-pi:dth:(pi);
%
% for i=1:length(theta)
% Coords1(:,i)=radius*[cos(theta(i)) sin(theta(i))]';
% end
%
% SolNum_ew=10;
%
% mi_kernel = mphinterp(md_ew,'mi_kernel','Complexfun','on','ComplexOut','on','CoordErr','on','SolNum',SolNum_ew,'coord',Coords1);
%
% egd_int=mpheval(md_ew,'egd_int','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
%
% Nrf_ew=sqrt(1/abs(egd_int.d1(1,1)/2));
%
% mi_kernel=Nrf_ew*mi_kernel;
%
% figure; hold on
% plot(theta, real(mi_kernel), 'r');
% plot(theta, imag(mi_kernel), 'g');
coupcoef_data=zeros(Num_ew,3); % array to save the opto-acoustic coupling coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f_count=0;
tic
% for SolNum_ew=1001:2:2000
% for SolNum_ew=1:Num_ew
% for SolNum_ew=1:1
for SolNum_ew=1:2:Num_ew
% for SolNum_ew=5
data=mpheval(md_ew,'freq','Complexfun','on','ComplexOut','on','Solnum',SolNum_ew);
freq=abs(real(data.d1(:,1)))
if abs(freq)>1e8
f_count=f_count+1
rp_coef=mpheval(md_ew,'rp_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
% rp_coef=0;
% es_coef=mpheval(md_ew,'es_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
egd_int=mpheval(md_ew,'egd_int','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
es_coef=mpheval(md_ew,'es_coef','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
% es_coef_2=mpheval(md_ew,'es_coef_2','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
% es_coef_3=mpheval(md_ew,'es_coef_3','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
% es_coef_4=mpheval(md_ew,'es_coef_4','Complexfun','on','ComplexOut','on','solnum',SolNum_ew);
Nrf_ew=sqrt(1/abs(egd_int.d1(1,1)/2));
Imb=Nrf_ew*rp_coef.d1(1,1) %pay attention to the proportional factor
% Imb=0;
Ies=Nrf_ew*es_coef.d1(1,1)
% Ies_1=Nrf_ew*es_coef_1.d1(1,1)
% Ies_2=Nrf_ew*es_coef_2.d1(1,1)
% Ies_3=Nrf_ew*es_coef_3.d1(1,1)
% Ies_4=Nrf_ew*es_coef_4.d1(1,1)
% Ies=0;
% coupcoef_data(f_count,:)=[freq Imb Ies];
coupcoef_data(f_count,:)=[freq Imb Ies];
end
end
toc
%generating the phonon spectral data:
nf=0;
nQf=320;
ndf=800;
for i=1:f_count
freq0=abs(coupcoef_data(i,1));
eg_freq=0;
freq1=freq0-nQf*freq0/Qf;
freq2=freq0+nQf*freq0/Qf;
df=(freq2-freq1)/ndf;
for j=1:(ndf+1)
freq=freq1+(j-1)*df;
eg_freq=0;
for k=1:f_count
% eg_freq=eg_freq+abs(coupcoef_data(k,2)+coupcoef_data(k,3))^2*freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; %derivation???
% % eg_freq=eg_freq+abs(coupcoef_data(k,2))^2*freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2;
eg_freq=eg_freq+abs(coupcoef_data(k,2)+coupcoef_data(k,3))^2*coupcoef_data(k,1)^2/(coupcoef_data(k,1)^2+4*Qf^2*(freq-coupcoef_data(k,1))^2);
% coupcoef_data(k,1)^2/(coupcoef_data(k,1)^2+4*Qf^2*(freq-coupcoef_data(k,1))^2);
% freq^4/abs(freq^2-(coupcoef_data(k,1)*(1+sqrt(-1)/(2*Qf)))^2)^2; %derivation???
end
nf=nf+1;
eg_freq_data(nf,1)=abs(freq);
eg_freq_data(nf,2)=eg_freq/4;
end
end
% eg_freq_data(:,2)=2*2*pi*freq_opt*Qf*eg_freq_data(:,2);
eg_freq_data(:,2)=2*pi*freq_opt*Qf*eg_freq_data(:,2);
[a,b]=sort(eg_freq_data(:,1));
eg_freq_data(:,1)=eg_freq_data(b,1);
eg_freq_data(:,2)=eg_freq_data(b,2);
figure;hold on
filename='d2h_fsbs_si_opt1_Rev.xlsx';
xlswrite(filename,eg_freq_data);
plot(eg_freq_data(:,1),log10(eg_freq_data(:,2)),'r')
% [a,b]=sort(eg_freq_data(:,1));
% figure;hold on
% plot(eg_freq_data(b,1),log10(eg_freq_data(b,2)),'r')
% save data
toc
”根据上述建议修改后的代码如上,通过Matlab运行程序后的图片坐标值相同,纵坐标数值相差较大,图线的形状相似,请结合相关的物理知识解释一下,其中存在的原因,详细具体的解释一下。