function Process_Current_Speed_RunAver_3;

clear all;
set(0,'DefaultFigureColor','w',...
      'DefaultAxesColor','w',...
      'DefaultAxesXColor','k',...
      'DefaultAxesYColor','k',...
      'DefaultAxesZColor','k',...
      'DefaultTextColor','k',...
      'DefaultLineColor','k');
close all;

sourcePath  = 'E:\707756POSBuoyStudy\EnvironmentData\NDBC_41012\current\';
fname{1}    = '4101212008.txt';
fname{2}    = '4101232008.txt';
fname{3}    = '4101242008.txt';
fname{4}    = '4101252008.txt';
fname{5}    = '4101262008.txt';
imagePath   = sourcePath;

%fileselect = [2 3 4 5];      %-- if gluing all data sets together
                              %-- January cannot be glued because of
                              %-- shorter record (20 dys), and a hole in February
                              
fileselect = [5];             %-- If selecting a single data set                        
deg6     = [];
win6     = [];
deg20    = [];
win20    = [];
time     = [];
for ifile = fileselect
   name = [sourcePath fname{ifile}]
   p2 = findstr('.', name);
   p1 = findstr('\', name);
   p1 = p1(length(p1));
   title = name(p1+1:p2-1);
   aa = importdata(name);
   if(isstruct(aa))
      header1 = aa.textdata{1};
      header2 = aa.textdata{2};
      dd = aa.data;
   else
      dd = aa;
   end;
   %----- unpack the imported data
   if fileselect == 1       %-- NOTE: January field has different format
      deg6read    = dd(:, 7);
      win6read    = dd(:, 8);
      deg20read   = dd(:, 28);
      win20read   = dd(:, 29);
   else
      deg6read    = dd(:, 10);
      win6read    = dd(:, 11);
      deg20read   = dd(:, 31);
      win20read   = dd(:, 32);
   end;
   timeread    = datenum(dd(:, 1), dd(:, 2), dd(:, 3), dd(:, 4),dd(:, 5), 0);
   timeread    = timeread - datenum(dd(1,1), 1, 1);     %-- Days from the beginning of the year
   %----- stuff into a continuous data stream
   deg6        = [deg6, deg6read'];
   deg20       = [deg20, deg20read'];
   win6        = [win6, win6read'];
   win20       = [win20, win20read'];
   time        = [time, timeread'];
end;

deg6        = deg6';
deg20       = deg20';
win6        = win6';
win20       = win20';
time        = time'; 
len         = length(time);
%----- Convert to vectors------------
u6    = win6 .* sin(pi*deg6/180);
v6    = win6 .* cos(pi*deg6/180);
u20   = win20 .* sin(pi*deg20/180);
v20   = win20 .* cos(pi*deg20/180);   

%----- Data Processing ----------------------------------------------------
%----- Select a set of symmetric RunAver windows
averWindow = [1,3,5,7,9,11, 13,15,17,19,21,25,31,41,51,61,71,81, 91, 121, 141];
averTime   = averWindow * (time(length(time)) - time(1))/len;

load jetplus3;
colormap(jetplus3);
noaver = length(averWindow);
for run = 1: noaver
   cc(run, 1:3) = jetplus3(4*run, 1:3);
end;
sample_limit = round((averWindow(length(averWindow))-1)/2);
istart = sample_limit;
iend   = len - sample_limit;

day_start = round(time(istart));
day_end   = round(time(iend));
title       = ['2008 Buoy data, days ', num2str(day_start), '-', num2str(day_end)];
imageTitle  = ['2008_Buoy_data_days_', num2str(day_start), '-', num2str(day_end) '_'];
imageName   = [imagePath, imageTitle]; 


%---- Process speed data for different lenght of the averaging windows
for run = 1: length(averWindow)
   ww = averWindow(run);
   w = round((ww - 1)/2);
   
   %----- Assess the running averages ---------
   for i =1:len
      i1 = i - w;
      i2 = i + w;
      if i1 < 1, i1 = 1; end;
      if i2 > len, i2 = len; end;
      
      u6aver = mean(u6(i1:i2));
      v6aver = mean(v6(i1:i2));
      w6aver = sqrt(u6aver .^2 + v6aver.^2);
      u6averdata(run, i) = u6aver;
      v6averdata(run, i) = v6aver;
      w6averdata(run, i) = w6aver;
      
      u20aver = mean(u20(i1:i2));
      v20aver = mean(v20(i1:i2));
      w20aver = sqrt(u20aver .^2 + v20aver.^2);
      u20averdata(run, i) = u20aver;
      v20averdata(run, i) = v20aver;
      w20averdata(run, i) = w20aver;
   end;   
end;

%----- Limit data to a window where all averages are well defined
time = time(istart: iend);
win6 = win6(istart: iend);
win20 = win20(istart: iend);
w6averdata  = w6averdata(:, istart: iend);
w20averdata = w20averdata(:, istart: iend);
len = length(time);

%---- make 2 3D plot witj all running averages
figav = figure('name', 'Running averages plot', 'position',[0,140,600,1000]);
set(figav, 'PaperPositionMode', 'auto');          
set(figav,'units','pixels', 'color', 'w');
up_panel = uipanel('position',[0.02, 0.52, 0.96, 0.46], 'parent', figav);
down_panel = uipanel('position',[0.02, 0.02, 0.96, 0.46], 'parent', figav);

val = zeros(noaver, len);
val(1, 1:len) = win6(1:len);
val(2:noaver+1, 1:len) = w6averdata(1:noaver, 1:len);
taver = [0, averTime];
[X, Y] = meshgrid(taver, time);
h = surf(X, Y, val');
set(h, 'edgecolor', 'none');
shading interp;
az = 66;
el = 32; 
view(70, 20);
ax = get(gca, 'xlabel');
set(ax, 'string', 'Aver. Time','fontsize', 12);
ay = get(gca, 'ylabel');
set(ay, 'string', 'Time (Year-day)','fontsize', 12);
az = get(gca, 'zlabel');
set(az, 'string', 'Current (cm/s)','fontsize', 12);
at = get(gca, 'title');
set(at, 'string', [title ' at 6m'], 'color', 'k', 'fontsize', 12);
set(gca, 'parent', up_panel);
load jetplus3;
colormap(jetplus3);

ff = figure;
set(ff, 'PaperPositionMode', 'auto');          
set(ff,'units','pixels', 'color', 'w');
val = zeros(noaver, len);
val(1, 1:len) = win20(1:len);
val(2:noaver+1, 1:len) = w20averdata(1:noaver, 1:len);
taver = [0, averTime];
[X, Y] = meshgrid(taver, time);
h = surf(X, Y, val');
set(h, 'edgecolor', 'none');
shading interp;
az = 66;
el = 32; 
view(70, 20);
ax = get(gca, 'xlabel');
set(ax, 'string', 'Aver. Time','fontsize', 12);
ay = get(gca, 'ylabel');
set(ay, 'string', 'Time (Year-day)','fontsize', 12);
az = get(gca, 'zlabel');
set(az, 'string', 'Current (cm/s)','fontsize', 12);
at = get(gca, 'title');
set(at, 'string', [title ' at 20m'], 'color', 'k', 'fontsize', 12);
set(gca, 'parent', down_panel);

close gcf;

eval(['print -r100 -djpeg ' imageName '3Dplots',  '_r1']);
eval(['print -r300 -djpeg ' imageName '3Dplots'  '_r3']);



%----- Process statistics
for run = 1: length(averWindow)
   dd = squeeze(w6averdata(run, :));
   dc = win6;
   dc = dc(averWindow:len);
   dd = dd(averWindow:len);
   diff = dc - dd'; 
   %diff = -diff;
   rat = dd' ./ (dc + 0.0001);
   diff_mean6(run) = mean(diff);
   diff_std6(run) = std(diff);
   ratio_mean6(run) = mean(rat);
   ratio_std6(run) = std(rat);
   cor      = corrcoef(dd', dc);
   cor6(run) = cor(1,2);
   
   dd = squeeze(w20averdata(run, :));
   dc = win20;
   dc = dc(averWindow:len);
   dd = dd(averWindow:len);
   diff = dc - dd'; 
   %diff = -diff;
   rat = dd' ./ (dc + 0.0001);
   
   diff_mean20(run) = mean(diff);
   diff_std20(run) = std(diff);
   ratio_mean20(run) = mean(rat); 
   ratio_std20(run) = std(rat);
   cor      = corrcoef(dd', dc);
   cor20(run) = cor(1,2);
end;

%------ Setup figure with two panels ----
%------ Ratio plots
f = figure('name', 'Ratio plots', 'position', [150,100,600, 1000]);
set(f, 'PaperPositionMode', 'auto');          
set(f,'units','pixels', 'color', 'w');
up_panel = uipanel('position',[0.02, 0.52, 0.96, 0.46], 'parent', f);
down_panel = uipanel('position',[0.02, 0.02, 0.96, 0.46], 'parent', f);

%----- Make the first plot
plot_up = plot(averTime, ratio_mean6, 'o', 'color', 'k');
hold on;
plot_up = errorbar(averTime, ratio_mean6, ratio_std6, '-', 'color', 'r');
ax = get(gca, 'xlabel');
set(ax, 'string', 'Averaging time (days)','fontsize', 12);
ay = get(gca, 'ylabel');
set(ay, 'string', 'Curr. Aver./Current Instant.','fontsize', 12);
at = get(gca, 'title');
set(at, 'string', [title ' at 6m; Averaged/Instantaneous'], 'color', 'k', 'fontsize', 12);
set(gca, 'parent', up_panel);

figure;  %-- temporary figure
plot_down = plot(averTime, ratio_mean20, 'o', 'color', 'k');
hold on;
plot_up = errorbar(averTime, ratio_mean20, ratio_std20, '-', 'color', 'r');
ax = get(gca, 'xlabel');
set(ax, 'string', 'Averaging time (days)','fontsize', 12);
ay = get(gca, 'ylabel');
set(ay, 'string', 'Curr. Aver./Current Instant.','fontsize', 12);
at = get(gca, 'title');
set(at, 'string', [title ' at 20m; Averaged/Instant.'], 'color', 'k', 'fontsize', 12);
set(gca, 'parent', down_panel);  %-- move plot to lower lanel
close gcf;  %-- close temporary figure
eval(['print -r100 -djpeg ' imageName 'ratio_error'  '_r1']);
eval(['print -r300 -djpeg ' imageName 'ratio_error'  '_r3']);


%------ Setup figure with two panels ----
%------ Difference plots
f = figure('name', 'Difference plots', 'position', [300,70, 600, 1000]);
set(f, 'PaperPositionMode', 'auto');          
set(f,'units','pixels', 'color', 'w');
up_panel = uipanel('position',[0.02, 0.52, 0.96, 0.46], 'parent', f);
down_panel = uipanel('position',[0.02, 0.02, 0.96, 0.46], 'parent', f);

%----- Make the first plot
plot_up = plot(averTime, diff_mean6, 'o', 'color', 'k');
hold on;
plot_up = errorbar(averTime, diff_mean6, diff_std6, '-', 'color', 'r');
ax = get(gca, 'xlabel');
set(ax, 'string', 'Averaging time (days)','fontsize', 12);
ay = get(gca, 'ylabel');
set(ay, 'string', 'Curr. Aver. - Curr. Inst. Decorrelation','fontsize', 12);
at = get(gca, 'title');
set(at, 'string', [title ' at 6m; Instantaneous - Averaged'], 'color', 'k', 'fontsize', 12);
set(gca, 'parent', up_panel);

figure;     %-- temporary figure
plot_down = plot(averTime, diff_mean20, 'o', 'color', 'k');
hold on;
plot_up = errorbar(averTime, diff_mean20, diff_std20, '-', 'color', 'r');
ax = get(gca, 'xlabel');
set(ax, 'string', 'Averaging time (days)','fontsize', 12);
ay = get(gca, 'ylabel');
set(ay, 'string', 'Curr. Aver. - Curr. Inst. Decorrelation','fontsize', 12);
at = get(gca, 'title');
set(at, 'string', [title ' at 20m; Instantaneous - Averaged'], 'color', 'k', 'fontsize', 12);
set(gca, 'parent', down_panel);   %-- move plot to lower lanel
close gcf;  %-- close temporary figure
eval(['print -r100 -djpeg ' imageName 'diff_error'  '_r1']);
eval(['print -r300 -djpeg ' imageName 'diff_error'  '_r3']);

                                         %-- close temporary figure

%---- Decorrelation plot --------------------
f = figure('name', 'Decorrelation plots', 'position', [450, 40, 600, 1000]);
set(f, 'PaperPositionMode', 'auto');          
set(f,'units','pixels', 'color', 'w');
up_panel = uipanel('position',[0.02, 0.52, 0.96, 0.46], 'parent', f);
down_panel = uipanel('position',[0.02, 0.02, 0.96, 0.46], 'parent', f);

%----- Make the first plot
plot_up = plot(averTime, cor6, 'o', 'color', 'k');
set(ax, 'string', 'Averaging time (days)','fontsize', 12);
ay = get(gca, 'ylabel');
set(ay, 'string', 'Correlation coeff','fontsize', 12);
at = get(gca, 'title');
set(at, 'string', [title ' at 6m; Instantaneous vs Averaged Decorrelation'], 'color', 'k', 'fontsize', 12);
set(gca, 'parent', up_panel);

figure;  %-- temporary figure                                        
plot_down = plot(averTime, cor20, 'o', 'color', 'k');
ax = get(gca, 'xlabel');
set(ax, 'string', 'Averaging time (days)','fontsize', 12);
ay = get(gca, 'ylabel');
set(ay, 'string', 'Correlation coeff','fontsize', 12);
at = get(gca, 'title');
set(at, 'string', [title ' at 20m; Instantaneous vs Averaged Decorrelation'], 'color', 'k', 'fontsize', 12);
set(gca, 'parent', down_panel);  %-- move plot to lower lanel
close gcf;   %-- close temporary figure                                   
colormap(jetplus3);

eval(['print -r100 -djpeg ' imageName 'corr1'  '_r1']);
eval(['print -r300 -djpeg ' imageName 'corr1'  '_r3']);


i = 1;