mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-29 20:51:05 +00:00
First commit
This commit is contained in:
Joost van den Born committed
commit
10bf0fa17d
7 files changed
+14613
No files matched your search
+179
@@ -0,0 +1,179 @@
|
||||
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
|
||||
Binary file not shown.
File diff suppressed because it is too large.
Load diff
@@ -0,0 +1,139 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"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."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"#!/usr/bin/env Python\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from atmosphericModel import *\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# Parameters at Cerro Armazones\n",
|
||||
"T = 279.65 # k\n",
|
||||
"P = 71200 # Pa\n",
|
||||
"H = 0.22\n",
|
||||
"xc = 450 # ppm\n",
|
||||
"lat = -24.5983 # degrees\n",
|
||||
"h = 3064 # m\n",
|
||||
"l1 = 1.49 # micron\n",
|
||||
"l2 = 1.78 # micron\n",
|
||||
"\n",
|
||||
"# Initializing dispersion model\n",
|
||||
"at = Ciddor1996(p=P, T=T, H=H, xc=xc)\n",
|
||||
"\n",
|
||||
"# Calculating indices of refraction for l1 and l2\n",
|
||||
"n1 = at.n_tph(l=l1, T=T, p=P, RH=H, xc=xc)\n",
|
||||
"n2 = at.n_tph(l=l2, T=T, p=P, RH=H, xc=xc)\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can also calculate the density of the air, which will be useful when calculating the scale height of the atmosphere, necessary for the refraction calculations."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"# Density of the atmosphere (following CIPM-81/91 equations)\n",
|
||||
"rho = at.rho(p=P, T=T, RH=H, xc=xc)\n",
|
||||
"\n",
|
||||
"# Initializing refraction model\n",
|
||||
"ref = refraction(lat, h, P, rho)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The refraction for 1.49 um light at a zenith angle of 45 degrees, as calculated by the Cassini model, is found to be:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"The refraction is 40.697 arc seconds\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Calculating refraction for edge FoV\n",
|
||||
"zenith_obs = 45\n",
|
||||
"refraction_obs = ref.cassini(n1, zenith_obs) * 3.6e3\n",
|
||||
"print ('The refraction is %.03f arc seconds' %(refraction_obs))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Dispersion is simply the refraction differential between two wavelengths. This is calculated with this module by"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"The dispersion is 31.019 milli arc seconds\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Initializing dispersion object\n",
|
||||
"disp = dispersion(lat, h, P, rho)\n",
|
||||
"dispersion_obs = disp.cassini(n1, n2, zenith=zenith_obs)*3.6e6\n",
|
||||
"print ('The dispersion is %.03f milli arc seconds' %(dispersion_obs))"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.7.1"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
@@ -0,0 +1,808 @@
|
||||
#!/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 Ciddor1996:
|
||||
def __init__(self, p=101325, T=288.15, H=0, xc=450):
|
||||
self.p = p # Pressure in Pa
|
||||
self.T = T # Temperature in Kelvin
|
||||
self.H = H # Relative humidity on a scale from 0 to 1
|
||||
self.xc = xc # ppm of CO2
|
||||
pass
|
||||
|
||||
|
||||
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=1/1.7564043320974863, xc=xc) # Moist air component at standard 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) )
|
||||
|
||||
|
||||
class error(Ciddor1996):
|
||||
def __init__(self, p=101325, T=288.15, H=0, xc=450, dT=0.2, dP=20, dH=0.02, dCO2=20, dl=0.001):
|
||||
Ciddor1996.__init__(self, p, T, H)
|
||||
self.dT = 0#dT # in K
|
||||
self.dP = 0#dP # in Pa
|
||||
self.dH = 0#dH # relative humidity between 0 and 1
|
||||
self.dCO2 = 0#dCO2 # in ppm
|
||||
self.dl = 0#dl # in micron
|
||||
|
||||
|
||||
def dn_as(self, l, dl=None):
|
||||
'''
|
||||
Formula that calculate the error on the refractive index for standard
|
||||
atmospheric conditions, i.e. 15 C, 101,325 Pa, 0% RH, 450 ppm CO2.
|
||||
It returns dn_as, not d(n_as - 1) * 1e8!
|
||||
'''
|
||||
if not dl:
|
||||
dl = self.dl
|
||||
|
||||
_k0 = 238.0185 # um^-2
|
||||
_k1 = 5792105 # um^-2
|
||||
_k2 = 57.362 # um^-2
|
||||
_k3 = 167917 # um^-2
|
||||
_dn_as = np.sqrt((( - (2 * _k1 * l ) / ((_k0 * l**2 - 1)**2) - (2 * _k3 * l) / ((_k2 * l**2 - 1)**2) )*dl)**2)
|
||||
return _dn_as / 1e8
|
||||
|
||||
|
||||
def dn_axs(self, l, xc, dCO2=None, dl=None):
|
||||
'''
|
||||
Formula that calculates the error on the refractive index for standard
|
||||
atmospheric conditions, i.e. 15 C, 101,325 Pa, 0% RH, but with xc ppm of CO2.
|
||||
It returns dn_axs, not d(n_axs - 1)
|
||||
'''
|
||||
if not dCO2:
|
||||
dCO2 = self.dCO2
|
||||
if not dl:
|
||||
dl = self.dl
|
||||
_dn_as = self.dn_as(l, dl)
|
||||
_n_as = self.n_as(l)
|
||||
_dn_axs = np.sqrt( ( (_n_as - 1) * 0.534e-6 * dCO2 )**2 + ( (0.534e-6 * (xc - 450) + 1) * _dn_as )**2 )
|
||||
return _dn_axs
|
||||
|
||||
|
||||
def dn_ws(self, l, dl=None):
|
||||
'''
|
||||
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
|
||||
'''
|
||||
if not dl:
|
||||
dl = self.dl
|
||||
_cf = 1.022
|
||||
_w0 = 295.235 # um^-2
|
||||
_w1 = 2.6422 # um^-2
|
||||
_w2 = -0.032380 # um^-4
|
||||
_w3 = 0.004028 # um^-6
|
||||
|
||||
_dndl = _cf * (-1 * _w1 * 2/l**3 - _w2 * 4/l**5 - _w3 * 6/l**7)
|
||||
return np.sqrt( (_dndl * dl)**2 ) / 1e8
|
||||
|
||||
|
||||
def df(self, p, T, dP=None, dT=None):
|
||||
'''
|
||||
Uncertainty in enhancement factor of
|
||||
water vapour in air.
|
||||
t in degrees celsius
|
||||
'''
|
||||
if not dP:
|
||||
dP = self.dP
|
||||
if not dT:
|
||||
dT = self.dT
|
||||
|
||||
_t = T - 273.15
|
||||
_b = 3.14e-8 # Pa^-1
|
||||
_c = 5.6e-7 # deg C^-2
|
||||
return np.sqrt( (_b*dP)**2 + (2*_c*_t*dT)**2 )
|
||||
|
||||
|
||||
def dsvp(self, T, dT=None):
|
||||
'''
|
||||
Uncertainty in saturation vapour pressure in Pa.
|
||||
'''
|
||||
if not dT:
|
||||
dT = self.dT
|
||||
_A = 1.2378847e-5 # K^-2
|
||||
_B = -1.9121316e-2 # K^-2
|
||||
_C = 33.93711047
|
||||
_D = -6.3431645e3 # K^-2
|
||||
_p1 = 2*_A*T + _B - _D/(T**2)
|
||||
_p2 = np.exp(_A*T**2 + _B*T + _C + _D/T)
|
||||
return np.sqrt( (_p1 * _p2 * dT)**2 )
|
||||
|
||||
|
||||
def dxw(self, p, T, RH, dP=None, dT=None, dRH=None):
|
||||
'''
|
||||
Uncertainty in molar fraction of water vapour in moist air.
|
||||
'''
|
||||
if not dP:
|
||||
dP = self.dP
|
||||
if not dT:
|
||||
dT = self.dT
|
||||
if not dRH:
|
||||
dRH = self.dH
|
||||
|
||||
_f = self.f(p, T)
|
||||
_svp = self.svp(T)
|
||||
|
||||
_df = self.df(p=p, T=T, dP=dP, dT=dT)
|
||||
_dsvp = self.dsvp(T=T, dT=dT)
|
||||
_p1 = (RH * _svp / p) * _df
|
||||
_p2 = (_f * _svp / p) * dRH
|
||||
_p3 = (_f * RH / p) * _dsvp
|
||||
_p4 = (-1 * _f * RH * _svp / p**2) * dP
|
||||
return np.sqrt( _p1**2 + _p2**2 + _p3**2 + _p4**2)
|
||||
|
||||
|
||||
def dZ(self, p, T, RH, xw='', dP=None, dT=None, dRH=None):
|
||||
'''
|
||||
Uncertainty in compressibility of moist air
|
||||
'''
|
||||
if not dP:
|
||||
dP = self.dP
|
||||
if not dT:
|
||||
dT = self.dT
|
||||
if not dRH:
|
||||
dRH = self.dH
|
||||
_dt = dT
|
||||
_dxw = self.dxw(p=p, T=T, RH=RH, dP=dP, dT=dT, dRH=dRH)
|
||||
|
||||
_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
|
||||
|
||||
_dZdp = -1/T * ( _a0 + _a1*_t + _a2*_t**2 + (_b0 + _b1*_t)*_xw + (_c0 + _c1*_t)*_xw**2) + (2*p/T**2) * (_d + _e*_xw**2)
|
||||
_dZdT = p/T**2 * ( _a0 + _a1*_t + _a2*_t**2 + (_b0 + _b1*_t)*_xw + (_c0 + _c1*_t)*_xw**2) - 2*p**2/T**3 * (_d + _e*_xw**2)
|
||||
_dZdt = -1*p/T * (_a1 + 2*_a2*_t + _b1*_xw + _c1*_xw**2)
|
||||
_dZdxw = -1*p/T * ( (_b0 + _b1*_t) + 2*(_c0 + _c1*_t)*_xw ) + 2 * (p/T)**2 * _e * _xw
|
||||
_dZ = np.sqrt( (_dZdp * dP)**2 + (_dZdT * dT)**2 + (_dZdt * _dt)**2 + (_dZdxw * _dxw)**2 )
|
||||
return _dZ
|
||||
|
||||
|
||||
def drho(self, p, T, RH, xc, xw='', dP=None, dT=None, dRH=None, dCO2=None):
|
||||
'''
|
||||
Calculates the uncertainty in the density of air for the given parameters.
|
||||
T in Kelvin
|
||||
p in Pascal
|
||||
RH between 0 and 1
|
||||
xc in ppm
|
||||
'''
|
||||
if not dP:
|
||||
dP = self.dP
|
||||
if not dT:
|
||||
dT = self.dT
|
||||
if not dRH:
|
||||
dRH = self.dH
|
||||
if not dCO2:
|
||||
dCO2 = self.dCO2
|
||||
|
||||
_Mw = 0.018015 # kg/mol
|
||||
_R = 8.314510 # J mol^-1 K^-1, gas constant
|
||||
_Ma = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) ) # Molar mass in kg/mol of dry air containing xc ppm of CO2
|
||||
if xw == '':
|
||||
_xw = self.xw(p, T, RH)
|
||||
else:
|
||||
_xw = xw
|
||||
_Z = self.Z(p=p, T=T, RH=RH, xw=_xw)
|
||||
|
||||
_dMa = np.sqrt( (12.011e-9*dCO2)**2 )
|
||||
_dZ = self.dZ(p=p, T=T, RH=RH, xw=_xw, dP=dP, dT=dT, dRH=dRH)
|
||||
_dxw = self.dxw(p=p, T=T, RH=RH, dP=dP, dT=dT, dRH=dRH)
|
||||
|
||||
_drhodp = _Ma / (_Z * _R * T) * (1 - _xw * (1 - _Mw / _Ma))
|
||||
_drhodMa = p / (_Z * _R * T) * (1 - _xw)
|
||||
_drhodZ = -1 * (p * _Ma) / (_Z**2 * _R * T) * (1 - _xw * (1 - _Mw / _Ma))
|
||||
_drhodT = -1 * (p * _Ma) / (_Z * _R * T**2) * (1 - _xw * (1 - _Mw / _Ma))
|
||||
_drhodxw = p / (_Z * _R * T) * (_Mw - _Ma)
|
||||
|
||||
return np.sqrt( (_drhodp * dP)**2 + (_drhodMa * _dMa)**2 + (_drhodZ * _dZ)**2 + (_drhodT * dT)**2 + (_drhodxw * _dxw)**2 )
|
||||
|
||||
|
||||
def dn_tph(self, l, T, p, RH, xc, dl=None, dT=None, dP=None, dRH=None, dCO2=None):
|
||||
'''
|
||||
Calculates the uncertainty in the refraction of the atmosphere for any conditions.
|
||||
'''
|
||||
if not dP:
|
||||
dP = self.dP
|
||||
if not dT:
|
||||
dT = self.dT
|
||||
if not dRH:
|
||||
dRH = self.dH
|
||||
if not dCO2:
|
||||
dCO2 = self.dCO2
|
||||
|
||||
_n_a = self.n_axs(l, xc)
|
||||
_n_w = self.n_ws(l)
|
||||
_dn_a = self.dn_axs(l=l, xc=xc, dl=dl, dCO2=dCO2)
|
||||
_dn_w = self.dn_ws(l=l, dl=dl)
|
||||
|
||||
_La = (_n_a**2 - 1) / (_n_a**2 + 2)
|
||||
_Lw = (_n_w**2 - 1) / (_n_w**2 + 2)
|
||||
_dLa = np.sqrt( ( (6*_n_a) / (_n_a**2+2)**2 * _dn_a)**2 )
|
||||
_dLw = np.sqrt( ( (6*_n_w) / (_n_w**2+2)**2 * _dn_w)**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
|
||||
_pw = RH*self.svp(T)
|
||||
_pa = p - _pw
|
||||
_rho_a = self.rho(p=_pa, 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.
|
||||
|
||||
_dpa = np.sqrt( (dP)**2 + (-1*self.svp(T) * dRH)**2 + (-1 * RH * self.dsvp(T=T, dT=dT) )**2 )
|
||||
_dpw = np.sqrt( (RH * self.dsvp(T=T, dT=dT))**2 + (self.svp(T) * dRH)**2 )
|
||||
_drho_a = self.drho(p=_pa, T=T, RH=RH, xc=xc, dP=_dpa, dT=dT, dRH=dRH, dCO2=dCO2)
|
||||
_drho_axs= 0
|
||||
|
||||
if RH == 0:
|
||||
# nan is returned for _rho_w if we put in a (partial) pressure of 0
|
||||
_rho_w = 0
|
||||
_drho_w = 0
|
||||
else:
|
||||
# Partial pressure of water vapour is p_water = h*svp
|
||||
_rho_w = self.rho(p=_pw, T=T, RH=RH, xc=xc) # Moist air component for the actual conditions.
|
||||
_drho_w = self.drho(p=_pw, T=T, RH=RH, xc=xc, dP=_dpw, dT=dT, dRH=dRH, dCO2=dCO2)
|
||||
_rho_ws = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc) # Moist air component at standard conditions.
|
||||
_drho_ws = 0
|
||||
|
||||
|
||||
_dLdL_a = _rho_a / _rho_axs
|
||||
_dLdL_w = _rho_w / _rho_ws
|
||||
_dLdrho_a = _La / _rho_axs
|
||||
_dLdrho_axs = - (_rho_a * _La) / (_rho_axs**2)
|
||||
_dLdrho_w = _Lw / _rho_ws
|
||||
_dLdrho_ws = - (_rho_w * _Lw) / _rho_ws**2
|
||||
|
||||
_Ltotal = ( _rho_a / _rho_axs ) * _La + ( _rho_w / _rho_ws ) * _Lw
|
||||
_dLtotal = np.sqrt( (_dLdL_a * _dLa)**2 + (_dLdL_w * _dLw)**2 + (_dLdrho_a * _drho_a)**2
|
||||
+ (_dLdrho_axs * _drho_axs)**2 + (_dLdrho_w * _drho_w)**2 + (_dLdrho_ws * _drho_ws)**2 )
|
||||
|
||||
_nLL = np.sqrt( (1 + 2 * _Ltotal) / (1 - _Ltotal) )
|
||||
_dnLL = np.sqrt( ( (3*np.sqrt(1 - _Ltotal )) / ( 2 * (1 - _Ltotal)**2 * np.sqrt(1 + 2*_Ltotal) ) * _dLtotal )**2 )
|
||||
|
||||
return (_nLL, _dnLL)
|
||||
|
||||
|
||||
class dispersion:
|
||||
def __init__(self, lat, h, P, rho):
|
||||
self.h = h # height above sea level in meters
|
||||
self.lat = lat # Latitude of the observer in degrees
|
||||
self.rc = self.rc() * 1000 # Radius of curvature of the earth at the observer
|
||||
self.H = self.H(P, rho) # Reduced height of the atmosphere assuming ideal gas law
|
||||
# self.H2 = self.H_isotherm(T) # Reduced height of isotherm atmosphere
|
||||
|
||||
|
||||
|
||||
def planeParallel(self, n1, n2, zenith):
|
||||
'''
|
||||
Refraction based on a flat atmosphere, easily derived from Snell's law.
|
||||
'''
|
||||
# _theta1 = np.degrees( np.arcsin( np.sin( np.deg2rad(zenith) / n1) ) )
|
||||
# _theta2 = np.degrees( np.arcsin( np.sin( np.deg2rad(zenith) / n2) ) )
|
||||
_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
|
||||
'''
|
||||
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]:
|
||||
|
||||
_refractionIntegralList = []
|
||||
_nList = []
|
||||
|
||||
for _i in range(len(_h)):
|
||||
_R = R0 + _h[_i]*1000
|
||||
_n_h = 1 + (_n - 1) / _Dens[0] * _Dens[_i] # Gladstone-Dale relation between refractive index and density
|
||||
_nList.append(_n_h)
|
||||
_refractionIntegralList.append( (_n * R0 * np.sin(z0))/(_n_h*np.sqrt((_n_h*_R)**2 - (_n * R0 * np.sin(z0))**2 )) )
|
||||
_dn = (_nList[0] - _nList[-1]) / ((len(_h)-1))
|
||||
|
||||
_rList.append(np.sum(np.degrees(np.asarray(_refractionIntegralList) * _dn)))
|
||||
|
||||
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 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 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)
|
||||
|
||||
class refraction:
|
||||
def __init__(self, lat, h, P, rho):
|
||||
self.h = h # height above sea level in meters
|
||||
self.lat = lat # Latitude of the observer in degrees
|
||||
self.rc = self.rc() * 1000 # Radius of curvature of the earth at the observer
|
||||
self.H = self.H(P, rho) # Reduced height of the atmosphere assuming ideal gas law
|
||||
# self.H2 = self.H_isotherm(T) # Reduced height of isotherm atmosphere
|
||||
|
||||
|
||||
|
||||
def planeParallel(self, n, zenith):
|
||||
'''
|
||||
Refraction based on a flat atmosphere, easily derived from Snell's law.
|
||||
'''
|
||||
# _theta1 = np.degrees( np.arcsin( np.sin( np.deg2rad(zenith) / n) ) )
|
||||
_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
|
||||
|
||||
_refractionIntegralList = []
|
||||
_nList = []
|
||||
z0 = np.deg2rad(z0)
|
||||
|
||||
# Integration of the refraction integral over the atmosphere layers.
|
||||
for _i in range(len(_h)):
|
||||
_R = R0 + _h[_i]*1000
|
||||
_n_h = 1 + (n0 - 1) / _Dens[0] * _Dens[_i] # Gladstone-Dale relation between refractive index and density
|
||||
_nList.append(_n_h)
|
||||
_refractionIntegralList.append( (n0 * R0 * np.sin(z0))/(_n_h*np.sqrt((_n_h*_R)**2 - (n0 * R0 * np.sin(z0))**2 )) )
|
||||
_dn = (_nList[0] - _nList[-1]) / ((len(_h)-1))
|
||||
|
||||
return np.sum(np.degrees(np.asarray(_refractionIntegralList) * _dn))
|
||||
|
||||
|
||||
|
||||
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 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 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)
|
||||
|
||||
BIN
Binary file not shown.
@@ -0,0 +1,213 @@
|
||||
#!/usr/bin/env Python
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.axes_grid1 import host_subplot
|
||||
from atmosphericModel import *
|
||||
from pdb import set_trace
|
||||
import pandas as pd
|
||||
|
||||
'''
|
||||
@ Author: Joost van den Born
|
||||
@ Contact: born@astron.nl
|
||||
@ Description: Modeling the atmospheric refraction using multiple
|
||||
models of the atmosphere.
|
||||
|
||||
Models included:
|
||||
- Plane-parallel atmosphere
|
||||
- Cassini's model
|
||||
- Tan5 model (A*tan(z) + B*tan^3(z) + C*tan^5(z))
|
||||
- Oriani's theorem (A*tan(z) + B*tan^3(z))
|
||||
- Erf model (from Corbard et al., 2019)
|
||||
- I do not yet consider full integration of the refraction integral.
|
||||
'''
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
# Conditions at Cerro Armazones
|
||||
# T = 279.65 # k
|
||||
# P = 71200 # Pa
|
||||
# H = 0.22
|
||||
# xc = 450 # ppm
|
||||
T = 288.15 # k
|
||||
P = 101325 # Pa
|
||||
H = 0.0
|
||||
xc = 314 # ppm
|
||||
lat = -24.5983 # degrees
|
||||
h = 3064 # m
|
||||
l1 = 1.49 # micron
|
||||
l2 = 1.78 # micron
|
||||
zenith = 30 # degrees
|
||||
|
||||
# Calculate the index of refraction
|
||||
at = Ciddor1996(p=P, T=T, H=H, xc=xc)
|
||||
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)
|
||||
rho = at.rho(p=P, T=T, RH=H, xc=xc)
|
||||
|
||||
refractionAnalysis = True
|
||||
dispersionAnalysis = True
|
||||
tableCompare = True
|
||||
if dispersionAnalysis:
|
||||
# Initiate the refraction class
|
||||
ref = dispersion(lat, h, P, rho)
|
||||
|
||||
# Compare the refraction models to check if the output is okay
|
||||
print ('Plane - parallel atmosphere dispersion: %s mas' %(ref.planeParallel(n1, n2, zenith) * 3.6e6))
|
||||
print ('Cassini dispersion: %s mas' %(ref.cassini(n1, n2, zenith) * 3.6e6))
|
||||
print ('Tan5 dispersion: %s mas' %(ref.tan5(n1, n2, zenith) * 3.6e6))
|
||||
print ('Oriani dispersion: %s mas' %(ref.oriani(n1, n2, zenith) * 3.6e6))
|
||||
print ('Corbard dispersion: %s mas' %(ref.corbard(n1, n2, zenith)*3.6e6))
|
||||
print ('Full integration: %s mas' %(ref.refractionIntegral(n1, n2, zenith)*3.6e6))
|
||||
|
||||
|
||||
# Now we calculate the dispersion over a wide range of zenith angles
|
||||
# for ideal gas atmosphere
|
||||
zenithRange = np.linspace(5, 85, 1000)
|
||||
ppList = ref.planeParallel(n1, n2, zenithRange) * 3.6e6
|
||||
cassiniList = ref.cassini(n1, n2, zenithRange) * 3.6e6
|
||||
tan5List = ref.tan5(n1, n2, zenithRange) * 3.6e6
|
||||
orianiList = ref.oriani(n1, n2, zenithRange) * 3.6e6
|
||||
corbardList = [ref.corbard(n1, n2, z) * 3.6e6 for z in zenithRange]
|
||||
fullIntegrationList = [ref.refractionIntegral(n1, n2, z) * 3.6e6 for z in zenithRange]
|
||||
|
||||
# for isothermal atmosphere
|
||||
cassiniList_isotherm = ref.cassini(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
|
||||
tan5List_isotherm = ref.tan5(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
|
||||
orianiList_isotherm = ref.oriani(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6
|
||||
corbardList_isotherm = [ref.corbard(n1, n2, z, reduced_height=ref.H_isotherm(T)) * 3.6e6 for z in zenithRange]
|
||||
|
||||
# check when the variation between the different models pass different values
|
||||
dispersionWidth = []
|
||||
widthList = []
|
||||
for z in range(len(zenithRange)):
|
||||
compList = [cassiniList[z], tan5List[z], orianiList[z], corbardList[z], cassiniList_isotherm[z], tan5List_isotherm[z], orianiList_isotherm[z], corbardList_isotherm[z]]#, fullIntegrationList[z]]
|
||||
width = np.max(compList) - np.min(compList)
|
||||
widthList.append(width)
|
||||
dispersionWidth.append((zenithRange[z], width))
|
||||
|
||||
masBorder = np.argmax(np.asarray(widthList)>1)
|
||||
p1masBorder = np.argmax(np.asarray(widthList)>0.1)
|
||||
p01masBorder = np.argmax(np.asarray(widthList)>0.01)
|
||||
p001masBorder = np.argmax(np.asarray(widthList)>0.001)
|
||||
|
||||
|
||||
# Plotting it all, with vertical lines where the variation between the models exceed 1, 0.1, 0.01 and 0.001 mas.
|
||||
fig = plt.figure(figsize=(11,6))
|
||||
ax = fig.add_subplot(111)
|
||||
ax.set_xlim([30,85])
|
||||
ax.set_xlabel('Zenith (degrees)')
|
||||
ax.set_ylabel('Dispersion (mas)')
|
||||
ax.grid(True, alpha=0.3)
|
||||
ax.axvline(dispersionWidth[masBorder][0], linestyle='--', color='k', linewidth=0.5)
|
||||
ax.text(dispersionWidth[masBorder][0]-0.5, 402, '1 mas', ha='center', va='bottom', rotation='vertical')
|
||||
ax.axvline(dispersionWidth[p1masBorder][0], linestyle='--', color='k', linewidth=0.5)
|
||||
ax.text(dispersionWidth[p1masBorder][0]-0.5, 402, '0.1 mas', ha='center', va='bottom', rotation='vertical')
|
||||
ax.axvline(dispersionWidth[p01masBorder][0], linestyle='--', color='k', linewidth=0.5)
|
||||
ax.text(dispersionWidth[p01masBorder][0]-0.5, 402, '10 $\mu$as', ha='center', va='bottom', rotation='vertical')
|
||||
ax.axvline(dispersionWidth[p001masBorder][0], linestyle='--', color='k', linewidth=0.5)
|
||||
ax.text(dispersionWidth[p001masBorder][0]-0.5, 402, '1 $\mu$as', ha='center', va='bottom', rotation='vertical')
|
||||
|
||||
ax.plot(zenithRange, ppList, 'k-', label='Plane-parallel', alpha=0.2)
|
||||
ax.plot(zenithRange, cassiniList, 'g-', label='Cassini')
|
||||
ax.plot(zenithRange, tan5List, 'b-', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$')
|
||||
ax.plot(zenithRange, orianiList, 'y-', label='$Atan(z) - Btan^3(z)$')
|
||||
ax.plot(zenithRange, corbardList, 'r-', label='Error function')
|
||||
ax.plot(zenithRange, fullIntegrationList, 'c-', label='Full integration')
|
||||
|
||||
ax.plot(zenithRange, cassiniList_isotherm, 'g--', label='Cassini (isotherm)')
|
||||
ax.plot(zenithRange, tan5List_isotherm, 'b--', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$ (isotherm)')
|
||||
ax.plot(zenithRange, orianiList_isotherm, 'y--', label='$Atan(z) - Btan^3(z)$ (isotherm)')
|
||||
ax.plot(zenithRange, corbardList_isotherm, 'r--', label='Error function (isotherm)')
|
||||
|
||||
ax.legend()
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
|
||||
if refractionAnalysis:
|
||||
# Initiate the refraction class
|
||||
ref = refraction(lat, h, P, rho)
|
||||
|
||||
# Compare the refraction models to check if the output is okay
|
||||
print ('Plane - parallel atmosphere refraction: %s arcsec' %(ref.planeParallel(n1, zenith) * 3.6e3))
|
||||
print ('Cassini refraction: %s arcsec' %(ref.cassini(n1, zenith) * 3.6e3))
|
||||
print ('Tan5 refraction: %s arcsec' %(ref.tan5(n1, zenith) * 3.6e3))
|
||||
print ('Oriani refraction: %s arcsec' %(ref.oriani(n1, zenith) * 3.6e3))
|
||||
print ('Corbard refraction: %s arcsec' %(ref.corbard(n1, zenith)*3.6e3))
|
||||
print ('Full integration: %s mas' %(ref.refractionIntegral(n1, zenith)*3.6e3))
|
||||
|
||||
|
||||
# Now we calculate the dispersion over a wide range of zenith angles
|
||||
# for ideal gas atmosphere
|
||||
zenithRange = np.linspace(0, 88, 1000)
|
||||
ppList = ref.planeParallel(n1, zenithRange) * 3.6e3
|
||||
cassiniList = ref.cassini(n1, zenithRange) * 3.6e3
|
||||
tan5List = ref.tan5(n1, zenithRange) * 3.6e3
|
||||
orianiList = ref.oriani(n1, zenithRange) * 3.6e3
|
||||
corbardList = [ref.corbard(n1, z) * 3.6e3 for z in zenithRange]
|
||||
fullIntegrationList = [ref.refractionIntegral(n1, z) * 3.6e3 for z in zenithRange]
|
||||
|
||||
# for isothermal atmosphere
|
||||
cassiniList_isotherm = ref.cassini(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3
|
||||
tan5List_isotherm = ref.tan5(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3
|
||||
orianiList_isotherm = ref.oriani(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3
|
||||
corbardList_isotherm = [ref.corbard(n1, z, reduced_height=ref.H_isotherm(T)) * 3.6e3 for z in zenithRange]
|
||||
|
||||
# Plotting it all, with vertical lines where the variation between the models exceed 1, 0.1, 0.01 and 0.001 arcsec.
|
||||
fig = plt.figure(figsize=(11,6))
|
||||
ax = fig.add_subplot(111)
|
||||
ax.set_xlim([80,88])
|
||||
ax.set_xlabel('Zenith (degrees)')
|
||||
ax.set_ylabel('Refraction (arcsec)')
|
||||
ax.grid(True, alpha=0.3)
|
||||
|
||||
ax.plot(zenithRange, ppList, 'k-', label='Plane-parallel', alpha=0.2)
|
||||
ax.plot(zenithRange, cassiniList, 'g-', label='Cassini')
|
||||
ax.plot(zenithRange, tan5List, 'b-', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$')
|
||||
ax.plot(zenithRange, orianiList, 'y-', label='$Atan(z) - Btan^3(z)$')
|
||||
ax.plot(zenithRange, corbardList, 'r-', label='Error function')
|
||||
ax.plot(zenithRange, fullIntegrationList, 'c-', label='Full integration')
|
||||
|
||||
ax.plot(zenithRange, cassiniList_isotherm, 'g--', label='Cassini (isotherm)')
|
||||
ax.plot(zenithRange, tan5List_isotherm, 'b--', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$ (isotherm)')
|
||||
ax.plot(zenithRange, orianiList_isotherm, 'y--', label='$Atan(z) - Btan^3(z)$ (isotherm)')
|
||||
ax.plot(zenithRange, corbardList_isotherm, 'r--', label='Error function (isotherm)')
|
||||
|
||||
ax.legend()
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
|
||||
if tableCompare:
|
||||
refTableData = []
|
||||
dispTableData = []
|
||||
ref = refraction(lat, h, P, rho)
|
||||
disp = dispersion(lat, h, P, rho)
|
||||
zenithRange = [10, 30, 50, 60, 70, 80, 85]
|
||||
headerList = np.array(['z', 'Plane-parallel', 'Cassini', 'tan3', 'tan5', 'Error function', 'Full integration'])
|
||||
refTableData.append(headerList)
|
||||
dispTableData.append(headerList)
|
||||
|
||||
for z in zenithRange:
|
||||
ppValue = round(ref.planeParallel(n1, z) * 3.6e3, 3)
|
||||
cassiniValue = round(ref.cassini(n1, z) * 3.6e3, 3)
|
||||
tan3Value = round(ref.oriani(n1, z) * 3.6e3, 3)
|
||||
tan5Value = round(ref.tan5(n1, z) * 3.6e3, 3)
|
||||
erfValue = round(ref.corbard(n1, z) * 3.6e3, 3)
|
||||
fullIntValue = round(ref.refractionIntegral(n1, z) * 3.6e3, 3)
|
||||
refTableData.append([z, ppValue, cassiniValue, tan3Value, tan5Value, erfValue, fullIntValue])
|
||||
|
||||
ppValue = round(disp.planeParallel(n1, n2, z) * 3.6e6, 3)
|
||||
cassiniValue = round(disp.cassini(n1, n2, z) * 3.6e6, 3)
|
||||
tan3Value = round(disp.oriani(n1, n2, z) * 3.6e6, 3)
|
||||
tan5Value = round(disp.tan5(n1, n2, z) * 3.6e6, 3)
|
||||
erfValue = round(disp.corbard(n1, n2, z) * 3.6e6, 3)
|
||||
fullIntValue = round(disp.refractionIntegral(n1, n2, z) * 3.6e6, 3)
|
||||
dispTableData.append([z, ppValue, cassiniValue, tan3Value, tan5Value, erfValue, fullIntValue])
|
||||
|
||||
dfr = pd.DataFrame(np.asarray(refTableData))
|
||||
dfd = pd.DataFrame(np.asarray(dispTableData))
|
||||
print(dfr.to_string(index=False))
|
||||
print(dfd.to_string(index=False))
|
||||
|
||||
|
||||
Reference in new issue
Block a user