function [] = lightning_and_iar_functions(data_folder, datestring, randstring, stationstring, fstring, starthour, stophour, npt, direction, filterchoice, utc_or_mlt, winuniquestr, spikeloc, local, nonlocal, subsize, heightfact, heightfact2, strikepeakdt, peakwint, IARwin)

%% Why the name iar_lightning_6plot.m?
% Well, once upon a time Charles wrote a function called
% iar_lightning_4plot.m that plotted four subplots in one figure, one for the
% PSD, one for the time series, and two maps of where lightning occured.
% Then another function was added to plot two more graphs and the function
% was updated to be called iar_lightning_6plot.m. Maybe it is time to
% rename it as the lightning_and_iar_functions.m.
%
% This is another top level function for plotting data for studying the
% lightning excitation of the IAR. Rough guess is that 80% of the code is
% plotting related, the other 20% is data processing.
%
% By: Charles Nokes
% Latest revision: January 2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Inputs:
% data_folder: see lightning_and_iar.m
% datestring: 'YYYYMMDD'
% randstring: see lightning_and_iar.m
% stationstring: four character string denoting the CARISMA station. E.g.
% for Ministik Lake it is 'MSTK'
% fstring: string indicating the sample rate, '20' or '100'
% starthour: 0-23, initial time for the window being generated (min 0)
% stophour: 0-23, end time for window being generated (max 23)
% npt:  power of two points for fft calculation, the PSD window size,
% see overlap example below
% direction: mag direction to analyze, 'H', 'D', 'H\_band', or 'D\_band', see lightning_and_iar.m
% filterchoice: 'yes' or 'no', toggle filtering the time series to remove DC component
% utc_or_mlt: 'UTC' or 'MLT'; toggle analysis in utc time or mlt time
% winuniquestr: 'unique' or ''; used in plot_strikes_peaks_coin.m and
% menmednormpol.m, for setting only superposition windows with one
% strike/peak coincidence.
% spikeloc: see lightning_and_iar.m, no longer really used
% local: Define lightning local to a CARISMA station such as within 1000
% km. Or set to 'outside' to only look at lightning a minimum distance
% away.
% nonlocal: Define lightning nonlocal to a CARISMA station as between 1000 km and 2000 km
% subsize: plotting subset variable (remove)
% heightfact: quantile for B field peak threshold, between
% 0 and 1. 0 sets a low threshold, 1 sets the highest threshold. See figure
% 3.1 of Charles Nokes' MSc. Thesis, Lightning Excitation of the IAR
% heightfact2: lightning energy threshold, quantile based
% strikepeakdt: propagation time limits
% peakwint: superposed window size
% IARwin: start and end hours of IAR and bin size data dependent
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Functions (indented functions are used inside the function):
% carisma_icm_psd.m
    %  importCARISMAICMFolder.m
        % importCARISMAICM.m
% plot_psd.m
% winrms.m
% IAR_quiet_wins.m
% winlnum.m
% winpeaksrms.m
% matchedtimes.m
% rand_peak_times.m
% plot_strikes_peaks_coin.m
% plot_icm_lightning.m
    % samplewindow.m
% plot_lightning_map.m
% dt_hist_code.m
% menmednormpol.m
    % singsidefft.m
% plot_density_scatter.m
    % dscatter.m
    % medians_B_d.m
% plot_property_vs_winltimes
% create_hist_vec.m
% plot_polarity_latlon.m
% singsidefft.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Important Outputs:
% time: this vector contains the times of each of the mag data points, in
% datetime format
% B_output: this is the mag date, together with the time vector they form
% the time series. B_output is defined by the user selected direction ('H'
% or 'D') and whether or not to filter out the DC.
% fs = 20 or 100; sample frequency, no need to declare here, it is an
% output of function carisma_icm_psd.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Global plotting variables
global fntsize res xsize ysize lmarg bmarg widthv heightv xstr titlestr fstr

%% Global analysis variables
global fs overlap

overlap=round(0.3*npt); % Overlap for psd. Derived from npt. In order
% to produce a coherent PSD it is useful to take
% FFT's of overlapped windows of time series
% datapoints, not sequenced windows. Here's a
% visual emplanation. Imagine this is for 100 samp/sec
% data and each 'window' has 2^14 time series
% points.  Sequenced windows 'look' like this (a different
% bracket is used for each window in the example because
% it makes the overlap case below clearer):
% |_______________|[_______________]{_______________}...
% Overlapped windows, overlapping 20% (3/15 underscores
% in this example) look like this:
% |____________[___|_________{___]____________}...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Figure properties

fntsize=15;
res=170;
xsize=6*res;
ysize=2.5*res;
lmarg=0.5*res;
bmarg=0.5*res;
widthv=5.1*res;
heightv=1.7*res;
fstr = 'Frequency (Hz)';
if strcmp(utc_or_mlt,'UTC')
    xstr = 'Universal Time (hh:mm)';
elseif strcmp(utc_or_mlt,'MLT')
    xstr = 'Mean Local Time (hh:mm)';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Data import
STN_folder = strcat(data_folder, datestring, randstring, '/', stationstring, '/', fstring); % directory needed

% Lightning folder
lfilename = strcat(data_folder, datestring, randstring, '/AE',datestring,'.mat'); % Filename for lightning data

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

if strcmp(utc_or_mlt,'UTC')
    % Call function to import ICM data and calculate PSD for UTC
    [time, B_output, x, f, t, spectraltime_axis, fs, titlestr, siteData] = carisma_icm_psd(STN_folder, starthour, stophour, npt, direction, filterchoice, utc_or_mlt);
elseif strcmp(utc_or_mlt,'MLT')
    % Call function to import ICM data and calculate PSD for MLT; utc_or_mlt must be set to 'MLT'; needs fixing
    [time, B_output, x, f, t, spectraltime_axis, fs, titlestr, siteData, mlt_time, mlt_starttime, xspectrstart, spectraltime_axis_mlt] = carisma_icm_psd(folder, starthour, stophour, npt, direction, filterchoice, utc_or_mlt);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Figure characteristics
% This section is to set up some figure and axis parameters. CHANGE AT
% WILL! If an error is thrown that 's' is not defined or an axis is not
% defined, check here, and fix as needed. Also, feel free to change the
% definition of each axis. If you want a 3x1 figure or a 2x3 figure, choose
% what is needed. This code below just makes sure the rest of the code will
% work without errors until you are comfortable making your own changes.

s = get(0, 'ScreenSize');
ufi0 = input('New figure fit to screen? (yes/no): ','s');
if strcmp(ufi0,'yes')
    f1 = figure('Position', [0 0 s(3) s(4)]); % New figure
end

% Axes handles
ufi1 = input('Are you running plot_psd to make a PSD? (yes/no): ','s');
if strcmp(ufi1,'yes')
    ax1 = subplot(3,1,1); %subplot(3,2,1); %axes('Position',[0.08 0.7 0.4 0.2]); % Left half, top
    ax2 = subplot(3,1,2); %subplot(2,2,2); % Right half bottom box
else
    ufi2 = input('Are you running plot_strikes_peaks_coin to make a plot of mag pulses, lightning strikes, and matches? (yes/no): ','s');
    if strcmp(ufi2,'yes')
        ax2 = subplot(3,1,2); %subplot(2,2,2); % Right half bottom box
    end
end

ufi3 = input('Are you running plot_lightning_map to make a plot of lightning strike locations? (yes/no): ','s');
if strcmp(ufi3,'yes')
    ax3 = subplot(3,1,3); %subplot(3,2,3); %axes('Position',[0.08 0.4 0.4 0.2]); % Left half, middle
    ax5 = 'off'; %subplot(3,2,5); %axes('Position',[0.55 0.5 0.42 0.5],'Box','on'); % Right half top upper box
end

ufi4 = input('Are you running plot_property_vs_winltimes? (yes/no): ','s');
if strcmp(ufi4,'yes')
    ax4 = subplot(2,2,4);
end

ufi5 = input('Are you running plot_icm_lightning to colour code the time series based on lightning? (yes/no): ','s');
if strcmp(ufi5,'yes')
    ax5 = subplot(3,2,5); %axes('Position',[0.55 0.5 0.42 0.5],'Box','on'); % Right half top upper box
    ax6 = subplot(3,2,4);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function calls to process data, calculate values of interest, and plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot power spectral density plot (PSD)
ufu0 = input('Generate plot of PSD? (yes/no): ','s');
if strcmp(ufu0,'yes')
    if strcmp(utc_or_mlt,'UTC')
        % Plot psd for UTC
        plot_psd(x, f, t, time,  spectraltime_axis, B_output, f1, ax1, ax2, 'on', 'off', siteData.unit);
    elseif strcmp(utc_or_mlt,'MLT')
        % Plot PSD for MLT
        plot_psd(x, f, t, mlt_time,  spectraltime_axis_mlt, B_output, f1, ax1, ax2, 'on', 'off', siteData.unit);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Calculate RMS of B for each psd window
[Brms, Brmssmooth, Brmsquiet] = winrms(npt,overlap,spectraltime_axis,B_output);
Brms = rms(Brms);
quietwin_index = round(Brmsquiet(1)*(npt-overlap)*1/IARwin(3)/(3600*fs));
if quietwin_index > 48
    quietwin_index = 48;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Calculate windows of quiet period and IAR periods
[quietwin, IARwins, lwinsvec, midwinpoints] = IAR_quiet_wins(Brmsquiet(1),IARwin,length(time),npt);
lwinsvec = [lwinsvec length(time)];
midwinpoints = [midwinpoints lwinsvec(end-1)+IARwin(3)*3600*fs/2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check if input dt for matching allows for lightning anywhere in the world
max_Earth_radius = 6378; % Greatest Earth radius in km
if strikepeakdt(2) < pi * max_Earth_radius / 300000
    max_arc_len_c = strikepeakdt(2) * 300000;  % max distance a matching strike can be away from the ICM station in km
else
    max_arc_len_c = pi * max_Earth_radius;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Calculate number of lightning strikes within the defined zone for each window
[winltimes, winlnums, winllatlon, winlEnergy] = winlnum(lfilename, siteData, local, nonlocal, time, lwinsvec, midwinpoints, npt, overlap, peakwint, max_arc_len_c);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Calculate number of magnetic peaks above threshold in each window
[winpeaksrmstime, winpeaksrmsnumt, winpeaksrmsfract, minpeakheight] = winpeaksrms(npt, overlap, lwinsvec, midwinpoints, time, B_output, heightfact, Brmsquiet(2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Find the matches between lightning strikes and magnetic peaks above the threshold
[wintlims, ltimes, ptimes, ptimesnum, llatlon, ldist, lEnergy, matched_dt, wintr3, wintc3] = matchedtimes(winltimes, winllatlon, winlEnergy, winpeaksrmstime, strikepeakdt, peakwint, time, B_output, siteData, heightfact2);
for i = 1:length(winltimes)
    matched_B{i,1} = B_output(ptimes{i,1});
    matched_B{i,2} = B_output(ptimes{i,1}(ptimes{i,13}));
end

ufu1 = input('Match for random case (random case needs fixing)? (yes/no): ','s');
if strcmp(ufu1,'yes')
    [rand_ltimes, rand_ptimes, rand_ptimesnum, rand_llatlon, rand_lEnergy, rand_vecl, rand_vec] = rand_peak_times(winltimes, winllatlon, winlEnergy, winpeaksrmstime, strikepeakdt, peakwint, B_output, siteData, wintr3, wintc3);
    rand_matched_B{i} = B_output(rand_ptimes{i,1});
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Saving outputs to .mat file for live analysis or debugging
saveoutputs = input('Save matched times, lat, long, N, matched_dt, B and El? (yes/no): ','s');
if strcmp(saveoutputs,'yes')
    if strcmp(local,'outside')
        titlestr{1}(titlestr{1} == ' ')= '-';
        titlestr{1}(24) = '';
        if strcmp(ufu1,'no')
            save(strcat(titlestr{1},'_matched_times_lat_lon_N_',local,'_',num2str(nonlocal),'_',num2str(strikepeakdt),'.mat'),'titlestr','IARwin','quietwin_index','llatlon','local','nonlocal','matched_dt','matched_B','lEnergy','wintr3','wintc3');
        elseif strcmp(ufu1,'yes')
            save(strcat(titlestr{1},'_matched_times_lat_lon_N_',local,'_',num2str(nonlocal),'_',num2str(strikepeakdt),'.mat'),'titlestr','IARwin','quietwin_index','llatlon','local','nonlocal','matched_dt','matched_B','lEnergy','wintr3','wintc3','rand_vecl','rand_vec','rand_llatlon','rand_lEnergy','rand_matched_B');
        end
        x_lpm = time(midwinpoints);
        y_l = winlnums;
        y_p = winpeaksrmsnumt;
        y_m = ptimesnum(1,:);
        save(strcat(titlestr{1},'_threshold_matching_statistics_',local,'_',num2str(nonlocal),'_t',num2str(strikepeakdt(1)),'_t',num2str(strikepeakdt(2)),'_quartile_',num2str(heightfact),'.mat'),'x_lpm','y_l','y_p','y_m');
    else
        titlestr{1}(titlestr{1} == ' ')= '-';
        titlestr{1}(24) = '';
        if strcmp(ufu1,'no')
            save(strcat(titlestr{1},'_matched_times_lat_lon_N_',num2str(local),'_',num2str(nonlocal),'_',num2str(strikepeakdt),'.mat'),'titlestr','IARwin','quietwin_index','llatlon','local','nonlocal','matched_dt','matched_B','lEnergy','wintr3','wintc3');
        elseif strcmp(ufu1,'yes')
            save(strcat(titlestr{1},'_matched_times_lat_lon_N_',num2str(local),'_',num2str(nonlocal),'_',num2str(strikepeakdt),'.mat'),'titlestr','IARwin','quietwin_index','llatlon','local','nonlocal','matched_dt','matched_B','lEnergy','wintr3','wintc3','rand_vecl','rand_vec','rand_llatlon','rand_lEnergy','rand_matched_B');
        end
        x_lpm = time(midwinpoints);
        y_l = winlnums;
        y_p = winpeaksrmsnumt;
        y_m = ptimesnum(1,:);
        save(strcat(titlestr{1},'_threshold_matching_statistics_',num2str(local),'_',num2str(nonlocal),'_t',num2str(strikepeakdt(1)),'_t',num2str(strikepeakdt(2)),'_quartile_',num2str(heightfact),'.mat'),'x_lpm','y_l','y_p','y_m');
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot of the total number of mag pulses above the threshold in each time bin, the total number of lightning strikes within the defined region within each time bin, and the number of cases of lightning associated mag pulses in each time bin
ufu2 = input('Plot mag pulses, lightning strikes, and matched cases? (yes/no): ','s');
if strcmp(ufu2,'yes')
    if strcmp(ufu1,'no')
        plot_strikes_peaks_coin(ax2, winuniquestr,midwinpoints, time, winlnums, winpeaksrmsnumt, winpeaksrmsfract, ptimesnum, strikepeakdt,datestring, stationstring);
    elseif strcmp(ufu1,'yes') %random case below
        plot_strikes_peaks_coin(ax3, midwinpoints, time, winlnums, winpeaksrmsnumt, winpeaksrmsfract, rand_ptimesnum, strikepeakdt);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot ICM time series window and associated FFT
ufu3 = input('Plot just ICM time series and its FFT in separate figure? (yes/no): ','s');
if strcmp(ufu3,'yes')
    figure;
    ax5a = subplot(2,1,1);
    ax6a = subplot(2,1,2);
    [sampwin, samplloc, sampllocl, sampllocp] = plot_icm_lightning(lfilename, siteData, time, B_output, npt, ax5a, ax6a, spikeloc, local, nonlocal, subsize, 9, 'log', 'smooth');
    
    % Mark which window in the PSD is being examined with a asterisk
    figure(f1);
    axes(ax1);
    hold on;
    plot(datenum(sampwin{1}(length(sampwin{1})/2)),0,'k*','MarkerSize',10);
    hold off;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot the lightning map
ufu4 = input('Plot map with colour coded locations of lightning strikes? (yes/no): ','s');
if strcmp(ufu4,'yes')
    samplloc = [];
    sampllocl = [];
    sampllocp = [];
    if strcmp(ufu1,'no')
        plot_lightning_map(lfilename, time, 'off', ax5, ax3, siteData, samplloc, sampllocl, sampllocp, local, nonlocal,ltimes,llatlon,IARwin);
    elseif strcmp(ufu1,'yes') % random case below
        plot_lightning_map(lfilename, time, 'off', ax5, ax2, siteData, samplloc, sampllocl, sampllocp, local, nonlocal,rand_ltimes,rand_llatlon,IARwin);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot B magnitude or lightning energy vs. distance or matched dt
ufu5 = input('Plot histogram of magnitude of B pulses or lightning energy vs. distance or matched dt? (yes/no): ','s');
if strcmp(ufu5,'yes')
    if strcmp(ufu1,'no')
        dt_hist_code(IARwin, quietwin_index, llatlon, matched_dt, matched_B, lEnergy);
    elseif strcmp(ufu1,'yes') % Random version
        dt_hist_code(IARwin, quietwin_index, rand_llatlon, rand_matched_dt, rand_matched_B, rand_lEnergy);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot peak shape and FFT resulting from Superposed Epoch Analysis
ufu6 = input('Calculate and plot superposed epoch analysis time series and FFT? (yes/no): ','s');
% String inputs for menmenormpol.m are:
% ('mean' or 'quartiles' or 'deciles'), ('unnormalized' or 'normalized')
% ('normalpolarity' or 'positivepolarity'), (rmsthreshold or Brmsquiet(2)),
% filter by lEnergy ('on' or '')
%
% Include only superposition windows with one strike/peak coincidence by
% setting the winunique variable to 'unique', otherwise ''.
if strcmp(ufu6,'yes')
    % Define variable to make guessing where echo pulse is more accurate
    delayguess = 3/10; %5/8
    if strcmp(ufu1,'no')
        [R, C, dt_reflection] = menmednormpol(quietwin, IARwins, IARwin, 'mean', 'unnormalized', 'positivepolarity','','', ptimes, B_output, peakwint,s,'',minpeakheight,Brmsquiet,datestring,stationstring, wintlims, delayguess);
    elseif strcmp(ufu1,'yes') % Random version
        [rand_R, rand_C] = menmednormpol(quietwin, IARwins, IARwin, 'quartiles', 'unnormalized', 'positivepolarity', rand_ptimes, B_output, peakwint,s,'',rmsthreshold); % ('mean' or 'quartiles' or 'deciles') ('unnormalized' or 'normalized') ('normalpolarity' or 'positivepolarity') (rmsthreshold or Brmsquiet(2))
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Identify the echo pulse lag times for each time bin
ufu7 = input('Identify the echo pulse lag times for each time bin in order to plot over the PSD in next function? (yes/no): ','s');
if strcmp(ufu7,'yes')
    figure;
    plot(dt_reflection(:,1),'*');
    hold on
    plot(dt_reflection(:,2),'x');
    plot(dt_reflection(:,3),'o');
    disp(dt_reflection);
    dt_reflection_locs = input('Actual reflection indices [i,j]?');
    dt_reflection_actual = dt_reflection(sub2ind(size(dt_reflection),dt_reflection_locs(:,1),dt_reflection_locs(:,2)));
    dt_reflection_wins = IARwins(dt_reflection_locs(:,1))+(IARwin(3)/2*3600*fs);
    titlestr{1}(titlestr{1} == ' ')= '-';
    titlestr{1}(24) = '';
    save(strcat(titlestr{1},'_dt_reflection_',randstring,num2str(local),'_',num2str(nonlocal),'_t',num2str(strikepeakdt(1)),'_t',num2str(strikepeakdt(2)),'_quartile_',num2str(heightfact),'.txt'), 'dt_reflection_actual','dt_reflection_wins','-ascii');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot reciprocals of echo pulse lag time for each time bin on the psd
ufu8 = input('Plot reciprocals of echo pulse lag time for each time bin on the psd? (yes/no): ','s');
if strcmp(ufu8,'yes')
    axes(ax1);
    hold on
    plot(ax1,datenum(time(dt_reflection_wins)),1./(dt_reflection_actual),'m-o');
    plot(ax1,datenum(time(dt_reflection_wins)),2./(dt_reflection_actual),'m-o');
    plot(ax1,datenum(time(dt_reflection_wins)),3./(dt_reflection_actual),'m-o');
    plot(ax1,datenum(time(dt_reflection_wins)),4./(dt_reflection_actual),'m-o');
    plot(ax1,datenum(time(dt_reflection_wins)),5./(dt_reflection_actual),'m-o');
    plot(ax1,datenum(time(dt_reflection_wins)),6./(dt_reflection_actual),'m-o');
    plot(ax1,datenum(time(dt_reflection_wins)),7./(dt_reflection_actual),'m-o');
    plot(ax1,datenum(time(dt_reflection_wins)),8./(dt_reflection_actual),'m-o');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Scatter plot of distance vs. lag time dt between strike and mag peak,
% Options for density highlighting with colour, contour, or surface plot
ufu9 = input('Scatter plot of distance vs. lag time dt between strike and mag pulse? (yes/no): ','s');
if strcmp(ufu9,'yes')
    [Centroids] = plot_density_scatter(IARwin, quietwin_index, matched_dt(:,2),'abs', '', ldist(:,2),'','','Matched dt (s)','Distance (km)',s,datestring,stationstring);
% Uncomment below for other scenarios
    %[Centroids] = plot_density_scatter(IARwin, quietwin_index, matched_dt,'abs','', lEnergy,'','','Matched dt (s)','Lightning Energy (J)',s,datestring,stationstring);
    %[Centroids] = plot_density_scatter(IARwin, quietwin_index, matched_dt,'abs','', matched_B,'abs','log','Matched dt (s)','Matched B (nT)',s,datestring,stationstring);
    %[Centroids] = plot_density_scatter(IARwin, quietwin_index, lEnergy,'', ldist,'','Lightning Energy (J)','Distance (km)',s,datestring,stationstring);
    %[Centroids] = plot_density_scatter(IARwin, quietwin_index, ldist(:,2),'', '', matched_B(:,2),'abs','log','Distance (km)','Matched B (nT)',s,datestring, stationstring);
    %[Centroids] = plot_density_scatter(IARwin, quietwin_index, matched_B,'abs','',lEnergy,'','log','Matched B','lEnergy',s,datestring,stationstring);
    savecentroids = input('Save centroid datapoints to .txt file? (yes/no): ','s');
    if strcmp(savecentroids,'yes')
        Ctosave = [];
        for i=2:length(Centroids)
            Ctosave = [Ctosave; Centroids{i}(:,1),Centroids{i}(:,2)];
        end
        titlestr{1}(titlestr{1} == ' ')= '-';
        titlestr{1}(24) = '';
        save(strcat(titlestr{1},'_B_vs_dt_Maximas_all_day',randstring,num2str(local),'_',num2str(nonlocal),'_t',num2str(strikepeakdt(1)),'_t',num2str(strikepeakdt(2)),'_quartile_',num2str(heightfact),'.txt'), 'Ctosave','-ascii');
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot lEnergy quantiles for each 30 minute window
ufu10 = input('Plot lEnergy quantiles for each 30 minute window? (yes/no): ','s');
if strcmp(ufu10,'yes')
plot_property_vs_winltimes(ax4, midwinpoints, time, lEnergy)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot histograms of N vs. B above threshold and N vs. matched B
ufu11 = input('Plot histograms of N vs. B above threshold and N vs. matched B? (yes/no): ','s');
if strcmp(ufu11,'yes')
create_hist_vec('unmatched', 'intervals', winltimes, ptimes, B_output, quietwin_index, IARwin, lwinsvec, winpeaksrmstime, s, quietwin, IARwins);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Save mean of matched B in each window to compare with other stations
ufu12 = input('Save mean of matched B in each window to compare with other stations? (yes/no): ','s');
if strcmp(ufu12,'yes')
for i = 1:length(matched_B)
y_B(i) = mean(abs(matched_B{i}));
end
x_B = time(midwinpoints);
titlestr{1}(titlestr{1} == ' ')= '-';
titlestr{1}(24) = '';
save(strcat(titlestr{1},'_mean_matched_B_',num2str(local),'_',num2str(nonlocal),'_t',num2str(strikepeakdt(1)),'_t',num2str(strikepeakdt(2)),'_quartile_',num2str(heightfact),'.mat'), 'y_B', 'x_B');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plot peak polarity as function of lightning strike lat lon position in quadrants
ufu13 = input('Plot peak polarity as function of lightning strike lat lon position in quadrants? (yes/no): ','s');
if strcmp(ufu13,'yes')
plot_polarity_latlon(llatlon, siteData, quietwin, IARwin, 'frames', R, C, s); % 'day' or 'frames'
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end