空气折射率公式有很多，比如：
目录Birch and Downs 1994提出的Updated Edlén Equation for the Refractive Index of Air （n 0.35-0.65 μm）（Ciddor 1996: n 0.23-1.690 μm）MATLABPython
Birch and Downs 1994提出的Updated Edlén Equation for the Refractive Index of Air （n 0.35-0.65 μm）
《光在湍流大气中的传播》1.2节给出的公式是错误的，我下面给出的py代码是正确的。具体可以和【原论文Correction to the Updated Edlén Equation for the Refractive Index of Air】比对
# -*- coding: utf-8 -*-
# Author: aspen138
# Last modified: 2024-02-21
# Original data: Birch and Downs 1994, https://doi.org/10.1088/0026-1394/31/4/006
###############################################################################
import numpy as np
import matplotlib.pyplot as plt
π = np.pi
def n(λ, t, p, h, xc):
# λ: wavelength, 0.3 to 1.69 μm
# t: temperature, -40 to +100 °C
# p: pressure, 80000 to 120000 Pa
# h: fractional humidity, 0 to 1
# xc: CO2 concentration, 0 to 2000 ppm
esT= 6.112*1e-3*np.exp(17.62*t /(243.12+t)) # the saturation vapor pressure, using the Magnus-Tetens approximation
pv= esT*h # pressure of water vapour
A = 8342.54
B = 2406147
C = 15998
D = 96095.43
E = 0.601
F = 0.00972
G = 0.003661
ns = 1 + 1e-8 * (A + B / (130 - λ**(-2) ) + C / (38.9 - λ**(-2)))
ntp=1+p*(ns-1)*(1+1e-8*(E-F*t)*p)/(1+G*t)/D
n=ntp-pv*(3.73345-0.0401/ (λ)**2 )*1e-10
return n
# output - modify code below the line to match your needs
##############################################################################
# use this to calculate n at particular conditions
print("n =",n(0.633,20,101325,0,450))
#plot n vs μm
λ = np.arange(0.3, 1.691, 0.01)
n1 = n(λ,15,101325,0,450) #dry air, 15 °C, 450 ppm
n2 = n(λ,15,101325,0.5,450) #50% humidity, 15 °C, 450 ppm
n3 = n(λ,15,101325,0,370) #dry air, 15 °C, 370 ppm
n4 = n(λ,26.85,101325,0,450) #dry air, 300K, 450 ppm
plt.rc('font', family='Arial', size='14')
plt.figure(1)
# plt.plot(λ, n1-1, label="dry air, 15 °C, 101325 Pa, 450 ppm CO2")
plt.plot(λ, n2-1, label="50% humidity")
plt.plot(λ, n3-1, label="370 ppm CO2")
plt.plot(λ, n4-1, label="300K (26.75 °C)")
plt.xlabel('Wavelength (μm)')
plt.ylabel('n-1')
plt.legend()
t = np.arange(-40, 100.1, 1)
n5 = [None] * len(t)
for i in range(0, len(t)):
n5[i] = n(0.6328,t[i],101325,0,450)-1 #dry air, 450 ppm @ HeNe wavelength
plt.figure(2)
plt.plot(t, n5, label="dry air, 101325 Pa, 450 ppm CO2, 632.8 nm")
plt.xlabel('Temperature (°C)')
plt.ylabel('n-1')
plt.legend()
p = np.arange(80000, 120001, 250)
n6 = n(0.6328,15,p,0,450) #dry air, 15 °, 450 ppm @ HeNe wavelength
plt.figure(3)
plt.plot(p, n6-1, label="dry air, 15 °, 450 ppm CO2, 632.8 nm")
plt.xlabel('Pressure (Pa)')
plt.ylabel('n-1')
plt.legend()
h = np.arange(0, 1.001, 0.01)
n7 = n(0.6328, 15, 101325, h, 450) #dry air, 15 °, 450 ppm @ HeNe wavelength
plt.figure(4)
plt.plot(h*100, n7-1, label="15 °, 101325 Pa, 450 ppm CO2, 632.8 nm")
plt.xlabel('Humidity (%)')
plt.ylabel('n-1')
plt.legend()
# --- the formula doesn't depend on xc: CO2 concentration, so I comment the code
# xc = np.arange(0, 2001, 100)
# n8 = n(0.6328, 15, 101325, 0, xc) #dry air, 15 °, 450 ppm @ HeNe wavelength
# plt.figure(5)
# plt.plot(xc, n8-1, label="dry air, 15 °, 101325 Pa, 632.8 nm")
# plt.xlabel('CO2 concentration (ppm)')
# plt.ylabel('n-1')
# plt.legend()
plt.show()
（Ciddor 1996: n 0.23-1.690 μm）
（Ciddor 1996: n 0.23-1.690 μm） 论文是这篇 Ciddor 1996, https://doi.org/10.1364/AO.35.001566
MATLAB
代码从fileExchange下的。
% 下面这些入参的含义：
% 此颜色的光在真空中波长是633$\mu m$，$CO_2$的含量为$450 \mu mol/mol$即450ppm ， 温度是20℃，大气压强是101325Pa，相对湿度RH(%)是0
>> format LongG
>> ciddor(0.633, 450, 20, 101325, 0)
ans =
1.00027179983163
function [n] = ciddor(w, c, t, p, h)
% ciddor: Calculates refractive index using ciddor equation
% Dr Jody Muelaner 2015 jody@muelaner.com
% Department of Mechanical Engineering, The University of Bath
%
% [n] = ciddor(w, c, t, p, h)
% OUTPUTS:
% * n = refractive index (dimensionless)s
% INPUTS:
% * w = vacuum wavelength of laser (micron) valid from 0.3 to 1.7
% * c = CO2 Level (ppm mole) valid from 0 to 2000, typically 450
% * t = Air Temperature (C) valid from -40 to 100
% * p = Air Pressure (Pa) valid from 1E4 to 14E4 typically 101325
% * h = Relative Humidity (%) valid from 0 to 100
% * All inputs may be either scalars or column vectors
% * All inputs must have the same dimensions
%
% Sub-functions:
% * humidity_as_mole_fraction function
% * Saturation_Vapor_Pressure
% * Saturation_Vapor_Pressure_ice
%
% Validaton:
% Validated against the NIST Engineering Metrology Toolbox
% to 9 decimal places i.e n=1.000296813
% emtoolbox.nist.gov/
%% Validate inputs (this section can be removed if not required)
%Find size of w
sizew = size(w); rowsw = sizew(1,1); colsw = sizew(1,2);
%Find size of c
sizec = size(c); rowsc = sizec(1,1); colsc = sizec(1,2);
%Find size of t
sizet = size(t); rowst = sizet(1,1); colst = sizet(1,2);
%Find size of p
sizep = size(p); rowsp = sizep(1,1); colsp = sizep(1,2);
%Find size of h
sizeh = size(h); rowsh = sizeh(1,1); colsh = sizeh(1,2);
%Check all inputs are scalar or column vectors
if colsw == 1 && colsc == 1 && colst == 1 && colsp == 1 && colsh == 1
%Do nothing
else
errordlg('Error in ciddor function to calculate refractive index: All inputs must be either scalars or column vectors')
RETURN
end
%Check all inputs are same size
if rowsc==rowsw && rowst==rowsw && rowsp==rowsw && rowsh==rowsw
%Do nothing
else
errordlg('Error in ciddor function to calculate refractive index: All inputs must be the same size')
return
end
%Check wavelength is within valid range
lessthan = w <= 1.7;
lessthan = sum(lessthan);
greaterthan = w >= 0.3;
greaterthan = sum(greaterthan);
if greaterthan == rowsw && lessthan == rowsw
%Do nothing
else
errordlg('Refractive index not accurate: Ciddor function requires that wavelength (w) must be between 0.3 and 1.7 (microns)')
end
%Check CO2 level is within valid range
lessthan = c <= 2000;
lessthan = sum(lessthan);
greaterthan = c >= 0;
greaterthan = sum(greaterthan);
if greaterthan == rowsc && lessthan == rowsc
%Do nothing
else
errordlg('Refractive index not accurate: Ciddor function requires that CO2 level (c) must be between 0 and 2000 (ppm mole)')
end
%Check temperature is within valid range
lessthan = t <= 100;
lessthan = sum(lessthan);
greaterthan = t >= -40;
greaterthan = sum(greaterthan);
if greaterthan == rowst && lessthan == rowst
%Do nothing
else
errordlg('Refractive index not accurate: Ciddor function requires that temperature (t) must be between -40 and 100 (deg C)')
end
%Check pressure is within valid range
lessthan = p <= 140000;
lessthan = sum(lessthan);
greaterthan = p >= 10000;
greaterthan = sum(greaterthan);
if greaterthan == rowsp && lessthan == rowsp
%Do nothing
else
errordlg('Refractive index not accurate: Ciddor function requires that Pressure (p) must be between 1E4 to 14E4 (pa)')
end
%Check Relative Humidity is within valid range
lessthan = h <= 100;
lessthan = sum(lessthan);
greaterthan = h >= 0;
greaterthan = sum(greaterthan);
if greaterthan == rowsh && lessthan == rowsh
%Do nothing
else
errordlg('Refractive index not accurate: Ciddor function requires that Relative Humidity (h) must be between 0 and 100 (%)')
end
%% Assign Values to Constants
w0 = 295.235;
w1 = 2.6422;
w2 = -0.03238;
w3 = 0.004028;
k0 = 238.0185;
k1 = 5792105;
k2 = 57.362;
k3 = 167917;
a0 = 1.58123 * (10 ^ (-6));
a1 = -2.9331 * (10 ^ (-8));
a2 = 1.1043 * (10 ^ (-10));
b0 = 5.707 * (10 ^ (-6));
b1 = -2.051 * (10 ^ (-8));
c0 = 1.9898 * (10 ^ (-4));
c1 = -2.376 * (10 ^ (-6));
d = 1.83 * (10 ^ (-11));
e = -0.765 * (10 ^ (-8));
pr1 = 101325;
tr1 = 288.15;
za = 0.9995922115;
r = 8.314472;
mv = 0.018015;
%% Perform Intermediate Calculations
s = 1 ./ (w .^ 2);
ras = (10 ^ (-8)) .* ((k1 ./ (k0 - s)) + (k3 ./ (k2 - s)));
rvs = 1.022 * (10 ^ (-8)) * (w0 + (w1 .* s) + (w2 .* (s .^ 2)) + (w3 .* (s .^ 3)));
ma = 0.0289635 + (1.2011 * (10 ^ (-8)) .* (c - 400));
raxs = ras .* (1 + (5.34 * (10 ^ (-7)) .* (c - 450)));
tk = t + 273.15;
% Calls sub-function to calculate relative humidity as mole fraction
x = humidity_as_mole_fraction(t, p, h);
zm1 = (p ./ tk);
zm2 = a0 + (a1 .* t) + (a2 .* (t .^ 2));
zm3 = (b0 + (b1 .* t)) .* x;
zm4 = (c0 + (c1 .* t)) .* (x .^ 2);
zm5 = ((p ./ tk) .^ 2) .* (d + (e .* (x .^ 2)));
zm = 1 - (zm1 .* (zm2 + zm3 + zm4)) + zm5;
% Calculate Zw
pref2 = 1333;
tref2 = 293.15;
zw1 = (pref2 / tref2);
zw2 = a0 + (a1 * (tref2 - 273.15)) + (a2 * ((tref2 - 273.15) ^ 2));
zw3 = (b0 + (b1 * (tref2 - 273.15))) * 1;
zw4 = (c0 + (c1 * (tref2 - 273.15))) * (1 ^ 2);
zw5 = ((pref2 / tref2) ^ 2) * (d + (e * (1 ^ 2)));
zw = 1 - (zw1 * (zw2 + zw3 + zw4)) + zw5;
rhoaxs = (pr1 .* ma) ./ (za * r * tr1);
rhov = (x .* p .* ma) ./ (zm * r .* tk) .* (1 - 1 .* (1 - mv ./ ma));
rhoa = ((1 - x) .* p .* ma) ./ (zm * r .* tk);
rhovs = (pref2 .* ma ./ (zw * r * tref2)) .* (1 - 1 .* (1 - mv ./ ma));
%% Final Result (the refractive index)
n = 1 + ((rhoa ./ rhoaxs) .* raxs) + ((rhov ./ rhovs) .* rvs);
end
%% Sub Funcition used to cacluate the humidity as a mole fraction
function [hmol] = humidity_as_mole_fraction(t, p, h)
% Calculates humidity_as_mole_fraction
% Dr Jody Muelaner 2015, University of Bath, jody@muelaner.com
% [hmol] = humidity_as_mole_fraction(t, p, h)
% OUTPUTS:
% * hmol = humidity as mole fraction (% mole)
% INPUTS:
% * t = Air Temperature (C)
% * p = Air Pressure (Pa)
% * h = Relative Humidity (%)
%
% Requires sub-functions:
% * Saturation_Vapor_Pressure
% * Saturation_Vapor_Pressure_ice
%% Assign Values to Constants
alpha = 1.00062;
beta = 3.14 * (10 ^ (-8));
gama = 5.6 * (10 ^ (-7));
%% Intermediate Calculations
f = alpha + (beta .* p) + (gama .* (t .^ 2));
%% Determine whether over water or ice and call relevant fuction to calculate the saturation vapor preassure
% Note: This section not vectorized
psv = zeros(max(size(t)),1); %preallocate psv before looping
for i=1:size(t)
if t > 0
psv(i) = Saturation_Vapor_Pressure(t(i));
else
psv(i) = Saturation_Vapor_Pressure_ice(t(i));
end %End if else statement
end %End loop through each temperature value
%% Final Result (the humidity as mole fraction)
hmol = ((h ./ 100) .* f .* psv) ./ p;
end
%% Sub Funcition used to cacluate the saturation vapour pressure over water
function [svp] = Saturation_Vapor_Pressure(t)
% Calculates the saturation vapour pressure over water using IAPWS formula
% Dr Jody Muelaner 2015, University of Bath, jody@muelaner.com
% [svp] = Saturation_Vapor_Pressure_ice(t)
% OUTPUTS:
% * svp = saturation vapour pressure (Pa)
% INPUTS:
% * t = Air Temperature (C)
%% Assign Values to Constants
k1 = 1.16705214528 * (10 ^ 3);
k2 = -7.24213167032 * (10 ^ 5);
k3 = -1.70738469401 * 10;
k4 = 1.20208247025 * (10 ^ 4);
k5 = -3.23255503223 * (10 ^ 6);
k6 = 1.49151086135 * 10;
k7 = -4.82326573616 * (10 ^ 3);
k8 = 4.05113405421 * (10 ^ 5);
k9 = -2.38555575678 * (10 ^ -1);
k10 = 6.50175348448 * (10 ^ 2);
%% Intermediate Calculations
tk = t + 273.15;
omega = tk + (k9 / (tk - k10));
a = (omega ^ 2) + (k1 * omega) + k2;
b = (k3 * (omega ^ 2)) + (k4 * omega) + k5;
c = (k6 * (omega ^ 2)) + (k7 * omega) + k8;
x = -b + sqrt((b ^ 2) - (4 * a * c));
%% Final Result (the saturation vapour pressure)
svp = (10 ^ 6) * (((2 * c) / x) ^ 4);
end
%% Sub Funcition used to cacluate the saturation vapour pressure over ice
function [svp] = Saturation_Vapor_Pressure_ice(t)
% Calculates the saturation vapour pressure over ice using IAPWS formula
% Dr Jody Muelaner 2015, University of Bath, jody@muelaner.com
% [svp] = Saturation_Vapor_Pressure_ice(t)
% OUTPUTS:
% * svp = saturation vapour pressure (?)
% INPUTS:
% * t = Air Temperature (C)
%% Assign Values to Constants
a1 = -13.928169;
a2 = 34.7078238;
%%Intermediate Calculations
tk = t + 273.15;
thi = tk / 273.16;
y = (a1 * (1 - (thi ^ (-1.5)))) + (a2 * (1 - (thi ^ (-1.25))));
%% Final Result (the saturation vapour pressure)
svp = 611.657 * ((exp(1)) ^ y);
end
Python
从refractiveindex.info-scripts仓库复制的。
import numpy as np
# import matplotlib.pyplot as plt
π = np.pi
def Z(T, p, xw): # compressibility
t = T - 273.15
a0 = 1.58123e-6 # K·Pa^-1
a1 = -2.9331e-8 # Pa^-1
a2 = 1.1043e-10 # K^-1·Pa^-1
b0 = 5.707e-6 # K·Pa^-1
b1 = -2.051e-8 # Pa^-1
c0 = 1.9898e-4 # K·Pa^-1
c1 = -2.376e-6 # Pa^-1
d = 1.83e-11 # K^2·Pa^-2
e = -0.765e-8 # K^2·Pa^-2
return 1 - (p / T) * (a0 + a1 * t + a2 * t ** 2 + (b0 + b1 * t) * xw + (c0 + c1 * t) * xw ** 2) + (p / T) ** 2 * (
d + e * xw ** 2)
def n(λ, t, p, h, xc):
# λ: wavelength, 0.3 to 1.69 μm
# t: temperature, -40 to +100 °C
# p: pressure, 80000 to 120000 Pa
# h: fractional humidity, 0 to 1
# xc: CO2 concentration, 0 to 2000 ppm
σ = 1 / λ # μm^-1
T = t + 273.15 # Temperature °C -> K
R = 8.314510 # gas constant, J/(mol·K)
k0 = 238.0185 # μm^-2
k1 = 5792105 # μm^-2
k2 = 57.362 # μm^-2
k3 = 167917 # μm^-2
w0 = 295.235 # μm^-2
w1 = 2.6422 # μm^-2
w2 = -0.032380 # μm^-4
w3 = 0.004028 # μm^-6
A = 1.2378847e-5 # K^-2
B = -1.9121316e-2 # K^-1
C = 33.93711047
D = -6.3431645e3 # K
α = 1.00062
β = 3.14e-8 # Pa^-1,
γ = 5.6e-7 # °C^-2
# saturation vapor pressure of water vapor in air at temperature T
if (t >= 0):
svp = np.exp(A * T ** 2 + B * T + C + D / T) # Pa
else:
svp = 10 ** (-2663.5 / T + 12.537)
# enhancement factor of water vapor in air
f = α + β * p + γ * t ** 2
# molar fraction of water vapor in moist air
xw = f * h * svp / p
# refractive index of standard air at 15 °C, 101325 Pa, 0% humidity, 450 ppm CO2
nas = 1 + (k1 / (k0 - σ ** 2) + k3 / (k2 - σ ** 2)) * 1e-8
# refractive index of standard air at 15 °C, 101325 Pa, 0% humidity, xc ppm CO2
naxs = 1 + (nas - 1) * (1 + 0.534e-6 * (xc - 450))
# refractive index of water vapor at standard conditions (20 °C, 1333 Pa)
nws = 1 + 1.022 * (w0 + w1 * σ ** 2 + w2 * σ ** 4 + w3 * σ ** 6) * 1e-8
Ma = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400)) # molar mass of dry air, kg/mol
Mw = 0.018015 # molar mass of water vapor, kg/mol
Za = Z(288.15, 101325, 0) # compressibility of dry air
Zw = Z(293.15, 1333, 1) # compressibility of pure water vapor
# Eq.4 with (T,P,xw) = (288.15, 101325, 0)
ρaxs = 101325 * Ma / (Za * R * 288.15) # density of standard air
# Eq 4 with (T,P,xw) = (293.15, 1333, 1)
ρws = 1333 * Mw / (Zw * R * 293.15) # density of standard water vapor
# two parts of Eq.4: ρ=ρa+ρw
ρa = p * Ma / (Z(T, p, xw) * R * T) * (1 - xw) # density of the dry component of the moist air
ρw = p * Mw / (Z(T, p, xw) * R * T) * xw # density of the water vapor component
nprop = 1 + (ρa / ρaxs) * (naxs - 1) + (ρw / ρws) * (nws - 1)
return nprop
print("n =",n(λ=0.633, t=20, p=101325, h=0, xc=450))
# n = 1.0002717998316348