addpath(genpath('W:\Cassini\Titan\valerio\Anthony_backscatter'));
addpath(genpath('W:\Anthony\Plains Backscatter'));
%% SITE 1 DEFINITION

%Load files and coregister
sar_filenames = ["BISQH69N225_D210_T064S01_V02.mat"; "BISQH09N141_D166_T043S01_V02.mat"; "BISQH61N226_D243_T086S01_V02.mat"];
% inc_filenames = ["BIEQH69N225_D210_T064S01_V02.mat"; "BIEQH09N141_D166_T043S01_V02.mat"; "BIEQH61N226_D243_T086S01_V02.mat"];

%Add altimetry
% csv = readmatrix("ABDR_SUMMARY_04_D131_V03.csv");

% Find the indices of the rows of the altimetry bursts we want
% ind = (csv(:,1) >= 131033440 & csv(:,1) <= 131033459)...
%     | (csv(:,1) >= 131033532 & csv(:,1) <= 131033551);

%new altimetry because the old one was bad
% csv = readmatrix("ABDR_SUMMARY_04_D127_V03.csv");
% ind = (csv(:,1) >= 127017856 & csv(:,1) <= 127017869);

%that one was also bad

% csv = readmatrix("ABDR_SUMMARY_04_D285_P1_V02.csv");
% ind = (csv(:,1) >= 285038053 & csv(:,1) <= 285038112);

load T126_burstID_sigma0.mat
ind = (burst_id >= 285038053 & burst_id <= 285038112);


site = 'Site1';


%% SITE 2 DEFINITION NOT USING. ALTIMETER SHOWS IT IS A HILL AND ALSO BETWEEN LAKES
%AND WHEN ALTIMETER IS CALIBRATED, -5 to 0 DB BACKSCATTER. STRANGE!
%Load files and coregister
sar_filenames = ["BISQH33N028_D120_T025S01_V02.mat", "BISQH65N015_D127_T029S01_V02.mat", "BISQH69N225_D210_T064S01_V02.mat", "BISQG79N344_D250_T092S02_V02.mat"];
% inc_filenames = ["BIEQH33N028_D120_T025S01_V02.mat", "BIEQH65N015_D127_T029S01_V02.mat", "BIEQH69N225_D210_T064S01_V02.mat", "BIEQG79N344_D250_T092S02_V02.mat"];
%Add altimetry
csv = readmatrix("ABDR_SUMMARY_04_D265_P1_V02.csv");

% Find the indices of the rows of the altimetry bursts we want
ind = (csv(:,1) >= 265025481 & csv(:,1) <= 265025511);

site = 'Site2';

%% SITE 3 DEFINITION

sar_filenames = ["BISQH72S104_D149_T036S03_V02.mat","BISQH51S121_D177_T049S01_V02.mat", "BISQG70S033_D181_T050S02_V02.mat"];
% inc_filenames = ["BIEQH72S104_D149_T036S03_V02.mat","BIEQH51S121_D177_T049S01_V02.mat", "BIEQG70S033_D181_T050S02_V02.mat"];

csv = readmatrix("ABDR_SUMMARY_04_D181_V03.csv");

% Find the indices of the rows of the altimetry bursts we want
ind = (csv(:,1) >= 181035022 & csv(:,1) <= 181035028)...
    | (csv(:,1) >= 181035048 & csv(:,1) <= 181035053);

site = 'Site3';

%% SITE ? DEFINITION
sar_filenames = ["BISQH52N010_D126_T028S01_V02.mat", "T95NoCorr.mat", "BISQH64N243_D285_T126S01_V02.mat"];
% inc_filenames = ["BIEQH52N010_D126_T028S01_V02.mat", "T95InA.mat", "BIEQH64N243_D285_T126S01_V02.mat"];

csv = readmatrix("ABDR_SUMMARY_04_D248_P1_V03.csv");

% Find the indices of the rows of the altimetry bursts we want
 ind = (csv(:,1) >= 248037248 & csv(:,1) <= 248037261);
   

site = 'Site?';

%% SITE 4 DEFINITION TODO NOT USING. AROUND ONTARIO LACUS
%BUT I WANT TO DO IT ANYWAYS, THERE ARE 1 MILLION TRACKS. ;(

sar_filenames = ["BISQH62S167_D200_T058S01_V02.mat","BISQH78S004_D211_T065S01_V02.mat","BISQI46S161_D199_T057S01_V02.mat"];
% inc_filenames = [""];

csv = readmatrix("ABDR_SUMMARY_04_D181_V03.csv");

% Find the indices of the rows of the altimetry bursts we want
ind = (csv(:,1) >= 181035022 & csv(:,1) <= 181035028)...
    | (csv(:,1) >= 181035048 & csv(:,1) <= 181035053);

site = 'Site4';
%% SITE 5 DEFINITION 

sar_filenames = ["BISQH80N051_D087_T016S01_V02.mat", "BISQH69N342_D100_T019S01_V02.mat", "BISQG40N156_D276_T120S05_V02.mat"];
% inc_filenames = ["BIEQH80N051_D087_T016S01_V02.mat", "BIEQH69N342_D100_T019S01_V02.mat", "BIEQH43N144_D100_T019S02_V02.mat", "BIEQG40N156_D276_T120S05_V02.mat"];

csv = readmatrix("ABDR_SUMMARY_04_D195_V03.csv");

% Find the indices of the rows of the altimetry bursts we want
ind = csv(:,1) >= 31212 & csv(:,1) <= 31296;

site = 'Site5';

%% SITE 6 DEFINITION 
%NOTE: No Altimetry. Also, might not be plains?
sar_filenames = ["BISQH36S166_D195_T056S01_V02.mat","BISQH35S109_D161_T041S01_V02.mat"];
csv = [];

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site6';

%% SITE 7 DEFINITION 
%Alt: T30 
%Note: The T30 echoes on PDS node are analyzed pulse-limited
%but I am using Valerio's data which is alalyze beam-limited because T30 is
%so low altitude. Using the (linear units) power from the new EchoesT30...
%file with ind = 1972-2007, the datapoint is 10log10(mean(sigma_0(1972:2007)))

%Note: T108 has very low incidence angle. Focus on that one, it is the most
%important. 
sar_filenames = ["BISQH26N009_D111_T023S01_V02.mat", "BISQH36N017_D265_T108S01_V02.mat", "BISQH49N282_D240_T084S01_V02.mat"];

% csv = readmatrix("ABDR_SUMMARY_04_D131_V03.csv"); 
% % Find the indices of the rows of the altimetry bursts we want

ind = csv(:,1) >= 131027153 & csv(:,1) <= 131027188;


site = 'Site7';

%% Site 8 Definition
%Alt: no
sar_filenames = ["BISQH09N141_D166_T043S01_V02.mat", "BISQH30N213_D166_T043S03_V02.mat", "BISQH69N225_D210_T064S01_V02.mat"];
csv = [];

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site8';

%% Site 9 Definition
%Alt: no
sar_filenames = ["BISQI36N017_D265_T108S01_V02.mat"];
%Site 9a
csv = []

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site9';
%% Site 9a Definition
%Alt: no
sar_filenames = ["BISQI36N017_D265_T108S01_V02.mat"];
csv = [];

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site9a';

%% Site 10 Definition
%Alt: no
sar_filenames = ["BISQH03N158_D167_T044S01_V02.mat", "BISQH01N196_D149_T036S04_V02.mat","BISQH09N141_D166_T043S01_V02.mat","BISQH75N177_D239_T083S01_V02.mat"];
csv = [];

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site10';
 
%% Site 11 Definition NOT USING. SOUTHERN PALEO-BASIN? OUTLIER, MUCH DARKER THAN PLAINS 
% %Alt: no
% sar_filenames = ["BISQH62S167_D200_T058S01_V02.mat","BISQH69S314_D220_T071S01_V02.mat","BISQH46S161_D199_T057S01_V02.mat"];
% csv = [];
% 
% % Find the indices of the rows of the altimetry bursts we want
% ind = [];
% 
% site = 'Site11';
% 
%% Site 12 Definition
%Alt: no
sar_filenames = ["BISQH61N250_D131_T030S01_V02.mat", "BISQH05N175_D174_T048S03_V02.mat"];
csv = [];

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site12';

%% Site 13 Definition 
%Alt: no
sar_filenames = ["BISQH24S001_D261_T104S03_V02.mat"];
csv = [];

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site13';

%% Site 14 Definition NOT USING - MORPHOLOGICALLY TOO COMPLEX, TOO BRIGHT OUTLIER
% %Alt: no
% sar_filenames = ["BISQH68N341_D243_T086S03_V03.mat", "BISQH66N015_D098_T018S01_V02.mat"];
% csv = [];
% 
% % Find the indices of the rows of the altimetry bursts we want
% ind = [];
% 
% site = 'Site14';

%% Site 15 Definition
%Alt: two swaths. I want to vertcat them after doing the indices so I can
%use both swaths for this alt datapoint
sar_filenames = ["BISQH88S198_D157_T039S01_V02.mat", "BISQH72S104_D149_T036S03_V02.mat"];
csv = readmatrix("ABDR_SUMMARY_04_D157_V03.csv");
% Find the indices of the rows of the altimetry bursts we want
ind = (csv(:,1) >=157026803 & csv(:,1) <= 157026812);

site = 'Site15';

%% Site 15a Definition
%TODO: decide if this site is worth using - it is very RADAR bright. Can be
%seen even by eye on the arcmap. At all inc angles, so high dialectric? Hm.
%Not sure it's fair to call it an outlier, but it widens the spread of
%points a whole lot. 
%Using T007 BIE *V01*. Is that ok? Should be.
sar_filenames = ["BISQH51S012_D059_T007S01_V03.mat","BISQH30S030_D086_T015S02_V02.mat","BISQH29S032_D149_T036S01_V02.mat","BISQH29S025_D149_T036S02_V02.mat", ...
    "BISQH88S198_D157_T039S01_V02.mat", "BISQH36S029_D157_T039S02_V02.mat", "BISQH38S011_D157_T039S07_V02.mat"];



csv = readmatrix("ABDR_SUMMARY_04_D261_P1_V03.csv");
ind = (csv(:,1) >=261034452  & csv(:,1) <= 261034505);

site = 'Site15a';

%% Site 16 Definition 
%Alt: 
sar_filenames = ["BISQH69N342_D100_T019S01_V02.mat", "BISQH22N352_D261_T104S01_V02.mat", "BISQG54N343_D240_T084S02_V02.mat"];
csv = [];

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site16';

%% Site 17 Definition 
%Alt: 
sar_filenames = ["BISQH62S167_D200_T058S01_V02.mat", "BISQH17S147_D200_T058S02_V03.mat", "BISQG18S132_D229_T077S02_V02.mat"];
csv = [];

% Find the indices of the rows of the altimetry bursts we want
ind = [];

site = 'Site17';

%% Site 18 Definition NOT USING, PROBABLY NOT PLAINS
%Alt: 
sar_filenames = ["BISQH61N250_D131_T030S01_V02.mat", "T95S2"];

load  %Valerio reproc 291
ind = (burst_id >= 248037488  & burst_id <= 248037505); % Find the indices of the rows of the altimetry bursts we want

site = 'Site18';

%% Site 19 Definition TODO 
%Alt: 
sar_filenames = ["BISQH71N100_D248_T091S04_V02.mat", "BISQH56N094_D111_T023S02_V02.mat"];
csv = readmatrix("ABDR_SUMMARY_04_D111_V03.csv");

% Find the indices of the rows of the altimetry bursts we want
ind = (csv(:,1) >=111018640  & csv(:,1) <= 111018705);

site = 'Site19';

%% Now, analyze the selected site
[co_sar, co_inc] = Coreg(sar_filenames);

%%
%Mask out plains region.
%mask_sar, mask_inc probably unused but might be useful for visualization?
[mask_sar, mask_inc, mask] = Mask(co_sar, co_inc);

close
Plotsite(mask_sar,mask_inc);

%%
figure
hold on

% Define colors
C = distinguishable_colors(numel(sar_filenames));

% Create histograms with percentile indicators
for i = 1:size(mask_inc, 3)
    % Convert data to dB
    data_db = 10*log10(mask_sar(:,:,i));
    
    % Calculate percentiles
    median_val = median(data_db(:));
    percentile_25 = prctile(data_db(:), 25);
    percentile_75 = prctile(data_db(:), 75);
    
    % Plot histogram
    h = histogram(data_db, 100, 'FaceColor', C(i,:), 'FaceAlpha', 0.3, 'EdgeColor', 'none');
    
    % Get the maximum height of the histogram for this dataset
    max_height = max(h.Values);
    
    % Plot vertical lines for percentiles
    line([median_val median_val], [0 max_height], 'Color', C(i,:), 'LineWidth', 2, 'LineStyle', '-');
    line([percentile_25 percentile_25], [0 max_height], 'Color', C(i,:), 'LineWidth', 1.5, 'LineStyle', '--');
    line([percentile_75 percentile_75], [0 max_height], 'Color', C(i,:), 'LineWidth', 1.5, 'LineStyle', '--');
    
    % Add text labels for percentiles
    text(median_val, max_height*1.05, 'Median', 'Color', C(i,:), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 8);
    text(percentile_25, max_height*1.02, '25th', 'Color', C(i,:), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 8);
    text(percentile_75, max_height*1.02, '75th', 'Color', C(i,:), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 8);
end

% Customize the plot
xlabel('Backscatter (dB)');
ylabel('Frequency');
title('Histogram of Backscatter with Percentile Indicators');
legend(cellstr(num2str((1:size(mask_inc, 3))')), 'Location', 'best');
grid on

% Adjust axis limits to accommodate labels
ylim_current = ylim;
ylim([ylim_current(1), ylim_current(2)*1.1]);  % Extend the upper y-limit by 10%

hold off
%% Extract sigma

for ii=1:size(co_sar,3)
    tempsar = co_sar(:,:,ii);
    sigma0(:,ii) = tempsar(mask);
    
    tempinc = co_inc(:,:,ii);
    inc(:,ii) = tempinc(mask);
end

%% EXTRACT ALTIMETRY MEASUREMENTS FROM CSV
%USE "IND" DEFINED IN SITE DEFINITIONS

%COL 1:SAB COUNTER
%COL 5:LON
%COL 6:LAT
%COL 14:INCIDENCE ANGLE
%COL 15: SIGMA0

%Compile data into columns for fitting, including and not including
%altimetry
%NOTE: from here, all angles are in radians'

if ~isempty(csv)

    % Extract the values in the 14th and 15th columns for those rows
    alt_sigma0 = db2pow(csv(ind,15));
    alt_inc = csv(ind,14);

    data.alt.sigma0 = [reshape(sigma0,[],1); alt_sigma0];
    data.alt.inc = deg2rad([reshape(inc,[], 1); alt_inc]);
else
    
    alt_sigma0 = [];
    alt_inc = [];

    data.alt.sigma0 = reshape(sigma0,[],1);
    data.alt.inc = deg2rad(reshape(inc,[], 1));

end

data.noalt.sigma0 = reshape(sigma0,[],1);
data.noalt.inc = deg2rad(reshape(inc,[], 1));

%% Fit the individual site data

fits = fitBackscatter(data, 1);

% %% 
% fitsfieldnames = fieldnames(fits);
% coeffs = zeros(3,1);
% for kk = 1:numel(fitsfieldnames)
%         cfitarray = fits.(fitsfieldnames{kk});
%         cfitarray = cfitarray(1);
%         cfitarray = cfitarray{1};
%         coeffs = coeffs + coeffvalues(cfitarray)';
% end
% coeffs = coeffs ./3;
% 
% 
% fitsfieldnames = fieldnames(fits);
% nocoeffs = zeros(3,1);
% for kk = 1:numel(fitsfieldnames)
%         cfitarray = fits.(fitsfieldnames{kk});
%         cfitarray = cfitarray(2);
%         cfitarray = cfitarray{1};
%         nocoeffs = nocoeffs + coeffvalues(cfitarray)';
% end
% nocoeffs = nocoeffs ./3;
% 
% fits.avg.withalt = coeffs;
% fits.avg.withoutalt = nocoeffs;
% 
% 
% %% SAVE FITS DATA CORRECTLY
% switch site
%     case 'Site1'
%         Site1fits = fits;
%     case 'Site2'
%         Site2fits = fits;
%     case 'Site3'
%         Site3fits = fits;
%     case 'Site5'
%         Site5fits = fits;
% end
