clear all

% Fit measured MHO ACF's to Brownian theory
date = 'Apr_2019';
%date = 'Nov_2013';
pos = '3';

pp =2; % Plot index

ion_comp = 1;   % 1 for MH Chapman profile, 2 for IRI value

plot_height = [100, 300];
plot_var =2; % 1 for range, 2 for altitude, 3 for aspect angle
plot_angles = [0, 1, 2, 3, 4, 5];



filename = [date, '/acf_pos', pos, '_', date, '.h5'];
infoname = [date, '/looker', pos, '_', date, '.csv'];
IRIname = [date, '/iri', pos, '_', date, '.txt'];
MSISname = [date, '/msis_', date, '.txt'];

frad = 440*10^6;
B = 5*10^-5;

%% Read in Data
if strcmp(date, 'Feb_2019')
    % Feb. 22, 2013, Az -20 or -40, El < 10, 6 km range resolution
    acf = double(h5read(filename, '/acf')');
    lags = double(h5read(filename, '/lag_us')); % In microseconds
    range_h5 = double(h5read(filename, '/range_km')); % In km
elseif strcmp(date, 'Apr_2019')
    % Feb. 22, 2013, Az -20 or -40, El < 10, 6 km range resolution
    acf = double(h5read(filename, '/acf')');
    lags = double(h5read(filename, '/lag_us')); % In microseconds
    range_h5 = double(h5read(filename, '/range_km')); % In km
elseif strcmp(date, 'Nov_2013')
    % Nov. 11, 2013, Az -12, El 4.4, 15 km range resolution
    acf = double(h5read(filename, '/acfreal')');
    lags = double(h5read(filename, '/lags')); % In microseconds
    range_h5 = double(h5read(filename, '/ranges')); % In km
end

fv = (1:(length(lags))/2)/length(lags)/(lags(2)*10^-6);
fv = [-fv(end:-1:1), 0, fv];

for i = 1:length(range_h5)
    spec(i,:) = abs(fftshift(fft(acf(i,:), length(acf(i,:)))));
    spec_norm(i,:) = spec(i,:)/max(spec(i,:));
    spec_r2(i,:) = spec(i,:)*range_h5(i)^2;
    acf_norm(i, :) = acf(i,:)/acf(i,2);
end

if size(fv,2) ~= size(spec,2)
   fv = fv(2:end); 
end

% Read in geometry parameters
formatSpec = '%f%f%f%f%f%f%f%f%f%f%f%f%f%f%[^\n\r]';
fileID = fopen(infoname,'r','n','UTF-8');
fseek(fileID, 3, 'bof');
dataArray = textscan(fileID, formatSpec, 'Delimiter', ',', 'TextType', 'string',  'ReturnOnError', false);
fclose(fileID);

% 1) AZ, 2) EL, 3) RANGE, 4) LAT, 5) LON, 6) ALT, 7) AP_LAT, 8) AP_LON, 9) L_SHELL,
% 10) INV_LAT, 11) ASPECT, 12) CXR, 13) CYR, 14) CZR

range = dataArray{:,3};
alt = dataArray{:,6};
asp = dataArray{:,11}-90;

% Constrain altitude range
nind = [find(abs(alt-plot_height(1))==min(abs(alt-plot_height(1)))),...
    find(abs(alt-plot_height(2))==min(abs(alt-plot_height(2))))];
%nind(1) = 1;


%% Read in IRI parameters

delimiter = ' ';
formatSpec = '%f%f%f%f%f%f%f%[^\n\r]';
fileID = fopen(IRIname,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string',  'ReturnOnError', false);
fclose(fileID);
iri = [dataArray{1:end-1}];
clearvars filename delimiter formatSpec fileID dataArray ans;

% Interpolate to measurement altitudes
h_iri = iri(:,1);
ne_iri = interp1(h_iri, iri(:,2), alt);
Ti_iri = interp1(h_iri, iri(:,3), alt);
Te_iri = interp1(h_iri, iri(:,4), alt);
nOr_iri = interp1(h_iri, iri(:,5), alt)/100; % Fraction of O+
nO2_iri = interp1(h_iri, iri(:,6), alt)/100; % Fraction of O2+ or NO+
nNO_iri = interp1(h_iri, iri(:,7), alt)/100; % Fraction of NO+ or O2+


% Millstone Hill Chapman profile
zz1 = min(-(alt-120)/40, 50);
HS = 10-6*exp(zz1);
zz2 = min(-(alt-180)./HS, 50);
nOr_chap = 1 - (1 - 2./(1 + sqrt(1 + 8*exp(zz2))));

figure(10); clf
semilogx(ne_iri,alt)
hold on
plot(ne_iri.*nOr_iri,alt)
plot(ne_iri.*nOr_chap,alt)
plot(ne_iri.*(nNO_iri+nO2_iri),alt);
legend('ne', 'nO IRI', 'nO Chap', 'nNO + NO2')
axis([10^9 10^12 100 400])

if ion_comp == 1
    nOr = nOr_chap;
elseif ion_comp == 2
    nOr = nOr_iri;
end





%% Plots


% Select a few aspect angles
for i= 1:length(plot_angles)
    aind(i) = find(abs(asp-plot_angles(i))==min(abs(asp-plot_angles(i))),1);
    if aind(i) < nind(1)
        aind(i) = nind(1);
    elseif aind(i) > nind(2)
        aind(i) = nind(2);
    end
    alegend{i} = [num2str(plot_angles(i)), '^o'];
end

if plot_var == 1
    yvar = range;
    ytit = 'Range (km)';
elseif plot_var == 2
    yvar = alt;
    ytit = 'Altitude (km)';
elseif plot_var == 3
    yvar = asp;
    ytit = 'Aspect Angle (deg)';
end


figure(1); clf
subplot(1,2,1)
p = pcolor(lags, yvar(nind(1):nind(2)), acf_norm(nind(1):nind(2),:));
caxis([-1 1])
set(p,'EdgeColor','none');
colormap jet
colorbar
xlabel('Time Lag (\mus)', 'FontSize', 16);
ylabel(ytit, 'FontSize', 16);
title('ACF, Normalized', 'FontSize', 20);
set(gca, 'FontSize', 16, 'LineWidth', 1.8);
axis([0 lags(end) yvar(nind(1)) yvar(nind(2))]);

subplot(1,2,2)
p = pcolor(lags, yvar(nind(1):nind(2)), (acf(nind(1):nind(2),:)));
%caxis([-1 1])
set(p,'EdgeColor','none');
colormap jet
colorbar
xlabel('Time Lag (\mus)', 'FontSize', 16);
ylabel(ytit, 'FontSize', 16);
title('ACF, Measured', 'FontSize', 20);
set(gca, 'FontSize', 16, 'LineWidth', 1.8);
axis([0 lags(end) yvar(nind(1)) yvar(nind(2))]);


figure(2); clf
subplot(1,2,1)
p = pcolor(fv, yvar(nind(1):nind(2)), (spec_norm(nind(1):nind(2),:)));
caxis([0 1])
set(p,'EdgeColor','none');
colormap jet
colorbar
xlabel('Frequency (\mus)', 'FontSize', 16);
ylabel(ytit, 'FontSize', 16);
title('Spectra, Normalized', 'FontSize', 20);
set(gca, 'FontSize', 16, 'LineWidth', 1.8);
axis([fv(1) fv(end) yvar(nind(1)) yvar(nind(2))]);

subplot(1,2,2)
p = pcolor(fv, yvar(nind(1):nind(2)), (spec(nind(1):nind(2),:)));
caxis([0 5*10^-4])
set(p,'EdgeColor','none');
colormap jet
colorbar
xlabel('Frequency (\mus)', 'FontSize', 16);
ylabel(ytit, 'FontSize', 16);
title('Spectra, Measured', 'FontSize', 20);
set(gca, 'FontSize', 16, 'LineWidth', 1.8);
axis([fv(1) fv(end) yvar(nind(1)) yvar(nind(2))]);



figure(3); 
subplot(1,2,pp)
p = pcolor(fv, yvar(nind(1):nind(2)), (spec_r2(nind(1):nind(2),:)));
caxis([0 300])
set(p,'EdgeColor','none');
colormap jet
colorbar
xlabel('Frequency (\mus)', 'FontSize', 16);
ylabel(ytit, 'FontSize', 16);
title(['Spectra, Pos.', pos], 'FontSize', 20);
set(gca, 'FontSize', 16, 'LineWidth', 1.8);
axis([fv(1) fv(end) yvar(nind(1)) yvar(nind(2))]);



figure(4); %clf
hold on
plot(sum(spec_r2(nind(1):nind(2),:),2),  yvar(nind(1):nind(2)))
xlabel('Ion line power', 'FontSize', 16)
ylabel(ytit, 'FontSize', 16)
set(gca, 'FontSize', 16, 'LineWidth', 1.8);