Sept 1, 2022
Bought 4 seat licenses for MODTRAN6.
...
http://cds-espri.ipsl.fr/tapas/project?methodName=home_en
RUNNING MODTRAN6
Feb 20, 2023
Go to
/usr/local/bin/MODTRAN6.0/bin/macos
...
Clicking on Edit function brings up a range of panels
Unit conversions
It wants Ozone content as gm/cm^2 and we have it in Dobson units. Wikipedia says 1 DU is 2.69×1020 molecules per meter squared.
Molecular weight of O_3^16 is 47.992 gm/mole. One mole is 6.023E23 molecules. One square meter is 1E4 cm^2. So:
1 DU = 2.69E20 molecules/m^2 * (1m^2/1e4 cm^2) * ( 47.992 gm/mole) * (1 mole/6.023E23 molecules) = 2.1434e-06 gm/cm^2
so to get O3 in gm/cm^2 do (DU*2.1434e-06)
250 DU = 5.3586e-04 gm/cm^2
PWV is in mm, collapsing entire column density into some depth of liquid. MODTRAN6 wants gm/cm^2 of H2O. Density of water is 1 gm/cc, so 1 cm of water is 1 gm/cc.
So for this the conversion is 10 mm PWV → 1 gm/cc. Typical Pachon PWV is a few mm. Let's use a default 5mm PWV which corresponds to 0.5 gm/cm^2.
Set up point-to-point geometry from LSST elevation to 80 km, straight up.
Data go to ~/MODTRAN6
rename .tp7 output files with X.XX***.dat where X.XX is airmass to 3 sig figures
Matlab code ReadModtran6.m will then read in and act on these files. Need at least 4 of them.
The ***.psc output file is really simple; nm, T in two columns.
Default Atmosphere at zenith
Used MODTRAN6 to generate transmission data for Rubin elevation up to 100 km. Used 5mm PWV and 250 Dobson units. Used ReadModtran6.m MATLAB code to convert from cm^-1 to nm and extract total transmission, ozone, and water.
Stable Transmission is Transmission/(water*Ozone), and is basically the Rayleigh scattering portion in blue.
The quadnotch filter avoids the O2 lines, doesn't go that far into the red.
...
Expected range of PWV is given by The GPS water vapor monitor and thermal astronomy at Gemini South (note this is an estimate from Paranal). So 5mm is a sensible default value.
Targets
Calspec standards that are both bright and rather AB-flat.
From https://www.eso.org/sci/observing/tools/standards/spectra/stanlis.html
...
Here's a list of bright A0 stars from Simbad in the right dec band:
Polar A0 stars:
Exposure times
HD 185975 is V=8.1
In 30 sec integration it gets hologram peak of 2500, sum of 60K ADU at peak, in typical seeing. A mag 6 star is 7X brighter so we'd get peak of 17K ADU, not saturated. So 30 sec seems good for V=6-8.
March 5 2023- synthetic photometry
Downloaded MuCol CalSpec data files from STSCI. https://archive.stsci.edu/hlsps/reference-atlases/cdbs/current_calspec/
data files are mucol_mod_004.fits and mucol_stis_006.fits:
mucol_stis_006.fits
mucol_mod_004.fits
...
Here is a plot of (lambda/8000)*F_lambda, which is photon spectrum for this source:
LATISS throughput estimates, ab initio
data files are in laptop directory /Users/cstubbs/data/AuxTel/throughput
...
band | Mu Col calculated extinctions mag/airmass | Mu Col magnitude extrapolated to zero airmass | photon-flat source, mag/airmass | photon-flat source extrapolated to zero airmass |
---|---|---|---|---|
1 | 0.3892 | 20.5042 | 0.3857 | 22.4303 |
2 | 0.1899 | 20.4675 | 0.1866 | 21.9940 |
3 | 0.1171 | 21.9831 | 0.1169 | 23.1427 |
4 | 0.0970 | 21.5152 | 0.0963 | 22.3052 |
Color-airmass terms:
For a color difference of Mu Col -flat = delta(B1-B2) = (22.4303 - 21.9940) - (20.5042 - 20.4675) = 0.4363 - 0.0367 = +0.3996, we get an extinction difference of 3.5 mmag/airmass.
For a color difference of Mu Col -flat = delta(B3-B4) = (21.9831-21.5152) - (23.1427-22.3052) = 0.4679 - 0.8375 = -0.3696 we get an extinction difference of 0.7 mmag/airmass
...
Compare to photon-flat source:
Quick-look comparison with observations from March 2, 2023
SEE BELOW FOR ERROR FOUND AND CORRECTED IN BAND 4
...
color | predicted | observed | Observed-predicted color difference |
---|---|---|---|
B1-B2 | 0.0367 | -0.0349 | -0.0716 |
B1-B3 | -1.4789 | -1.5849 | -0.1059 |
B1-B4 | -1.0110 | -1.0864 | -0.0753 |
A toy example of determining excess extinction.
March 13 2023
Determined effective band center for observations of Mu Col by computing, band per band, lambda_eff=sum(lambda*detected photon spectrom on ground)/sum(photon specrum on ground) for atmosphere at 1.4 airmasses and Mu Col photon spectrum and estimated throughput.
Resulting band centers are
...
define lambdanorm=lambda/500.
March 19 2023, wavelength solution for quadnotch + 300 l/mm disperser.
Images from March 16 2023 seq num 477 has a star with nice stellar atmosphere features. HD 73495 = Eta Pyxidis HR 3420. HD 73495. HIP 42334 is an A0V star.
Spectrum from RubinTV is
Copied images from 20230316 to local disk on laptop. Need to include bias frames as well as images of interest. Note that dispersion depends critically on disperser-to-CCD spacing so we should solve for it each time.
Note apparent m=0 stellar contamination at blue end of band3. We need to either subtract those out or median-filter with rotations.
Balmer lines are at
486.135
434.047
410.173
397.007
388.906
383.540
Downloaded seqnums
477
745-755 bias frames
Ran InjestAndAnalyze.m to create bias frames and full frame debiased images.
Ran Specexam2.m on frame AT_O_20230316_000477.full.debias.fits. m=0 star centroid is row 300.4 and col 1737.1
zoom on spectrum at absorption lines
March 21, 2023
...
403, 404, 405, 436, 437, 438, 444, 445, 446, 477, 478, 479, 485, 486, 487, 512, 513, 514, 536, 537, 538, 560, 561, 562, 595, 596, 597, 695, 696, 697
...
November 23, 2023 - revisions for quadband paper 1.
We've done an analysis of observations from Oct 10, 2023 on Aux Tel with masked quadband filter. There is about a 1% discrepancy between photometry and MODTRAN predictions.
...
do geographical search on latitude -30:14:40.68 longitude -70:44:57.90; = -30.245, -70.75
Well, that didn't work so well. Try this:
https://www.esrl.noaa.gov/gmd/grad/neubrew/SatO3DataTimeSeries.jsp
That works! Can get plots as well as CSV data files:
For our site:
Ozone data file for Rubin site: Ozone2023.csv
...
angle-corrected quadband throughput: Quadband.dat
Ozone for Oct 10 2023 is 300 283.7 Dobson units, interpolated to our site location.
Try to get barometric pressure right. Used https://weatherspark.com/h/d/25822/2023/10/10/Historical-Weather-on-Tuesday-October-10-2023-in-La-Serena-Chile#Figures-Pressure for barometric pressure at La Serena airport.
On Oct 10 2023 at 10 pm local the pressure at airport was 30.06 inches of mercury which is 1018 mbar. On Nov 23 2023 at 10 pm local it was the same value (precision is 0.01 inches). And at the summit we had 744.35 at that same time.
So a good pressure value to use for Oct 10 2023 is 0.74435 mbar.
Let's explore sensitivity to MODTRAN parameter choices.
Over the course of a year, barometric pressure at La Serena airport ranged from 29.8 in to 30.3 in of mercury. That's less than +- 1% variation.
PWV varies (very conservatively) from 0- 10 mm
Ozone varies from (see plot above) 250 to 300 Dobson units
Stellar colors go from -1 to 1.
So introduce perturbations that amount to mean-to-peak excursions, i.e. half the peak-to-peak value. This will show peak extinction excursions about the mean.
Let's explore sensitivity to MODTRAN parameter choices.
MATLAB program takes Star temperature, PWV, barometric pressure, and Ozone as inputs.
Pressure first:
T in 1000K. Dobson. PWV (mm). P(mbar) m1-m4. E1. E2. E3. E4. E14. E24. E34
10.5000 298.0000 5.0000 0.7327 0.0036 0.3892 0.1887 0.1194 0.1030 0.2862 0.0857 0.0164
10.5000 298.0000 5.0000 0.7400 0.0068 0.3930 0.1906 0.1205 0.1037 0.2894 0.0869 0.0168
10.5000 298.0000 10.0000 0.7400 0.0067 0.3930 0.1906 0.1205 0.1038 0.2892 0.0868 0.0167
10.5000 275.0000 5.0000 0.7400 0.0096 0.3930 0.1905 0.1195 0.1008 0.2922 0.0896 0.0187
10.5000 250.0000 10.0000 0.7400 0.0125 0.3930 0.1904 0.1185 0.0980 0.2950 0.0924 0.0206
6.0000 250.0000 10.0000 0.7400 1.0046 0.3872 0.1875 0.1184 0.0976 0.2896 0.0900 0.0208
E1 is extinction in band 1 in mag per airmass, bluest band. E14=E1-E4, etc.
We see color-extinction changes of
Final selection for Oct 10 2023:
which gives:
10.5000 283.7000 5.0000 0.7443 0.0105 0.3953 0.1916 0.1205 0.1023 0.2930 0.0893 0.0182
5.0000 283.7000 5.0000 0.7443 1.4842 0.3868 0.1875 0.1203 0.1018 0.2850 0.0857 0.0185
15.0000 283.7000 5.0000 0.7443 -0.3492 0.3975 0.1927 0.1206 0.1024 0.2951 0.0902 0.0181
8.0000 283.7000 5.0000 0.7443 0.4131 0.3929 0.1905 0.1205 0.1022 0.2908 0.0883 0.0183
Except OOPS we are using c34 as the definition of color, to reduce airmass sensitivity, recreate table with c34 in final column:
5.0000 283.7000 5.0000 0.7443 1.4842 0.3868 0.1875 0.1203 0.1018 0.2850 0.0857 0.0185 1.0369
10.5000 283.7000 5.0000 0.7443 0.0105 0.3953 0.1916 0.1205 0.1023 0.2930 0.0893 0.0182 0.6828
15.0000 283.7000 5.0000 0.7443 -0.3492 0.3975 0.1927 0.1206 0.1024 0.2951 0.0902 0.0181 0.5998
Fitting to C34 color (typcially around 0.5) gives
E14= 0.3090 -0.0231*C34
E24= 0.0964 -0.0103*C34
E34= 0.0176 +0.0009*C34
and for an A star we get E14=0.2930, E24=0.0893, E34=0.0182.
MATLAB program is here: MATLAB MODTRAN code
Nov 25, 2023.
Looked up some papers on short term ozone fluctuations Typical RMS daily flluctuations at -30 latitude is around 5%.
Also looked at ESA Sentinel satellite data for Ozone. See https://dataspace.copernicus.eu/browser/?zoom=6&lat=-29.84835&lng=-71.2672&themeId=DEFAULT-THEME&visualizationUrl=https%3A%2F%2Fsh.dataspace.copernicus.eu%2Fogc%2Fwms%2F0b0f5a61-f3d1-4c6e-8d11-4e58e2d454ef&datasetId=S5_O3_CDAS&fromTime=2023-11-24T00%3A00%3A00.000Z&toTime=2023-11-24T23%3A59%3A59.999Z&layerId=O3_VISUALIZED&demSource3D=%22MAPZEN%22&cloudCoverage=30
stubbs@g.harvard.edu, cat4Tropomi!
Sentinel 5 has ozone data with 3 day sliding averages. Can place a pin and download data file.
they use different units, moles/m^2. Value for 10 Oct 2023 was 0.13.
So the conversation factor is 2.2379e+03. The Sentinel5 ozone value for Oct 10 2023 is therefore 0.1248* 2.2379e+03 = 279.3 Dobson units. NASA said 284.
The conversion factor is 2.2379e+03 to go from Sentinal 5 units to Dobson units.
For Oct 10 2023 the Sentinel 5 value was 0.1248 (better precision obtained by downloading CSV file) which corresponds to 279.3 Dobson units. This is to be compared to NASA OMI value of 283.7.
November 26, 2023
Trying to understand how MODTRAN makes excess extinction prediction. An excellent paper is Patat et al about extinction above Paranal. They also see, even after accounting for aerosols, and excess extinction predicted by LBLRTM (line by line radiative transfer model).
Their Figure 4:
If observed - model is negative, then model extinction is > observed. Difference between 570 nm and 380 nm is 0.75 percent. We see about 1%
So this is indeed possibly a real discrepancy, and our measurement might not be the culprit.
How about aerosols? Take our band centers as 373, 445, 508, 575. Patat aerosol estimate from 2011 would imply (see figure below)
band | k | difference |
---|---|---|
B1 | 0.055 | |
B2 | 0.040 | |
B3 | 0.035 | |
B4 | 0.030 | |
E14 | 0.025 | |
E24 | 0.01 | |
E34 | 0.005 |