function [winltimes, winlnums, winllatlon, winlEnergy] = winlnum(lfilename, siteData, local, proximal, time, lwinsvec, midwinpoints, npt, overlap, peakwint, max_arc_len_c)
ldata = importdata(lfilename); % Load WWLLN "data" cell

Latitude=ldata(:,7); % latitude vector: Lat, Lon in fractional degrees (geographic coordinated)
Longitude=ldata(:,8); % longitude vector

% 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') & arclen < km2deg(max_arc_len_c,'earth'),1:6),'TimeZone','UTC'); % Date/time in UTC, datetime vector
    llatlons = ldata(arclen > km2deg(proximal,'earth') & arclen < km2deg(max_arc_len_c,'earth'),7:8);
    lEnergy = ldata(arclen > km2deg(proximal,'earth') & arclen < km2deg(max_arc_len_c,'earth'),11);
    ltimessize=size(ltimes);
else
    ltimes = datetime(ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth') & arclen < km2deg(max_arc_len_c,'earth'),1:6),'TimeZone','UTC'); % Date/time in UTC, datetime vector
    llatlons = ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth') & arclen < km2deg(max_arc_len_c,'earth'),7:8);
    lEnergy = ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth') & arclen < km2deg(max_arc_len_c,'earth'),11);
    ltimessize=size(ltimes);
    %ltimesl = datetime(ldata(arclen < km2deg(local,'earth'),1:6),'TimeZone','UTC'); % Date/time in UTC, datetime vector
    %ltimesp = datetime(ldata(arclen > km2deg(local,'earth') & arclen < km2deg(proximal,'earth'),1:6),'TimeZone','UTC'); % Date/time in UTC, datetime vector
end

% Convert winstep into points
% winltimes{1} = ltimes(ltimes > time(1) & ltimes < time(npt));
% checkmindiff = 3600*hour(winltimes{1}) + 60*minute(winltimes{1}) + second(winltimes{1});
% checkmindiff = diff(checkmindiff);
% mindiffloc = find(checkmindiff < 0.6 * peakwint);
% if isempty(mindiffloc)
%     winltimes{1} = winltimes{1};
% else
%     mindiffloc2 = mindiffloc+1;
%     winltimes{1}([mindiffloc, mindiffloc2]) = [];
% end
% winlnums(1) = length(winltimes{1}); %FIXME?
%winltimesp(1) = length(ltimesp(ltimesp > time(1) & ltimesp < time(npt)));
for i = 1:length(midwinpoints)
    [winltimes{i}, I] = sort(ltimes(ltimes > time(lwinsvec(i)) & ltimes < time(lwinsvec(i+1))));
    winllatlon{i} = llatlons(ltimes > time(lwinsvec(i)) & ltimes < time(lwinsvec(i+1)),:);
    winllatlon{i} = winllatlon{i}(I,:);
    winlEnergy{i} = lEnergy(ltimes > time(lwinsvec(i)) & ltimes < time(lwinsvec(i+1)),:);
    winlEnergy{i} = winlEnergy{i}(I);
    winlnums(i) = length(winltimes{i}); %FIXME?
%     checkmindiffa = 3600*hour(winltimes{i}(1:end-1)) + 60*minute(winltimes{i}(1:end-1)) + second(winltimes{i}(1:end-1));
%     if length(checkmindiffa) < 7
%         continue
%     end
%     checkmindiffb = 3600*hour(winltimes{i}(2:end)) + 60*minute(winltimes{i}(2:end)) + second(winltimes{i}(2:end));
%     [arclenb,azb] = distance(STN_Lat,STN_Long,winllatlon{i}(:,1),winllatlon{i}(:,2));
%     checkmindiffa_dist = deg2km(arclenb(1:end-1));
%     checkmindiffb_dist = deg2km(arclenb(2:end));
%     diffab = checkmindiffb - checkmindiffa';
%     diffab_dist = checkmindiffb_dist - checkmindiffa_dist';
%     figure;
%     subplot(2,1,1);
%     dscatter(diffab(diffab>0 & diffab_dist>0),diffab_dist(diffab>0 & diffab_dist>0),'PLOTTYPE','contour');
%     subplot(2,1,2);
%     dscatter(diffab(diffab>0 & diffab_dist>0),diffab_dist(diffab>0 & diffab_dist>0),'PLOTTYPE','surf');
%     view([0 0 1]);
%     shading flat;
% %     checkmindiff = diff(checkmindiff);
% %     mindiffloc = find(checkmindiff < 0.4 * peakwint);
% %     if isempty(mindiffloc)
% %         winltimes{i} = winltimes{i};
% %         winllatlon{i} = winllatlon{i};
% %     else
% %         %     mindiffloc2 = mindiffloc+1;
% %         %     winltimes{i}([mindiffloc, mindiffloc2]) = [];
% %         winltimes{i}(mindiffloc+1) = [];
% %         winllatlon{i}(mindiffloc+1,:) = [];
% %         mindifflocsize1(i) = length(mindiffloc);
% %     end
% %     
% %     checkmindiff = 3600*hour(winltimes{i}) + 60*minute(winltimes{i}) + second(winltimes{i});
% %     checkmindiff = diff(checkmindiff);
% %     mindiffloc = find(checkmindiff < 0.6 * peakwint);
% %     if isempty(mindiffloc)
% %         winltimes{i} = winltimes{i};
% %         winllatlon{i} = winllatlon{i};
% %     else
% %         %     mindiffloc2 = mindiffloc+1;
% %         %     winltimes{i}([mindiffloc, mindiffloc2]) = [];
% %         winltimes{i}(mindiffloc) = [];
% %         winllatlon{i}(mindiffloc+1,:) = [];
% %         mindifflocsize2(i) = length(mindiffloc);
% %     end
    
    
end
% ltimessize(1)-sum(mindifflocsize1)-sum(mindifflocsize2)
% sum(winlnums)

% Add capability to differentiate between local, between local and
% proximal, and proximal

end