
loss_tan_range = [.0005 0.015];

%Calculate range of alpha from paillou paper
% tholins_alpha = 2.*pi./(0.0217).*sqrt((1.336./2).*(sqrt(1+.00282.^2)-1));
% 2.*pi./(0.0217).*sqrt((1.336./2).*(sqrt(1+.01084.^2)-1));

eps2_range = [1.5, 4];  % Range for real_eps2
% alpha_range = [0.1, 2.5];  % Range for alpha from above calc1
alpha_range = 2.*pi./(0.0217).*sqrt((1.336./2).*(sqrt(1+loss_tan_range.^2)-1));


[h_grid, eps2_grid, alpha_grid, loss_tan_grid, p2] = calculate_plot( ...
    fits.simpleTwoLayer{1}, eps2_range, loss_tan_range);

h_grid(h_grid<0) = NaN; %get rid of everything less than a few wavelengths

average_h_grid = mean(h_grid,3);
% average_h_grid = h_grid(:,:,3);


%Convert "alpha" to "electrical skin depth" via formula from Alice Le Gall
%Mare

skin_depth_grid = (1./alpha_grid)./2;
%%
fig1 = figure;


axes1 = axes('Parent',fig1,'YScale','log','XScale','linear','Layer','top');

% set(gca, 'YTick', logspace(log10(alpha_range(1)), log10(alpha_range(2)),11))
%colorbar;
% yyaxis right
contours = [0.02, 0.05, 0.1, 0.2, 0.4, 0.8, 1.6];
[C,H] = contourf(eps2_grid,loss_tan_grid,average_h_grid, ...
    contours, ...
    "LabelFormat","%0.2f m", ...
    "FaceAlpha",0.35, ...
    "LineWidth", 2);
clabel([],H,contours, 'LabelSpacing', 1000, 'FontSize', 13, 'FontWeight', 'bold')
yscale("log")
ylabel("Surface Loss Tangent tan(\delta)");
xlabel("Basal Layer Dielectric Constant \epsilon_2");
% title("Two-Layer Phase Space")
colormap linspecer
xlim([2.7 4])
% ylim([0.1 2.5])

xline(3.1,'-.',{'Water Ice Subsurface'}, ...
    'FontSize', 13, "LineWidth", 1.5, ...
    'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'bottom');
yline(.00282,'-.',{'Powdered Tholins Surface'}, ...
    'FontSize', 13, "LineWidth", 1.5, ...
    'LabelHorizontalAlignment', 'Right', 'LabelVerticalAlignment', 'top');


% yyaxis right
% set(gca,'ydir','reverse')
% yscale('log')
% ylim([skin_depth_grid(1,end) skin_depth_grid(1,1)]);
% ytick = [0.2 0.5 1 2 3 4 5 6 7];
% yticks(ytick)
% ylabel("Surface Electromagnetic Skin Depth (m)")
% 
% ax = gca;
% ax.YAxis(1).Color = 'k';
% ax.YAxis(2).Color = 'k';

fixfig


%%

% disp()

e_grid = skin_depth_tb(fits.simpleTwoLayer{1}.eps1, eps2_grid, skin_depth_grid, h_grid(:,:,1));

figure;
contourf(eps2_grid,skin_depth_grid,e_grid,...
    "ShowText",true, ...
    "FaceAlpha",0.25, ...
    "LineWidth", 2);
set(gca,'YScale','log')
ylabel("Electrical Skin Depth of Top Layer (m)");
xlabel("Effective Dielectric of the Subsurface");
title("Emissivity Phase Space. Color = emissivity")
colormap jet
xlim([2.5 4])
% ylim([0.1 2.5])

% xline(3.1,'--',{'Solid Water Ice','Baselayer'}, ...
%     'FontSize', 15, "LineWidth", 2);
% yline(1/tholins_alpha/2,'--',{'Powdered Tholins', 'Toplayer'}, ...
%     'FontSize', 15, "LineWidth", 2);

fixfig

function [h_grid, eps2_grid, alpha_grid, loss_tan_grid, p2] = calculate_plot(simple_two_layer_fit,eps2_range, loss_tan_range)
% For an output of (so called) eps2 from simple_two_layer_fit, calculate the "power coeff"
% p2. From p2, create a phase diagram (matrix) which holds possible
% combinations of {eps2, h, alpha} which could produce that power
% coefficient. eps2 is the x axis, alpha is the y axis, and h is the color.

eps1 = simple_two_layer_fit.eps1;
fake_eps2 = simple_two_layer_fit.eps2;


x = deg2rad(40);


p2 =((1 - sqrt(fake_eps2))/(1 + sqrt(fake_eps2)))^2;
disp(p2)

% %FOR EXAMPLE. IN THE FUTURE, A WHOLE MATRIX OF REAL_EPS2, ALPHA
% real_eps2 = 2.9;
% alpha = 0.1;



% Define the range and resolution for real_eps2 and alpha

eps2_resolution = 500;  % Number of points for real_eps2
loss_tan_resolution = 1000;  % Number of points for alpha

real_eps2_axis = linspace(eps2_range(1), eps2_range(2), eps2_resolution);
loss_tan_axis = logspace(log10(loss_tan_range(1)), log10(loss_tan_range(2)), loss_tan_resolution);


x_resolution = 5;
x_axis = linspace(deg2rad(10), deg2rad(40), x_resolution);

h_grid = zeros([eps2_resolution loss_tan_resolution x_resolution]);
[eps2_grid,loss_tan_grid] = ndgrid(real_eps2_axis, loss_tan_axis);

alpha_grid = 2.*pi./(0.0217).*sqrt((1.336./2).*(sqrt(1+loss_tan_grid.^2)-1));


for k = 1:x_resolution
for i = 1:eps2_resolution
    for j = 1:loss_tan_resolution

        real_eps2 = real_eps2_axis(i);
        alpha = alpha_grid(1,j);

        real_ref = ((1 - sqrt(real_eps2))/(1 + sqrt(real_eps2)))^2;

        % Define the equation to solve
        equation_to_solve = @(h) coeff(eps1, real_eps2, alpha, h, x_axis(k)) .* real_ref - p2;
        
        % Initial guess for h
        initial_guess = 1; % Your initial guess for h
        
        % Solve for h
        h_solution = fzero(equation_to_solve, initial_guess);
        h_grid(i,j,k) = h_solution;
    end
    if mod(i, fix(eps2_resolution/5)) == 0
        disp(i)
    end
end
end


end


function AttRe = coeff(eps1, eps2, alpha, h, x)
%Att = attenuation coeff through top layer
%reflectivity of sub relative to its free space value
%T2 = Two-way Power Transmission Coefficient, VV polarizatoin


n1 = sqrt(eps1); %refractive index of surf
n2 = sqrt(eps2); %refractive index of sub

Att = exp(-4.*alpha.*h./cos(x)); %attenuation coeff through top layer

Re = (((n1 - n2)/(n1 + n2))./((1 - n2)/(1 + n2)))^2; %reflectivity of sub relative to its free space value

AttRe = Att.*Re;

end

function e = skin_depth_tb(eps1, eps2, skin_depth, h)
%returns a theoretical brightness temperature based on Le Gall Composition,
% seasonal change, and bathymetry of Ligeia Mare, Titan, derived from its
%microwave thermal emission
%applied to the plains

%inc and pol are angles in rad
%lat in deg
%t in cassini standard ephem time, seconds since J2000


G1 = fresn(1, eps1, 0, 0); %Fresnel air to medium
G2 = fresn(eps1, eps2, 0, 0); %Fresnel medium to baselayer

inc_prime = asin(sin(0)./sqrt(eps1));

% skin_depth = (0.0217./4.*pi)./sqrt((eps1./2).*(sqrt(1+loss_tan.^2)-1));
L1 = exp(-h./(skin_depth.*cos(inc_prime)));

e =  ((1-G1)./(1-G1.*G2.*L1.^2)).*(((1+G2.*L1).*(1-L1))+((1-G2).*L1));


end



