% function lma_channels_13_07_08_20_33_02
clear
close all
addpath('/home/jeff/Research/20130708/');
addpath('/home/jeff/MATLAB/Research_code/')
addpath('/home/jeff/Documents/MATLAB/StormAnalysis-2.1/');

%% Load data 
flshdate = '20130708';
time = '2033';
% LMA
fname = 'LYLOUT_130708_203302_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};
lmatime=t1-3600*str2double(time(1:2))-60*str2double(time(3:4));
lat=C1{3};
% size(lat)
lon=C1{4};
alt=C1{5}/1000;
% % 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)/1000;
ylma=(UTMyLMA-y0)/1000;

%Load E-field data
t1=2.3;
t2=3.5;
[E3,Etime3]=load_E('20130708','2033','lefa3','5',t1,t2);
E3=E3/1000;
Edata = E3;

%% Filter
minX = -15;
maxX = 2;
minY = -4;
maxY = 5;
maxZ = 9;
idx=find((xlma>minX & xlma<maxX) & (ylma>minY & ylma<maxY) & alt<maxZ);
lmatime=lmatime(idx);
lon=lon(idx);
lat=lat(idx);
alt=alt(idx);
xlma=xlma(idx);
ylma=ylma(idx);
%%
%             1st RS, 
RSbegintimes = [2.4751,2.6109,2.7202,2.7971,2.884,2.9485,3.042,3.1254,3.2416,3.3339];
RSendtimes =   [2.4770,2.6120,2.7230,2.8800,2.892,3.0050,3.043,3.2200,3.2450,3.3900];

chan_growth(lmatime,lon,lat,alt,RSbegintimes,RSendtimes);

% chan_growth(lmatime,lon,lat,alt,RSbegintimes-0.01,RSbegintimes);
return
%% Make Video of Flash
fid = fopen(strcat('growth_rate_',flshdate,'_',time,'.txt'));
C2 = textscan(fid,'%f%f');
divtime = C2{1};
edgediff = C2{2};
growthfig = figure('Units','inches','Position',[2,2,18,8]);
minE = -50;
maxE=0;
minGrow=0;
maxGrow=15;

%Plot Growth rate
growax = subplot('Position',[0.1,0.5,0.4,0.4]);
plot(divtime,edgediff,'LineWidth',2)
set(gca,'XGrid','on');
xlim([t1 t2])
ylim([minGrow,maxGrow])
ylabel('Channel Growth (km)','FontSize',24)
set(gca,'XTickLabel','')
set(gca,'FontSize',24)
hold on
growln = plot(growax,[divtime(1) divtime(1)],[minGrow maxGrow],'Color','r','LineWidth',2);

%Plot E-field
Eax = subplot('Position',[0.1,0.1,0.4,0.4]);
plot(Etime3,E3,'LineWidth',2)
set(gca,'XGrid','on');
xlim([t1 t2])
ylim([minE,maxE])
ylabel('E-field (kV/m)','FontSize',24)
xlabel(['Time in seconds after ',time(1:2),':',time(3:4),' UTC'],'FontSize',24)
set(gca,'FontSize',24)
% set(gca,'YTickLabel',{'-4','-2','0',''})
% Eln = line([divtime(1) divtime(1)],[min(E3) max(E3)],'Color','r');
hold on
Eln = plot(Eax,[divtime(1) divtime(1)],[minE maxE],'Color','r','LineWidth',2);

%Plot Pulse Map
pulseax = subplot('Position',[0.56,0.1,0.39,0.8]);
xlabel('Longitude (km)','FontSize',24)
ylabel('Latitude (km)','FontSize',24)
set(gca,'FontSize',24)
hold on
EdgeList_old = [];
m=1;
for i = 2:4%length(divtime)
    %Select new data
    Ifirst = find(lmatime>0,1,'first');
    Ilast = find(lmatime<divtime(i),1,'last');
    lmatimetemp=lmatime(Ifirst:Ilast);
    lontemp=lon(Ifirst:Ilast);
    lattemp=lat(Ifirst:Ilast);
    xlmatemp=xlma(Ifirst:Ilast);
    ylmatemp=ylma(Ifirst:Ilast);
    alttemp=alt(Ifirst:Ilast);
    lmatimeSAM = lmatimetemp + 3600*str2double(time(1:2))-60*str2double(time(3:4));
    Data = [lmatimeSAM lattemp lontemp alttemp];
    [EdgeList, ~, ~] = PulseGraph_new ('input',Data);
    %Plot new data
    nedges = length(EdgeList);
    for j = 1:nedges
        k = EdgeList (1, j) ;
        l = EdgeList (2, j) ;
        x1 = xlmatemp (k) ;
        x2 = xlmatemp (l) ;
        X = [ x1 x2 ] ;
        y1 = ylmatemp (k) ;
        y2 = ylmatemp (l) ;
        Y = [ y1 y2 ] ;
        plot (pulseax,X, Y, '-r', 'LineWidth', 1) ;
    end
    %Select old data
    Ifirst = find(lmatime>0,1,'first');
    Ilast = find(lmatime<divtime(i-1),1,'last');
    lmatimetemp=lmatime(Ifirst:Ilast);
    lontemp=lon(Ifirst:Ilast);
    lattemp=lat(Ifirst:Ilast);
    xlmatemp=xlma(Ifirst:Ilast);
    ylmatemp=ylma(Ifirst:Ilast);
    alttemp=alt(Ifirst:Ilast);
    nedges = length(EdgeList_old);
    %Plot old data
    for j = 1:nedges
        k = EdgeList_old (1, j) ;
        l = EdgeList_old (2, j) ;
        x1 = xlma (k) ;
        x2 = xlma (l) ;
        X = [ x1 x2 ] ;
        y1 = ylma (k) ;
        y2 = ylma (l) ;
        Y = [ y1 y2 ] ;
        plot (pulseax,X, Y, '-k', 'LineWidth', 1) ;
    end
    xlim([minX maxX]);
    ylim([minY,maxY]);
    EdgeList_old = EdgeList;
    %Delete and create new lines
    delete(growln);
    delete(Eln);
    growln = plot(growax,[divtime(i) divtime(i)],[minGrow maxGrow],'Color','r','LineWidth',2);
    Eln = plot(Eax,[divtime(i) divtime(i)],[minE maxE],'Color','r','LineWidth',2);
    savefig(['growth_pics_',flshdate,'_',time,'/growth_rate_temp',num2str(m)],growthfig,'png')
    im = imread(['growth_pics_',flshdate,'_',time,'/growth_rate_temp',num2str(m),'.png']);
    if m == 1
        [mov(:,:,1,m), map] = rgb2ind(im, 256, 'nodither');
    else
        mov(:,:,1,m) = rgb2ind(im, map, 'nodither');
    end
    m=m+1;
end
hold off
%Write gif
% imwrite(mov, map, strcat('growth_rate_',date,'_',time,'.gif'), 'DelayTime', 0.5, 'LoopCount', inf);
%Write avi
makevideo(flshdate,time)
return

%% Find channel growth vs time
Ifirst = find(lmatime>t1,1,'first');
Ilast = find(lmatime<t2,1,'last');

[divtime edgediff] = chan_growth_rate(time,lmatime(Ifirst:Ilast),lon(Ifirst:Ilast),lat(Ifirst:Ilast),alt(Ifirst:Ilast));

fid = fopen(strcat('growth_rate_',flshdate,'_',time,'.txt'),'w');
y = [divtime; edgediff];
fprintf(fid,'%12.8f %12.8f\n',y);
fclose(fid);

edgefig = figure();
subplot(2,1,1)
xlim([t1 t2])
plot(divtime,edgediff)
set(gca,'XGrid','on');
title(strcat('Channel growth rate 20130708 - ',time));
subplot(2,1,2)
xlim([t1 t2])
plot(Etime3,E3)
set(gca,'XGrid','on');
saveas(edgefig,strcat('growth_rate_20130708-',time,'.fig'))
% 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.355 0.493 0.55 0.66 0.79 0.9904 0.8935];
% endtimes = [0.359 0.4945 0.565 0.68 0.88 0.993 0.8955];
% mindist = 100;
% % for m = 1:length(begintimes)
% div = 400;
% totalt = 1.1;
% tgrowth = totalt / div
% dpairs = [];
% for m = 1:div
%     m
% %     Ifirst = find(lmatime>begintimes(m) +lmatime(1),1,'first');
% %     Ilast = find(lmatime<endtimes(m) + lmatime(1),1,'last');
%     Ifirst = find(lmatime>lmatime(1) + tgrowth*(m-1),1,'first');
%     Ilast = find(lmatime>lmatime(1) + tgrowth*(m),1,'first');
%     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);
%         if d > mindist
% %             disp('No pair')
%         else
%             dpairs(itr) = [d];
% %             Ipairs(itr,:) = [j I];
% %             xyz(itr,:) = [[xlma(j) ylma(j) alt(j)] [xlma(I) ylma(I) alt(I)]];
%         end
%     end
%     d1 = sum(dpairs);
%     dpairs = [];
%     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);
%         if d > mindist
%     %         disp('No pair')
%         else
%             dpairs(itr) = [d];
% %             Ipairs(itr,:) = [j I];
% %             xyz(itr,:) = [[xlma(j) ylma(j) alt(j)] [xlma(I) ylma(I) alt(I)]];
%         end
%     end
%     d2 = sum(dpairs);
%     ddiff(m) = d2-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
% 
% end
% divtime = linspace(0,totalt,div) + lmatime(1);
% figure
% subplot(2,1,1)
% plot(divtime,ddiff)
% subplot(2,1,2)
% plot(Etime3,E3)
% xlim([25 26.2])
    
    
    