function [sampwin, samplloc, sampllocl, sampllocp] = plot_icm_lightning(lfilename, siteData, time, B_output, npt, ax2n1, ax3n2, spikeloc, local, proximal, subsize, smoothspan, logging, smoothing)
% function to take in the sample set time series ICM data for a window in
% of time in the PSD and the lightning local and proximal to that ICM data
% in that time frame, and plot them together using plotyy.

global fntsize res xsize ysize lmarg bmarg widthv heightv xstr titlestr fstr pstr

% Global analysis variables
global fs overlap

ldata = importdata(lfilename); % Load WWLLN "data" cell

% Assign data to vectors
%Time = datetime(ldata(:,1:6)); % Date/time in UTC, datetime vector
Latitude=ldata(:,7); % latitude vector: Lat, Lon in fractional degrees (geographic coordinated)
Longitude=ldata(:,8); % longitude vector
Resid=ldata(:,9); % resid is the residual fit error in Microseconds (always <30)
nstn=ldata(:,10); % nstn is the number of WWLLN stations participarting in the location fit
Energy=ldata(:,11); % energy is the RMS energy (in Joules) of the stroke (from 1.3 ms waveform sampling between 7 and 18 kHz)
E_uncertainty=ldata(:,12);% energy uncertainty (energy error of the fit in Joules)
nstn_energy=ldata(:,13);% nstn_eneregy is the subset of nstn stations between 1000 and 8000 km from the stroke whose energy data were used in the energy estimate.

% Map plotting of WWLLN Data between 5 and 12 UTC
STN_Lat=siteData.geoLat; STN_Long=siteData.geoLon-360;

% Find data points inside and outside the circle
[arclen,az] = distance(STN_Lat,STN_Long,Latitude,Longitude);

if strcmp(local,'outside')
    ltimes = datetime(ldata(arclen > km2deg(proximal,'earth'),1:6),'TimeZone','UTC'); % Date/time in UTC, datetime vector;
    lenergy = ldata(arclen > km2deg(proximal,'earth'),11); % Energy withint this radius;
    llat = ldata(arclen > km2deg(proximal,'earth'),7);
    llng = ldata(arclen > km2deg(proximal,'earth'),8);
else
    % Create vectors of times when lightning strikes occur locally and
    % proximally of associated station
%     ltimesl = datetime(ldata(arclen < km2deg(local,'earth'),1:6),'TimeZone','UTC'); % Date/time in UTC, datetime vector
%     lenergyl = ldata(arclen < km2deg(local,'earth'),11); % Energy withint this radius
%     ltimesp = datetime(ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),1:6),'TimeZone','UTC'); % Date/time in UTC, datetime vector
%     lenergyp = ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),11); % Energy between these radii
%     llatl = ldata(arclen < km2deg(local,'earth'),7);
%     llngl = ldata(arclen < km2deg(local,'earth'),8);
%     llatp = ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),7);
%     llngp = ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),8);
    ltimes = datetime(ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),1:6),'TimeZone','UTC');
    lenergy = ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),11);
    llat = ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),7);
    llng = ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),8);
end

% sampwin is cell with datetime time vector and data vector for ICM time series data for one psd window; set1 is a cell of the first half set; set2 is a cell of the second half set ; subset1 is same as set 1 but in 10ths ; subset2 is same as set 2 but in 10ths;
%samplewindow(time, B_output, spikeloc, '', 9,'log', siteData.unit); %
[sampwin Ntitlestr] = samplewindow(time, B_output,...
    spikeloc, 'smooth', 9,'log', siteData.unit, npt, 'off'); % final input value is the smoothing span

% Identify subset of lightning strikes inside the geographic region that
% happen within the sampwin time frame
if strcmp(local,'outside')
    sampltimes = ltimes(ltimes > sampwin{1}(1) & ltimes < sampwin{1}(end));
    samplenergy = lenergy(ltimes > sampwin{1}(1) & ltimes < sampwin{1}(end));
else
%     sampltimesl = ltimesl(ltimesl > sampwin{1}(1) & ltimesl < sampwin{1}(end));
%     samplenergyl = lenergyl(ltimesl > sampwin{1}(1) & ltimesl < sampwin{1}(end));
    %sampllocl = {llatl(ltimesl > sampwin{1}(1) & ltimesl < sampwin{1}(end)) llngl(ltimesl > sampwin{1}(1) & ltimesl < sampwin{1}(end))};
    % yvecl = max(sampwin{2})*ones(1,length(sampltimesl));
%     sampltimesp = ltimesp(ltimesp > sampwin{1}(1) & ltimesp < sampwin{1}(end));
%     samplenergyp = lenergyp(ltimesp > sampwin{1}(1) & ltimesp < sampwin{1}(end));
    %sampllocp = {llatp(ltimesp > sampwin{1}(1) & ltimesp < sampwin{1}(end)) llngp(ltimesp > sampwin{1}(1) & ltimesp < sampwin{1}(end))};
    % yvecp = max(sampwin{2})*ones(1,length(sampltimesp));
    sampltimes = ltimes(ltimes > sampwin{1}(1) & ltimes < sampwin{1}(end));
    samplenergy = lenergy(ltimes > sampwin{1}(1) & ltimes < sampwin{1}(end));
end

if strcmp(ax2n1,'off')
    samplloc = [];
        sampllocl = [];
        sampllocp = [];
else
% Plot lightning strike times over time series ICM in window
% plot(sampltimesl,yvecl,'mo',sampltimesp,yvecp,'go')
colpl1 = {[1 0 1] [0.5 0 0] [1 0 0] [1 0.9 0] [0.4 0.8 0] [1 0.5 0] [0 0 1] [0.5 0.5 0.5]}; % Predefined colour pattern for set 1
axes(ax2n1);

% if strcmp(local,'outside')
    if isempty(sampltimes)
        for i = 1:subsize
            h1(i) = plot(sampwin{1}((i-1)*length(sampwin{1})/subsize+1:i*length(sampwin{1})/subsize),sampwin{2}((i-1)*length(sampwin{2})/subsize+1:i*length(sampwin{2})/subsize),'color',colpl1{i});
        end
        
        %legend([h1],'ICM','Location','southoutside','Orientation','horizontal');
        %ylim([round(min(sampwin{2})) round(max(sampwin{2}))]);
        ylim([-75 75]);
        yticks([-75 0 75]);
        xlabel('Universal Time (hh:mm:ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
        ylabel(strcat('B (',siteData.unit,')'),'fontsize', fntsize, 'fontname', 'helvetica');
        samplloc = [];
        sampllocl = [];
        sampllocp = [];
    else
        [ax h1 h2] = plotyy(sampwin{1}, sampwin{2},sampltimes,samplenergy,'plot','plot');
        axes(ax(1));
        hold(ax(1));
        for i = 1:subsize
            h8(i) = plot(sampwin{1}((i-1)*length(sampwin{1})/subsize+1:i*length(sampwin{1})/subsize),sampwin{2}((i-1)*length(sampwin{2})/subsize+1:i*length(sampwin{2})/subsize),'color',colpl1{i});
        end
        xlim([sampwin{1}(1) sampwin{1}(end)]);
        %ylim([round(min(sampwin{2})) round(max(sampwin{2}))]);
        ylim([-75 75]);
        yticks([-75 0 75]);
        xlabel('Universal Time (hh:mm:ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
        ylabel(strcat('B (',siteData.unit,')'),'fontsize', fntsize, 'fontname', 'helvetica');
        hold off
        
        axes(ax(2));
        set(h2,'Marker','*', 'MarkerSize', 10, 'Color','k', 'LineStyle','none');
        set(ax(2),'YColor','k');
%         hold(ax(2));
%         h3 = plot(sampltimes,samplenergy,'ko','MarkerSize',10);
        %xlim([sampwin{1}(1) sampwin{1}(end)]);
        %ylim([min(samplenergy) max(samplenergy)]);
        %yticks(linspace(round(min(samplenergy)), floor(max(samplenergy)),3));
        legend([h1,h2],'ICM','Lightning','Location','southoutside','Orientation','horizontal');
        ylabel('Lightning Energy (J)','fontsize', fntsize, 'fontname', 'helvetica');
        hold off
        samplloc = {llat(ltimes > sampwin{1}(1) & ltimes < sampwin{1}(end)) llng(ltimes > sampwin{1}(1) & ltimes < sampwin{1}(end))};
        sampllocl = [];
        sampllocp = [];
    end
% else
%     if isempty(sampltimesl) && isempty(sampltimesp)
%         for i = 1:subsize
%             h1(i) = plot(sampwin{1}((i-1)*length(sampwin{1})/subsize+1:i*length(sampwin{1})/subsize),sampwin{2}((i-1)*length(sampwin{2})/subsize+1:i*length(sampwin{2})/subsize),'color',colpl1{i});
%         end
%         
%         %legend([h1],'ICM','Location','southoutside','Orientation','horizontal');
%         %ylim([round(min(sampwin{2})) round(max(sampwin{2}))]);
%         ylim([-75 75]);
%         yticks([-75 0 75]);
%         xlabel('Universal Time (hh:mm:ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
%         ylabel(strcat('B (',siteData.unit,')'),'fontsize', fntsize, 'fontname', 'helvetica');
%         sampllocl = [];
%         sampllocp = [];
%         samplloc = [];
%         
%     elseif isempty(sampltimesl) && not(isempty(sampltimesp))
%         [ax h1 h2] = plotyy(sampwin{1}, sampwin{2},sampltimesp,samplenergyp,'plot','plot');
%         axes(ax(1));
%         hold(ax(1));
%         for i = 1:subsize
%             h8(i) = plot(sampwin{1}((i-1)*length(sampwin{1})/subsize+1:i*length(sampwin{1})/subsize),sampwin{2}((i-1)*length(sampwin{2})/subsize+1:i*length(sampwin{2})/subsize),'color',colpl1{i});
%         end
%         xlim([sampwin{1}(1) sampwin{1}(end)]);
%         %ylim([round(min(sampwin{2})) round(max(sampwin{2}))]);
%         ylim([-75 75]);
%         yticks([-75 0 75]);
%         xlabel('Universal Time (hh:mm:ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
%         ylabel(strcat('B (',siteData.unit,')'),'fontsize', fntsize, 'fontname', 'helvetica');
%         hold off
%         
%         axes(ax(2));
%         set(h2,'Marker','o', 'MarkerSize', 10, 'Color','k', 'LineStyle','none');
%         set(ax(2),'YColor','k');
%         xlim([sampwin{1}(1) sampwin{1}(end)]);
%         ylabel('Lightning Energy (J)','fontsize', fntsize, 'fontname', 'helvetica');
%         legend([h1,h2],'ICM','Proximal Lightning','Location','southoutside','Orientation','horizontal');
%         samplloc = [];
%         sampllocl = [];
%         sampllocp = {llatp(ltimesp > sampwin{1}(1) & ltimesp < sampwin{1}(end)) llngp(ltimesp > sampwin{1}(1) & ltimesp < sampwin{1}(end))};
%         
%     elseif not(isempty(sampltimesl)) && isempty(sampltimesp)
%         [ax h1 h2] = plotyy(sampwin{1}, sampwin{2},sampltimesl,samplenergyl,'plot','plot');
%         axes(ax(1));
%         hold(ax(1));
%         for i = 1:subsize
%             h8(i) = plot(sampwin{1}((i-1)*length(sampwin{1})/subsize+1:i*length(sampwin{1})/subsize),sampwin{2}((i-1)*length(sampwin{2})/subsize+1:i*length(sampwin{2})/subsize),'color',colpl1{i});
%         end
%         xlim([sampwin{1}(1) sampwin{1}(end)]);
%         %ylim([round(min(sampwin{2})) round(max(sampwin{2}))]);
%         ylim([-75 75]);
%         yticks([-75 0 75]);
%         xlabel('Universal Time (hh:mm:ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
%         ylabel(strcat('B (',siteData.unit,')'),'fontsize', fntsize, 'fontname', 'helvetica');
%         hold off
%         
%         axes(ax(2));
%         set(h2,'Marker','*', 'MarkerSize', 10, 'Color','k', 'LineStyle','none');
%         set(ax(2),'YColor','k');
%         xlim([sampwin{1}(1) sampwin{1}(end)]);
%         ylabel('Lightning Energy (J)','fontsize', fntsize, 'fontname', 'helvetica');
%         legend([h1,h2],'ICM','Local Lightning','Location','southoutside','Orientation','horizontal');
%         samplloc = [];
%         sampllocl = {llatl(ltimesl > sampwin{1}(1) & ltimesl < sampwin{1}(end)) llngl(ltimesl > sampwin{1}(1) & ltimesl < sampwin{1}(end))};
%         sampllocp = [];
%         
%     else
%         [ax h1 h2] = plotyy(sampwin{1}, sampwin{2},sampltimesl,samplenergyl,'plot','plot');
%         axes(ax(1));
%         hold(ax(1));
%         for i = 1:subsize
%             h8(i) = plot(sampwin{1}((i-1)*length(sampwin{1})/subsize+1:i*length(sampwin{1})/subsize),sampwin{2}((i-1)*length(sampwin{2})/subsize+1:i*length(sampwin{2})/subsize),'color',colpl1{i});
%         end
%         xlim([sampwin{1}(1) sampwin{1}(end)]);
%         %ylim([round(min(sampwin{2})) round(max(sampwin{2}))]);
%         ylim([-75 75]);
%         yticks([-75 0 75]);
%         xlabel('Universal Time (hh:mm:ss)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
%         ylabel(strcat('B (',siteData.unit,')'),'fontsize', fntsize, 'fontname', 'helvetica');
%         hold off
%         
%         axes(ax(2));
%         set(h2,'Marker','*', 'MarkerSize', 10, 'Color','k', 'LineStyle','none');
%         set(ax(2),'YColor','k');
%         hold(ax(2));
%         h3 = plot(sampltimesp,samplenergyp,'ko','MarkerSize',10);
%         xlim([sampwin{1}(1) sampwin{1}(end)]);
%         ylim([min([samplenergyl; samplenergyp]) max([samplenergyl; samplenergyp])])
%         yticks(linspace(round(min([samplenergyl; samplenergyp])), floor(max([samplenergyl; samplenergyp])),3));
%         legend([h1,h2,h3],'ICM','Local Lightning','Proximal Lightning','Location','southoutside','Orientation','horizontal');
%         ylabel('Lightning Energy (J)','fontsize', fntsize, 'fontname', 'helvetica');
%         hold off
%         samplloc = [];
%         sampllocl = {llatl(ltimesl > sampwin{1}(1) & ltimesl < sampwin{1}(end)) llngl(ltimesl > sampwin{1}(1) & ltimesl < sampwin{1}(end))};
%         sampllocp = {llatp(ltimesp > sampwin{1}(1) & ltimesp < sampwin{1}(end)) llngp(ltimesp > sampwin{1}(1) & ltimesp < sampwin{1}(end))};
%         
%     end
% end
end

if strcmp(ax3n2,'off')
    samplloc = [];
        sampllocl = [];
        sampllocp = [];
else
axes(ax3n2)
hold on;
for i = 1:subsize
    subset{i,1} = sampwin{1}((i-1)*length(sampwin{1})/subsize+1:i*length(sampwin{1})/subsize); % time vector
    subset{i,2} = sampwin{2}((i-1)*length(sampwin{2})/subsize+1:i*length(sampwin{2})/subsize); % time series y vector
    
    L1 = length(subset{i,1}); % length of subset
    subset{i,3} = fs*(0:(L1/2))/L1; % frequency vector
    sampfft1 = fft(subset{i,2},length(subset{i,2})); % fft of time series y values
    P2samp1 = abs(sampfft1/L1); % two sided spectrum of the fft
    P1samp1 = P2samp1(1:L1/2+1); % single sided spectrum based on P2samp1
    P1samp1(2:end-1) = 2*P1samp1(2:end-1); % single sided spectrum based on P2samp1
    if strcmp(logging,'log')
        subset{i,4} = log(P1samp1); % take log of single sided spectrum and save in the cell
        pstr = 'log|P1(f)|';
    else
        subset{i,4} = P1samp1; % save in the cell
        pstr = '|P1(f)|';
    end
    subset{i,5} = smooth(subset{i,4},smoothspan,'moving'); % smoothing fft of time series y values with a moving average
    
    if strcmp(smoothing,'smooth')
        plot(subset{i,3},subset{i,5},'color',colpl1{i});
        title(strcat('Smoothed single sided amplitude spectrum of B, subset 1-',num2str(subsize)), 'fontsize', fntsize, 'fontname', 'helvetica');
    else
        plot(subset{i,3},subset{i,4},'color',colpl1{i});
        title(strcat('Single sided amplitude spectrum of B, subset 1-',num2str(subsize)), 'fontsize', fntsize, 'fontname', 'helvetica');
    end
    
end
axes(ax3n2);
axis tight;
xlabel(fstr, 'fontsize', fntsize, 'fontname', 'helvetica');
ylabel(pstr, 'fontsize', fntsize, 'fontname', 'helvetica');
hold off
end
end