MATLAB MODTRAN code

%% MakeAtmosphereExtinction.m computes atmospheric transmission vs. zenith angle, barometric pressure, PWV and Ozone content

% C. Stubbs Feb 21, 2023

% Zenith-looking transmission at nominal values of 250 Dobson units of

% ozone, 5mm of precipitable water, US 1976 standard atmosphere,

% observatory elevation of 2.647 km.

% Note - no aerosols included here. That is done separately

% To convert to a different airmass and conditions, what matters is line

% integral of scatterers. That depends on airmass which in turn depends on

% zenith (or elevation) angle.

clear all

close all

ZenithAngle=0; % in degrees

ThisPWV=5; % in mm

ThisOzone=283.7; % in Dobson units

ThisPressure=0.74435;

StarTemp=10500; % stellar temp in K

% these are the values used for baseline MODTRAN modeling. DON'T CHANGE!

PWVDefault=5; % mm of water

OzoneDefault=250 ; % default Ozone column integral in Dobson units

PressureDefault=0.73269 ; % default barometric pressure in millibars

%% Import data from text file

% Script for importing data from the following text file:

%

% filename: /Users/cstubbs/MODTRAN6/1.00z_Rubin_5mmH2O.dat

%

% Auto-generated by MATLAB on 21-Feb-2023 06:45:59

%% Set up the Import Options and import the data

opts = delimitedTextImportOptions("NumVariables", 4);

% Specify range and delimiter

opts.DataLines = [2, Inf];

opts.Delimiter = ",";

% Specify column names and types

opts.VariableNames = ["Nm", "Trans", "Ozone1", "Water"];

opts.VariableTypes = ["double", "double", "double", "double"];

% Specify file level properties

opts.ExtraColumnsRule = "ignore";

opts.EmptyLineRule = "read";

% Import the data

tbl = readtable("/Users/cstubbs/MODTRAN6/1.00z_Rubin_5mmH2O.dat", opts);

%% Convert to output type, for deault values

nm = tbl.Nm;

Trans = tbl.Trans;

Ozone = tbl.Ozone1;

Water = tbl.Water;

%% determine non-ozone and non-water component

StableTrans=Trans./Ozone;

StableTrans=StableTrans./Water;

%% Clear temporary variables

clear opts tbl

%% sanity plot. These are the baseline values

figure(10)

plot(nm,Trans,'ko-')

hold on

plot(nm,Ozone,'g-','LineWidth',2)

plot(nm,Water,'b-','LineWidth',2)

plot(nm,StableTrans,'c-','LineWidth',2)

hold off

legend('Total','Qzone','Water','Stable')

grid on

shg

%% compute transmission for given zenith angle, PWV and Ozone

% If we add more material, need to use the product of transmissions. So compute a dimensionless quantity EAM, effective airmass, where EAM=(X)*(P/0.73269) and X is computed from zenith angle using formulation given above.Â

% Also compute OEAM, Ozone effective airmass, as OEAM=(X)*(O3/250) where O3 is in Dobson units. This has no pressure dependence since Ozone is already total column density

% Also compute WEAM, Water effective airmass, as WEAM=sqrt(X*PWV/5).Â

% compute effective airmass terms from zenith angle, pressure, ozone, and

% PWV

X=secd(ZenithAngle)*(1-0.0012*(secd(ZenithAngle)^2-1));

EAM=X*(ThisPressure/PressureDefault);

OEAM=X*(ThisOzone/OzoneDefault);

WEAM=sqrt(X*ThisPWV/PWVDefault);

ThisStableTrans=StableTrans.^EAM;

ThisOzoneTrans=Ozone.^OEAM;

ThisWaterTrans=Water.^WEAM;

ThisTrans=ThisStableTrans.*ThisOzoneTrans.*ThisWaterTrans;

%% read in transmission inforamtion

filename="/Users/cstubbs/Desktop/OldDesktop/projects/lsst/AuxTel/QuadBand/Quadband.dat";

[Quadnm,Quadband]=importQuad(filename);

filename="/Users/cstubbs/Desktop/OldDesktop/projects/lsst/AuxTel/QuadBand/TotalAuxtelThroughput.dat";

[AuxTelnm,AuxTelT]=importAuxTel(filename);

% note this has quadband in there already!

% spline

Quad=spline(Quadnm,Quadband,nm);

Quad(nm>650)=0;

AuxTel=spline(AuxTelnm,AuxTelT,nm);

%% make blackbody star

StarSpec=calculatePhotonSpectrum(StarTemp,nm);

%% compute extinctions at one airmass

% band edges are [335,405] [405,490] [490,534] [534,625]. Find indices in

% nm vector:

TOA=(StarSpec.*AuxTel);

OneAirmass=TOA.*ThisTrans;

[foo,index1]=min(abs(nm-335));

[foo,index2]=min(abs(nm-405));

[foo,index3]=min(abs(nm-490));

[foo,index4]=min(abs(nm-534));

[foo,index5]=min(abs(nm-625));

m1=-2.5*log10(sum(OneAirmass(index1:index2)));

m2=-2.5*log10(sum(OneAirmass(index2:index3)));

m3=-2.5*log10(sum(OneAirmass(index3:index4)));

m4=-2.5*log10(sum(OneAirmass(index4:index5)));

E1=-2.5*log10(sum(OneAirmass(index1:index2))/sum(TOA(index1:index2)));

E2=-2.5*log10(sum(OneAirmass(index2:index3))/sum(TOA(index2:index3)));

E3=-2.5*log10(sum(OneAirmass(index3:index4))/sum(TOA(index3:index4)));

E4=-2.5*log10(sum(OneAirmass(index4:index5))/sum(TOA(index4:index5)));

E14=E1-E4;

E24=E2-E4;

E34=E3-E4;

%%

% print out parameters, C14 color, extinctions, extinction differences, C34

% color

out=[StarTemp/1000 ThisOzone ThisPWV ThisPressure m1-m4 E1 E2 E3 E4 E14 E24 E34 m3-m4]

%%

figure(20)

plot(nm,Trans,'ko-')

hold on

plot(nm,ThisTrans,'r-','LineWidth',2)

plot(nm,Quad.*AuxTel,'b-','LineWidth',2)

legend('Default','This Transmission','Instrument')

hold off

grid on

shg

%%

function photon_spectrum = calculatePhotonSpectrum(temperature, wavelengths_nm)

% Constants

h = 6.62607004e-34; % Planck constant (JÂ·s)

c = 2.99792458e8; % Speed of light in vacuum (m/s)

k = 1.38064852e-23; % Boltzmann constant (J/K)

% Convert wavelengths from nm to m

wavelengths_m = wavelengths_nm * 1e-9; % nm to m

% Calculate the spectrum using Planck's law

spectrum = arrayfun(@(wavelength) planckLaw(wavelength, temperature, h, c, k), wavelengths_m);

spectrum=spectrum.*wavelengths_nm; % convert to photon spectrum, arb normalization

% Normalize the spectrum to its value at 510 nm

normalization_factor = planckLaw(510e-9, temperature, h, c, k);

photon_spectrum = spectrum / normalization_factor;

end

function spectral_radiance = planckLaw(wavelength, temperature, h, c, k) % energy per wavelength

% Planck's law to calculate spectral radiance

spectral_radiance = (2 * h * c^2) ./ (wavelength.^5 .* (exp((h * c) ./ (wavelength * k * temperature)) - 1));

end

%% function to set up throughput

function [Quadnm, Quad] = importQuad(filename)

%IMPORTFILE1 Import data from a text file

% [QUADNM, QUAD] = IMPORTFILE1(FILENAME) reads data from text file

% FILENAME for the default selection. Returns the data as column

% vectors.

%

% [QUADNM, QUAD] = IMPORTFILE1(FILE, DATALINES) reads data for the

% specified row interval(s) of text file FILENAME. Specify DATALINES as

% a positive scalar integer or a N-by-2 array of positive scalar

% integers for dis-contiguous row intervals.

%

% Example:

% filemane="/Users/cstubbs/Desktop/OldDesktop/projects/lsst/AuxTel/QuadBand/Quadband.dat";

%

% See also READTABLE.

%

% Auto-generated by MATLAB on 23-Nov-2023 13:06:39

%% Input handling

% If dataLines is not specified, define defaults

if nargin < 2

dataLines = [1, Inf];

end

%% Set up the Import Options and import the data

opts = delimitedTextImportOptions("NumVariables", 2);

% Specify range and delimiter

opts.DataLines = dataLines;

opts.Delimiter = "\t";

% Specify column names and types

opts.VariableNames = ["Quadnm", "Quad"];

opts.VariableTypes = ["double", "double"];

% Specify file level properties

opts.ExtraColumnsRule = "ignore";

opts.EmptyLineRule = "read";

% Import the data

tbl = readtable(filename, opts);

%% Convert to output type

Quadnm = tbl.Quadnm;

Quad = tbl.Quad;

end

function [AuxTelnm, AuxTel] = importAuxTel(filename)

%IMPORTFILE1 Import data from a text file

% [AUXTELNM, AUXTEL] = IMPORTFILE1(FILENAME) reads data from text file

% FILENAME for the default selection. Returns the data as column

% vectors.

%

% [AUXTELNM, AUXTEL] = IMPORTFILE1(FILE, DATALINES) reads data for the

% specified row interval(s) of text file FILENAME. Specify DATALINES as

% a positive scalar integer or a N-by-2 array of positive scalar

% integers for dis-contiguous row intervals.

%

% Example:

% filename="/Users/cstubbs/Desktop/OldDesktop/projects/lsst/AuxTel/QuadBand/TotalAuxtelThroughput.dat"

%

% See also READTABLE.

%

% Auto-generated by MATLAB on 23-Nov-2023 13:08:07

%% Input handling

% If dataLines is not specified, define defaults

if nargin < 2

dataLines = [1, Inf];

end

%% Set up the Import Options and import the data

opts = delimitedTextImportOptions("NumVariables", 2);

% Specify range and delimiter

opts.DataLines = dataLines;

opts.Delimiter = ",";

% Specify column names and types

opts.VariableNames = ["AuxTelnm", "AuxTel"];

opts.VariableTypes = ["double", "double"];

% Specify file level properties

opts.ExtraColumnsRule = "ignore";

opts.EmptyLineRule = "read";

% Import the data

tbl = readtable(filename, opts);

%% Convert to output type

AuxTelnm = tbl.AuxTelnm;

AuxTel = tbl.AuxTel;

end