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

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')