function [] = plot_lightning_map(lfilename, time, tltng, ax3, ax4, siteData, samplloc, sampllocl, sampllocp, local, proximal,ltimes,llatlon,IARwin)

global fntsize

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.


% generate time-axis for time series
sz=size(Time); time_axis=linspace(hour(Time(1)), hour(Time(end))+minute(Time(end))/60+second(Time(end))/3600, sz(1));

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

% Define plot properties
%fntsize=15; % fontsize

% Plot lightning on map of Canada with various stuff
axes(ax4) % Right, center box
geoshow(ax4,'landareas.shp', 'FaceColor', [0.85 0.85 0.85]); % generate world map
%ylim(ax4,[24 75]) % zoom in to Canada
%xlim(ax4,[-170 -52])
xlabel(ax4,'Longitude (degrees)   ','fontsize',15); % axes labels
ylabel(ax4,'Latitude (degrees)   ','fontsize',15);
hold on
% Highlight lightning occurring during IAR times
% Step 1: find lightning that occurs within a certain geographic zone "local" lightning
%find_WWLLN_MSTK_Lightning = find(ldata(:,7)>=35.5 & ldata(:,7)<=71.5 & ldata(:,8)>=-143 & ldata(:,8)<=-83);

% Step 2: Create a vector of just "local" lighting from Step 1
%WWLLN_MSTK_lightning = ldata(find_WWLLN_MSTK_Lightning,:);

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

% Plot on map where lightning is inside or outside circles
%plot(ax4,ldata(:,8),ldata(:,7),'r.'); % plot all lightning data for that day on map
if strcmp(local,'outside')
    %p41 = plot(ax4,ldata(arclen > km2deg(proximal,'earth'),8),ldata(arclen > km2deg(proximal,'earth'),7),'b.'); % plot lightning outside 2000 km
    %p42 = plot(ax4,ldata(arclen < km2deg(proximal,'earth'),8),ldata(arclen < km2deg(proximal,'earth'),7),'k.'); % plot lightning inside 2000 km
    p41a = plot(ax4,ldata(arclen < km2deg(local,'earth'),8),ldata(arclen < km2deg(local,'earth'),7),'k.'); % plot lightning outside 2000 km
    p41b = plot(ax4,ldata(arclen > km2deg(proximal,'earth'),8),ldata(arclen > km2deg(proximal,'earth'),7),'k.');
    p42 = plot(ax4,ldata(arclen > km2deg(proximal,'earth'),8),ldata(arclen > km2deg(proximal,'earth'),7),'b.');
    [latc,longc] = scircle1(STN_Lat,STN_Long,km2deg(proximal,'earth')); % 2000 km circle around station
    plot(ax4,longc,latc) % plot circle
    radstr = strcat(num2str(proximal),' km'); % string to label the circle
    text(STN_Long-km2deg(2*proximal,'earth'),STN_Lat+km2deg(0.5*proximal,'earth'),radstr); % label the circle on the map
else
    p41a = plot(ax4,ldata(arclen < km2deg(local,'earth'),8),ldata(arclen < km2deg(local,'earth'),7),'k.'); % plot lightning outside 2000 km
    p41b = plot(ax4,ldata(arclen > km2deg(proximal,'earth'),8),ldata(arclen > km2deg(proximal,'earth'),7),'k.'); % plot lightning outside 2000 km
    p42 = plot(ax4,ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),8),ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),7),'b.'); % plot lightning between 1000 and 2000 km
    %p43 = plot(ax4,ldata(arclen < km2deg(local,'earth'),8),ldata(arclen < km2deg(local,'earth'),7),'m.'); % plot lightning inside 1000 km
    
% Map features
%for i =1:2 %fixmeeeee
    %[latc,longc] = scircle1(STN_Lat,STN_Long,km2deg(i*local,'earth')); % 1000-2000 km circle around station
    [latc1,longc1] = scircle1(STN_Lat,STN_Long,km2deg(local,'earth'));
    [latc2,longc2] = scircle1(STN_Lat,STN_Long,km2deg(proximal,'earth'));
    plot(ax4,longc1,latc1) % plot circle
    plot(ax4,longc2,latc2)
    %radstr = strcat(num2str(i*local),' km'); % string to label the circle
    %text(STN_Long-km2deg(2*i*local,'earth'),STN_Lat+km2deg(0.5*i*local,'earth'),radstr); % label the circle on the map
%end
end
plot(STN_Long,STN_Lat,'ks','markersize',5); % Insert a marker for MSTK on the map
text(STN_Long,STN_Lat,siteData.siteID); % Insert text for MSTK on the map

% Lightning matched to B
WWLLN_matched = [];
for i = 1:length(llatlon)
    WWLLN_matched = [WWLLN_matched; llatlon{i}(:,1:2)];
end
p43 = plot(ax4, WWLLN_matched(:,2),WWLLN_matched(:,1),'r.'); % plot the IAR time lightning over the existing lightning on the map

% Find lightning that occurs within the time period where IAR is visible in
% the ICM data and matched to B peaks
%WWLLN_IAR_time = ldata(ldata(:,4)>=IARwin(1) & ldata(:,4)<=IARwin(2),:); 
WWLLN_IAR_matched_time = [];
WWLLN_IAR_matched_datetime = [];
for i = 2*IARwin(1)+1:2*IARwin(2)
    WWLLN_IAR_matched_time = [WWLLN_IAR_matched_time; llatlon{i}(:,1:2)];
    WWLLN_IAR_matched_datetime = [WWLLN_IAR_matched_datetime; ltimes{i}];
end
p44 = plot(ax4, WWLLN_IAR_matched_time(:,2),WWLLN_IAR_matched_time(:,1),'g.','MarkerSize',8); % plot the IAR time lightning over the existing lightning on the map
save(strcat('WWLLN Lightning IAR time matched',' ', datestr(time(1),'yyyy-mm-dd'),'-',siteData.siteID,'-',num2str(IARwin(1)),'-',num2str(IARwin(2)),'.mat'),'siteData','STN_Lat','STN_Long','WWLLN_IAR_matched_datetime','WWLLN_IAR_matched_time');

% Plot format features
title(strcat('WWLLN Lightning',{' '}, datestr(time(1),'yyyy-mm-dd')), 'fontsize', 12, 'fontname', 'helvetica','fontweight','bold');
% Plot the ICM psd sample window local and proximal lightning strikes
%if strcmp(local,'outside')
    %if isempty(samplloc)
       %legend([p42 p41],{strcat(['lightning at <',num2str(proximal),' km']),strcat(['lightning at >',num2str(proximal),' km'])},'Fontsize',12,'Location','southwest'); 
       legend([p42 p43 p44 p41a],{strcat(['lightning within ',num2str(local),'-',num2str(proximal),' km']),'B matched',strcat('B matched - ',num2str(IARwin(1)),'-',num2str(IARwin(2)),' UTC'),'Other Lightning'},'Fontsize',12,'Location','southwest');

       %     else
%         p43 = plot(samplloc{2},samplloc{1},'k*','MarkerSize',10);
%     legend([p42 p41 p43],{strcat(['lightning at <',num2str(proximal),' km']),strcat(['lightning at >',num2str(proximal),' km']), 'Distant lightning'},'Fontsize',12,'Location','southwest');
    %end
% else
% if isempty(sampllocl) && isempty(sampllocp)
%     legend([p43 p42 p41],{strcat(['lightning at <',num2str(local),' km']),strcat(['lightning at ',num2str(local),'-',num2str(proximal),' km']),strcat(['lightning at >',num2str(proximal),' km'])},'Fontsize',12,'Location','southwest');
% elseif isempty(sampllocl) && not(isempty(sampllocp))
%     p45 = plot(sampllocp{2},sampllocp{1},'ko','MarkerSize',10);
%     legend([p43 p42 p41 p45],{strcat(['lightning at <',num2str(local),' km']),strcat(['lightning at ',num2str(local),'-',num2str(proximal),' km']),strcat(['lightning at >',num2str(proximal),' km']), 'Proximal lightning'},'Fontsize',12,'Location','southwest');
%     
% elseif not(isempty(sampllocl)) && isempty(sampllocp)
%     p44 = plot(sampllocl{2},sampllocl{1},'k*','MarkerSize',10);
%     legend([p43 p42 p41 p44],{strcat(['lightning at <',num2str(local),' km']),strcat(['lightning at ',num2str(local),'-',num2str(proximal),' km']),strcat(['lightning at >',num2str(proximal),' km']), 'Local lightning'},'Fontsize',12,'Location','southwest');
%     
% else
%     p44 = plot(sampllocl{2},sampllocl{1},'k*','MarkerSize',10);
%     p45 = plot(sampllocp{2},sampllocp{1},'ko','MarkerSize',10);
%     legend([p43 p42 p41 p44 p45],{strcat(['lightning at <',num2str(local),' km']),strcat(['lightning at ',num2str(local),'-',num2str(proximal),' km']),strcat(['lightning at >',num2str(proximal),' km']), 'Local lightning', 'Proximal lightning'},'Fontsize',12,'Location','southwest');
%     
% end
% end
hold off

% Plot lightning on map of Canada with various stuff
% axes(ax5) % Right, center box
% geoshow(ax5,'landareas.shp', 'FaceColor', [0.85 0.85 0.85]); % generate world map
% ylim(ax5,[24 75]) % zoom in to Canada
% xlim(ax5,[-170 -52])
% xlabel(ax5,'Longitude (degrees)   ','fontsize',15); % axes labels
% ylabel(ax5,'Latitude (degrees)   ','fontsize',15);
% hold on
%
% % Find data points inside and outside the circle
% [arclen,az] = distance(STN_Lat,STN_Long,Latitude,Longitude);
%
% % Plot on map where lightning is inside or outside circles
% %plot(ax4,ldata(:,8),ldata(:,7),'r.'); % plot all lightning data for that day on map
% p41 = plot(ax5,ldata(arclen > km2deg(2000,'earth'),8),ldata(arclen > km2deg(2000,'earth'),7),'k.'); % plot lightning outside 2000 km
% p42 = plot(ax5,ldata(arclen > km2deg(1000,'earth') & arclen < km2deg(2000,'earth'),8),ldata(arclen > km2deg(1000,'earth') & arclen < km2deg(2000,'earth'),7),'g.'); % plot lightning between 1000 and 2000 km
% p43 = plot(ax5,ldata(arclen < km2deg(1000,'earth'),8),ldata(arclen < km2deg(1000,'earth'),7),'m.'); % plot lightning inside 1000 km
%
% % Map features
% for i =1:2
% [latc,longc] = scircle1(STN_Lat,STN_Long,km2deg(i*1000,'earth')); % 1000-2000 km circle around station
% plot(ax5,longc,latc) % plot circle
% radstr = strcat(num2str(i*1000),' km'); % string to label the circle
% text(STN_Long-km2deg(2*i*1000,'earth'),STN_Lat+km2deg(0.5*i*1000,'earth'),radstr); % label the circle on the map
% end
%
% % Find lightning that occurs within the time period where IAR is visible in the ICM data
% WWLLN_IAR_time = ldata(ldata(:,4)>=1 & ldata(:,4)<=11,:); % FIXMEEEEE
% p44 = plot(ax5, WWLLN_IAR_time(:,8),WWLLN_IAR_time(:,7),'b.'); % plot the IAR time lightning over the existing lightning on the map
%
% plot(STN_Long,STN_Lat,'ko','markersize',5); % Insert a marker for MSTK on the map
% text(STN_Long,STN_Lat,siteData.siteID); % Insert text for MSTK on the map
%
% % Plot format features
% title(strcat('WWLLN Lightning',{' '}, datestr(time(1),'yyyy-mm-dd')), 'fontsize', 12, 'fontname', 'helvetica','fontweight','bold');
% legend([p43 p42 p41 p44],{'lightning at <1000 km','lightning at 1000-2000 km','lightning at >2000 km','lightning during IAR'},'Fontsize',12,'Location','southwest');


% Plot lightning time series (fix for distant lightning mods January 24,
% 2018)
if strcmp(tltng,'on')
    axes(ax3) % Left side, bottom
    %stem(time_axis,Energy,'Marker','none'); % stem plot instead of regular plot
    %plot(ax3,time_axis,Energy); % plot of WWLLN lightning energy vs. time
    p31 = plot(ax3,time_axis(arclen > km2deg(proximal,'earth')),Energy(arclen > km2deg(proximal,'earth')),'k');
    set(gca,'YScale','log'); % change to log scale
    xlim([time_axis(1) ceil(time_axis(end))]); % set xlimits of plot
    set(gca,'XTick',[0 3 6 9 12 15 18 21 24]); % set tick labels
    xlabel(ax3,'Universal Time (h)', 'fontsize', fntsize, 'fontname', 'helvetica'); % label the axes
    ylabel(ax3,'Lightning Energy (J)','fontsize', fntsize, 'fontname', 'helvetica');
    title(strcat('WWLLN Lightning Energy (J)',{' '}, datestr(time(1),'yyyy-mm-dd')), 'fontsize', 12, 'fontname', 'helvetica','fontweight','bold');
    hold on
    
    p32 = plot(ax3,time_axis(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth')),Energy(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth')),'g');
    p33 = plot(ax3,time_axis(arclen < km2deg(local,'earth')),Energy(arclen < km2deg(local,'earth')),'m');
    
    legend([p33 p32 p31],strcat(['lightning at <',num2str(local),' km']),strcat(['lightning at ',num2str(local),'-',num2str(proximal),' km']),strcat(['lightning at >',num2str(proximal),' km']))
end
end