function [R3var,R4var,R5var,PssPs3var,PssPs4var,PssPs5var,epsvar,stratvar] = strat_vary_analysis(eps_profile, strat, numruns)
%Reads in a single permittivity profile and strat profile, then varies them
%numruns times and calculates the PssPs for each variant at each frequency. 

%eps_ice             = 3.15 + i*0.0006;
eps_ice             = matzler(4e6,175);
eps_base            = 8.8 - 1i*((2e-2)*8.8);
eps_spld            = maxwell_garnett(0.10,eps_ice,eps_base);
%w                   = 2 * pi * (4e6);
%c                   = 3e8;
[rs4, ts]           = fresnel(1, real(eps_spld), 1, imag(eps_spld), 4e6, 0, 'normal');
[rs3, ts]           = fresnel(1, real(eps_spld), 1, imag(eps_spld), 3e6, 0, 'normal');
[rs5, ts]           = fresnel(1, real(eps_spld), 1, imag(eps_spld), 5e6, 0, 'normal');
Rs4                 = abs(rs4)^2;
Ts4                 = 1-Rs4;
Rs3                 = abs(rs3)^2;
Ts3                 = 1-Rs3;
Rs5                 = abs(rs5)^2;
Ts5                 = 1-Rs5;

epsvar(:,1) = eps_profile;
stratvar(:,1) = strat;
[R3var(1),R4var(1),R5var(1)] = marsis_reflectivity(stratvar(:,1),epsvar(:,1));
parfor ii=2:numruns+1
    [epsvar(:,ii), stratvar(:,ii)] = vary_layers(eps_profile,strat);
    [R3var(ii),R4var(ii),R5var(ii)] = marsis_reflectivity(stratvar(:,ii),epsvar(:,ii));
%     [rho4(ii), R4(ii),T4(ii),A4(ii)] = multi_layer_rta(ones(size(epsvar,1)),[eps_ice; conj(epsvar(:,ii)); eps_base],[1; ones(size(epsvar,1),1); 1],4e6,0,0);
%     [rho3(ii), R3(ii),T3(ii),A3(ii)] = multi_layer_rta(ones(size(epsvar,1)),[eps_ice; conj(epsvar(:,ii)); eps_base],[1; ones(size(epsvar,1),1); 1],3e6,0,0);
%     [rho5(ii), R5(ii),T5(ii),A5(ii)] = multi_layer_rta(ones(size(epsvar,1)),[eps_ice; conj(epsvar(:,ii)); eps_base],[1; ones(size(epsvar,1),1); 1],5e6,0,0);
end

attenuation4 = exp(-(2 * pi * (4e6)) * (17e-6) * abs(imag(eps_spld)/real(eps_spld)) ); %from Grima 2021
attenuation3 = exp(-(2 * pi * (3e6)) * (17e-6) * abs(imag(eps_spld)/real(eps_spld)) );
attenuation5 = exp(-(2 * pi * (5e6)) * (17e-6) * abs(imag(eps_spld)/real(eps_spld)) );

R4lin = 10.^(R4var./10);
PssPs4var = (R4lin./Rs4) .* (Ts4.^2) .* attenuation4; %from Grima 2021
R3lin = 10.^(R3var./10);
PssPs3var = (R3lin./Rs3) .* (Ts3.^2) .* attenuation3;
R5lin = 10.^(R5var./10);
PssPs5var = (R5lin./Rs5) .* (Ts5.^2) .* attenuation5;

% median(10.*log10(PssPs3))
% median(10.*log10(PssPs4))
% median(10.*log10(PssPs5))
% 
% iqr(10.*log10(PssPs3))
% iqr(10.*log10(PssPs4))
% iqr(10.*log10(PssPs5))

% subplot(1,3,1)
% histogram(PssPs3,20)
% histogram(10.*log10(PssPs3),20)
% subplot(1,3,2)
% histogram(10.*log10(PssPs4),20)
% subplot(1,3,3)
% histogram(10.*log10(PssPs5),20)


end