【物理应用】大气吸收损耗附matlab代码

发布于:2023-02-02 ⋅ 阅读:(427) ⋅ 点赞:(0)

 1 内容介绍

大气对太赫兹辐射传输存在一定的非协作性,即存在吸收衰减.为了实现太赫兹辐射的有效应用必 须细致地了解太赫兹辐射大气传输的窗口位置,宽度及大气透过率.选取了处理大气非均匀路径,吸收带重叠等大气辐射传输问题的最精确方法——逐线积分法,发 展计算程序,并基于HITRAN分子谱线数据,对水汽,氧气,臭氧,氮气,二氧化碳等单组分气体分子对太赫兹辐射传输的吸收衰减情况进行了计算与分析,并 给出了在太赫兹电磁波大气传输衰减中占主要因素的水汽和氧气的衰减峰位置.​

2 仿真代码
 

function gam=gamw(f,rp,rt,rho);

eta1=0.955*rp*rt^0.68+0.006*rho;
eta2=0.735*rp*rt^0.5+0.0353*rt^4*rho;

gam=1e-4*rho*rt^2.5*f^2*(...
     3.98*eta1*exp(2.23*(1-rt))/((f-22.235)^2+9.42*eta1^2)*gfun(f,22)+...
    11.96*eta1*exp( 0.7*(1-rt))/((f-183.31)^2+11.14*eta1^2)+...
    0.081*eta1*exp(6.44*(1-rt))/((f-321.226)^2+6.29*eta1^2)+...
     3.66*eta1*exp(1.60*(1-rt))/((f-325.153)^2+9.22*eta1^2)+...
    25.37*eta1*exp(1.09*(1-rt))/((f-380)^2)+...
    17.40*eta1*exp(1.46*(1-rt))/((f-448)^2)+...
    844.6*eta1*exp(0.17*(1-rt))/((f-557)^2)*gfun(f,557)+...
    290.0*eta1*exp(0.41*(1-rt))/((f-752)^2)*gfun(f,752)+...
    83328*eta2*exp(0.99*(1-rt))/((f-1780)^2)*gfun(f,1780));

   

% ITU-R P.676-9, Annex 2 method for computing atmospheric attenuation
close all; clear all;

p=1013; % pressure in hPa (1 atm=1013 hPa)
t=15;    % atmospheric temp in C, determine from maps in P.1510 if not known
rho=7.5  % water vapor density (g/m^3)

rp=p/1013;
rt=288/(273+t);

% Compute and plot specific attenuation
i1=1;
for f=1:350
    ff(i1)=f;
    gamdry(i1)=gamo(f,rp,rt);
    gamwat(i1)=gamw(f,rp,rt,rho);
    i1=i1+1;
end
figure
set(gca,'Fontsize',14)
loglog(ff,gamdry,'linewidth',3)
hold on
loglog(ff,gamwat,'r--','linewidth',3)
loglog(ff,gamdry+gamwat,'ko','markersize',8)

xlabel('Frequency (GHz)')
ylabel('Specific attenuation (dB/km)')
grid on
axis([1 350 1e-3 1e2])
set(gca,'Xtick',[1:10 20:10:100 200 350])
set(gca,'Xticklabel',{'1';'2';'';'';'5';'';'';'';'';'10';'20';'';'';'50';'';'';'';'';'100';'200';'350'})

set(gca,'Ytick',[0.001 0.01 0.1 1 10 100])
set(gca,'Yticklabel',['0.001';' 0.01';' 0.1 ';'  1  ';'  10 ';' 100 '])
legend('Oxygen','Water Vapor','Total')
title('1 atm, 15^\circ C, \rho=7.5 g/m^3')

% Compute and plot zenith attenuation

i1=1;
for f=1:350
  ff(i1)=f;

  t1=4.64/(1+0.066*rp^-2.3)*exp(-((f-59.7)/(2.87+12.4*exp(-7.9*rp)))^2);
  t2=0.14*exp(2.12*rp)/((f-118.75)^2+0.031*exp(2.2*rp));
  t3=0.0114/(1+0.14*rp^-2.6)*f*(-0.0247+0.0001*f+1.61e-6*f^2)/(1-0.0169*f+4.1e-5*f^2+3.2e-7*f^3);
  ho=6.1/(1+0.17*rp^-1.1)*(1+t1+t2+t3);
  if (f<70)&(ho>10.7*rp^0.3)
     ho=10.7*rp^0.3;
     display('Violated condition!'); drawnow;
  end;

  sigw=1.013/(1+exp(-8.6*(rp-0.57)));
  hw=1.66*(1+1.39*sigw/((f-22.235)^2+2.56*sigw)+3.37*sigw/((f-183.31)^2+4.69*sigw)+1.58*sigw/((f-325.1)^2+2.89*sigw));
  
  attendry(i1)=ho*gamo(f,rp,rt);
  attenwat(i1)=hw*gamw(f,rp,rt,rho);
  i1=i1+1;
end
figure
set(gca,'Fontsize',14)
loglog(ff,attendry,'b','linewidth',3)
hold on
loglog(ff,attenwat,'r--','linewidth',3)
loglog(ff,attendry+attenwat,'ko','markersize',8)

xlabel('工作频率 (GHz)')
ylabel('大气气体吸收衰减 (dB)')
grid on
axis([1 350 1e-3 1e3])
set(gca,'Xtick',[1:10 20:10:100 200 350])
set(gca,'Xticklabel',{'1';'2';'';'';'5';'';'';'';'';'10';'20';'';'';'50';'';'';'';'';'100';'200';'350'})

set(gca,'Ytick',[0.001 0.01 0.1 1 10 100 1000])
set(gca,'Yticklabel',['0.001';' 0.01';' 0.1 ';'  1  ';'  10 ';' 100 ';' 1000'])
legend('氧气','水蒸气','总吸收损耗')
title('1 atm, 15^\circ C, \rho=7.5 g/m^3')

 1 内容介绍

大气对太赫兹辐射传输存在一定的非协作性,即存在吸收衰减.为了实现太赫兹辐射的有效应用必 须细致地了解太赫兹辐射大气传输的窗口位置,宽度及大气透过率.选取了处理大气非均匀路径,吸收带重叠等大气辐射传输问题的最精确方法——逐线积分法,发 展计算程序,并基于HITRAN分子谱线数据,对水汽,氧气,臭氧,氮气,二氧化碳等单组分气体分子对太赫兹辐射传输的吸收衰减情况进行了计算与分析,并 给出了在太赫兹电磁波大气传输衰减中占主要因素的水汽和氧气的衰减峰位置.​

2 仿真代码
 

function gam=gamw(f,rp,rt,rho);

eta1=0.955*rp*rt^0.68+0.006*rho;
eta2=0.735*rp*rt^0.5+0.0353*rt^4*rho;

gam=1e-4*rho*rt^2.5*f^2*(...
     3.98*eta1*exp(2.23*(1-rt))/((f-22.235)^2+9.42*eta1^2)*gfun(f,22)+...
    11.96*eta1*exp( 0.7*(1-rt))/((f-183.31)^2+11.14*eta1^2)+...
    0.081*eta1*exp(6.44*(1-rt))/((f-321.226)^2+6.29*eta1^2)+...
     3.66*eta1*exp(1.60*(1-rt))/((f-325.153)^2+9.22*eta1^2)+...
    25.37*eta1*exp(1.09*(1-rt))/((f-380)^2)+...
    17.40*eta1*exp(1.46*(1-rt))/((f-448)^2)+...
    844.6*eta1*exp(0.17*(1-rt))/((f-557)^2)*gfun(f,557)+...
    290.0*eta1*exp(0.41*(1-rt))/((f-752)^2)*gfun(f,752)+...
    83328*eta2*exp(0.99*(1-rt))/((f-1780)^2)*gfun(f,1780));

   

% ITU-R P.676-9, Annex 2 method for computing atmospheric attenuation
close all; clear all;

p=1013; % pressure in hPa (1 atm=1013 hPa)
t=15;    % atmospheric temp in C, determine from maps in P.1510 if not known
rho=7.5  % water vapor density (g/m^3)

rp=p/1013;
rt=288/(273+t);

% Compute and plot specific attenuation
i1=1;
for f=1:350
    ff(i1)=f;
    gamdry(i1)=gamo(f,rp,rt);
    gamwat(i1)=gamw(f,rp,rt,rho);
    i1=i1+1;
end
figure
set(gca,'Fontsize',14)
loglog(ff,gamdry,'linewidth',3)
hold on
loglog(ff,gamwat,'r--','linewidth',3)
loglog(ff,gamdry+gamwat,'ko','markersize',8)

xlabel('Frequency (GHz)')
ylabel('Specific attenuation (dB/km)')
grid on
axis([1 350 1e-3 1e2])
set(gca,'Xtick',[1:10 20:10:100 200 350])
set(gca,'Xticklabel',{'1';'2';'';'';'5';'';'';'';'';'10';'20';'';'';'50';'';'';'';'';'100';'200';'350'})

set(gca,'Ytick',[0.001 0.01 0.1 1 10 100])
set(gca,'Yticklabel',['0.001';' 0.01';' 0.1 ';'  1  ';'  10 ';' 100 '])
legend('Oxygen','Water Vapor','Total')
title('1 atm, 15^\circ C, \rho=7.5 g/m^3')

% Compute and plot zenith attenuation

i1=1;
for f=1:350
  ff(i1)=f;

  t1=4.64/(1+0.066*rp^-2.3)*exp(-((f-59.7)/(2.87+12.4*exp(-7.9*rp)))^2);
  t2=0.14*exp(2.12*rp)/((f-118.75)^2+0.031*exp(2.2*rp));
  t3=0.0114/(1+0.14*rp^-2.6)*f*(-0.0247+0.0001*f+1.61e-6*f^2)/(1-0.0169*f+4.1e-5*f^2+3.2e-7*f^3);
  ho=6.1/(1+0.17*rp^-1.1)*(1+t1+t2+t3);
  if (f<70)&(ho>10.7*rp^0.3)
     ho=10.7*rp^0.3;
     display('Violated condition!'); drawnow;
  end;

  sigw=1.013/(1+exp(-8.6*(rp-0.57)));
  hw=1.66*(1+1.39*sigw/((f-22.235)^2+2.56*sigw)+3.37*sigw/((f-183.31)^2+4.69*sigw)+1.58*sigw/((f-325.1)^2+2.89*sigw));
  
  attendry(i1)=ho*gamo(f,rp,rt);
  attenwat(i1)=hw*gamw(f,rp,rt,rho);
  i1=i1+1;
end
figure
set(gca,'Fontsize',14)
loglog(ff,attendry,'b','linewidth',3)
hold on
loglog(ff,attenwat,'r--','linewidth',3)
loglog(ff,attendry+attenwat,'ko','markersize',8)

xlabel('工作频率 (GHz)')
ylabel('大气气体吸收衰减 (dB)')
grid on
axis([1 350 1e-3 1e3])
set(gca,'Xtick',[1:10 20:10:100 200 350])
set(gca,'Xticklabel',{'1';'2';'';'';'5';'';'';'';'';'10';'20';'';'';'50';'';'';'';'';'100';'200';'350'})

set(gca,'Ytick',[0.001 0.01 0.1 1 10 100 1000])
set(gca,'Yticklabel',['0.001';' 0.01';' 0.1 ';'  1  ';'  10 ';' 100 ';' 1000'])
legend('氧气','水蒸气','总吸收损耗')
title('1 atm, 15^\circ C, \rho=7.5 g/m^3')

 

 
     

3 运行结果

4 参考文献

[1]李瀚宇等. "太赫兹电磁波大气吸收衰减逐线积分计算." 强激光与粒子束 25.6(2013):5.

[2]库流杰. 窄带RCS测量大气吸收衰减修正方法研究[D]. 西安电子科技大学.

[3]贺明超, 禹胜林, 翟光贤,等. 基于MATLAB的毫米波降雨衰减特性研究[J].  2021.

博主简介:擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真,相关matlab代码问题可私信交流。

部分理论引用网络文献,若有侵权联系博主删除。

 
     

3 运行结果

4 参考文献

[1]李瀚宇等. "太赫兹电磁波大气吸收衰减逐线积分计算." 强激光与粒子束 25.6(2013):5.

[2]库流杰. 窄带RCS测量大气吸收衰减修正方法研究[D]. 西安电子科技大学.

[3]贺明超, 禹胜林, 翟光贤,等. 基于MATLAB的毫米波降雨衰减特性研究[J].  2021.

博主简介:擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真,相关matlab代码问题可私信交流。

部分理论引用网络文献,若有侵权联系博主删除。