% function lma_channels_12_09_07_00_49_35
clear
close all
% addpath('/home/jeff/Research/20120814/');
addpath('/home/jeff/MATLAB/Research/')
addpath('/home/jeff/Research/20120907/lma/80us/');
addpath('/home/jeff/Documents/MATLAB/StormAnalysis-2.1/');
% addpath('/home/jeff/Research/20120814/LMA/20us/');
% addpath('/home/jeff/Research/20120814/LMA/10us/');
% addpath('/home/jeff/Research/20120814/Interferometer/');

%% Load data 
% LMA
fname = 'LYLOUT_120907_004935_0002_no_header_10_stations_1_chi2.flt'
fid1=fopen(fname);
% fid2=fopen('LYLOUT_120814_214225_0001_12_stations_1_chi2.flt');
% Load 80us data
C1=textscan(fid1,'%n%n%n%n%n%n%n%s',-1);
fclose(fid1)
t1=C1{2};
% flashtime=24.93117+0.08;
lmatime_all=t1-3600*00-60*49;
lat=C1{3};
% size(lat)
lon=C1{4};
alt=C1{5};
% % Load 20us data
% C2=textscan(fid2,'%n%n%n%n%n%n%n%s',-1);
% t2lma=C2{2};
% lat2=C2{3};
% size(lat2)
% lon2=C2{4};
% alt2=C2{5}/1000;

% % Find overlap of two files
% I_patch1 = find(t1lma<t2lma(1),1,'last')
% I_patch2 = find(t1lma>t2lma(end),1)

% % Patch in 20 us data
% t1 = [t1lma(1:I_patch1); t2lma; t1lma(I_patch2:end)];
% lat = [lat1(1:I_patch1); lat2; lat1(I_patch2:end)];
% lon = [lon1(1:I_patch1); lon2; lon1(I_patch2:end)];
% alt = [alt1(1:I_patch1); alt2; alt1(I_patch2:end)];
% size(lat)

% Convert from lat-long to x-y
[x0,y0]=deg2utm(33.97533222,-107.18100389);
[UTMxLMA,UTMyLMA]=deg2utm(lat,lon);
xlma=(UTMxLMA-x0);
ylma=(UTMyLMA-y0);


%% Filter
lmatime = lmatime_all;
% tlimits = [0.5 0.95];
% xlimits = [0 13000];
% ylimits = [-17000 -2000];
% zlimits = [5000 9000];
% % lmatime(i)>tlimits(1) & lmatime(i)<tlimits(2) & 
% size(lmatime)
% for i = 1:length(lmatime)
%     if all(xlma(i)>xlimits(1) & xlma(i)<xlimits(2) & ylma(i)>ylimits(1) & ylma(i)<ylimits(2) & alt(i)>zlimits(1) & alt(i)<zlimits(2))
%     else
% %         disp('Filtered')
%         lmatime(i)=NaN;
%         xlma(i)=NaN;
%         ylma(i)=NaN;
%         alt(i)=NaN;
%         lat(i) = NaN;
%         lon(i) = NaN;
%     end
% end
% % Remove NaNs
% lmatime = lmatime(~any(isnan(lmatime),2),:);
% xlma = xlma(~any(isnan(xlma),2),:);
% ylma = ylma(~any(isnan(ylma),2),:);
% lat = lat(~any(isnan(lat),2),:);
% lon = lon(~any(isnan(lon),2),:);
% alt = alt(~any(isnan(alt),2),:);
% % size(lmatime)
%%

begintimes = [35.855 35.965 36.135];
endtimes = [35.866 35.995 36.21];
for i = 1:3
    %Find EdgeLength up to flash
    Ifirst1 = find(lmatime>0,1,'first');
    lmatime_all(Ifirst1);
    Ilast1 = find(lmatime<begintimes(i),1,'last');
    lmatime_all(Ilast1);
%     cmap = jet(Ilast1-Ifirst1);
%     plot(lmatime(Ifirst1:Ilast1),alt(Ifirst1:Ilast1),'.')
    Data = [];
%     for k = 2:7
%         Data = [Data C1{k}(Ifirst1:Ilast1)];
%     end
%     Data = cell2mat(C1{:}(Ifirst1:Ilast1));
    save('Data','Data')
    lmatime_new = lmatime(Ifirst1:Ilast1)+3600*21+60*42;%+flashtime;
    Data = [lmatime_new lat(Ifirst1:Ilast1) lon(Ifirst1:Ilast1) alt(Ifirst1:Ilast1)];
    fileID = fopen('Data.txt','w');
    for j=1:length(Data)
        fprintf(fileID,'%d %d %d %d\n',Data(j,1),Data(j,2),Data(j,3),Data(j,4));
    end
    fclose(fileID);
    fid2=fopen('Data.txt');
    Data_new = textscan(fid2,'%d %d %d %d\n');
    xlma_new = Data_new{2};
    ylma_new = Data_new{3};
    plot(xlma_new,ylma_new,'.')
    [EdgeList, Vertex1, Vertex2] = PulseGraph ('Data.txt','quiet');
    dist = sqrt((Vertex2(1,:) - Vertex1(1,:)).^2 + (Vertex2(2,:) - Vertex1(2,:)).^2 + (Vertex2(3,:) - Vertex1(3,:)).^2);
    edgelength1 = sum(dist);

    %Find EdgeLength after CC
    Ifirst2 = find(lmatime>0,1,'first');
    % lmatime(Ifirst1)
    Ilast2 = find(lmatime<endtimes(i),1,'last');
    % lmatime(Ilast1)
    Data = [];
    % for k = 2:7
    %     Data = [Data C1{k}(Ifirst1:Ilast1)];
    % end
    % Data = cell2mat(C1{:}(Ifirst1:Ilast1));
    % save('Data','Data')
    lmatime_new = lmatime(Ifirst2:Ilast2)+3600*21+60*42;%+flashtime;
    Data = [lmatime_new lat(Ifirst2:Ilast2) lon(Ifirst2:Ilast2) alt(Ifirst2:Ilast2)];
    fileID = fopen('Data2.txt','w');
    for j=1:length(Data)
        fprintf(fileID,'%d %d %d %d\n',Data(j,1),Data(j,2),Data(j,3),Data(j,4));
    end
    fclose(fileID);
    [EdgeListf, Vertex1f, Vertex2f] = PulseGraph ('Data2.txt','quiet');
    distf = sqrt((Vertex2f(1,:) - Vertex1f(1,:)).^2 + (Vertex2f(2,:) - Vertex1f(2,:)).^2 + (Vertex2f(3,:) - Vertex1f(3,:)).^2);
    edgelength2 = sum(distf);
    edgediff = edgelength2 - edgelength1
    % Plot lines
%     figure(i+1);
%     cmap = jet(length(Vertex1f));
%     patch([Vertex1f(1,:)' Vertex2f(1,:)'],[Vertex1f(2,:)' Vertex2f(2,:)'],cmap)
%     line([Vertex1f(1,:) Vertex2f(1,:)],[Vertex1f(2,:) Vertex2f(2,:)], 'LineWidth', 2,'Color',cmap);
%     ah = axes('Parent', i+1);
%     Vertex1f(1:2,1)
%     Vertex2f(1:2,1)
%     plot(ah, Vertex1f(1:2,1),Vertex2f(1:2,1), 'LineWidth', 2);
% %     return
%     for n = 2:length(Vertex2f)
%         line(ah, Vertex1f(1:2,n),Vertex2f(1:2,n), 'LineWidth', 2);
%     end
%     return
end
return
%% Find distance between points (distance between consecutive points)
% dist = zeros(length(lmatime)-1);
% for j = 2:length(lmatime)
%     dist(j) = sqrt((xlma(j) - xlma(j-1))^2 + (ylma(j) - ylma(j-1))^2 + (alt(j) - alt(j-1))^2);
% end
% Ifirst1 = find(lmatime>0.53,1,'first');
% Ilast1 = find(lmatime<0.545,1,'last');
% dist1 = sum(dist((Ifirst1:Ilast1)))
% 
% Ifirst2 = find(lmatime>0.64,1,'first');
% Ilast2 = find(lmatime<0.655,1,'last');
% dist2 = sum(dist((Ifirst2:Ilast2)))
% 
% Ifirst3 = find(lmatime>0.77,1,'first');
% Ilast3 = find(lmatime<0.785,1,'last');
% dist3 = sum(dist((Ifirst3:Ilast3)))


%% Find branching
% k = 1;
% branch = 0;
% while k<length(lmatime) 
%     dist = sqrt((xlma(k+5) - xlma(k+4))^2 + (ylma(k+5) - ylma(k+4))^2 + (alt(k+5) - alt(k+4))^2);
%     if dist>200
%         branch = branch + 1;
%         k = k+5;
%     end
%         
%     k = k+1;
% end
% branch
%% Find distance between points (distance between closest neighbour)
% begintimes = [0.55 0.66 0.79];
% endtimes = [0.565 0.68 0.88];
% % for m = 1:3
% Ifirst = find(lmatime>begintimes(1),1,'first')
% Ilast = find(lmatime<endtimes(1),1,'last')
% % Ilast-Ifirst
% clear dpairs Ipairs xyz
% for j = Ifirst:Ilast-1
%     dist = zeros(Ilast-j-1,1); %clear dist for every lma point
%     for k = j+1:Ilast
%         dist(k-j+1,1) = sqrt((xlma(k) - xlma(j))^2 + (ylma(k) - ylma(j))^2 + (alt(k) - alt(j))^2);
%     end
%     itr = j-Ifirst+1;
%     [d I] = min(dist);
%     dpairs(itr) = [d];
%     Ipairs(itr,:) = [j I];
%     xyz(itr,:) = [[xlma(j) ylma(j) alt(j)] [xlma(I) ylma(I) alt(I)]];
% end
% sum(dpairs)
% clear dpairs Ipairs xyz
% for j = 1:Ifirst-1
%     dist = zeros(Ifirst-j-1,1); %clear dist for every lma point
%     for k = j+1:Ifirst
%         dist(k-j,1) = sqrt((xlma(k) - xlma(j))^2 + (ylma(k) - ylma(j))^2 + (alt(k) - alt(j))^2);
%     end
%     itr = j;
%     [d I] = min(dist);
%     dpairs(itr) = [d];
%     Ipairs(itr,:) = [j I];
%     xyz(itr,:) = [[xlma(j) ylma(j) alt(j)] [xlma(I) ylma(I) alt(I)]];
% end
% d1 = sum(dpairs);
% clear dpairs Ipairs xyz
% for j = 1:Ilast-1
%     dist = zeros(Ilast-j-1,1); %clear dist for every lma point
%     for k = j+1:Ilast
%         dist(k-j,1) = sqrt((xlma(k) - xlma(j))^2 + (ylma(k) - ylma(j))^2 + (alt(k) - alt(j))^2);
%     end
%     itr = j;
%     [d I] = min(dist);
%     dpairs(itr) = [d];
%     Ipairs(itr,:) = [j I];
%     xyz(itr,:) = [[xlma(j) ylma(j) alt(j)] [xlma(I) ylma(I) alt(I)]];
% end
% sum(dpairs)-d1
% fig = figure;
% ah = axes('Parent', fig);
% plot(ah, xyz(1,[1 4]),xyz(1,[2 5]), 'LineWidth', 2);
% [col row] = size(xyz)
% for n = 2:col
%   line(xyz(n,[1 4]), xyz(n,[2 5]), 'LineWidth', 2);
% end
% return
% end
 
    
    
    