yosata/AtmosphericModel
1
0 Fork 0
mirror of https://gitlab.astro-wise.org/micado/atmosphericmodel.git synced 2026-10-29 20:51:05 +00:00
Code Issues Projects Releases Wiki Activity

Added installation functionality and renamed the package

Browse Source
This commit is contained in:
Joost van den Born committed 2020-05-19 12:07:29 +02:00
1 parent 9e0b09787c
commit f66d26c649
16 files changed
+40 -1700
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

.gitignore
+1 -1
View File
@@ -1,3 +1,3 @@
.DS_Store
__pychache__/
.ipynb_checkpoints/
atmospheric refraction example.html
+19 -2
View File
@@ -13077,8 +13077,25 @@ div#notebook {
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="Example:-Calculating-atmospheric-dispersion-for-H-band">Example: Calculating atmospheric dispersion for H-band<a class="anchor-link" href="#Example:-Calculating-atmospheric-dispersion-for-H-band">&#182;</a></h1><p>Using atmosphericModel.py is made as simple as importing the python file and calling its functions. Several parameters are needed, such as location, temperature, pressure and relative humidity. Then the dispersion model object is initiated, allowing us to calculate the index of refraction for any set of TPH and CO2 content values.</p>
<h1 id="Example:-Calculating-atmospheric-dispersion-for-H-band">Example: Calculating atmospheric dispersion for H-band<a class="anchor-link" href="#Example:-Calculating-atmospheric-dispersion-for-H-band">&#182;</a></h1><p>Using atmosphericModel.py is made as simple as importing the python package and calling its functions. Several parameters are needed, such as location, temperature, pressure and relative humidity. Then the dispersion model object is initiated, allowing us to calculate the index of refraction for any set of TPH and CO2 content values.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="Installation">Installation<a class="anchor-link" href="#Installation">&#182;</a></h3><p>We install the package using the <code>setup.py</code> file and running the following command:
<code>python setup.py install</code></p>
<p>Now the package can be imported using <code>import AstroAtmosphere</code> in Python.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="How-it-works">How it works<a class="anchor-link" href="#How-it-works">&#182;</a></h3>
</div>
</div>
</div>
@@ -13090,7 +13107,7 @@ div#notebook {
<div class=" highlight hl-ipython3"><pre><span></span><span class="ch">#!/usr/bin/env Python</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">from</span> <span class="nn">atmosphericModel</span> <span class="k">import</span> <span class="o">*</span>
<span class="kn">from</span> <span class="nn">AstroAtmosphere</span> <span class="k">import</span> <span class="o">*</span>
<span class="c1"># Parameters at Cerro Armazones</span>
atmospheric refraction example.ipynb
+20 -2
View File
@@ -5,7 +5,25 @@
"metadata": {},
"source": [
"# Example: Calculating atmospheric dispersion for H-band\n",
"Using atmosphericModel.py is made as simple as importing the python file and calling its functions. Several parameters are needed, such as location, temperature, pressure and relative humidity. Then the dispersion model object is initiated, allowing us to calculate the index of refraction for any set of TPH and CO2 content values."
"Using atmosphericModel.py is made as simple as importing the python package and calling its functions. Several parameters are needed, such as location, temperature, pressure and relative humidity. Then the dispersion model object is initiated, allowing us to calculate the index of refraction for any set of TPH and CO2 content values."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Installation\n",
"We install the package using the `setup.py` file and running the following command:\n",
"`python setup.py install`\n",
"\n",
"Now the package can be imported using `import AstroAtmosphere` in Python."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### How it works"
]
},
{
@@ -17,7 +35,7 @@
"#!/usr/bin/env Python\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from atmosphericModel import *\n",
"from AstroAtmosphere import *\n",
"\n",
"\n",
"# Parameters at Cerro Armazones\n",
atmosphericModel/1976USSA.txt
-179
View File
@@ -1,179 +0,0 @@
alt sigma delta theta temp press dens a visc k.visc
Km K N/sq.m kg/cu.m m/sec kg/m-s sq.m/s
-2.0 1.2067E+00 1.2611E+00 1.0451 301.2 1.27783E+05 1.478E+00 347.9 18.51 1.25E-05
-1.5 1.1522E+00 1.1912E+00 1.0338 297.9 1.20696E+05 1.411E+00 346.0 18.36 1.30E-05
-1.0 1.0996E+00 1.1244E+00 1.0226 294.7 1.13931E+05 1.347E+00 344.1 18.21 1.35E-05
-0.5 1.0489E+00 1.0607E+00 1.0113 291.4 1.07478E+05 1.285E+00 342.2 18.05 1.40E-05
0.0 1.0000E+00 1.0000E+00 1.0000 288.1 1.01325E+05 1.225E+00 340.3 17.89 1.46E-05
0.5 9.5288E-01 9.4213E-01 0.9887 284.9 9.54613E+04 1.167E+00 338.4 17.74 1.52E-05
1.0 9.0748E-01 8.8701E-01 0.9774 281.7 8.98763E+04 1.112E+00 336.4 17.58 1.58E-05
1.5 8.6376E-01 8.3454E-01 0.9662 278.4 8.45597E+04 1.058E+00 334.5 17.42 1.65E-05
2.0 8.2168E-01 7.8462E-01 0.9549 275.2 7.95014E+04 1.007E+00 332.5 17.26 1.71E-05
2.5 7.8119E-01 7.3715E-01 0.9436 271.9 7.46917E+04 9.570E-01 330.6 17.10 1.79E-05
3.0 7.4225E-01 6.9204E-01 0.9324 268.7 7.01211E+04 9.093E-01 328.6 16.94 1.86E-05
3.5 7.0482E-01 6.4920E-01 0.9211 265.4 6.57804E+04 8.634E-01 326.6 16.78 1.94E-05
4.0 6.6885E-01 6.0854E-01 0.9098 262.2 6.16604E+04 8.193E-01 324.6 16.61 2.03E-05
4.5 6.3432E-01 5.6997E-01 0.8986 258.9 5.77525E+04 7.770E-01 322.6 16.45 2.12E-05
5.0 6.0117E-01 5.3341E-01 0.8873 255.7 5.40482E+04 7.364E-01 320.5 16.28 2.21E-05
5.5 5.6936E-01 4.9878E-01 0.8760 252.4 5.05393E+04 6.975E-01 318.5 16.12 2.31E-05
6.0 5.3887E-01 4.6600E-01 0.8648 249.2 4.72176E+04 6.601E-01 316.5 15.95 2.42E-05
6.5 5.0964E-01 4.3499E-01 0.8535 245.9 4.40754E+04 6.243E-01 314.4 15.78 2.53E-05
7.0 4.8165E-01 4.0568E-01 0.8423 242.7 4.11052E+04 5.900E-01 312.3 15.61 2.65E-05
7.5 4.5485E-01 3.7799E-01 0.8310 239.5 3.82996E+04 5.572E-01 310.2 15.44 2.77E-05
8.0 4.2921E-01 3.5185E-01 0.8198 236.2 3.56515E+04 5.258E-01 308.1 15.27 2.90E-05
8.5 4.0470E-01 3.2721E-01 0.8085 233.0 3.31541E+04 4.958E-01 306.0 15.10 3.05E-05
9.0 3.8128E-01 3.0398E-01 0.7973 229.7 3.08006E+04 4.671E-01 303.8 14.93 3.20E-05
9.5 3.5891E-01 2.8211E-01 0.7860 226.5 2.85846E+04 4.397E-01 301.7 14.75 3.36E-05
10.0 3.3756E-01 2.6153E-01 0.7748 223.3 2.64998E+04 4.135E-01 299.5 14.58 3.53E-05
10.5 3.1720E-01 2.4219E-01 0.7635 220.0 2.45401E+04 3.886E-01 297.4 14.40 3.71E-05
11.0 2.9780E-01 2.2403E-01 0.7523 216.8 2.26998E+04 3.648E-01 295.2 14.22 3.90E-05
11.5 2.7545E-01 2.0710E-01 0.7519 216.7 2.09847E+04 3.374E-01 295.1 14.22 4.21E-05
12.0 2.5464E-01 1.9146E-01 0.7519 216.7 1.93993E+04 3.119E-01 295.1 14.22 4.56E-05
12.5 2.3541E-01 1.7699E-01 0.7519 216.7 1.79340E+04 2.884E-01 295.1 14.22 4.93E-05
13.0 2.1763E-01 1.6363E-01 0.7519 216.7 1.65795E+04 2.666E-01 295.1 14.22 5.33E-05
13.5 2.0119E-01 1.5127E-01 0.7519 216.7 1.53275E+04 2.465E-01 295.1 14.22 5.77E-05
14.0 1.8600E-01 1.3985E-01 0.7519 216.7 1.41702E+04 2.279E-01 295.1 14.22 6.24E-05
14.5 1.7196E-01 1.2929E-01 0.7519 216.7 1.31005E+04 2.107E-01 295.1 14.22 6.75E-05
15.0 1.5898E-01 1.1953E-01 0.7519 216.7 1.21117E+04 1.948E-01 295.1 14.22 7.30E-05
15.5 1.4698E-01 1.1051E-01 0.7519 216.7 1.11976E+04 1.801E-01 295.1 14.22 7.90E-05
16.0 1.3589E-01 1.0217E-01 0.7519 216.7 1.03527E+04 1.665E-01 295.1 14.22 8.54E-05
16.5 1.2564E-01 9.4465E-02 0.7519 216.7 9.57164E+03 1.539E-01 295.1 14.22 9.24E-05
17.0 1.1616E-01 8.7339E-02 0.7519 216.7 8.84961E+03 1.423E-01 295.1 14.22 9.99E-05
17.5 1.0740E-01 8.0752E-02 0.7519 216.7 8.18215E+03 1.316E-01 295.1 14.22 1.08E-04
18.0 9.9302E-02 7.4662E-02 0.7519 216.7 7.56512E+03 1.216E-01 295.1 14.22 1.17E-04
18.5 9.1815E-02 6.9032E-02 0.7519 216.7 6.99471E+03 1.125E-01 295.1 14.22 1.26E-04
19.0 8.4893E-02 6.3828E-02 0.7519 216.7 6.46738E+03 1.040E-01 295.1 14.22 1.37E-04
19.5 7.8494E-02 5.9017E-02 0.7519 216.7 5.97989E+03 9.616E-02 295.1 14.22 1.48E-04
20.0 7.2578E-02 5.4569E-02 0.7519 216.7 5.52921E+03 8.891E-02 295.1 14.22 1.60E-04
20.5 6.6979E-02 5.0460E-02 0.7534 217.1 5.11291E+03 8.205E-02 295.4 14.24 1.74E-04
21.0 6.1807E-02 4.6670E-02 0.7551 217.6 4.72884E+03 7.571E-02 295.7 14.27 1.88E-04
21.5 5.7045E-02 4.3173E-02 0.7568 218.1 4.37446E+03 6.988E-02 296.0 14.29 2.05E-04
22.0 5.2660E-02 3.9945E-02 0.7585 218.6 4.04741E+03 6.451E-02 296.4 14.32 2.22E-04
22.5 4.8622E-02 3.6965E-02 0.7603 219.1 3.74551E+03 5.956E-02 296.7 14.35 2.41E-04
23.0 4.4902E-02 3.4214E-02 0.7620 219.6 3.46678E+03 5.500E-02 297.0 14.38 2.61E-04
23.5 4.1474E-02 3.1674E-02 0.7637 220.1 3.20939E+03 5.081E-02 297.4 14.40 2.83E-04
24.0 3.8316E-02 2.9328E-02 0.7654 220.6 2.97167E+03 4.694E-02 297.7 14.43 3.07E-04
24.5 3.5404E-02 2.7161E-02 0.7672 221.1 2.75206E+03 4.337E-02 298.1 14.46 3.33E-04
25.0 3.2721E-02 2.5158E-02 0.7689 221.6 2.54915E+03 4.008E-02 298.4 14.48 3.61E-04
25.5 3.0246E-02 2.3308E-02 0.7706 222.0 2.36163E+03 3.705E-02 298.7 14.51 3.92E-04
26.0 2.7964E-02 2.1597E-02 0.7723 222.5 2.18831E+03 3.426E-02 299.1 14.54 4.24E-04
26.5 2.5858E-02 2.0016E-02 0.7740 223.0 2.02807E+03 3.168E-02 299.4 14.57 4.60E-04
27.0 2.3916E-02 1.8553E-02 0.7758 223.5 1.87991E+03 2.930E-02 299.7 14.59 4.98E-04
27.5 2.2124E-02 1.7201E-02 0.7775 224.0 1.74289E+03 2.710E-02 300.1 14.62 5.39E-04
28.0 2.0470E-02 1.5950E-02 0.7792 224.5 1.61614E+03 2.508E-02 300.4 14.65 5.84E-04
28.5 1.8943E-02 1.4793E-02 0.7809 225.0 1.49888E+03 2.320E-02 300.7 14.67 6.32E-04
29.0 1.7533E-02 1.3722E-02 0.7826 225.5 1.39037E+03 2.148E-02 301.0 14.70 6.84E-04
29.5 1.6231E-02 1.2731E-02 0.7844 226.0 1.28994E+03 1.988E-02 301.4 14.73 7.41E-04
30.0 1.5028E-02 1.1813E-02 0.7861 226.5 1.19698E+03 1.841E-02 301.7 14.75 8.01E-04
30.5 1.3917E-02 1.0964E-02 0.7878 227.0 1.11092E+03 1.705E-02 302.0 14.78 8.67E-04
31.0 1.2891E-02 1.0177E-02 0.7895 227.5 1.03122E+03 1.579E-02 302.4 14.81 9.38E-04
31.5 1.1942E-02 9.4488E-03 0.7912 228.0 9.57402E+02 1.463E-02 302.7 14.83 1.01E-03
32.0 1.1065E-02 8.7740E-03 0.7930 228.5 8.89024E+02 1.355E-02 303.0 14.86 1.10E-03
32.5 1.0228E-02 8.1493E-03 0.7968 229.6 8.25726E+02 1.253E-02 303.8 14.92 1.19E-03
33.0 9.4469E-03 7.5724E-03 0.8016 231.0 7.67272E+02 1.157E-02 304.7 14.99 1.30E-03
33.5 8.7297E-03 7.0395E-03 0.8064 232.4 7.13279E+02 1.069E-02 305.6 15.07 1.41E-03
34.0 8.0709E-03 6.5470E-03 0.8112 233.7 6.63380E+02 9.887E-03 306.5 15.14 1.53E-03
34.5 7.4654E-03 6.0917E-03 0.8160 235.1 6.17243E+02 9.145E-03 307.4 15.21 1.66E-03
35.0 6.9085E-03 5.6705E-03 0.8208 236.5 5.74564E+02 8.463E-03 308.3 15.29 1.81E-03
35.5 6.3961E-03 5.2807E-03 0.8256 237.9 5.35066E+02 7.835E-03 309.2 15.36 1.96E-03
36.0 5.9245E-03 4.9198E-03 0.8304 239.3 4.98495E+02 7.257E-03 310.1 15.43 2.13E-03
36.5 5.4901E-03 4.5854E-03 0.8352 240.7 4.64618E+02 6.725E-03 311.0 15.51 2.31E-03
37.0 5.0899E-03 4.2756E-03 0.8400 242.1 4.33223E+02 6.235E-03 311.9 15.58 2.50E-03
37.5 4.7209E-03 3.9883E-03 0.8448 243.4 4.04115E+02 5.783E-03 312.8 15.65 2.71E-03
38.0 4.3806E-03 3.7218E-03 0.8496 244.8 3.77116E+02 5.366E-03 313.7 15.72 2.93E-03
38.5 4.0666E-03 3.4746E-03 0.8544 246.2 3.52061E+02 4.982E-03 314.6 15.79 3.17E-03
39.0 3.7767E-03 3.2450E-03 0.8592 247.6 3.28801E+02 4.626E-03 315.4 15.87 3.43E-03
39.5 3.5089E-03 3.0318E-03 0.8640 249.0 3.07198E+02 4.298E-03 316.3 15.94 3.71E-03
40.0 3.2616E-03 2.8337E-03 0.8688 250.4 2.87125E+02 3.995E-03 317.2 16.01 4.01E-03
40.5 3.0329E-03 2.6496E-03 0.8736 251.7 2.68467E+02 3.715E-03 318.1 16.08 4.33E-03
41.0 2.8214E-03 2.4783E-03 0.8784 253.1 2.51116E+02 3.456E-03 318.9 16.15 4.67E-03
41.5 2.6257E-03 2.3190E-03 0.8832 254.5 2.34974E+02 3.216E-03 319.8 16.22 5.04E-03
42.0 2.4445E-03 2.1708E-03 0.8880 255.9 2.19952E+02 2.995E-03 320.7 16.29 5.44E-03
42.5 2.2768E-03 2.0327E-03 0.8928 257.3 2.05965E+02 2.789E-03 321.5 16.36 5.87E-03
43.0 2.1214E-03 1.9041E-03 0.8976 258.6 1.92937E+02 2.599E-03 322.4 16.43 6.32E-03
43.5 1.9774E-03 1.7843E-03 0.9024 260.0 1.80799E+02 2.422E-03 323.3 16.50 6.81E-03
44.0 1.8438E-03 1.6727E-03 0.9072 261.4 1.69484E+02 2.259E-03 324.1 16.57 7.34E-03
44.5 1.7199E-03 1.5685E-03 0.9120 262.8 1.58932E+02 2.107E-03 325.0 16.64 7.90E-03
45.0 1.6050E-03 1.4714E-03 0.9168 264.2 1.49090E+02 1.966E-03 325.8 16.71 8.50E-03
45.5 1.4983E-03 1.3808E-03 0.9216 265.5 1.39905E+02 1.835E-03 326.7 16.78 9.14E-03
46.0 1.3992E-03 1.2961E-03 0.9263 266.9 1.31330E+02 1.714E-03 327.5 16.85 9.83E-03
46.5 1.3071E-03 1.2171E-03 0.9311 268.3 1.23322E+02 1.601E-03 328.4 16.92 1.06E-02
47.0 1.2215E-03 1.1433E-03 0.9359 269.7 1.15841E+02 1.496E-03 329.2 16.99 1.14E-02
47.5 1.1437E-03 1.0743E-03 0.9393 270.6 1.08849E+02 1.401E-03 329.8 17.04 1.22E-02
48.0 1.0748E-03 1.0095E-03 0.9393 270.6 1.02287E+02 1.317E-03 329.8 17.04 1.29E-02
48.5 1.0100E-03 9.4865E-04 0.9393 270.6 9.61215E+01 1.237E-03 329.8 17.04 1.38E-02
49.0 9.4912E-04 8.9147E-04 0.9393 270.6 9.03286E+01 1.163E-03 329.8 17.04 1.47E-02
49.5 8.9193E-04 8.3776E-04 0.9393 270.6 8.48858E+01 1.093E-03 329.8 17.04 1.56E-02
50.0 8.3819E-04 7.8728E-04 0.9393 270.6 7.97716E+01 1.027E-03 329.8 17.04 1.66E-02
50.5 7.8770E-04 7.3986E-04 0.9393 270.6 7.49663E+01 9.649E-04 329.8 17.04 1.77E-02
51.0 7.4026E-04 6.9530E-04 0.9393 270.6 7.04511E+01 9.068E-04 329.8 17.04 1.88E-02
51.5 6.9630E-04 6.5342E-04 0.9384 270.4 6.62082E+01 8.530E-04 329.7 17.02 2.00E-02
52.0 6.5758E-04 6.1395E-04 0.9336 269.0 6.22084E+01 8.055E-04 328.8 16.96 2.10E-02
52.5 6.2085E-04 5.7668E-04 0.9289 267.7 5.84321E+01 7.605E-04 328.0 16.89 2.22E-02
53.0 5.8599E-04 5.4150E-04 0.9241 266.3 5.48678E+01 7.178E-04 327.1 16.82 2.34E-02
53.5 5.5293E-04 5.0831E-04 0.9193 264.9 5.15047E+01 6.773E-04 326.3 16.75 2.47E-02
54.0 5.2158E-04 4.7700E-04 0.9145 263.5 4.83322E+01 6.389E-04 325.4 16.68 2.61E-02
54.5 4.9187E-04 4.4748E-04 0.9098 262.1 4.53405E+01 6.025E-04 324.6 16.61 2.76E-02
55.0 4.6370E-04 4.1964E-04 0.9050 260.8 4.25201E+01 5.680E-04 323.7 16.54 2.91E-02
55.5 4.3702E-04 3.9341E-04 0.9002 259.4 3.98620E+01 5.354E-04 322.9 16.47 3.08E-02
56.0 4.1175E-04 3.6869E-04 0.8954 258.0 3.73576E+01 5.044E-04 322.0 16.40 3.25E-02
56.5 3.8782E-04 3.4541E-04 0.8907 256.6 3.49988E+01 4.751E-04 321.2 16.33 3.44E-02
57.0 3.6516E-04 3.2349E-04 0.8859 255.3 3.27778E+01 4.473E-04 320.3 16.26 3.64E-02
57.5 3.4373E-04 3.0286E-04 0.8811 253.9 3.06871E+01 4.211E-04 319.4 16.19 3.85E-02
58.0 3.2344E-04 2.8344E-04 0.8763 252.5 2.87198E+01 3.962E-04 318.6 16.12 4.07E-02
58.5 3.0426E-04 2.6518E-04 0.8716 251.1 2.68692E+01 3.727E-04 317.7 16.05 4.31E-02
59.0 2.8612E-04 2.4800E-04 0.8668 249.8 2.51289E+01 3.505E-04 316.8 15.98 4.56E-02
59.5 2.6897E-04 2.3186E-04 0.8620 248.4 2.34929E+01 3.295E-04 315.9 15.91 4.83E-02
60.0 2.5276E-04 2.1668E-04 0.8573 247.0 2.19554E+01 3.096E-04 315.1 15.84 5.11E-02
60.5 2.3746E-04 2.0243E-04 0.8525 245.6 2.05110E+01 2.909E-04 314.2 15.77 5.42E-02
61.0 2.2300E-04 1.8904E-04 0.8477 244.3 1.91546E+01 2.732E-04 313.3 15.69 5.75E-02
61.5 2.0935E-04 1.7647E-04 0.8430 242.9 1.78811E+01 2.565E-04 312.4 15.62 6.09E-02
62.0 1.9647E-04 1.6468E-04 0.8382 241.5 1.66860E+01 2.407E-04 311.5 15.55 6.46E-02
62.5 1.8431E-04 1.5361E-04 0.8334 240.2 1.55648E+01 2.258E-04 310.7 15.48 6.86E-02
63.0 1.7285E-04 1.4323E-04 0.8287 238.8 1.45132E+01 2.117E-04 309.8 15.41 7.28E-02
63.5 1.6204E-04 1.3351E-04 0.8239 237.4 1.35275E+01 1.985E-04 308.9 15.33 7.72E-02
64.0 1.5185E-04 1.2439E-04 0.8191 236.0 1.26036E+01 1.860E-04 308.0 15.26 8.20E-02
64.5 1.4225E-04 1.1585E-04 0.8144 234.7 1.17382E+01 1.743E-04 307.1 15.19 8.72E-02
65.0 1.3321E-04 1.0785E-04 0.8096 233.3 1.09277E+01 1.632E-04 306.2 15.12 9.26E-02
65.5 1.2470E-04 1.0036E-04 0.8048 231.9 1.01690E+01 1.528E-04 305.3 15.04 9.85E-02
66.0 1.1668E-04 9.3354E-05 0.8001 230.5 9.45909E+00 1.429E-04 304.4 14.97 1.05E-01
66.5 1.0914E-04 8.6800E-05 0.7953 229.2 8.79502E+00 1.337E-04 303.5 14.90 1.11E-01
67.0 1.0204E-04 8.0672E-05 0.7906 227.8 8.17408E+00 1.250E-04 302.6 14.82 1.19E-01
67.5 9.5371E-05 7.4944E-05 0.7858 226.4 7.59372E+00 1.168E-04 301.7 14.75 1.26E-01
68.0 8.9101E-05 6.9593E-05 0.7811 225.1 7.05149E+00 1.091E-04 300.7 14.67 1.34E-01
68.5 8.3209E-05 6.4595E-05 0.7763 223.7 6.54509E+00 1.019E-04 299.8 14.60 1.43E-01
69.0 7.7675E-05 5.9929E-05 0.7715 222.3 6.07234E+00 9.515E-05 298.9 14.53 1.53E-01
69.5 7.2478E-05 5.5576E-05 0.7668 221.0 5.63120E+00 8.879E-05 298.0 14.45 1.63E-01
70.0 6.7601E-05 5.1515E-05 0.7620 219.6 5.21972E+00 8.281E-05 297.1 14.38 1.74E-01
70.5 6.3026E-05 4.7728E-05 0.7573 218.2 4.83607E+00 7.721E-05 296.1 14.30 1.85E-01
71.0 5.8734E-05 4.4199E-05 0.7525 216.8 4.47851E+00 7.195E-05 295.2 14.23 1.98E-01
71.5 5.4712E-05 4.0912E-05 0.7478 215.5 4.14543E+00 6.702E-05 294.3 14.15 2.11E-01
72.0 5.0905E-05 3.7852E-05 0.7436 214.3 3.83531E+00 6.236E-05 293.4 14.08 2.26E-01
72.5 4.7296E-05 3.5007E-05 0.7402 213.3 3.54709E+00 5.794E-05 292.8 14.03 2.42E-01
73.0 4.3927E-05 3.2365E-05 0.7368 212.3 3.27939E+00 5.381E-05 292.1 13.98 2.60E-01
73.5 4.0786E-05 2.9912E-05 0.7334 211.3 3.03083E+00 4.996E-05 291.4 13.92 2.79E-01
74.0 3.7856E-05 2.7635E-05 0.7300 210.3 2.80012E+00 4.637E-05 290.7 13.87 2.99E-01
74.5 3.5125E-05 2.5522E-05 0.7266 209.4 2.58605E+00 4.303E-05 290.1 13.81 3.21E-01
75.0 3.2580E-05 2.3563E-05 0.7232 208.4 2.38749E+00 3.991E-05 289.4 13.76 3.45E-01
75.5 3.0209E-05 2.1746E-05 0.7198 207.4 2.20338E+00 3.701E-05 288.7 13.70 3.70E-01
76.0 2.8001E-05 2.0061E-05 0.7164 206.4 2.03272E+00 3.430E-05 288.0 13.65 3.98E-01
76.5 2.5946E-05 1.8501E-05 0.7131 205.5 1.87458E+00 3.178E-05 287.4 13.59 4.28E-01
77.0 2.4033E-05 1.7055E-05 0.7097 204.5 1.72810E+00 2.944E-05 286.7 13.54 4.60E-01
77.5 2.2253E-05 1.5716E-05 0.7063 203.5 1.59247E+00 2.726E-05 286.0 13.48 4.95E-01
78.0 2.0597E-05 1.4477E-05 0.7029 202.5 1.46692E+00 2.523E-05 285.3 13.43 5.32E-01
78.5 1.9058E-05 1.3331E-05 0.6995 201.6 1.35076E+00 2.335E-05 284.6 13.37 5.73E-01
79.0 1.7627E-05 1.2270E-05 0.6961 200.6 1.24331E+00 2.159E-05 283.9 13.32 6.17E-01
79.5 1.6298E-05 1.1290E-05 0.6927 199.6 1.14396E+00 1.996E-05 283.2 13.26 6.64E-01
80.0 1.5063E-05 1.0384E-05 0.6893 198.6 1.05213E+00 1.845E-05 282.5 13.21 7.16E-01
80.5 1.3917E-05 9.5464E-06 0.6860 197.7 9.67291E-01 1.705E-05 281.8 13.15 7.71E-01
81.0 1.2853E-05 8.7731E-06 0.6826 196.7 8.88933E-01 1.574E-05 281.1 13.10 8.32E-01
81.5 1.1866E-05 8.0591E-06 0.6792 195.7 8.16590E-01 1.454E-05 280.4 13.04 8.97E-01
82.0 1.0950E-05 7.4002E-06 0.6758 194.7 7.49827E-01 1.341E-05 279.7 12.98 9.68E-01
82.5 1.0101E-05 6.7924E-06 0.6724 193.8 6.88237E-01 1.237E-05 279.0 12.93 1.04E+00
83.0 9.3145E-06 6.2318E-06 0.6690 192.8 6.31441E-01 1.141E-05 278.3 12.87 1.13E+00
83.5 8.5856E-06 5.7151E-06 0.6657 191.8 5.79087E-01 1.052E-05 277.6 12.82 1.22E+00
84.0 7.9106E-06 5.2391E-06 0.6623 190.8 5.30848E-01 9.690E-06 276.9 12.76 1.32E+00
84.5 7.2857E-06 4.8006E-06 0.6589 189.9 4.86416E-01 8.925E-06 276.2 12.70 1.42E+00
85.0 6.7074E-06 4.3968E-06 0.6555 188.9 4.45509E-01 8.217E-06 275.5 12.65 1.54E+00
85.5 6.1724E-06 4.0253E-06 0.6521 187.9 4.07863E-01 7.561E-06 274.8 12.59 1.67E+00
86.0 5.6777E-06 3.6835E-06 0.6488 186.9 3.73232E-01 6.955E-06 274.1 12.53 1.80E+00
atmosphericModel/__init__.py
-10
View File
@@ -1,10 +0,0 @@
#!/usr/bin/env Python3
import os
from .ciddorModel import *
from .dispersionModels import *
from .refractionModels import *
from .refractionModels import *
from .examples import *
from .refractivityModels import *
atmosphericModel/__pycache__/__init__.cpython-37.pyc
BIN
View File
Binary file not shown.
atmosphericModel/__pycache__/ciddorModel.cpython-37.pyc
BIN
View File
Binary file not shown.
atmosphericModel/__pycache__/dispersionModels.cpython-37.pyc
BIN
View File
Binary file not shown.
atmosphericModel/__pycache__/examples.cpython-37.pyc
BIN
View File
Binary file not shown.
atmosphericModel/__pycache__/refractionModels.cpython-37.pyc
BIN
View File
Binary file not shown.
atmosphericModel/__pycache__/refractivityModels.cpython-37.pyc
BIN
View File
Binary file not shown.
atmosphericModel/ciddorModel.py
-385
View File
@@ -1,385 +0,0 @@
#!/usr/bin/env Python3
import numpy as np
from scipy import integrate
'''
@ Author: Joost van den Born
@ Contact: born@astron.nl
@ Description: This script offers a complete model of the atmospheric refraction
and dispersion. It uses the Ciddor1996 paper to calculate the
index of refraction for a wide range of conditions, with corresponding
error propagation if needed.
The actual refraction is calculated by any of the following models:
- Plane parallel atmosphere
- Cassini spherical homogeneous atmosphere
- Other methods described in Corbard et al., 2019.
'''
class Observatory:
def __init__(self, dT=0.2, dP=20, dH=0.02, dCO2=20, dl=0.001):
self.dT = dT # in K
self.dP = dP # in Pa
self.dH = dH # relative humidity between 0 and 1
self.dCO2 = dCO2 # in ppm
self.dl = dl # in micron
def n_as(self, l):
'''
Formula that calculate the refractive index for standard atmospheric
conditions, i.e. 15 C, 101,325 Pa, 0% RH, 450 ppm CO2.
It returns n_as, not (n_as - 1) * 1e8!
'''
_s = 1 / l
_k0 = 238.0185 # um^-2
_k1 = 5792105 # um^-2
_k2 = 57.362 # um^-2
_k3 = 167917 # um^-2
return (_k1 / (_k0 - _s**2) + _k3 / (_k2 - _s**2)) / 1e8 + 1
def n_axs(self, l, xc):
'''
Formula that calculate the refractive index for standard atmospheric
conditions, i.e. 15 C, 101,325 Pa, 0% RH, but with xc ppm of CO2.
It returns n_axs, not (n_axs - 1)
'''
_n_as = self.n_as(l)
return ( (_n_as - 1) * (1 + 0.534e-6 * (xc - 450) ) ) + 1
def n_ws(self, l):
'''
Formula that calculates the refractivity of air for standard moist conditions
as defined by B&S (1939), i.e. 20 C, 1333 Pa
Returns n_ws, not (n_ws - 1) * 1e8
'''
_s = 1 / l
_cf = 1.022
_w0 = 295.235 # um^-2
_w1 = 2.6422 # um^-2
_w2 = -0.032380 # um^-4
_w3 = 0.004028 # um^-6
return (_cf * (_w0 + _w1 * _s**2 + _w2 * _s**4 + _w3 * _s**6) ) / 1e8 + 1
def f(self, p, T):
'''
Enhancement factor of water vapour in air
t in degrees celsius
p in Pascal (Pa)
'''
_t = T - 273.15
_a = 1.00062
_b = 3.14e-8 # Pa^-1
_c = 5.6e-7 # deg C^-2
return _a + _b * p + _c * _t**2
def svp(self, T):
'''
Saturation vapour pressure in Pa.
'''
_A = 1.2378847e-5 # K^-2
_B = -1.9121316e-2 # K^-2
_C = 33.93711047
_D = -6.3431645e3 # K^-2
return np.exp(_A*T**2 + _B*T + _C + _D/T)
def svp2(self, T):
'''
Calculation of svp from IAPWS formula, based on documentation by NIST toolbox.
'''
K1 = 1.16705214528e+03
K2 = -7.24213167032e+05
K3 = -1.70738469401e+01
K4 = 1.20208247025e+04
K5 = -3.23255503223e+06
K6 = 1.49151086135e+01
K7 = -4.82326573616e+03
K8 = 4.05113405421e+05
K9 = -2.38555575678e-01
K10 = 6.50175348448e+02
_W = T + K9 / (T - K10)
_A = _W**2 + K1*_W + K2
_B = K3*_W**2 + K4*_W + K5
_C = K6*_W**2 + K7*_W + K8
_X = -_B + np.sqrt(_B**2 - 4*_A*_C)
return (2*_C / _X)**4 * 1e6
def xw(self, p, T, RH):
'''
Molar fraction of water vapour in moist air.
'''
_f = self.f(p, T)
_svp = self.svp(T)
return _f * RH * _svp / p
def Z(self, p, T, RH, xw=''):
'''
Compressibility of moist air
'''
_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
_t = T - 273.15
if xw == '':
_xw = self.xw(p, T, RH)
else:
_xw = xw
Z = 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)
return Z
def rho(self, p, T, RH, xc, xw=''):
'''
Calculates the density of air for the given parameters.
T in Kelvin
p in Pascal
RH between 0 and 1
xc in ppm
'''
_Ma = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) ) # Molar mass in kg/mol of dry air containing xc ppm of CO2
_Mw = 0.018015 # kg/mol
_R = 8.314510 # J mol^-1 K^-1, gas constant
if xw == '':
_xw = self.xw(p, T, RH)
else:
_xw = xw
_Z = self.Z(p, T, RH, xw=_xw)
return ( ( p * _Ma ) / (_Z * _R * T) ) * ( 1 - _xw * (1 - _Mw/_Ma))
def n_tph(self, l, T, p, RH, xc):
'''
Calculates the refraction of the atmosphere for any conditions.
'''
_n_a = self.n_axs(l, xc)
_n_w = self.n_ws(l)
_La = (_n_a**2 - 1) / (_n_a**2 + 2)
_Lw = (_n_w**2 - 1) / (_n_w**2 + 2)
# Calculate the density of the dry air components, using the partial pressure of the dry air, i.e. p_dry = p_atm - p_water
_rho_a = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc) # Dry air component for the actual conditions.
_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc) # Dry air component at standard conditions.
if RH == 0:
# nan is returned for _rho_w if we put in a (partial) pressure of 0
_rho_w = 0
else:
# Partial pressure of water vapour is p_water = h*svp
_rho_w = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc) # Moist air component for the actual conditions.
_rho_ws = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc) # Moist air component at standard conditions.
_Ltotal = ( _rho_a / _rho_axs ) * _La + ( _rho_w / _rho_ws ) * _Lw
# print 'Dry air component = ', (_rho_a/_rho_axs)*(_n_a-1)*1e8
# print 'Wet component = ', ( _rho_w / _rho_ws ) * (_n_w - 1)*1e8
# print 'Sum of components = ', (_rho_a/_rho_axs)*(_n_a-1)*1e8+( _rho_w / _rho_ws ) * (_n_w - 1)*1e8
return np.sqrt( (1 + 2 * _Ltotal) / (1 - _Ltotal) )
def dn_tph(self, l, T, p, RH, xc, dl=None, dT=None, dP=None, dRH=None, dCO2=None):
'''
Calculates the uncertainty in the refractive index, following the Ciddor model,
by using fully propagated error analysis.
This involves almost endless use of the chain rule and evaluation of both the moist
air and dry air components. Note: dry air does not contain any water, i.e. H=0, x_w=0, etc.
Parameters:
l - Wavelength in micron (0.4 - 2.5 micron)
T - Temperature in Kelvin
p - Pressure in Pascals
RH - Relative humidity (no units, don't use percentage but value 0<x<1)
xc - CO2 density in ppm
dl - Uncertainty in the blue wavelength
dT - Uncertainty in the temperature, default: 0.2 K
dP - Uncertainty in the pressure, default: 20 Pa
dRH - Uncertainty in the relative humidity, default: 0.02
dCO2- Uncertainty in the CO2 density, default: 20
'''
if dP == None:
dP = self.dP
if dT == None:
dT = self.dT
if dRH == None:
dRH = self.dH
if dCO2 == None:
dCO2 = self.dCO2
# constants
_cf = 1.022
_w0 = 295.235 # um^-2
_w1 = 2.6422 # um^-2
_w2 = -0.032380 # um^-4
_w3 = 0.004028 # um^-6
_k0 = 238.0185 # um^-2
_k1 = 5792105 # um^-2
_k2 = 57.362 # um^-2
_k3 = 167917 # um^-2
_A = 1.2378847e-5 # K^-2
_B = -1.9121316e-2 # K^-2
_C = 33.93711047
_D = -6.3431645e3 # K^-2
_Mw = 0.018015 # kg/mol
_R = 8.314510 # J mol^-1 K^-1, gas constant
_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
_t = T - 273.15
_aa = 1.00062
_bb = 3.14e-8 # Pa^-1
_cc = 5.6e-7 # deg C^-2
_RH_a = 0
_RH_w = RH
# Calculated values
# _rho_obs = self.rho(p, T, RH, xc)
# _ref = refraction(lat, h, p, _rho_obs)
# _alpha = _ref.rc * np.sin(np.deg2rad(z)) / (_ref.rc + _ref.H)
_p_a = p - RH*self.svp(T)
_p_w = RH*self.svp(T)
_x_w_a = self.xw(_p_a, T, _RH_a)
_x_w_w = self.f(_p_w, T)
_Z_a = self.Z(p=_p_a, T=T, RH=_RH_a)
_Z_w = self.Z(p=_p_w, T=T, RH=_RH_w, xw=_x_w_w)
_Ma = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) ) # Molar mass in kg/mol of dry air containing xc ppm of CO2
_p_sv = self.svp(T)
_rho_a = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc)
_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc)
if RH == 0.0:
# nan is returned for _rho_w if we put in a (partial) pressure of 0
_rho_w = 0
else:
_rho_w = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc)
_rho_ws = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc)
_n_a1 = self.n_axs(l, xc)
_n_w1 = self.n_ws(l)
_L1_a = (_n_a1**2 - 1) / (_n_a1**2 + 2) # Dry air component of L
_L1_w = (_n_w1**2 - 1) / (_n_w1**2 + 2)
_L1 = ( _rho_a / _rho_axs ) * _L1_a + ( _rho_w / _rho_ws ) * _L1_w # Overall L
_n1 = np.sqrt( (1 + 2 * _L1) / (1 - _L1) )
# Partial derivatives for derivation to l1
_dn_wdl = _cf * 1e-8 * (-1 * _w1 * 2/l**3 - _w2 * 4/l**5 - _w3 * 6/l**7)
_dLdn_w1 = _rho_w / _rho_ws * 6 * _n_w1 / (_n_w1**2 + 2)**2
_dn_adl = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l ) / ((_k0 * l**2 - 1)**2) - (2 * _k3 * l) / ((_k2 * l**2 - 1)**2) )
_dLdn_a1 = _rho_a / _rho_axs * 6 * _n_a1 / (_n_a1**2 + 2)**2
_dndL1 = ( ( (1 + 2*_L1) / (1 - _L1)**2 ) + (2 / (1-_L1)) ) / (2 * np.sqrt( (1 + 2*_L1) / (1 - _L1) ) )
_dndl = _dndL1 * ( _dLdn_a1 * _dn_adl + _dLdn_w1 * _dn_wdl)
# # Partial derivatives for derivation to l2
# Partial derivatives for derivation to T
# Dry air parts:
_drho_adp_a = _Ma / (_Z_a * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
_dp_adp_sv = - RH
_dp_svdT = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
_drho_adTZ = - _p_a * _Ma / (_R * ( _Z_a * T)**2) * (1 - _x_w_a * (1 - _Mw / _Ma))
_dTZdT_a = _p_a / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_a + (_c0 + _c1*_t) * _x_w_a**2 ) - _p_a * (_a1 + 2*_a2*_t + _b1*_x_w_a + _c1*_x_w_a**2) - 2*(_p_a/T)**2 * (_d + _e * _x_w_a**2) + _Z_a
_dTZdxw_a = 2*_e*_p_a**2 / T * _x_w_a - _p_a * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_a)
_dxw_adf = _RH_a * _p_sv / _p_a
_dfdT = 2*_cc*_t
_dxw_adp_sv = self.f(_p_a, T) * _RH_a / _p_a
_drho_adxw_a = 0 # as rho_a has not water content, i.e. rho = p*M/ZRT
# Moist air parts:
_drho_wdp_w = _Ma / (_Z_w * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
_dp_wdp_sv = RH
_dp_svdT = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
_drho_wdTZ = - _p_w * _Ma / (_R * (_Z_w * T)**2) * (1 - _x_w_w * (1 - _Mw / _Ma))
_dTZdT_w = _p_w / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_w + (_c0 + _c1*_t) * _x_w_w**2 ) - _p_w * (_a1 + 2*_a2*_t + _b1*_x_w_w + _c1*_x_w_w**2) - 2*(_p_w/T)**2 * (_d + _e * _x_w_w**2) + _Z_w
_dTZdxw_w = 2*_e*_p_w**2 / T * _x_w_w - _p_w * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_w)
_dxw_wdf = 1
_dfdT = 2*_cc*_t
_dxw_wdp_sv = self.f(_p_w, T) / _p_sv
_drho_wdxw_w = _p_w * _Ma / (_Z_w * _R * T) * (_Mw / _Ma - 1)
_tA = _drho_adp_a * _dp_adp_sv * _dp_svdT + _drho_adTZ * (_dTZdT_a + _dTZdxw_a * (_dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT) ) + _drho_adxw_a * ( _dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT )
_tB = _drho_wdp_w * _dp_wdp_sv * _dp_svdT + _drho_wdTZ * (_dTZdT_w + _dTZdxw_w * (_dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT) ) + _drho_wdxw_w * ( _dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT )
_dL1drho_a = _L1_a / _rho_axs
_dL1drho_w = _L1_w / _rho_ws
_dndT = _tA * (_dndL1 * _dL1drho_a) + _tB * (_dndL1 * _dL1drho_w)
# Partial derivatives for derivation to p
# Dry air parts:
_drho_adZ = - _p_a * _Ma / ( _Z_a**2 * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
_dZdp_a = - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_a + (_c0 + _c1*_t)*_x_w_a**2 ) + 2*_p_a/(T**2) * (_d + _e * _x_w_a**2)
_dfdp = _bb
_dxw_adp_a = - self.f(_p_a, T) * _RH_a * _p_sv / (_p_a)**2
# Moist air parts
_dZdp_w = - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_w + (_c0 + _c1*_t)*_x_w_w**2 ) + 2*_p_w/(T**2) * (_d + _e * _x_w_w**2)
# Not necessary, since dpwdp=0
_pA = _drho_adp_a + _drho_adZ * _dZdp_a + _drho_adxw_a * (_dxw_adf * _dfdp + _dxw_adp_a)
_dndP = _pA * (_dndL1 * _dL1drho_a)
# Partial derivatives for derivation to H
# Dry air parts:
_dp_adH = - _p_sv
_dZdxw_a = (_p_a / T)**2 * 2 * _e * _x_w_a - (_p_a / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_a )
_dxw_adH = self.f(_p_a, T) * _p_sv / _p_a
# Moist air parts:
_dp_wdH = _p_sv
_dZdxw_w = (_p_w / T)**2 * 2 * _e * _x_w_w - (_p_w / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_w )
_dxw_wdH = 0
_dxw_wdp_w = 0
_drho_wdZ = - _p_w * _Ma / ( _Z_w**2 * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
_hA = _drho_adp_a * _dp_adH + _drho_adZ * ( _dZdp_a * _dp_adH + _dZdxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH) ) + _drho_adxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH)
_hB = _drho_wdp_w * _dp_wdH + _drho_wdZ * ( _dZdp_w * _dp_wdH + _dZdxw_w * ( _bb * _p_sv ) ) + _drho_wdxw_w * ( _bb * _p_sv )
_dndRH = _hA * (_dndL1 * _dL1drho_a) + _hB * (_dndL1 * _dL1drho_w)
# Partial derivatives for derivation to xc
# Dry air parts:
_drho_adMa = _p_a * (1 - _x_w_a) / (_Z_a * _R * T)
_dL1drho_axs = - (_rho_a * _L1_a)/(_rho_axs**2)
_drho_axsdMa = _drho_adMa
_dL1dL_a = _rho_a / _rho_axs
_dL_adn_a1 = 6 * _n_a1 / (_n_a1**2 + 2)**2
# Moist air parts:
_drho_wdMa = _p_w * (1 - _x_w_w) / (_Z_w * _R * T)
_dL1drho_ws = - (_rho_w * _L1_w)/(_rho_ws**2)
_drho_wsdMa = _drho_wdMa
# Other parts:
_dMadx_c = 12.011e-9
_dn_adx_c1 = (self.n_as(l) - 1) * 0.534e-6
_dL1dx_c = _dL1drho_a * _drho_adMa * _dMadx_c + _dL1drho_w * _drho_wdMa * _dMadx_c + _dL1drho_axs * _drho_axsdMa * _dMadx_c + _dL1drho_ws * _drho_wsdMa * _dMadx_c + _dL1dL_a * _dL_adn_a1 * _dn_adx_c1
_dndCO2 = _dndL1 * _dL1dx_c
# Final error:
_nError = np.sqrt( (_dndl * dl)**2 + (_dndT * dT)**2 + (_dndP * dP)**2 + (_dndRH * dRH)**2 + (_dndCO2 * dCO2)**2 )
return _nError
atmosphericModel/dispersionModels.py
-448
View File
@@ -1,448 +0,0 @@
#!/usr/bin/env Python3
import numpy as np
from scipy import integrate
from .ciddorModel import Observatory
'''
@ Author: Joost van den Born
@ Contact: born@astron.nl
@ Description: This script offers a complete model of the atmospheric refraction
and dispersion. It uses the Ciddor1996 paper to calculate the
index of refraction for a wide range of conditions, with corresponding
error propagation if needed.
The actual refraction is calculated by any of the following models:
- Plane parallel atmosphere
- Cassini spherical homogeneous atmosphere
- Other methods described in Corbard et al., 2019.
'''
class dispersion(Observatory):
def __init__(self, lat, h):
Observatory.__init__(self)
self.h = h # height above sea level in meters
self.lat = lat # Latitude of the observer in degrees
self.rc = self.set_rc() * 1000 # Radius of curvature of the earth at the observer
def setReducedHeight(self, P, rho):
self.H = self.set_H(P, rho) # Reduced height of the atmosphere assuming ideal gas law
def planeParallel(self, n1, n2, zenith):
'''
Refraction based on a flat atmosphere, easily derived from Snell's law.
'''
_R1 = np.arcsin(n1 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
_R2 = np.arcsin(n2 * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
return np.degrees(_R1) - np.degrees(_R2)
def cassini(self, n1, n2, zenith, reduced_height=None):
'''
Refraction of spherical atmosphere, derived from geometry using Sine law
'''
if reduced_height:
H = reduced_height
else:
H = self.H
_R1 = np.arcsin(n1 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )
_R2 = np.arcsin(n2 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )
return np.degrees(_R1) - np.degrees(_R2)
def tan5(self, n1, n2, zenith, reduced_height=None):
'''
Same as oriani's theorem, but including a higher order term.
Found in Corbard et al., 2019.
'''
if reduced_height:
H = reduced_height
else:
H = self.H
_rList = []
for _n in [n1, n2]:
_a = _n - 1
_b = H / self.rc
_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) ) - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3 + 3 * _a * (_b - _a/2)**2 * np.tan( np.deg2rad(zenith) )**5
_rList.append(np.degrees(_R))
return _rList[0] - _rList[1]
def oriani(self, n1, n2, zenith, reduced_height=None):
'''
Classic refraction formula with two tan(z) terms. This does not assume any
information about the structure of the atmosphere. Found with Oriani's theorem.
'''
if reduced_height:
H = reduced_height
else:
H = self.H
_rList = []
for _n in [n1, n2]:
_a = _n - 1
_b = H / self.rc
_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) ) - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3
_rList.append(np.degrees(_R))
return _rList[0] - _rList[1]
def corbard(self, n1, n2, zenith, reduced_height=None):
'''
Corbard et al., 2019, mention an additional formula based on the error function.
Oriani's theorem can be derived from this equation by 'keeping only the three first
terms of its asymptotic expansion'.
'''
if reduced_height:
H = reduced_height
else:
H = self.H
_rList = []
for _n in [n1, n2]:
_a = _n - 1
_b = H / self.rc
_R = _a * ( (2 - _a) / (np.sqrt( 2*_b - _a )) ) * np.sin( np.deg2rad(zenith) ) * self.psi( (np.cos(np.deg2rad(zenith))) / np.sqrt(2*_b - _a) )
_rList.append(np.degrees(_R))
return _rList[0] - _rList[1]
def refractionIntegral(self, n1, n2, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
'''
Using data that gives temperature, pressure and density as a function of height,
the gladstone-dale relation is invoked to calculate the path of a light ray
along the atmosphere. This is allows us to also include atmospheric activity in
different atmosphere layers.
Recommended use with the US Standard Atmosphere from 1976.
USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html
The method described in eq. 3 in Auer & Standish (2000) is used for the calculation.
'''
if not R0:
R0 = self.rc
if useStandardAtmosphere:
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.txt', unpack=True, skip_header=6)
else:
_Dens = rhoData
_h = heightData
_rList = []
z0 = np.deg2rad(z0)
# Integration of the refraction integral over the atmosphere layers.
for _n in [n1, n2]:
_Dens = np.asarray(_Dens)
_h = np.asarray(_h)
_R = R0 + _h * 1000
# Gladstone-Dale relation between refractive index and density
_n_h = 1 + (_n - 1) / _Dens[0] * _Dens
# Calculation of the zenith angle for given height and refractive index, from the refractive invariant
_z = np.arcsin( (_n * R0) / (_n_h * _R) * np.sin(z0) )
# Calculation of log space derivative of n and r
_lnn = np.log(_n_h)
_lnr = np.log(_R)
_dlnn = np.asarray([_lnn[i+1] - _lnn[i] for i in range(len(_lnn)-1)])
_dlnr = np.asarray([_lnr[i+1] - _lnr[i] for i in range(len(_lnr)-1)])
_dz = np.asarray([(_z[i] - _z[i+1]) for i in range(len(_z)-1)])
# Calculation of the integrand in eq. 3 in Auer & Standish
_Integrand = -1*(_dlnn/_dlnr) / (1 + _dlnn/_dlnr) * _dz
_rList.append(np.degrees(np.sum(_Integrand)))
return _rList[0] - _rList[1]
def psi(self, x):
_f = lambda a: np.exp(-a**2)
_int = integrate.quad(_f, x, np.inf)
return np.exp(x**2) * _int[0]
def set_rc(self):
'''
Returns the radius of the curvature at the observer in km,
taking into account the latitude of the observer and the corresponding
ellipsoidality of the earth.
'''
_A = 6378.137 # km
_B = 6356.752 # km
_phi = np.deg2rad(self.lat)
_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
return _rc0 + self.h/1000
def set_H(self, P, rho):
'''
Calculates the reduced height of the atmosphere, assuming ideal gas law.
'''
_phi = np.deg2rad(self.lat)
_c1 = 5.2790414e-3
_c2 = 2.32718e-5
_c3 = 1.262e-7
_c4 = 7e-10
_g0 = 9.780327 # m/s^2
_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
_g = _g0_local - (3.0877e-6 - 4.3e-9 * np.sin(_phi)**2) * self.h + 7.2e-13 * self.h**2
return P / (rho * _g)
def H_isotherm(self, T):
'''
Calculates the reduced height of the atmosphere, assuming an isotherm distribution.
'''
_kb = 1.38064852e-23 # kg m2 s-2 K-1
_m = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
_phi = np.deg2rad(self.lat)
_c1 = 5.2790414e-3
_c2 = 2.32718e-5
_c3 = 1.262e-7
_c4 = 7e-10
_g0 = 9.780327 # m/s^2
_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
return (_kb*T) / (_m * _g0_local)
def cassiniError(self, z, l1, l2, T, p, RH, xc, dl1=0, dl2=0, dT=0.2, dP=20, dRH=0.02, dCO2=20, dz=0, lat=None, h=None):
'''
Calculates the uncertainty in the atmospheric dispersion, following the Cassini model,
by using fully propagated error analysis.
This involves almost endless use of the chain rule and evaluation of both the moist
air and dry air components. Note: dry air does not contain any water, i.e. H=0, x_w=0, etc.
Not taken into account is the errror propagation in the scale height of the atmosphere or
the zenith error, from my results this should give a small or even negligible difference.
Parameters:
z - Zenith angle in degrees
l1 - Bluest wavelength in micron (0.4 - 2.5 micron)
l2 - Reddest wavelength in micron (0.4 - 2.5 micron)
T - Temperature in Kelvin
p - Pressure in Pascals
RH - Relative humidity (no units, don't use percentage but value 0<x<1)
xc - CO2 density in ppm
dl1 - Uncertainty in the blue wavelength
dl2 - Uncertainty in the red wavelength
dT - Uncertainty in the temperature, default: 0.2 K
dP - Uncertainty in the pressure, default: 20 Pa
dRH - Uncertainty in the relative humidity, default: 0.02
dCO2- Uncertainty in the CO2 density, default: 20
lat - [Optional] Latitude of the observatory in degrees
h - [Optional] Height in m of the observatory
'''
if lat == None:
lat = self.lat
else:
self.lat = lat
if h == None:
h = self.h
else:
self.h = h
# constants
_cf = 1.022
_w0 = 295.235 # um^-2
_w1 = 2.6422 # um^-2
_w2 = -0.032380 # um^-4
_w3 = 0.004028 # um^-6
_k0 = 238.0185 # um^-2
_k1 = 5792105 # um^-2
_k2 = 57.362 # um^-2
_k3 = 167917 # um^-2
_A = 1.2378847e-5 # K^-2
_B = -1.9121316e-2 # K^-2
_C = 33.93711047
_D = -6.3431645e3 # K^-2
_Mw = 0.018015 # kg/mol
_R = 8.314510 # J mol^-1 K^-1, gas constant
_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
_t = T - 273.15
_aa = 1.00062
_bb = 3.14e-8 # Pa^-1
_cc = 5.6e-7 # deg C^-2
_RH_a = 0
_RH_w = RH
# Calculated values
_rho_obs = self.rho(p, T, RH, xc)
self.rc = self.set_rc() * 1000
self.H = self.set_H(p, _rho_obs)
_alpha = self.rc * np.sin(np.deg2rad(z)) / (self.rc + self.H)
_p_a = p - RH*self.svp(T)
_p_w = RH*self.svp(T)
_x_w_a = self.xw(_p_a, T, _RH_a)
_x_w_w = self.f(_p_w, T)
_Z_a = self.Z(p=_p_a, T=T, RH=_RH_a)
_Z_w = self.Z(p=_p_w, T=T, RH=_RH_w, xw=_x_w_w)
_Ma = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) ) # Molar mass in kg/mol of dry air containing xc ppm of CO2
_p_sv = self.svp(T)
_rho_a = self.rho(p=p-RH*self.svp(T), T=T, RH=0, xc=xc)
_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc)
if RH == 0.0:
# nan is returned for _rho_w if we put in a (partial) pressure of 0
_rho_w = 0
else:
_rho_w = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc)
_rho_ws = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc)
_n_a1 = self.n_axs(l1, xc)
_n_w1 = self.n_ws(l1)
_n_a2 = self.n_axs(l2, xc)
_n_w2 = self.n_ws(l2)
_L1_a = (_n_a1**2 - 1) / (_n_a1**2 + 2) # Dry air component of L
_L1_w = (_n_w1**2 - 1) / (_n_w1**2 + 2)
_L1 = ( _rho_a / _rho_axs ) * _L1_a + ( _rho_w / _rho_ws ) * _L1_w # Overall L
_L2_a = (_n_a2**2 - 1) / (_n_a2**2 + 2) # Dry air component of L
_L2_w = (_n_w2**2 - 1) / (_n_w2**2 + 2)
_L2 = ( _rho_a / _rho_axs ) * _L2_a + ( _rho_w / _rho_ws ) * _L2_w # Overall L
_n1 = np.sqrt( (1 + 2 * _L1) / (1 - _L1) )
_n2 = np.sqrt( (1 + 2 * _L2) / (1 - _L2) )
# Partial derivatives for derivation to l1
_dn_wdl1 = _cf * 1e-8 * (-1 * _w1 * 2/l1**3 - _w2 * 4/l1**5 - _w3 * 6/l1**7)
_dLdn_w1 = _rho_w / _rho_ws * 6 * _n_w1 / (_n_w1**2 + 2)**2
_dn_adl1 = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l1 ) / ((_k0 * l1**2 - 1)**2) - (2 * _k3 * l1) / ((_k2 * l1**2 - 1)**2) )
_dLdn_a1 = _rho_a / _rho_axs * 6 * _n_a1 / (_n_a1**2 + 2)**2
_dndL1 = ( ( (1 + 2*_L1) / (1 - _L1)**2 ) + (2 / (1-_L1)) ) / (2 * np.sqrt( (1 + 2*_L1) / (1 - _L1) ) )
_dddn1 = _alpha / np.sqrt( 1 - _alpha**2 * _n1**2 )
_dddl1 = _dddn1 * _dndL1 * ( _dLdn_a1 * _dn_adl1 + _dLdn_w1 * _dn_wdl1)
# Partial derivatives for derivation to l2
_dn_wdl2 = _cf * 1e-8 * (-1 * _w1 * 2/l2**3 - _w2 * 4/l2**5 - _w3 * 6/l2**7)
_dLdn_w2 = _rho_w / _rho_ws * 6 * _n_w2 / (_n_w2**2 + 2)**2
_dn_adl2 = 1e-8 * (1 + 0.534e-6 * (xc - 450)) * ( - (2 * _k1 * l2 ) / ((_k0 * l2**2 - 1)**2) - (2 * _k3 * l2) / ((_k2 * l2**2 - 1)**2) )
_dLdn_a2 = _rho_a / _rho_axs * 6 * _n_a2 / (_n_a2**2 + 2)**2
_dndL2 = ( ( (1 + 2*_L2) / (1 - _L2)**2 ) + (2 / (1-_L2)) ) / (2 * np.sqrt( (1 + 2*_L2) / (1 - _L2) ) )
_dddn2 = _alpha / np.sqrt( 1 - _alpha**2 * _n2**2 )
_dddl2 = _dddn2 * _dndL2 * ( _dLdn_a2 * _dn_adl2 + _dLdn_w2 * _dn_wdl2)
# Partial derivatives for derivation to T
# Dry air parts:
_drho_adp_a = _Ma / (_Z_a * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
_dp_adp_sv = - RH
_dp_svdT = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
_drho_adTZ = - _p_a * _Ma / (_R * ( _Z_a * T)**2) * (1 - _x_w_a * (1 - _Mw / _Ma))
_dTZdT_a = _p_a / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_a + (_c0 + _c1*_t) * _x_w_a**2 ) - _p_a * (_a1 + 2*_a2*_t + _b1*_x_w_a + _c1*_x_w_a**2) - 2*(_p_a/T)**2 * (_d + _e * _x_w_a**2) + _Z_a
_dTZdxw_a = 2*_e*_p_a**2 / T * _x_w_a - _p_a * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_a)
_dxw_adf = _RH_a * _p_sv / _p_a
_dfdT = 2*_cc*_t
_dxw_adp_sv = self.f(_p_a, T) * _RH_a / _p_a
_drho_adxw_a = 0 # as rho_a has not water content, i.e. rho = p*M/ZRT
# Moist air parts:
_drho_wdp_w = _Ma / (_Z_w * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
_dp_wdp_sv = RH
_dp_svdT = (2*_A*T + _B - _D/(T**2)) * np.exp(_A * T**2 + _B * T + _C + _D/T )
_drho_wdTZ = - _p_w * _Ma / (_R * (_Z_w * T)**2) * (1 - _x_w_w * (1 - _Mw / _Ma))
_dTZdT_w = _p_w / T * (_a2 * _t**2 + _a1 * _t + _a0 + (_b1*_t + _b0) * _x_w_w + (_c0 + _c1*_t) * _x_w_w**2 ) - _p_w * (_a1 + 2*_a2*_t + _b1*_x_w_w + _c1*_x_w_w**2) - 2*(_p_w/T)**2 * (_d + _e * _x_w_w**2) + _Z_w
_dTZdxw_w = 2*_e*_p_w**2 / T * _x_w_w - _p_w * (_b1*_t + _b0 + 2 * (_c0 + _c1*_t) * _x_w_w)
_dxw_wdf = 1 # _RH_w * _p_sv / _p_w
_dfdT = 2*_cc*_t
_dxw_wdp_sv = self.f(_p_w, T) / _p_sv # * _RH_w / _p_w
_drho_wdxw_w = _p_w * _Ma / (_Z_w * _R * T) * (_Mw / _Ma - 1)
_tA = _drho_adp_a * _dp_adp_sv * _dp_svdT + _drho_adTZ * (_dTZdT_a + _dTZdxw_a * (_dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT) ) + _drho_adxw_a * ( _dxw_adf * _dfdT + _dxw_adp_sv * _dp_svdT )
_tB = _drho_wdp_w * _dp_wdp_sv * _dp_svdT + _drho_wdTZ * (_dTZdT_w + _dTZdxw_w * (_dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT) ) + _drho_wdxw_w * ( _dxw_wdf * _dfdT + _dxw_wdp_sv * _dp_svdT )
_dL1drho_a = _L1_a / _rho_axs
_dL2drho_a = _L2_a / _rho_axs
_dL1drho_w = _L1_w / _rho_ws
_dL2drho_w = _L2_w / _rho_ws
_dddT = _tA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a) + _tB * (_dddn1 * _dndL1 * _dL1drho_w - _dddn2 * _dndL2 * _dL2drho_w)
# Partial derivatives for derivation to p
# Dry air parts:
_drho_adZ = - _p_a * _Ma / ( _Z_a**2 * _R * T) * (1 - _x_w_a * (1 - _Mw / _Ma))
_dZdp_a = - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_a + (_c0 + _c1*_t)*_x_w_a**2 ) + 2*_p_a/(T**2) * (_d + _e * _x_w_a**2)
_dfdp = _bb
_dxw_adp_a = - self.f(_p_a, T) * _RH_a * _p_sv / (_p_a)**2
# Moist air parts
_dZdp_w = - 1 / T * (_a0 + _a1 * _t + _a2 * _t**2 + (_b0 + _b1*_t)*_x_w_w + (_c0 + _c1*_t)*_x_w_w**2 ) + 2*_p_w/(T**2) * (_d + _e * _x_w_w**2)
# Not necessary, since dpwdp=0
_pA = _drho_adp_a + _drho_adZ * _dZdp_a + _drho_adxw_a * (_dxw_adf * _dfdp + _dxw_adp_a)
_dddP = _pA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a)
# Partial derivatives for derivation to H
# Dry air parts:
_dp_adH = - _p_sv
_dZdxw_a = (_p_a / T)**2 * 2 * _e * _x_w_a - (_p_a / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_a )
_dxw_adH = self.f(_p_a, T) * _p_sv / _p_a
# Moist air parts:
_dp_wdH = _p_sv
_dZdxw_w = (_p_w / T)**2 * 2 * _e * _x_w_w - (_p_w / T) * ( (_b0 + _b1 * _t) + 2 * (_c0 + _c1 * _t) * _x_w_w )
_dxw_wdH = 0
_dxw_wdp_w = 0
_drho_wdZ = - _p_w * _Ma / ( _Z_w**2 * _R * T) * (1 - _x_w_w * (1 - _Mw / _Ma))
_hA = _drho_adp_a * _dp_adH + _drho_adZ * ( _dZdp_a * _dp_adH + _dZdxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH) ) + _drho_adxw_a * (_dxw_adH + _dxw_adp_a * _dp_adH + _dxw_adf * _dfdp * _dp_adH)
_hB = _drho_wdp_w * _dp_wdH + _drho_wdZ * ( _dZdp_w * _dp_wdH + _dZdxw_w * ( _bb * _p_sv ) ) + _drho_wdxw_w * ( _bb * _p_sv )
_dddRH = _hA * (_dddn1 * _dndL1 * _dL1drho_a - _dddn2 * _dndL2 * _dL2drho_a) + _hB * (_dddn1 * _dndL1 * _dL1drho_w - _dddn2 * _dndL2 * _dL2drho_w)
# Partial derivatives for derivation to xc
# Dry air parts:
_drho_adMa = _p_a * (1 - _x_w_a) / (_Z_a * _R * T)
_dL1drho_axs = - (_rho_a * _L1_a)/(_rho_axs**2)
_dL2drho_axs = - (_rho_a * _L2_a)/(_rho_axs**2)
_drho_axsdMa = _drho_adMa
_dL1dL_a = _rho_a / _rho_axs
_dL_adn_a1 = 6 * _n_a1 / (_n_a1**2 + 2)**2
_dL2dL_a = _dL1dL_a
_dL_adn_a2 = 6 * _n_a2 / (_n_a2**2 + 2)**2
# Moist air parts:
_drho_wdMa = _p_w * (1 - _x_w_w) / (_Z_w * _R * T)
_dL1drho_ws = - (_rho_w * _L1_w)/(_rho_ws**2)
_dL2drho_ws = - (_rho_w * _L2_w)/(_rho_ws**2)
_drho_wsdMa = _drho_wdMa
# Other parts:
_dMadx_c = 12.011e-9
_dn_adx_c1 = (self.n_as(l1) - 1) * 0.534e-6
_dn_adx_c2 = (self.n_as(l2) - 1) * 0.534e-6
_dL1dx_c = _dL1drho_a * _drho_adMa * _dMadx_c + _dL1drho_w * _drho_wdMa * _dMadx_c + _dL1drho_axs * _drho_axsdMa * _dMadx_c + _dL1drho_ws * _drho_wsdMa * _dMadx_c + _dL1dL_a * _dL_adn_a1 * _dn_adx_c1
_dL2dx_c = _dL2drho_a * _drho_adMa * _dMadx_c + _dL2drho_w * _drho_wdMa * _dMadx_c + _dL2drho_axs * _drho_axsdMa * _dMadx_c + _dL2drho_ws * _drho_wsdMa * _dMadx_c + _dL2dL_a * _dL_adn_a2 * _dn_adx_c2
_dddCO2 = _dddn1 * _dndL1 * _dL1dx_c - _dddn2 * _dndL2 * _dL2dx_c
# Partial derivatives for derivation to z
_zenith = np.deg2rad(z)
_dz = np.deg2rad(dz)
_R = self.rc
_H = self.H
_dddz = (_R * np.cos(_zenith)) / (_R + _H) * ( _n1 / np.sqrt(1-(_n1**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - _n2 / np.sqrt(1-(_n2**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )
# Final error:
_dispersionError = np.sqrt( (_dddl1 * dl1)**2 + (_dddl2 * dl2)**2 + (_dddT * dT)**2 + (_dddP * dP)**2 + (_dddRH * dRH)**2 + (_dddCO2 * dCO2)**2 + (_dddz * _dz)**2)
return np.degrees(_dispersionError)
atmosphericModel/examples.py
-106
View File
@@ -1,106 +0,0 @@
#!/usr/bin/env Python3
from .ciddorModel import Observatory
from .refractionModels import refraction
from .dispersionModels import dispersion
'''
Author: Joost van den Born
Some functions for quick testing.
'''
def quick_refractive_index(l, conditions='STANDARD'):
if conditions == 'STANDARD':
# Parameters at Cerro Armazones
T = 288.15 # k
P = 101325 # Pa
H = 0.0
xc = 450 # ppm
lat = 0 # degrees
h = 0 # m
elif conditions == 'CERRO_ARMAZONES':
# Parameters at Cerro Armazones
T = 279.65 # k
P = 71200 # Pa
H = 0.22
xc = 450 # ppm
lat = -24.5983 # degrees
h = 3064 # m
# Initializing dispersion model
at = Observatory()
# Calculating indices of refraction for l
n = at.n_tph(l=l, T=T, p=P, RH=H, xc=xc)
return n
def quick_refraction(l, z, conditions='STANDARD'):
if conditions == 'STANDARD':
# Parameters at Cerro Armazones
T = 288.15 # k
P = 101325 # Pa
H = 0.0
xc = 450 # ppm
lat = 0 # degrees
h = 0 # m
elif conditions == 'CERRO_ARMAZONES':
# Parameters at Cerro Armazones
T = 279.65 # k
P = 71200 # Pa
H = 0.22
xc = 450 # ppm
lat = -24.5983 # degrees
h = 3064 # m
# Initializing dispersion model
at = Observatory()
# Calculating indices of refraction for l
n = at.n_tph(l=l, T=T, p=P, RH=H, xc=xc)
# Density of the atmosphere (following CIPM-81/91 equations)
rho = at.rho(p=P, T=T, RH=H, xc=xc)
# Initializing refraction model and setting the reduced height
ref = refraction(lat, h)
ref.setReducedHeight(P, rho)
return ref.cassini(n, z)
def quick_dispersion(l1, l2, z, conditions='STANDARD'):
if conditions == 'STANDARD':
# Parameters at Cerro Armazones
T = 288.15 # k
P = 101325 # Pa
H = 0.0
xc = 450 # ppm
lat = 0 # degrees
h = 0 # m
elif conditions == 'CERRO_ARMAZONES':
# Parameters at Cerro Armazones
T = 279.65 # k
P = 71200 # Pa
H = 0.22
xc = 450 # ppm
lat = -24.5983 # degrees
h = 3064 # m
# Initializing dispersion model
at = Observatory()
# Calculating indices of refraction for l1 and l2
n1 = at.n_tph(l=l1, T=T, p=P, RH=H, xc=xc)
n2 = at.n_tph(l=l2, T=T, p=P, RH=H, xc=xc)
# Density of the atmosphere (following CIPM-81/91 equations)
rho = at.rho(p=P, T=T, RH=H, xc=xc)
# Initializing refraction model and setting the reduced height
disp = dispersion(lat, h)
disp.setReducedHeight(P, rho)
return disp.cassini(n1, n2, z)
atmosphericModel/refractionModels.py
-218
View File
@@ -1,218 +0,0 @@
#!/usr/bin/env Python3
import numpy as np
from scipy import integrate
from .ciddorModel import Observatory
'''
@ Author: Joost van den Born
@ Contact: born@astron.nl
@ Description: This script offers a complete model of the atmospheric refraction
and dispersion. It uses the Ciddor1996 paper to calculate the
index of refraction for a wide range of conditions, with corresponding
error propagation if needed.
The actual refraction is calculated by any of the following models:
- Plane parallel atmosphere
- Cassini spherical homogeneous atmosphere
- Other methods described in Corbard et al., 2019.
'''
class refraction(Observatory):
def __init__(self, lat, h):
Observatory.__init__(self)
self.h = h # height above sea level in meters
self.lat = lat # Latitude of the observer in degrees
self.rc = self.set_rc() * 1000 # Radius of curvature of the earth at the observer
# self.H2 = self.H_isotherm(T) # Reduced height of isotherm atmosphere
def setReducedHeight(self, P, rho):
self.H = self.set_H(P, rho) # Reduced height of the atmosphere assuming ideal gas law
def planeParallel(self, n, zenith):
'''
Refraction based on a flat atmosphere, easily derived from Snell's law.
'''
_R = np.arcsin(n * np.sin( np.deg2rad(zenith) )) - np.deg2rad(zenith)
return np.degrees(_R)
def cassini(self, n1, zenith, reduced_height=None):
'''
Refraction of spherical atmosphere, derived from geometry using Sine law
'''
if reduced_height:
H = reduced_height
else:
H = self.H
_R1 = np.arcsin(n1 * self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H)) - np.arcsin(self.rc * np.sin(np.deg2rad(zenith)) / (self.rc + H) )
return np.degrees(_R1)
def tan5(self, n1, zenith, reduced_height=None):
'''
Same as oriani's theorem, but including a higher order term.
Found in Corbard et al., 2019.
'''
if reduced_height:
H = reduced_height
else:
H = self.H
_a = n1 - 1
_b = H / self.rc
_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) ) - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3 + 3 * _a * (_b - _a/2)**2 * np.tan( np.deg2rad(zenith) )**5
return np.degrees(_R)
def oriani(self, n1, zenith, reduced_height=None):
'''
Classic refraction formula with two tan(z) terms. This does not assume any
information about the structure of the atmosphere. Found with Oriani's theorem.
'''
if reduced_height:
H = reduced_height
else:
H = self.H
_a = n1 - 1
_b = H / self.rc
_R = _a * (1 - _b) * np.tan( np.deg2rad(zenith) ) - _a * (_b - _a/2) * np.tan( np.deg2rad(zenith) )**3
return np.degrees(_R)
def corbard(self, n1, zenith, reduced_height=None):
'''
Corbard et al., 2019, mention an additional formula based on the error function.
Oriani's theorem can be derived from this equation by 'keeping only the three first
terms of its asymptotic expansion'.
'''
if reduced_height:
H = reduced_height
else:
H = self.H
_a = n1 - 1
_b = H / self.rc
_R = _a * ( (2 - _a) / (np.sqrt( 2*_b - _a )) ) * np.sin( np.deg2rad(zenith) ) * self.psi( (np.cos(np.deg2rad(zenith))) / np.sqrt(2*_b - _a) )
return np.degrees(_R)
def refractionIntegral(self, n0, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
'''
Using data that gives temperature, pressure and density as a function of height,
the gladstone-dale relation is invoked to calculate the path of a light ray
along the atmosphere. This is allows us to also include atmospheric activity in
different atmosphere layers.
Recommended use with the US Standard Atmosphere from 1976.
USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html
'''
if not R0:
R0 = self.rc
if useStandardAtmosphere:
_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.txt', unpack=True, skip_header=6)
else:
_Dens = rhoData
_h = heightData
z0 = np.deg2rad(z0)
_Dens = np.asarray(_Dens)
_h = np.asarray(_h)
_R = R0 + _h * 1000
# Gladstone-Dale relation between refractive index and density
_n_h = 1 + (n0 - 1) / _Dens[0] * _Dens
# Calculation of the zenith angle for given height and refractive index, from the refractive invariant
_z = np.arcsin( (n0 * R0) / (_n_h * _R) * np.sin(z0) )
# Calculation of log space derivative of n and r
_lnn = np.log(_n_h)
_lnr = np.log(_R)
_dlnn = np.asarray([_lnn[i+1] - _lnn[i] for i in range(len(_lnn)-1)])
_dlnr = np.asarray([_lnr[i+1] - _lnr[i] for i in range(len(_lnr)-1)])
_dz = np.asarray([(_z[i] - _z[i+1]) for i in range(len(_z)-1)])
# Calculation of the integrand in eq. 3 in Auer & Standish
_Integrand = -1*(_dlnn/_dlnr) / (1 + _dlnn/_dlnr) * _dz
return np.degrees(np.sum(_Integrand))
def psi(self, x):
_f = lambda a: np.exp(-a**2)
_int = integrate.quad(_f, x, np.inf)
return np.exp(x**2) * _int[0]
def set_rc(self):
'''
Returns the radius of the curvature at the observer in km,
taking into account the latitude of the observer and the corresponding
ellipsoidality of the earth.
'''
_A = 6378.137 # km
_B = 6356.752 # km
_phi = np.deg2rad(self.lat)
_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
return _rc0 + self.h/1000
def set_H(self, P, rho):
'''
Calculates the reduced height of the atmosphere, assuming ideal gas law.
'''
_phi = np.deg2rad(self.lat)
_c1 = 5.2790414e-3
_c2 = 2.32718e-5
_c3 = 1.262e-7
_c4 = 7e-10
_g0 = 9.780327 # m/s^2
_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
_g = _g0_local - (3.0877e-6 - 4.3e-9 * np.sin(_phi)**2) * self.h + 7.2e-13 * self.h**2
return P / (rho * _g)
def H_isotherm(self, T):
'''
Calculates the reduced height of the atmosphere, assuming an isotherm distribution.
'''
_kb = 1.38064852e-23 # kg m2 s-2 K-1
_m = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
_phi = np.deg2rad(self.lat)
_c1 = 5.2790414e-3
_c2 = 2.32718e-5
_c3 = 1.262e-7
_c4 = 7e-10
_g0 = 9.780327 # m/s^2
_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
return (_kb*T) / (_m * _g0_local)
def cassiniError(self, n, dn, zenith, dz=0, reduced_height=None):
'''
Calculates the uncertainty in the refraction angle in degrees.
'''
if reduced_height:
_H = reduced_height
else:
_H = self.H
_zenith = np.deg2rad(zenith)
_dz = np.deg2rad(dz)
_R = self.rc
_dRdn = (_R * np.sin(_zenith)) / ((_R + _H) * np.sqrt(1 - ((n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2)))
_dRdz = (_R * np.cos(_zenith)) / (_R + _H) * ( n / np.sqrt(1-(n**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - 1 / np.sqrt(1-(_R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )
_dR = np.sqrt( (_dRdn * dn)**2 + (_dRdz * _dz)**2)
return np.degrees(_dR)
atmosphericModel/refractivityModels.py
-349
View File
@@ -1,349 +0,0 @@
#!/usr/bin/env Python3
import numpy as np
'''
@ Author: Joost van den Born
@ Contact: born@astron.nl
@ Description: A comparison of various atmospheric refraction models.
Also these are used to calculate the difference in
atmospheric dispersion, given a set of conditions.
'''
def BarrellAndSears(l, T=288.15, p=101325, H=0):
'''
From Barrell & Sears, 1939
Equation 7.7, page 52
For atmospheric air (i.e. CO2 = 300 ppm)
This equation is applicable for the ranges T=10-30C and p=720-800 mmHg.
Input:
l - Wavelength in microns
T - Temperature in Kelvin
p - Atmospheric pressure in Pa
H - fractional humidity (0=<H=<1)
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
T = K2C(T)
P = hPA2mmHg(p/100)
R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * H / (1 + 0.003661*T)
n = R/1e6 + 1
return n
def Ciddor(l):
'''
From Ciddor, 1996
Equations 1
At default conditions only. For more flexibility use
Observatory.n_tph()
Input:
l - Wavelength in microns
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
s = 1 / l
R = 5792105 / (238.0185 - s**2) + 167917 / (57.362 - s**2)
n = R/1e8 + 1
return n
def Edlen1953(l):
'''
From Edlen, 1953
Equation 2
At standard conditions only
Input:
l - Wavelength in microns
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
s = 1 / l
R = 6432.8 + 2949810 / (146 - s**2) + 25540/(41-s**2)
n = R/1e8 + 1
return n
def Edlen1966(l):
'''
From Edlen, 1966
Equation 1
At standard conditions only
Input:
l - Wavelength in microns
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
s = 1 / l
R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2)
n = R/1e8 + 1
return n
def Owens(l, T=288.15, p=101325, H=0):
'''
From Owens, 1967
Equation 32
Note that the pressure P and H are in torricelli (Torr),
however, this is almost precisely 1 Torr = 1 mmHg
Input:
l - wavelength in micron
T - Temperature in Kelvin
p - Atmospheric pressure in Pa
H - Humdity (0<=H<=1)
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
P = hPA2mmHg(p/100)
T = K2C(T)
s = 1 / l
R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2) * (P/720.775) * ((1 + P * (0.817 - 0.0133) * 1e-6)/(1 + 0.0036610*T)) - H * (5.722 - 0.0457*s**2)
n = R/1e8 + 1
return n
def BonschPotulski(l, T=288.15, p=101325, CO2=450):
'''
From Bonsch & Potulski, 1998
Equations 6, 7, 8
Input:
l - Wavelength in microns
T - Temperature in Kelvin
p - Atmospheric pressure in Pa
CO2 - Particles per million of CO2
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
x = CO2 / 1e6
T = K2C(T)
s = 1 / l
# Rs = 8092.33 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
Rs = 8091.37 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
Rs = Rs / 1e8
Rx = Rs * (1 + 0.5327 * (x - 0.0004))
Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * T)*p) / (1 + 0.0036610 * T)
n = Rtp + 1
return n
def BirchDowns(l):
'''
From Birch & Downs, 1993 and 1994 (correction)
Equations 1 and 2
Only at standard conditions
Input:
l - wavelength in micron
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
s = 1 / l
R = 8342.54 + 2406147 / (130 - s**2) + 15998 / (38.9-s**2)
n = R/1e8 + 1
return n
def PeckReeder(l):
'''
From Peck & Reeder, 1972
Equation 3
Only at standard conditions
Input:
l - wavelength in micron
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
s = 1 / l
R = 8060.51 + 2480990 / (132.274 - s**2) + 17455.7 / (39.32957 - s**2)
n = R/1e8 + 1
return n
def slalib(l, T=288.15, p=101325, H=0, h=0, lat=-24, tlr=0.0065, eps=1e-10):
'''
SLALIB uses the same atmospheric model as ZEMAX
based on the works of Hohenkerk & Sinclair and Seidelmann.
They take the B&S refraction formula as a basis.
Requires the pySLALIB package to work.
Input:
l - Wavelength in microns
T - Temperature in Kelvin
p - Atmospheric pressure in Pa
H - fractional humidity (0=<H=<1)
h - Altitude of the observer (default=0 m)
lat - Latitude of the observer (default=-24)
tlr - temperature lapse rate (default=0.0065)
eps - precision at which the calculation is aborted (default=1e-10)
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
from pyslalib import slalib
zObs = np.deg2rad(45.06563)
p = p/100
R = []
if type(l) == float:
l = [l]
for wl in l:
R.append(slalib.sla_refro(zObs, h, T, p, H, wl, lat, tlr, eps))
if len(l) == 1:
return R[0]+1
else:
return np.asarray(R)+1
def Mathar(l, T=288.15, p=101325, H=0):
'''
Refractivity model by Mathar (2007), valid above 1.4 microns
Code taken and modified from:
https://github.com/polyanskiy/refractiveindex.info-scripts/blob/master/scripts/Mathar%202007%20-%20Air%201.3-2.5um.py
Input:
l - Wavelength in microns
T - Temperature in Kelvin
p - Atmospheric pressure in Pa
H - fractional humidity (0=<H=<1)
Returns:
n - Refractive index of atmospheric air at the given conditions
'''
if l < 1.3 or l > 24:
print('Using the Mathar model at this wavelength (l={:.3f}) is not recommended!'.format(l))
if l <= 2.65:
########################################################################
# model parameters for wavelengths between 1.3 and 2.5 microns
cref = [ 0.200192e-3, 0.113474e-9, -0.424595e-14, 0.100957e-16, -0.293315e-20, 0.307228e-24] # cm^j
cT = [ 0.588625e-1, -0.385766e-7, 0.888019e-10, -0.567650e-13, 0.166615e-16, -0.174845e-20] # cm^j · K
cTT = [-3.01513, 0.406167e-3, -0.514544e-6, 0.343161e-9, -0.101189e-12, 0.106749e-16] # cm^j · K^2
cH = [-0.103945e-7, 0.136858e-11, -0.171039e-14, 0.112908e-17, -0.329925e-21, 0.344747e-25] # cm^j · %^-1
cHH = [ 0.573256e-12, 0.186367e-16, -0.228150e-19, 0.150947e-22, -0.441214e-26, 0.461209e-30] # cm^j · %^-2
cp = [ 0.267085e-8, 0.135941e-14, 0.135295e-18, 0.818218e-23, -0.222957e-26, 0.249964e-30] # cm^j · Pa^-1
cpp = [ 0.609186e-17, 0.519024e-23, -0.419477e-27, 0.434120e-30, -0.122445e-33, 0.134816e-37] # cm^j · Pa^-2
cTH = [ 0.497859e-4, -0.661752e-8, 0.832034e-11, -0.551793e-14, 0.161899e-17, -0.169901e-21] # cm^j · K · %^-1
cTp = [ 0.779176e-6, 0.396499e-12, 0.395114e-16, 0.233587e-20, -0.636441e-24, 0.716868e-28] # cm^j · K · Pa^-1
cHp = [-0.206567e-15, 0.106141e-20, -0.149982e-23, 0.984046e-27, -0.288266e-30, 0.299105e-34] # cm^j · %^-1 · Pa^-1
sref = 1e4/2.25 # cm^−1, reference wavenumber
Tref = 273.15+17.5 # K, reference temperature
pref = 75000 # Pa, reference pressure
Href = 10 #%, reference humdity in percent
########################################################################
if l > 2.65 and l <= 4.275:
########################################################################
# model parameters for wavelengths between 2.8 and 4.2 microns
cref = [ 0.200049e-3, 0.145221e-9, 0.250951e-12, -0.745834e-15, -0.161432e-17, 0.352780e-20] # cm^j
cT = [ 0.588432e-1, -0.825182e-7, 0.137982e-9, 0.352420e-13, -0.730651e-15, -0.167911e-18] # cm^j · K
cTT = [-3.13579, 0.694124e-3, -0.500604e-6, -0.116668e-8, 0.209644e-11, 0.591037e-14] # cm^j · K^2
cH = [-0.108142e-7, 0.230102e-11, -0.154652e-14, -0.323014e-17, 0.630616e-20, 0.173880e-22] # cm^j · %^-1
cHH = [ 0.586812e-12, 0.312198e-16, -0.197792e-19, -0.461945e-22, 0.788398e-25, 0.245580e-27] # cm^j · %^-2
cp = [ 0.266900e-8, 0.168162e-14, 0.353075e-17, -0.963455e-20, -0.223079e-22, 0.453166e-25] # cm^j · Pa^-1
cpp = [ 0.608860e-17, 0.461560e-22, 0.184282e-24, -0.524471e-27, -0.121299e-29, 0.246512e-32] # cm^j · Pa^-2
cTH = [ 0.517962e-4, -0.112149e-7, 0.776507e-11, 0.172569e-13, -0.320582e-16, -0.899435e-19] # cm^j · K · %^-1
cTp = [ 0.778638e-6, 0.446396e-12, 0.784600e-15, -0.195151e-17, -0.542083e-20, 0.103530e-22] # cm^j · K · Pa^-1
cHp = [-0.217243e-15, 0.104747e-20, -0.523689e-23, 0.817386e-26, 0.309913e-28, -0.363491e-31] # cm^j · %^-1 · Pa^-1
sref = 1e4/3.4 # cm^−1, reference wavenumber
Tref = 273.15+17.5 # K, reference temperature
pref = 75000 # Pa, reference pressure
Href = 10 #%, reference humdity in percent
########################################################################
if l > 4.275 and l <= 6.35:
########################################################################
# model parameters for wavelengths between 4.35 and 5.2 microns
cref = [ 0.200020e-3, 0.275346e-9, 0.325702e-12, -0.693603e-14, 0.285610e-17, 0.338758e-18] # cm^j
cT = [ 0.590035e-1, -0.375764e-6, 0.134585e-9, 0.124316e-11, 0.508510e-13, -0.189245e-15] # cm^j · K
cTT = [-4.09830, 0.250037e-2, 0.275187e-6, -0.653398e-8, -0.310589e-9, 0.127747e-11] # cm^j · K^2
cH = [-0.140463e-7, 0.839350e-11, -0.190929e-14, -0.121399e-16, -0.898863e-18, 0.364662e-20] # cm^j · %^-1
cHH = [ 0.543605e-12, 0.112802e-15, -0.229979e-19, -0.191450e-21, -0.120352e-22, 0.500955e-25] # cm^j · %^-2
cp = [ 0.266898e-8, 0.273629e-14, 0.463466e-17, -0.916894e-23, 0.136685e-21, 0.413687e-23] # cm^j · Pa^-1
cpp = [ 0.610706e-17, 0.116620e-21, 0.244736e-24, -0.497682e-26, 0.742024e-29, 0.224625e-30] # cm^j · Pa^-2
cTH = [ 0.674488e-4, -0.406775e-7, 0.289063e-11, 0.819898e-13, 0.468386e-14, -0.191182e-16] # cm^j · K · %^-1
cTp = [ 0.778627e-6, 0.593296e-12, 0.145042e-14, 0.489815e-17, 0.327941e-19, 0.128020e-21] # cm^j · K · Pa^-1
cHp = [-0.211676e-15, 0.487921e-20, -0.682545e-23, 0.942802e-25, -0.946422e-27, -0.153682e-29] # cm^j · %^-1 · Pa^-1
σref = 1e4/4.8 # cm^−1, reference wavenumber
Tref = 273.15+17.5 # K, reference temperature
pref = 75000 # Pa, reference pressure
Href = 10 #%, reference humdity in percent
########################################################################
if l > 6.35 and l <= 15.05:
########################################################################
# model parameters for wavelengths between 7.5 and 14.1 microns
cref = [ 0.199885e-3, 0.344739e-9, -0.273714e-12, 0.393383e-15, -0.569488e-17, 0.164556e-19] # cm^j
cT = [ 0.593900e-1, -0.172226e-5, 0.237654e-8, -0.381812e-11, 0.305050e-14, -0.157464e-16] # cm^j · K
cTT = [-6.50355, 0.103830e-1, -0.139464e-4, 0.220077e-7, -0.272412e-10, 0.126364e-12] # cm^j · K^2
cH = [-0.221938e-7, 0.347377e-10, -0.465991e-13, 0.735848e-16, -0.897119e-19, 0.380817e-21] # cm^j · %^-1
cHH = [ 0.393524e-12, 0.464083e-15, -0.621764e-18, 0.981126e-21, -0.121384e-23, 0.515111e-26] # cm^j · %^-2
cp = [ 0.266809e-8, 0.695247e-15, 0.159070e-17, -0.303451e-20, -0.661489e-22, 0.178226e-24] # cm^j · Pa^-1
cpp = [ 0.610508e-17, 0.227694e-22, 0.786323e-25, -0.174448e-27, -0.359791e-29, 0.978307e-32] # cm^j · Pa^-2
cTH = [ 0.106776e-3, -0.168516e-6, 0.226201e-9, -0.356457e-12, 0.437980e-15, -0.194545e-17] # cm^j · K · %^-1
cTp = [ 0.77368e-6, 0.216404e-12, 0.581805e-15, -0.189618e-17, -0.198869e-19, 0.589381e-22] # cm^j · K · Pa^-1
cHp = [-0.206365e-15, 0.300234e-19, -0.426519e-22, 0.684306e-25, -0.467320e-29, 0.126117e-30] # cm^j · %^-1 · Pa^-1
σref = 1e4/10.1 # cm^−1, reference wavenumber
Tref = 273.15+17.5 # K, reference temperature
pref = 75000 # Pa, reference pressure
Href = 10 #%, reference humdity in percent
########################################################################
if l > 15.05:
########################################################################
# model parameters for wavelengths between 16 and 24 microns
cref = [ 0.199436e-3, 0.299123e-8, -0.214862e-10, 0.143338e-12, 0.122398e-14, -0.114628e-16] # cm^j
# something seems to be wrong with cT...
cT = [ 0.621723e-1, -0.177074e-4, 0.152213e-6, -0.954584-9, -0.996706e-11, 0.921476e-13] # cm^j · K
cTT = [-23.2409, 0.108557, -0.102439e-2, 0.634072e-5, 0.762517e-7, -0.675587e-9 ] # cm^j · K^2
cH = [-0.772707e-7, 0.347237e-9, -0.272675e-11, 0.170858e-13, 0.156889e-15, -0.150004e-17] # cm^j · %^-1
cHH = [-0.326604e-12, 0.463606e-14, -0.364272e-16, 0.228756e-18, 0.209502e-20, -0.200547e-22] # cm^j · %^-2
cp = [ 0.266827e-8, 0.120788e-14, 0.522646e-17, 0.783027e-19, 0.753235e-21, -0.228819e-24] # cm^j · Pa^-1
cpp = [ 0.613675e-17, 0.585494e-22, 0.286055e-24, 0.425193e-26, 0.413455e-28, -0.812941e-32] # cm^j · Pa^-2
cTH = [ 0.375974e-3, -0.171849e-5, 0.146704e-7, -0.917231e-10, -0.955922e-12, 0.880502e-14] # cm^j · K · %^-1
cTp = [ 0.778436e-6, 0.461840e-12, 0.306229e-14, -0.623183e-16, -0.161119e-18, 0.800756e-20] # cm^j · K · Pa^-1
cHp = [-0.272614e-15, 0.304662e-18, -0.239590e-20, 0.149285e-22, 0.136086e-24, -0.130999e-26] # cm^j · %^-1 · Pa^-1
σref = 1e4/20 # cm^−1, reference wavenumber
Tref = 273.15+17.5 # K, reference temperature
pref = 75000 # Pa, reference pressure
Href = 10 #%, reference humdity in percent
########################################################################
s = 1e4/l # cm^-1
n = 1
H = H*100
for j in range(0, 6):
n += ( cref[j] + cT[j]*(1/T-1/Tref) + cTT[j]*(1/T-1/Tref)**2
+ cH[j]*(H-Href) + cHH[j]*(H-Href)**2
+ cp[j]*(p-pref) + cpp[j]*(p-pref)**2
+ cTH[j]*(1/T-1/Tref)*(H-Href)
+ cTp[j]*(1/T-1/Tref)*(p-pref)
+ cHp[j]*(H-Href)*(p-pref) ) * (s-sref)**j
return n
def K2C(K):
'''
Kelvin to Celsius conversion
'''
return K - 273.15
def hPA2mmHg(hPa):
'''
hPa to mmHg conversion
'''
return 0.7500616827 * hPa