% MATLAB code comet_extinction.m to calculate the extinction corrected magnitude of a comet based on an observed comparison star
% Calculation is made according to the ICQ method from the July 1992 issue of *International Comet Quarterly*, Vol. 14, pages 55-59.  Copyright 1992.
% MAGNITUDE CORRECTIONS FOR ATMOSPHERIC EXTINCTION by Daniel W. E. Green
% http://www.icq.eps.harvard.edu/ICQExtinct.html
%
% Johan Warell, SAAF Solar System Section, Oct 20, 2024
% johan.warell at saaf.se, www.saaf.se

% Define extinction tables (Ia, Ib, Ic) for average, winter, and summer conditions
table_ia = [... % Zenith Distance (degrees), h=0, h=0.5, h=1, h=2, h=3
    1, 0.28, 0.24, 0.21, 0.16, 0.13;
    10, 0.29, 0.24, 0.21, 0.16, 0.13;
    20, 0.30, 0.25, 0.22, 0.17, 0.14;
    30, 0.32, 0.28, 0.24, 0.19, 0.15;
    40, 0.37, 0.31, 0.27, 0.21, 0.17;
    45, 0.40, 0.34, 0.29, 0.23, 0.19;
    50, 0.44, 0.37, 0.32, 0.25, 0.21;
    55, 0.49, 0.42, 0.36, 0.28, 0.23;
    60, 0.56, 0.48, 0.41, 0.32, 0.26;
    62, 0.60, 0.51, 0.44, 0.34, 0.28;
    64, 0.64, 0.54, 0.47, 0.37, 0.30;
    66, 0.69, 0.59, 0.51, 0.39, 0.32;
    68, 0.75, 0.64, 0.55, 0.43, 0.35;
    70, 0.82, 0.70, 0.60, 0.47, 0.39;
    71, 0.86, 0.73, 0.63, 0.49, 0.40;
    72, 0.91, 0.77, 0.66, 0.52, 0.43;
    73, 0.96, 0.81, 0.70, 0.55, 0.45;
    74, 1.02, 0.86, 0.74, 0.58, 0.48;
    75, 1.08, 0.92, 0.79, 0.62, 0.51;
    76, 1.15, 0.98, 0.84, 0.66, 0.54;
    77, 1.24, 1.05, 0.91, 0.71, 0.58;
    78, 1.34, 1.13, 0.98, 0.76, 0.63;
    79, 1.45, 1.23, 1.06, 0.83, 0.68;
    80, 1.59, 1.34, 1.16, 0.91, 0.74;
    81, 1.75, 1.48, 1.28, 1.00, 0.82;
    82, 1.94, 1.65, 1.42, 1.11, 0.91;
    83, 2.19, 1.86, 1.60, 1.25, 1.03;
    84, 2.50, 2.12, 1.83, 1.43, 1.17;
    85, 2.91, 2.46, 2.13, 1.66, 1.36;
    86, 3.45, 2.93, 2.53, 1.97, 1.62;
    87, 4.23, 3.59, 3.10, 2.42, 1.99;
    88, 5.41, 4.59, 3.96, 3.09, 2.54;
    89, 7.38, 6.26, 5.40, 4.22, 3.46;
    90, 11.24, 9.53, 8.23, 6.42, 5.28];

table_ib = [... % Zenith Distance (degrees), h=0, h=0.5, h=1, h=2, h=3 (Winter)
    1, 0.25, 0.21, 0.19, 0.15, 0.13;
    10, 0.25, 0.22, 0.19, 0.15, 0.13;
    20, 0.26, 0.23, 0.20, 0.16, 0.14;
    30, 0.28, 0.25, 0.22, 0.17, 0.15;
    40, 0.32, 0.28, 0.24, 0.20, 0.17;
    45, 0.35, 0.30, 0.26, 0.21, 0.18;
    50, 0.38, 0.33, 0.29, 0.24, 0.20;
    55, 0.43, 0.37, 0.33, 0.26, 0.22;
    60, 0.49, 0.42, 0.37, 0.30, 0.25;
    62, 0.52, 0.45, 0.40, 0.32, 0.27;
    64, 0.56, 0.48, 0.43, 0.34, 0.29;
    66, 0.60, 0.52, 0.46, 0.37, 0.31;
    68, 0.65, 0.57, 0.50, 0.40, 0.34;
    70, 0.72, 0.62, 0.55, 0.44, 0.37;
    71, 0.75, 0.65, 0.57, 0.46, 0.39;
    72, 0.79, 0.69, 0.60, 0.49, 0.41;
    73, 0.84, 0.72, 0.64, 0.52, 0.43;
    74, 0.89, 0.77, 0.68, 0.55, 0.46;
    75, 0.94, 0.82, 0.72, 0.58, 0.49;
    76, 1.01, 0.87, 0.77, 0.62, 0.52;
    77, 1.08, 0.94, 0.82, 0.67, 0.56;
    78, 1.16, 1.01, 0.89, 0.72, 0.60;
    79, 1.26, 1.10, 0.97, 0.78, 0.66;
    80, 1.38, 1.20, 1.06, 0.85, 0.72;
    81, 1.52, 1.32, 1.16, 0.94, 0.79;
    82, 1.70, 1.47, 1.29, 1.05, 0.88;
    83, 1.91, 1.65, 1.46, 1.18, 0.99;
    84, 2.18, 1.89, 1.66, 1.34, 1.13;
    85, 2.53, 2.20, 1.93, 1.56, 1.31;
    86, 3.01, 2.61, 2.30, 1.86, 1.56;
    87, 3.69, 3.20, 2.82, 2.28, 1.91;
    88, 4.72, 4.09, 3.60, 2.91, 2.45;
    89, 6.44, 5.58, 4.91, 3.97, 3.34;
    90, 9.80, 8.50, 7.49, 6.05, 5.08];

table_ic = [... % Zenith Distance (degrees), h=0, h=0.5, h=1, h=2, h=3 (Summer)
    1, 0.32, 0.26, 0.22, 0.17, 0.14;
    10, 0.32, 0.27, 0.23, 0.17, 0.14;
    20, 0.34, 0.28, 0.24, 0.18, 0.15;
    30, 0.37, 0.30, 0.26, 0.20, 0.16;
    40, 0.41, 0.34, 0.29, 0.22, 0.18;
    45, 0.45, 0.37, 0.32, 0.24, 0.19;
    50, 0.49, 0.41, 0.35, 0.26, 0.21;
    55, 0.55, 0.46, 0.39, 0.30, 0.24;
    60, 0.63, 0.53, 0.45, 0.34, 0.27;
    62, 0.68, 0.56, 0.48, 0.36, 0.29;
    64, 0.72, 0.60, 0.51, 0.39, 0.31;
    66, 0.78, 0.65, 0.55, 0.42, 0.34;
    68, 0.85, 0.70, 0.60, 0.45, 0.36;
    70, 0.93, 0.77, 0.65, 0.50, 0.40;
    71, 0.97, 0.81, 0.69, 0.52, 0.42;
    72, 1.02, 0.85, 0.72, 0.55, 0.44;
    73, 1.08, 0.90, 0.76, 0.58, 0.47;
    74, 1.15, 0.95, 0.81, 0.61, 0.49;
    75, 1.22, 1.01, 0.86, 0.65, 0.53;
    76, 1.30, 1.08, 0.92, 0.70, 0.56;
    77, 1.40, 1.16, 0.99, 0.75, 0.60;
    78, 1.51, 1.25, 1.07, 0.81, 0.65;
    79, 1.64, 1.36, 1.16, 0.88, 0.71;
    80, 1.79, 1.49, 1.26, 0.96, 0.77;
    81, 1.97, 1.64, 1.39, 1.06, 0.85;
    82, 2.19, 1.83, 1.55, 1.18, 0.95;
    83, 2.47, 2.06, 1.75, 1.32, 1.07;
    84, 2.82, 2.35, 1.99, 1.51, 1.22;
    85, 3.28, 2.73, 2.32, 1.76, 1.41;
    86, 3.90, 3.25, 2.75, 2.09, 1.68;
    87, 4.78, 3.98, 3.38, 2.56, 2.06;
    88, 6.11, 5.09, 4.32, 3.28, 2.63;
    89, 8.33, 6.93, 5.89, 4.47, 3.59;
    90, 12.68, 10.56, 8.97, 6.80, 5.47];

% Inputs from the user
delta_m = input('Observed magnitude difference between comet and comparison star (negative if comet is brighter): ');
m_star_catalog = input('Catalog magnitude of comparison star: ');
altitude_comet = input('Altitude of comet (in degrees): ');
altitude_star = input('Altitude of comparison star (in degrees): ');
observer_height_m = input('Observer height above sea level (meters): ');
extinction_type = input('Extinction coefficient type (1 for summer, 2 for average, 3 for winter): ');

% Convert observer height from meters to kilometers
observer_height_km = observer_height_m / 1000;

% Calculate zenith distance from altitude
zenith_comet = 90 - altitude_comet;
zenith_star = 90 - altitude_star;

% Select appropriate extinction table based on input type
if extinction_type == 1
    extinction_table = table_ic; % Summer
    fprintf('Selected summer extinction table\n');
elseif extinction_type == 2
    extinction_table = table_ia; % Average
    fprintf('Selected average extinction table\n');
elseif extinction_type == 3
    extinction_table = table_ib; % Winter
    fprintf('Selected winter extinction table\n');
else
    error('Invalid extinction type! Choose 1 for summer, 2 for average, or 3 for winter.');
end
fprintf('*****************************************\n');


% Set up the grid for interpolation (zenith distance and height)
zenith_distances = extinction_table(:, 1); % Zenith distances from the first column
heights = [0, 0.5, 1, 2, 3]; % Heights in km from the table headers

% Interpolate total extinction for comet and star using zenith distance and height
extinction_comet = interp1(zenith_distances, extinction_table(:, 2:end), zenith_comet, 'linear');
extinction_comet = interp1(heights, extinction_comet, observer_height_km, 'linear');

extinction_star = interp1(zenith_distances, extinction_table(:, 2:end), zenith_star, 'linear');
extinction_star = interp1(heights, extinction_star, observer_height_km, 'linear');

% Display interpolated extinction values for comet and star
fprintf('Interpolated extinction for comet: %.3f magnitudes\n', extinction_comet);
fprintf('Interpolated extinction for star: %.3f magnitudes\n', extinction_star);

% Calculate the corrected magnitude for the comparison star
m_star_corrected = m_star_catalog + extinction_star;
fprintf('Extinction corrected magnitude of comparison star: %.3f\n', m_star_corrected);

% Calculate the extinction corrected magnitude of the comet
m_comet_corrected = m_star_corrected + delta_m - extinction_comet;
fprintf('Extinction corrected magnitude of comet: %.3f\n', m_comet_corrected);

% Determine the ICQ note code based on the altitude of the comet and star
if altitude_comet < 10 || altitude_star < 10
    ICQ_code = '$';
    fprintf('ICQ Note Code: %s (Extinction corrections applied at altitudes <10 degrees)\n', ICQ_code);
elseif extinction_type == 1
    ICQ_code = 's';
    fprintf('ICQ Note Code: %s (Summer extinction correction applied)\n', ICQ_code);
elseif extinction_type == 2
    ICQ_code = 'a';
    fprintf('ICQ Note Code: %s (Average extinction correction applied)\n', ICQ_code);
elseif extinction_type == 3
    ICQ_code = 'w';
    fprintf('ICQ Note Code: %s (Winter extinction correction applied)\n', ICQ_code);
else
    fprintf('ICQ Note Code: None\n');
end