function DIF_struct = computeDIF( bldg_struct, HDF_path, out_opt_struct )

% Input arguments:

% bldg_struct           a structure array containing 
% HDF_path
%   
clr = get(0,'DefaultUicontrolBackgroundColor'); % Background color for figures and dialog boxes

% if input arguments are not specified, prompt user to select:
% Building input data:
if ~exist('bldg_struct') || isempty(bldg_struct)
    bldg_struct = readbldginfo( );
    if isfield(bldg_struct,'error')
        DIF_struct.error = 1; return;
    end
end
% Location of HDF files:
if ~exist('HDF_path') || isempty(HDF_path)
    % Load pathname for HDF files used in last analysis, if available, otherwise create empty string:
    if exist('tmp_HDF_path.mat','file'), load('tmp_HDF_path'); end
    if ~exist('HDF_path','var') || ~isstr(HDF_path), HDF_path = ''; end
    % Prompt use to select pressure database file:
    [HDF_filename HDF_path] = uigetfile({'*.hdf','HDF files (*.hdf)'},'Select an HDF pressure database file in the directory of interest',HDF_path);
    if any(HDF_filename~=0)
        save('tmp_HDF_path','HDF_path');
    else
        err = errordlg('HDF pressure database file selection cancelled', 'Analysis cancelled');
        uiwait(err); DIF_struct.error = 1; return;
    end
end
% Output options:
field_names = {'ts_plot','obs_pks','est_pks','mult_pks','pk_loads'};
if ~exist('out_opt_struct','var') || ~isstruct(out_opt_struct) || ~isempty(setdiff(field_names,fieldnames(out_opt_struct)))
    % Open dialog box to prompt user to select:
    [out_opt_struct, ok] = output_options;
    if ok~=1
        err = errordlg('Output options selection cancelled', 'Analysis cancelled');
        uiwait(err); DIF_struct.error = 1; return;
    end
end


% Extract building model information from 'bldg_struct':
d0 = bldg_struct.d0;
W0 = d0(1); L0 = d0(2); H0 = d0(3); R0 = d0(4);
length_units = bldg_struct.length_units;
attach_pts = bldg_struct.attach_pts;
frame_coords = bldg_struct.frame_coords;
n_f = size(frame_coords,1); % number of frames
resp_names = bldg_struct.resp_names; % descriptive names of responses of interest
resp_units = bldg_struct.resp_units; % units of responses of interest
N = bldg_struct.inf_coeff; % influence coefficients
[n_a, n_r] = size(N); % number of attachment points, number of response quantities
% Frame labels, for time series plots:
for i = 1:n_f
    frame_label{i} = ['Frame ' num2str(i) ': y = ' num2str(frame_coords(i,2)) ' ' length_units ];
end
% Response labels, for time series plots:
for i = 1:n_r
    resp_label{i} = ['Response ' num2str(i) ': ' resp_names{i}];
end

% Obtain list of which files in the pressure database folder to process:
dir_listing = dir(HDF_path);
file_list = dir_listing(find(cat(1,dir_listing.isdir)==0)); % eliminate directories from list
dir_ind = 1; % initialize index for counting the number of wind directions
for i = 1:length(file_list)
    filename = file_list(i).name; 
    hdf_ext_ind = strfind( filename, '.HDF' );
    if ~isempty(hdf_ext_ind)
        hdf_filename_root(dir_ind,:) = filename(1:hdf_ext_ind-5);
        theta_a_str(dir_ind,:) = filename(hdf_ext_ind-4:hdf_ext_ind-1);
        theta_a(dir_ind,1) = str2num( theta_a_str(dir_ind,:) )/10;
        dir_ind = dir_ind+1;
    end
end
hdf_filename_root = unique(hdf_filename_root,'rows');
if size(hdf_filename_root,1)~=1
    errordlg(['All HDF file names in the selected directory must be identical '...
        'except for the wind direction, specified by the last 4 characters before the file extension.'], ...
        'Error in pressure database folder contents');
    DIF_struct.error = 1; return;
end

% Display list of wind directions available in selected directory, allow user to select directions to use in analysis:
theta_a_cell = mat2cell(theta_a_str(:,1:3),ones(size(theta_a)));
[theta_a_sel,ok] = listdlg('ListString',theta_a_cell, 'InitialValue', [1:length(theta_a_cell)], ...
    'PromptString', {'Select which wind directions to','use in the analysis from the','following list for which pressure','databases are available:'}, ...
    'Name','Wind direction selection.');
if ok
    theta_a = theta_a(theta_a_sel); theta_a_str = theta_a_str(theta_a_sel,:); p = length(theta_a);
else
    errordlg('Wind direction selection cancelled by user.','Analysis cancelled');
    DIF_struct.error = 1; return;
end

theta_a(find(theta_a==360))=0; % ensure that >=0 & <360
theta_b = 360 - theta_a;
theta_b(find(theta_b==360))=0; % ensure that >=0 & <360
theta_c = 180 - theta_a;
theta_c(find(theta_c<0))=theta_c(find(theta_c<0))+360; % make >=0 & <360
theta_d = 180 + theta_a;
theta_d(find(theta_d>=360))=theta_d(find(theta_d>=360))-360; % make >=0 & <360

% Load pressure database model information from the first file in the list:
HDF_struct = load_hdf_coords( fullfile( HDF_path, [hdf_filename_root theta_a_str(1,:) '.HDF'] ) );

% extract full-scale dimensions of building:
W = HDF_struct.Building_Width_Ft; L = HDF_struct.Building_Length_Ft; H = HDF_struct.Building_Height_Ft;
R = HDF_struct.Roof_Slope_in_12*W/2; 
d = [W L H R];
if HDF_struct.A_Tap_Coordinates_3D_Units ~= 'FS_Feet'
    errordlg(['Unrecognized length units in HDF pressure database file:  >> ' HDF_struct.A_Tap_Coordinates_3D_Units ' << ("FS_Feet" expected)'],...
        'Error reading HDF file');
    DIF_struct.error = 1; return;
elseif ~strcmp(bldg_struct.terrain,'Directional') && ~strcmp(HDF_struct.Exposure_Name, bldg_struct.terrain)
    errordlg(['The terrain conditions specified in the HDF file ("' HDF_struct.Exposure_Name '") do not match the terrain' ...
        ' specified in the building input file ("' bldg_struct.terrain '").'], 'Error: terrain mismatch');
    DIF_struct.error = 1; return;
end
DIF_struct.terrain = HDF_struct.Exposure_Name; % 'Open_Country' or 'Suburban'

% store quantities for output in a structure array:
DIF_struct.resp_names = resp_names;
DIF_struct.resp_units = bldg_struct.resp_units;
% output dimensions in same units as building input file
DIF_struct.length_units = bldg_struct.length_units; 
DIF_struct.d0 = bldg_struct.d0;
DIF_struct.d = d;
DIF_struct.frame_coords = bldg_struct.frame_coords;
DIF_struct.attach_pts = bldg_struct.attach_pts;
DIF_struct.force_units = bldg_struct.force_units;
DIF_struct.ws_units = bldg_struct.ws_units;


% tol = 10^-3;
% % Compare dimensions in HDF file with dimensions in building input file and prompt user to confirm scaling of tap coordinates:
% button = 'Yes';
% if abs(H-H0)/H>tol || abs(W-W0)/W>tol || abs(L-L0)/L>tol || abs(R-R0)/R>tol
%     button = questdlg(['The building dimensions in the pressure database file ' ...
%         '(W=' num2str(W) ', L=' num2str(L) ', H=' num2str(H) ',  R=' num2str(R) ') ' ...
%         'do not match the dimensions specified in the building input file ' ...
%         '(W=' num2str(W0) ', L=' num2str(L0) ', H=' num2str(H0) ',  R=' num2str(R0) ') ' ...
%         'Proceed with scaling the tap coordinates to match the dimensions ' ...
%         'in the building input file?'], 'Mismatch of building dimensions: scale tap coordinates?','Yes','Cancel','Yes');
% end
% if ~strcmp(button,'Yes')
%     error('Mismatch of building dimensions: analysis cancelled by user.');
% end

% Assemble tributary matrices for specified frames:
Astruct = tributarymatrix( HDF_struct, frame_coords, attach_pts, 'plot_on', 0, ...
    'H0', H0, 'W0', W0, 'L0', L0, 'R0', R0);

ro = 0.0023769; % Air density (lb*s^2/ft^4) at sea level, 15 degrees C

% Initialize arrays for storing peak response quantities:
X_max_obs_a = zeros(n_r,p,n_f);
X_max_obs_b = zeros(n_r,p,n_f);
X_max_obs_c = zeros(n_r,p,n_f);
X_max_obs_d = zeros(n_r,p,n_f);
X_min_obs_a = zeros(n_r,p,n_f);
X_min_obs_b = zeros(n_r,p,n_f);
X_min_obs_c = zeros(n_r,p,n_f);
X_min_obs_d = zeros(n_r,p,n_f);
if out_opt_struct.est_pks==1
    X_max_est_a = zeros(n_r,p,n_f);
    X_max_est_b = zeros(n_r,p,n_f);
    X_max_est_c = zeros(n_r,p,n_f);
    X_max_est_d = zeros(n_r,p,n_f);
    X_min_est_a = zeros(n_r,p,n_f);
    X_min_est_b = zeros(n_r,p,n_f);
    X_min_est_c = zeros(n_r,p,n_f);
    X_min_est_d = zeros(n_r,p,n_f);
end    
if out_opt_struct.pk_loads == 1
    load_max_a = zeros(n_a, n_r, p, n_f);
    load_max_b = zeros(n_a, n_r, p, n_f);
    load_max_c = zeros(n_a, n_r, p, n_f);
    load_max_d = zeros(n_a, n_r, p, n_f);
    load_min_a = zeros(n_a, n_r, p, n_f);
    load_min_b = zeros(n_a, n_r, p, n_f);
    load_min_c = zeros(n_a, n_r, p, n_f);
    load_min_d = zeros(n_a, n_r, p, n_f);
end

h1 = waitbar(0, ['Please wait... Analyzing HDF file ' num2str(1) ' of ' num2str(p)],'Color',clr,'WindowStyle','modal',...
    'CreateCancelBtn',@WaitBarCancel,'Name','Analysis in progress...','CloseRequestFcn',@WaitBarCancel); % cancel button

% loop through list of wind directions
for i = 1:p
    
    HDF_filename_i = [hdf_filename_root theta_a_str(i,:) '.HDF'];
    for f = 1:n_f
        drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
        waitbar((i-1)/p, h1, ['Please wait... Analyzing HDF file ' num2str(i) ' of ' num2str(p)]); drawnow; % update waitbar

    	Cp = cpmatrix( fullfile(HDF_path, HDF_filename_i) , Astruct(f).tap_ind ); % Pressure coefficient time series
        TS_struct(f).n_s = size(Cp,2); % number of samples
        drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
        
        % Compute time series of response quantities of interest using influence vectors:
        TS_struct(f).ts_a = 0.5*ro*N'*Astruct(f).Aa*Cp;
        TS_struct(f).ts_b = 0.5*ro*N'*Astruct(f).Ab*Cp;
        TS_struct(f).ts_c = 0.5*ro*N'*Astruct(f).Ac*Cp; 
        TS_struct(f).ts_d = 0.5*ro*N'*Astruct(f).Ad*Cp;
    
        % Compute peak value (largest or smallest) of each internal force quantity;
        % Identify time sample at which peak occurs for each quantity and 
        % store the corresponding load distribution:
        [X_max_obs_a(:,i,f), TS_struct(f).max_sample_a] = max(TS_struct(f).ts_a,[],2);
        [X_max_obs_b(:,i,f), TS_struct(f).max_sample_b] = max(TS_struct(f).ts_b,[],2);
        [X_max_obs_c(:,i,f), TS_struct(f).max_sample_c] = max(TS_struct(f).ts_c,[],2);
        [X_max_obs_d(:,i,f), TS_struct(f).max_sample_d] = max(TS_struct(f).ts_d,[],2);

        [X_min_obs_a(:,i,f), TS_struct(f).min_sample_a] = min(TS_struct(f).ts_a,[],2);
        [X_min_obs_b(:,i,f), TS_struct(f).min_sample_b] = min(TS_struct(f).ts_b,[],2);
        [X_min_obs_c(:,i,f), TS_struct(f).min_sample_c] = min(TS_struct(f).ts_c,[],2);
        [X_min_obs_d(:,i,f), TS_struct(f).min_sample_d] = min(TS_struct(f).ts_d,[],2);

        if out_opt_struct.pk_loads == 1
            load_max_a( :, :, i, f ) = 0.5*ro*Astruct(f).Aa*Cp(:, TS_struct(f).max_sample_a );
            load_max_b( :, :, i, f ) = 0.5*ro*Astruct(f).Ab*Cp(:, TS_struct(f).max_sample_b );
            load_max_c( :, :, i, f ) = 0.5*ro*Astruct(f).Ac*Cp(:, TS_struct(f).max_sample_c );
            load_max_d( :, :, i, f ) = 0.5*ro*Astruct(f).Ad*Cp(:, TS_struct(f).max_sample_d );
            load_min_a( :, :, i, f ) = 0.5*ro*Astruct(f).Ad*Cp(:, TS_struct(f).min_sample_a );
            load_min_b( :, :, i, f ) = 0.5*ro*Astruct(f).Ab*Cp(:, TS_struct(f).min_sample_b );
            load_min_c( :, :, i, f ) = 0.5*ro*Astruct(f).Ac*Cp(:, TS_struct(f).min_sample_c );
            load_min_d( :, :, i, f ) = 0.5*ro*Astruct(f).Ad*Cp(:, TS_struct(f).min_sample_d );
        end
        drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
        if out_opt_struct.est_pks==1
            [X_max_est_a(:,i,f), X_min_est_a(:,i,f)] = maxminest( TS_struct(f).ts_a );
            [X_max_est_b(:,i,f), X_min_est_b(:,i,f)] = maxminest( TS_struct(f).ts_b );
            drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
            [X_max_est_c(:,i,f), X_min_est_c(:,i,f)] = maxminest( TS_struct(f).ts_c );
            [X_max_est_d(:,i,f), X_min_est_d(:,i,f)] = maxminest( TS_struct(f).ts_d );
        end
    end % end loop through frames
    if out_opt_struct.ts_plot
        drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
        set(h1,'Visible','off'); % hide waitbar
        h2 = TSplot;
        uiwait(h2);    
        if isfield(DIF_struct,'error'), return; end
        set(h1,'Visible','on');
    end % end if plot_on
end % end loop through wind directions

if exist('h1','var'), delete(h1); end % close waitbar

theta_all = [theta_a; theta_b; theta_c; theta_d];
theta = unique(theta_all);
n_q = length(theta);
DIF_struct.theta_all = theta_all(:);
DIF_struct.theta = theta(:);
if out_opt_struct.obs_pks
    X_max_obs_all = cat(2,X_max_obs_a, X_max_obs_b, X_max_obs_c, X_max_obs_d);
    X_min_obs_all = cat(2,X_min_obs_a, X_min_obs_b, X_min_obs_c, X_min_obs_d);
    X_max_obs = zeros(n_r,n_q,n_f);
    X_min_obs = zeros(n_r,n_q,n_f);
    for i = 1:n_q
        X_max_obs(:,i,:)=mean(X_max_obs_all(:,find(theta_all==theta(i)),:),2);
        X_min_obs(:,i,:)=mean(X_min_obs_all(:,find(theta_all==theta(i)),:),2);
    end
    if out_opt_struct.mult_pks
        DIF_struct.X_max_obs_all = X_max_obs_all;
        DIF_struct.X_min_obs_all = X_min_obs_all;
    end
    DIF_struct.X_max_obs = X_max_obs;
    DIF_struct.X_min_obs = X_min_obs;
end
if out_opt_struct.est_pks
    X_max_est_all = cat(2,X_max_est_a, X_max_est_b, X_max_est_c, X_max_est_d);
    X_min_est_all = cat(2,X_min_est_a, X_min_est_b, X_min_est_c, X_min_est_d);
    X_max_est = zeros(n_r,n_q,n_f);
    X_min_est = zeros(n_r,n_q,n_f);
    for i = 1:n_q
        X_max_est(:,i,:)=mean(X_max_est_all(:,find(theta_all==theta(i)),:),2);
        X_min_est(:,i,:)=mean(X_min_est_all(:,find(theta_all==theta(i)),:),2);
    end
    if out_opt_struct.mult_pks
        DIF_struct.X_max_est_all = X_max_est_all;
        DIF_struct.X_min_est_all = X_min_est_all;
    end
    DIF_struct.X_max_est = X_max_est;
    DIF_struct.X_min_est = X_min_est;
end
if out_opt_struct.pk_loads
    load_max_all = cat(3,load_max_a,load_max_b,load_max_c,load_max_d);
    load_min_all = cat(3,load_min_a,load_min_b,load_min_c,load_min_d);
    load_max = zeros(n_a,n_r,n_q,n_f);
    load_min = zeros(n_a,n_r,n_q,n_f);
    for i = 1:n_q
        load_max(:,:,i,:) = mean(load_max_all(:,:,find(theta_all==theta(i)),:),3);
        load_min(:,:,i,:) = mean(load_min_all(:,:,find(theta_all==theta(i)),:),3);
    end
    if out_opt_struct.mult_pks
        DIF_struct.load_max_all = load_max_all;
        DIF_struct.load_min_all = load_min_all;
    end
    DIF_struct.load_max = load_max;
    DIF_struct.load_min = load_min;
end


function WaitBarCancel(src, event)
    delete(h1);
    err = errordlg('Analysis cancelled by user', 'Analysis cancelled');
    uiwait(err); DIF_struct.error = 1;
end
% function WaitBarClose(src, event)
%     DIF_struct.error = 1;
%     delete(h1);
% end

function fig = TSplot

scrsz = get(0,'ScreenSize');
X0 = scrsz(3)*.05;
Y0 = scrsz(4)*.1;
dX = scrsz(3)*.9;
dY = scrsz(4)*.8;
fig = figure('Color',clr,'Position',[X0 Y0 dX dY],'NumberTitle','off','Name','Time Series Plots');

set( fig, 'Units', 'character');
Size = get( fig, 'Position');
X0 = Size(3);
Y0 = Size(4);
mw = 70; % menu width
Xp = X0-mw;
A = 5;
B = 13;
C = 7;

bx = 10; % plot border width in x-direction
by = 4; % plot border width in y-direction
plot_panel = uipanel('Parent',fig,'BackgroundColor',clr,'Units','character','Position',[0 0 Xp Y0],...
    'FontSize',12,'FontWeight','bold','TitlePosition','centertop','Title','Time Series Plots');
hA = axes('Parent', plot_panel, 'Units', 'character','Position',[bx Y0/2+by Xp/2-2*bx Y0/2-2*by]);
hB = axes('Parent', plot_panel, 'Units', 'character','Position',[bx by Xp/2-2*bx Y0/2-2*by]);
hC = axes('Parent', plot_panel, 'Units', 'character','Position',[Xp/2+bx Y0/2+by Xp/2-2*bx Y0/2-2*by]);
hD = axes('Parent', plot_panel, 'Units', 'character','Position',[Xp/2+bx by Xp/2-2*bx Y0/2-2*by]);

file_panel = uipanel('Parent',fig,'BackgroundColor',clr,'Units','character','Position',[X0-mw Y0-A mw-1 A],...
    'FontSize',12,'FontWeight','bold','TitlePosition','centertop','Title',['Current file (' num2str(i) ' of ' num2str(p) '):']);
select_panel = uipanel('Parent',fig,'BackgroundColor',clr,'Units','character','Position',[X0-mw Y0-(A+B) mw-1 B],...
    'FontSize',12,'FontWeight','bold','TitlePosition','centertop','Title','Select time series to plot (and save):');
proceed_panel = uipanel('Parent',fig,'BackgroundColor',clr,'Units','character','Position',[X0-mw Y0-(A+B+C) mw-1 C],...
    'FontSize',12,'FontWeight','bold','TitlePosition','centertop','Title','Analysis Options:');
%    'FontSize',12,'FontWeight','bold','BorderType','none');

% filename_txt = uicontrol(file_panel,'BackgroundColor',clr, 'Style', 'text',...
%     'String',['Current filename (file ' num2str(i) ' of ' num2str(p) '):'],...
%     'HorizontalAlignment','left','Units','characters','Position',[2 A-2 mw-4 1.5]);
filename_edit = uicontrol('Parent',file_panel,'BackgroundColor',clr,'Units','character','Position',[2 A-4 mw-6 1.5],...
    'HorizontalAlignment','left','Style','edit','Enable','inactive','String',HDF_filename_i);

frame_txt = uicontrol(select_panel,'BackgroundColor',clr, 'Style', 'text', 'String','Select frame:',...
    'HorizontalAlignment','left','Units','characters','Position',[2 B-4 mw-4 1.5]);
frame_popup = uicontrol(select_panel, 'Style', 'popupmenu', 'String',frame_label,'Value',1,...
    'Units','characters','Position',[2 B-5.5 mw-4 1.5],'BackgroundColor','w','Callback',@UpdatePlots );

resp_txt = uicontrol(select_panel,'BackgroundColor',clr, 'Style', 'text', 'String','Select response:',...
    'HorizontalAlignment','left','Units','characters','Position',[2 B-8 mw-4 1.5]);
resp_popup = uicontrol(select_panel, 'Style', 'popupmenu', 'String',resp_label,'Value',1,...
    'Units','characters','Position',[2 B-9.5 mw-4 1.5],'BackgroundColor','w','Callback',@UpdatePlots );
save_button = uicontrol(select_panel, 'Style','pushbutton','String','Save selected time series as...',...
    'Units','characters','Position',[2 B-12 45 1.5],'Callback',@SaveTS);

cont_button = uicontrol(proceed_panel, 'Style','pushbutton','String','Continue to next file',...
    'Units','characters','Position',[0.5*(mw-1)-25-2 C-4 25 1.5],'Callback',@ClosePlot);
plot_chkbox = uicontrol(proceed_panel, 'Style','checkbox','String','Suppress further time series plots',...
    'Value',0,'Units','characters','Position',[0.5*(mw-1)-25-2 C-6 40 1.5]);
cancel_button = uicontrol(proceed_panel, 'Style','pushbutton','String','Cancel analysis',...
    'Units','characters','Position',[0.5*(mw-1)+2 C-4 25 1.5],'Callback',@CancelAnalysis);

UpdatePlots;

function UpdatePlots( src, event )
    f = get(frame_popup,'Value');
    r = get(resp_popup,'Value');
    axes(hA);
    plot([1:TS_struct(f).n_s],TS_struct(f).ts_a(r,:),'k-',...
        TS_struct(f).max_sample_a(r), X_max_obs_a(r,i,f),'ro', TS_struct(f).min_sample_a(r), X_min_obs_a(r,i,f),'bo');
    if out_opt_struct.est_pks
        hold on; plot([1 TS_struct(f).n_s],[X_max_est_a(r,i,f) X_max_est_a(r,i,f)],'r-');
        plot([1 TS_struct(f).n_s],[X_min_est_a(r,i,f) X_min_est_a(r,i,f)],'b-'); hold off;
    end
    title(['Wind direction A = ' num2str(theta_a(i))]);
    xlabel('Sample number');
    ylabel([resp_names{r} ' (' resp_units{r} ')']);
    axes(hB);
    plot([1:TS_struct(f).n_s],TS_struct(f).ts_b(r,:),'k-',...
        TS_struct(f).max_sample_b(r), X_max_obs_b(r,i,f),'ro', TS_struct(f).min_sample_b(r), X_min_obs_b(r,i,f),'bo');
    if out_opt_struct.est_pks
        hold on; plot([1 TS_struct(f).n_s],[X_max_est_b(r,i,f) X_max_est_b(r,i,f)],'r-');
        plot([1 TS_struct(f).n_s],[X_min_est_b(r,i,f) X_min_est_b(r,i,f)],'b-'); hold off;
    end
    title(['Wind direction B = ' num2str(theta_b(i))]);
    xlabel('Sample number');
    ylabel([resp_names{r} ' (' resp_units{r} ')']);
    axes(hC);
    plot([1:TS_struct(f).n_s],TS_struct(f).ts_c(r,:),'k-',...
        TS_struct(f).max_sample_c(r), X_max_obs_c(r,i,f),'ro', TS_struct(f).min_sample_c(r), X_min_obs_c(r,i,f),'bo');
    if out_opt_struct.est_pks
        hold on; plot([1 TS_struct(f).n_s],[X_max_est_c(r,i,f) X_max_est_c(r,i,f)],'r-');
        plot([1 TS_struct(f).n_s],[X_min_est_c(r,i,f) X_min_est_c(r,i,f)],'b-'); hold off;
    end
    title(['Wind direction C = ' num2str(theta_c(i))]);
    xlabel('Sample number');
    ylabel([resp_names{r} ' (' resp_units{r} ')']);
    axes(hD);
    plot([1:TS_struct(f).n_s],TS_struct(f).ts_d(r,:),'k-',...
        TS_struct(f).max_sample_d(r), X_max_obs_d(r,i,f),'ro', TS_struct(f).min_sample_d(r), X_min_obs_d(r,i,f),'bo');
    if out_opt_struct.est_pks
        hold on; plot([1 TS_struct(f).n_s],[X_max_est_d(r,i,f) X_max_est_d(r,i,f)],'r-');
        plot([1 TS_struct(f).n_s],[X_min_est_d(r,i,f) X_min_est_d(r,i,f)],'b-'); hold off;
    end
    title(['Wind direction D = ' num2str(theta_d(i))]);
    xlabel('Sample number');
    ylabel([resp_names{r} ' (' resp_units{r} ')']);
end

function SaveTS( src, event )
    f = get(frame_popup,'Value');
    r = get(resp_popup,'Value');
    TS_filename = ['TS_' HDF_filename_i(1:end-4) '_f' num2str(f) 'r' num2str(r) '.csv']; 
    if exist('tmp_DIF_path.mat','file'), load('tmp_DIF_path'); end
    if ~exist('DIF_path','var') || ~isstr(DIF_path), DIF_path = ''; end
    [TS_filename, TS_path, filetype_ind]= uiputfile({'*.csv','Comma delimited (*.csv)';...
        '*.mat', 'Matlab MAT-files (*.mat)'},'Save time series as', fullfile(DIF_path,TS_filename));
    if filetype_ind==0
        figure(fig); return;
    end
    TS_full_filename = fullfile(TS_path, TS_filename);
    if filetype_ind==1 
        fw = fopen( TS_full_filename, 'wt');
        if fw==-1
            err = errordlg('Unable to open specified file for writing: make sure that the file is not currently open.',...
                'Error saving time series'); uiwait(err);
            figure(fig); return;
        end
        fprintf( fw, '%s\n',['HDF pressure file: ' HDF_filename_i]);
        fprintf( fw, '%s\n',frame_label{f});
        fprintf( fw, '%s\n',resp_label{r});
        fprintf( fw, '%s\n',['Number of samples: ' num2str(TS_struct(f).n_s)]);
        fprintf( fw, '%s\n','Wind Dir A,Wind Dir B,Wind Dir C,Wind Dir D');
        fprintf( fw, '%g,%g,%g,%g\n',[theta_a(i) theta_b(i) theta_c(i) theta_d(i)]);
        fprintf( fw, '%s\n','Time Series A,Time Series B,Time Series C,Time Series D');
        fprintf( fw, '%g,%g,%g,%g\n', [TS_struct(f).ts_a(r,:); TS_struct(f).ts_b(r,:); TS_struct(f).ts_c(r,:); TS_struct(f).ts_d(r,:)] );
        status = fclose(fw);
        if status==-1
            err = errordlg('Unable to close specified file after writing.', 'Error saving time series'); uiwait(err);
            figure(fig); return;
        end
    elseif filetype_ind==2
        TS_structure.HDF_filename = HDF_filename_i;
        TS_structure.frame_label = frame_label{f};
        TS_structure.resp_label = resp_label{r};
        TS_structure.n_s = TS_struct(f).n_s;
        TS_structure.theta_a = theta_a(i);
        TS_structure.theta_b = theta_b(i);
        TS_structure.theta_c = theta_c(i);
        TS_structure.theta_d = theta_d(i);
        TS_structure.ts_a = TS_struct(f).ts_a(r,:);
        TS_structure.ts_b = TS_struct(f).ts_b(r,:);
        TS_structure.ts_c = TS_struct(f).ts_c(r,:);
        TS_structure.ts_d = TS_struct(f).ts_d(r,:);
        try
            save(TS_full_filename, 'TS_structure' );
        catch
            msgstr = lasterr;
            err = errordlg(msgstr, 'Error saving time series'); uiwait(err);
            figure(fig); return;
        end
    end
    figure(fig);
end
        
function ClosePlot( src, event )
    if get(plot_chkbox,'Value')
        out_opt_struct.ts_plot = 0;
    end
    delete(fig);
end
function CancelAnalysis( src, event )
    DIF_struct.error = 1;
    delete(fig);
end

end
        
end

%ok = writeDIFfile( DIF_struct, DIF_full_filename )

% 
% h1 = figure;
% h2 = figure;
% for i = 1:n_r
%     figure(h1)
%     subplot(n_r,1,i)
%     plot(wind_dir_360, alpha_pk1_360(i,:),'ko'); hold on;
%     plot(wind_dir_list', abs(alpha_pk1(i,:)),'k+');
%     plot(wind_dir_list_op', abs(alpha_pk1_op(i,:)),'k.');
%     plot(wind_dir_360_S, alpha_pk1_360_S(i,:),'bo');
%     plot(wind_dir_list_S', abs(alpha_pk1_S(i,:)),'b+');
%     plot(wind_dir_list_S_op', abs(alpha_pk1_S_op(i,:)),'b.');
%     plot(wind_dir_360, alpha_pk1_360_ave(i,:),'m.-');
%     hold off;
%     figure(h2)
%     subplot(n_r,1,i)
%     plot(wind_dir_360, alpha_pk2_360(i,:),'ko'); hold on;
%     plot(wind_dir_list', abs(alpha_pk2(i,:)),'k+');
%     plot(wind_dir_list_op', abs(alpha_pk2_op(i,:)),'k.');
%     plot(wind_dir_360_S, alpha_pk2_360_S(i,:),'bo');
%     plot(wind_dir_list_S', abs(alpha_pk2_S(i,:)),'b+');
%     plot(wind_dir_list_S_op', abs(alpha_pk2_S_op(i,:)),'b.');
%     plot(wind_dir_360, alpha_pk2_360_ave(i,:),'m.-');
%     hold off;
% 
% end
% 



