function lma_channels_12_08_14_21_42_26
clear
close all

% addpath('/home/jeff/Research/20120814/');
% addpath('/home/jeff/MATLAB/Research_code/');
% addpath('/home/jeff/Research/20120814/LMA/10us/');
% 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 
flshdate = '20120814';
time = '2142';
% LMA
% fname = 'LYLOUT_120814_214225_0002_10_stations_1_chi2.flt';
fname = 'LYLOUT_120814_214224_0001_12_stations_0.4_chi2.flt';
% fname = 'LYLOUT_120814_214224_0001_12_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=24.95;
t2=26.1;
% [E3,Etime3]=load_fieldchange_20120814('load','lefa3','M');
[E3,Etime3]=load_E('20120814','2142','lefa3','4',t1,t2);
E3=E3/1000;
Edata = E3;

%% Filter
minX = -1;
maxX = 20;
minY = -25;
maxY = 0;
maxZ = 11;
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);
%%
%              %RS-1 BTW-1+2 RS-2  BTW-2+3 RS-3  BTW-3+4 RS-4  BTW-4+5 RS-5  BTW-5+6 RS-6   K-Last                   
% begintimes = [0.366 0.3700 0.5042 0.5070 0.5610 0.5820 0.6710 0.6960 0.8020 0.9520 1.0016 0.9980];
% endtimes =   [0.370 0.5042 0.5070 0.5610 0.5820 0.6710 0.6900 0.8020 0.9520 1.0015 1.0042 1.0003];
% 
% % #           Total   Leader K-1     RS-1   RS-2    RS-2+3  RS-3   BTW-3+4 RS-4   RS-5   RS-5+6 RS-6   K-8     BTW-2+3
% % tminlist = [25.0111,25.300,25.3007,25.366,25.5042,25.5042,25.5610,25.574,25.671,25.802,25.802,26.0016,25.3314,25.5070]
% % tmaxlist = [26.1100,25.365,25.3027,25.370,25.5070,25.5740,25.5820,25.671,25.686,25.952,26.004,26.0042,25.3322,25.5610]
% 
% 
               %RS-1  RS-2   RS-3   RS-4   RS-5   RS-6                      
% RSbegintimes = [0.366 0.5042 0.5610 0.6710 0.8020 1.0016];
% RSendtimes =   [0.370 0.5070 0.5820 0.6900 0.9520 1.0042];
% 
%              %BTW-1+2 BTW-2+3 BTW-3+4 BTW-4+5 BTW-5+6                    
% btwbegintimes = [0.3700  0.5070  0.5820  0.6900  0.9520];
% btwendtimes =   [0.5042  0.5610  0.6710  0.8020  1.0015];
% 
% chan_growth(lmatime,lon,lat,alt,RSbegintimes+25,RSendtimes+25);
% chan_growth(lmatime,lon,lat,alt,RSbegintimes+25,RSbegintimes+25+0.002);
% % % chan_growth(lmatime,lon,lat,alt,btwbegintimes+flashtime,btwendtimes+flasthime);
% return
%% Make Video of Flash
if 1
fid = fopen(['growth_rate_',flshdate,'_',time,'closer','.txt']);
C2 = textscan(fid,'%f%f');
divtime = C2{1};
edgediff = C2{2}/1000;
growthfig = figure('Units','inches','Position',[2,2,14,6]);
minE = -5.5;
maxE=0.2;
minGrow=0;
maxGrow=10;
% tflash1 = 25.5042; %Flash 2
% tflash2 = 25.5070;
% tflash1 = 25.561; %Flash 3
% tflash2 = 25.582;
% tflash1 = 25.671; %Flash 4
% tflash2 = 25.696;
tflash1 = 25.802; %Flash 5
tflash2 = 25.952;
% tflash2 = 26.111;

% for i = 1:floor((divtime(end)-divtime(1))*10)
% startt=25.365;
% tflash1 = startt + 0.1*(i-1)
% tflash2 = tflash1 + 0.1


picname = ['growth_pics_',flshdate,'_',time,'/growth_rate_temp'];

% Plot Growth rate
growax = subplot('Position',[0.06,0.69,0.4,0.27]);
plot(divtime,edgediff,'LineWidth',2,'Color','k')
Iflash1 = find(divtime>tflash1,1,'first');
Iflash2 = find(divtime<tflash2,1,'last');
% hold on
% plot(divtime(Iflash1:Iflash2),edgediff(Iflash1:Iflash2),'LineWidth',2,'Color','r')
% hold off
set(gca,'XGrid','on');
ylim([minGrow maxGrow])
xlim([t1 t2])
% ylim([minGrow,maxGrow])
% ylabel('(km)','FontSize',24)
text(25.37,9,'Chan growth/10 ms (km)','FontSize',30)
% title('Channal growth/10 \mus, LEFA3, and LMA vs Time','FontSize',30)
set(gca,'XTickLabel','')
set(gca,'FontSize',24)
hold on
% growln = plot(growax,[divtime(30) divtime(30)],[minGrow maxGrow],'Color','r','LineWidth',2);

%Plot E-field
Eax = subplot('Position',[0.06,0.42,0.4,0.27]);
plot(Etime3,E3,'LineWidth',2,'Color','k')
Iflash1 = find(Etime3>tflash1,1,'first');
Iflash2 = find(Etime3<tflash2,1,'last');
hold on
plot(Etime3(Iflash1:Iflash2),E3(Iflash1:Iflash2),'LineWidth',2,'Color','r')
hold off
set(gca,'XGrid','on');
xlim([t1 t2])
ylim([minE,maxE])
set(gca,'XTickLabel','')
% ylabel('(kV/m)','FontSize',24)
text(25.55,-0.5,'LEFA3 (kV/m)','FontSize',30);
% xlabel(['Time in seconds after ',time(1:2),':',time(3:4),' UTC'],'FontSize',24)
set(gca,'FontSize',24)
% set(gca,'YTickLabel',{'','-4','-2',''})
% hold on
% Eln = plot(Eax,[divtime(30) divtime(30)],[minE maxE],'Color','r','LineWidth',2);

%Plot LMA
LMAax = subplot('Position',[0.06,0.15,0.4,0.27]);
scatter(lmatime,alt,'.k')
Iflash1 = find(lmatime>tflash1,1,'first');
Iflash2 = find(lmatime<tflash2,1,'last');
hold on
scatter(lmatime(Iflash1:Iflash2),alt(Iflash1:Iflash2),'.r')
hold off
set(gca,'XGrid','on')
box on
set(gca,'FontSize',24)
xlim([t1 t2])
ylim([0,maxZ])
% ylabel('(km)','FontSize',24)
text(25.44,9.9,'LMA Altitude (km)','FontSize',29);
xlabel(['Time in seconds after ',time(1:2),':',time(3:4),' UTC'],'FontSize',24)
% hold on
% LMAln = plot(LMAax,[divtime(30) divtime(30)],[0 maxZ],'Color','r','LineWidth',2);

% xlmafig = figure('Units','inches','Position',[0.2,0.2,3,6]);
% XLMAax = subplot('Position',[0.15,0.5,0.8,0.4]);
% cxlma = jet(length(xlma));
% scatter(xlma,ylma,30,cxlma,'.')
% set(gca,'XTickLabel','')
% ylabel('Latitude (km)','FontSize',24)
% xlim([minX maxX]);
% ylim([minY,maxY]);
% text((maxX-minX)*0.48+minX,(maxY-minY)*0.95+minY,'LMA','FontSize',40);
% set(gca,'YTickLabel',{'','-20','-15','-10','-5','0'})

[EdgeList lmatimetemp xlmatemp ylmatemp alttemp] = growth_pic(time,tflash1,tflash2,minX,maxX,minY,maxY,lmatime,lon,lat,alt,xlma,ylma,divtime);
% lmaedgefname = '20120814-2142-lmaedge_old';
% lmaid1 = fopen([lmaedgefname,'.txt'],'w');
% fprintf(lmaid1,'%f %f %f %f\n',[lmatimetemp,xlmatemp,ylmatemp, alttemp]');
% fclose(lmaid1);
% edgefname = '20120814-2142-edgedata_old';
% lmaid2 = fopen([edgefname,'.txt'],'w');
% fprintf(lmaid2,'%f %f\n',[EdgeList]);
% fclose(lmaid2);
% k = waitforbuttonpress;
% end

% growth_video(flshdate,time,t1,t2,minX,maxX,minY,maxY,maxZ,minE,maxE,minGrow,maxGrow,lmatime,lon,lat,alt,xlma,ylma,divtime,growthfig,growax,Eax,LMAax,growln,Eln,LMAln);

return
end
%% 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(['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)
% plot(divtime,edgediff)
% set(gca,'XGrid','on');
% xlim([t1 t2])
% subplot(2,1,2)
% plot(Etime3,E3)
% set(gca,'XGrid','on');
% xlim([t1 t2])
% saveas(edgefig,['growth_rate_',date','_',time,'.fig'])

% dangfig = figure();
% subplot(2,1,1)
% plot(divtime,dangdiff)
% xlim([25 26.1])
% subplot(2,1,2)
% plot(Etime3,E3)
% xlim([25 26.1])
% saveas(dangfig,'dangle_rate_20120814-214225-100us.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])
    
    
end