function MRI_struct = computeMRIpeaks( bldg_struct, DIF_struct_cell, HURR_struct, MRI_list )

% check to make sure that bldg_struct and DIF_struct correspond
% if the terrain conditions are directional, then DIF_struct should be
% a structure array of length 2, with one "open" and one "suburban" 
% DIF_struct should be consistent


% PROCEDURE:


% for each building orientation:

% define list of wind directions relative to reference axis

% compute DIFs for each wind direction

err_dlgname = 'Error computing MRI peaks';

tol = 10^-3;
% Convert the input 'DIF_struct_cell' to a structure array 'DIF_struct' in which each element has the same fields:
if length(DIF_struct_cell)==1
    DIF_struct = DIF_struct_cell{1};
    n_r = length(DIF_struct(1).resp_names); % number of responses
    % make sure terrain in DIF results matches building input file:
    if ~ismember(bldg_struct.terrain, {DIF_struct.terrain,'Directional'})
        err = errordlg(['The terrain in the building input file (''' bldg_struct.terrain ...
            ''') does not match the terrain in the DIF results file (''' DIF_struct.terrain ''').'],err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    end
elseif length(DIF_struct_cell)==2
    % Make sure building input file specifies 'Directional' terrain:
    if ~strcmp(bldg_struct.terrain,'Directional')
        err = errordlg(['Only one DIF results file expected for the terrain specified in the building input file (''' bldg_struct.terrain ...
            '''). Multiple DIF results allowed for ''Directional'' terrain only.' ],err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    end
    fields1 = fieldnames(DIF_struct_cell{1});
    fields2 = fieldnames(DIF_struct_cell{2});
    fields = intersect(fields1,fields2);
    % Create a structure array containing only the common fields:
    for i=1:2
        for j=1:length(fields)
            DIF_struct(i).(fields{j}) = DIF_struct_cell{i}.(fields{j});
        end
    end
    % Make sure the terrain does not coincide for the two DIF results (one should be 'Open_Country' and one should be 'Suburban'):
    if strcmp(DIF_struct(1).terrain, DIF_struct(2).terrain)
        err = errordlg(['Both DIF files correspond to the same type of terrain (''' DIF_struct(1).terrain ...
            '''). DIF files for both ''Open_Country'' and ''Suburban'' terrain are required to compute MRIs by interpolation.'],err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    % make sure units, wind directions, frame locations, and building dimensions coincide:
    elseif ~strcmp(DIF_struct(1).force_units,DIF_struct(2).force_units)
        err = errordlg('The force units must coincide for all DIF results used in interpolation', err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    elseif ~strcmp(DIF_struct(1).length_units,DIF_struct(2).length_units)
        err = errordlg('The length units must coincide for all DIF results used in interpolation', err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    elseif ~strcmp(DIF_struct(1).ws_units,DIF_struct(2).ws_units)
        err = errordlg('The wind speed units must coincide for all DIF results used in interpolation', err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    elseif length(DIF_struct(1).theta)~=length(DIF_struct(i).theta) || any(DIF_struct(1).theta~=DIF_struct(i).theta)
        err = errordlg( 'The wind directions must coincide for all DIF results used in interpolation', err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    elseif any(size(DIF_struct(1).frame_coords)~=size(DIF_struct(2).frame_coords)) || max(max(abs(DIF_struct(1).frame_coords-DIF_struct(2).frame_coords)))>tol
        err = errordlg( 'The frame coordinates must coincide for all DIF results used in interpolation', err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    elseif max(max(abs(DIF_struct(1).d0-DIF_struct(2).d0)))>tol
        err = errordlg( 'The building dimensions must coincide for all DIF results used in interpolation',err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    elseif size(DIF_struct(1).d,1)~=size(DIF_struct(2).d,1) || any(any(DIF_struct(1).d~=DIF_struct(2).d))
        err = errordlg( 'The model dimensions must coincide for all DIF results used in interpolation.',err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    end
    % make sure the number of responses, the response names, and the response units coincide:
    if length(DIF_struct(1).resp_names)~=length(DIF_struct(2).resp_names)
        err = errordlg('The number of response quantities must coincide for all DIF results used in interpolation', err_dlgname);
        uiwait(err); MRI_struct.error = 1; return;
    end
    n_r = length(DIF_struct(1).resp_names);
    for j=1:n_r
        if ~strcmp(DIF_struct(1).resp_names(j), DIF_struct(2).resp_names(j))
            err = errordlg('The response names must coincide for all DIF results used in interpolation', err_dlgname);
            uiwait(err); MRI_struct.error = 1; return;
        elseif ~strcmp(DIF_struct(1).resp_units(j), DIF_struct(2).resp_units(j))
            err = errordlg('The response units must coincide for all DIF results used in interpolation', err_dlgname);
            uiwait(err); MRI_struct.error = 1; return;
        end
    end
else
    err = errordlg(['Unexpected number of DIF files (' num2str(length(DIF_struct_cell)) '; 1 or 2 expected.'], err_dlgname);
    uiwait(err); MRI_struct.error = 1; return;
end
MRI_struct.length_units = bldg_struct.length_units;
MRI_struct.force_units = bldg_struct.force_units;
MRI_struct.ws_units = bldg_struct.ws_units;
MRI_struct.d0 = bldg_struct.d0;
MRI_struct.d = DIF_struct(1).d;
MRI_struct.frame_coords = bldg_struct.frame_coords;
MRI_struct.terrain = bldg_struct.terrain;
MRI_struct.z0 = bldg_struct.z0;
MRI_struct.model_terrain = {DIF_struct.terrain};
MRI_struct.resp_names = DIF_struct(1).resp_names;
MRI_struct.resp_units = DIF_struct(1).resp_units;

% Identify which quantities are available to interpolate:
% initialize boolean values for which quantities to interpolate:
max_obs = 1;
min_obs = 1;
max_est = 1;
min_est = 1;

for i = 1:length(DIF_struct)
    if ~isfield(DIF_struct, 'X_max_obs') || isempty(DIF_struct(i).X_max_obs)
        max_obs = 0;
    elseif max_obs==1
        X_max_obs = DIF_struct(i).X_max_obs;
        n_f = size(X_max_obs,3);
        % 'wrap' DIFs above 360 and below 0 to facilitate resampling:
        DIF(i).X_max_obs_wrap = cat(2,X_max_obs,X_max_obs,X_max_obs);
    end

    if ~isfield(DIF_struct, 'X_min_obs') || isempty(DIF_struct(i).X_min_obs)
        min_obs = 0;
    elseif min_obs==1
        X_min_obs = DIF_struct(i).X_min_obs;
        n_f = size(X_min_obs,3);
        % 'wrap' DIFs above 360 and below 0 to facilitate resampling:
        DIF(i).X_min_obs_wrap = cat(2,X_min_obs,X_min_obs,X_min_obs);
    end

    if ~isfield(DIF_struct, 'X_max_est') || isempty(DIF_struct(i).X_max_est)
        max_est = 0;
    elseif max_est==1
        X_max_est = DIF_struct(i).X_max_est;
        n_f = size(X_max_est,3);
        % 'wrap' DIFs above 360 and below 0 to facilitate resampling:
        DIF(i).X_max_est_wrap = cat(2,X_max_est,X_max_est,X_max_est);
    end

    if ~isfield(DIF_struct, 'X_min_est') || isempty(DIF_struct(i).X_min_est)
        min_est = 0;
    elseif min_est==1
        X_min_est = DIF_struct(i).X_min_est;
        n_f = size(X_min_est,3);
        % 'wrap' DIFs above 360 and below 0 to facilitate resampling:
        DIF(i).X_min_est_wrap = cat(2,X_min_est,X_min_est,X_min_est);
    end
end
if ~max_obs & ~min_obs & ~max_est & ~min_est
    err = errordlg('Neither observed nor estimated DIF results are available for interpolation',err_dlgname);
    uiwait(err); MRI_struct.error = 1; return;
end

MRI_struct.max_obs = max_obs;
MRI_struct.min_obs = min_obs;
MRI_struct.max_est = max_est;
MRI_struct.min_est = min_est;

MRI_struct.DIF_struct = DIF_struct;

theta = DIF_struct(1).theta;

if any(theta~=sort(theta))
    error('DIF_struct.theta must be sorted in increasing order');
end    
% 'wrap' wind directions above 360 and below 0 to facilitate resampling:
theta_wrap = [theta(:)-360; theta(:); theta(:)+360];

% define upper and lower boundaries of sector for each wind direction:
theta_wrap_max = theta_wrap; theta_wrap_max(1:end-1) = theta_wrap(1:end-1)+ 0.5*diff(theta_wrap);
theta_wrap_min = theta_wrap; theta_wrap_min(2:end)   = theta_wrap(2:end)  - 0.5*diff(theta_wrap);

% define list of orientations of ridge axis of building:
d_a0 = 15; % increment in orientation (degrees)
alpha_0 = 0:d_a0:360-d_a0; % orientations (degrees clockwise from north)
n_o = length(alpha_0); % number of orientations considered

% Scale directional wind speeds to eave height over appropriate terrain:
Vsim = HURR_struct.Vsim; % matrix of simulated wind speeds (1-hour average in mph at 10 m elevation over open terrain)
% define list of wind directions corresponding to the columns of HURR_struct.Vsim:
if size(Vsim,2)==16
    alpha = 22.5:22.5:360;     % (degrees clockwise from north)
    delta_alpha = 22.5; % increment in wind direction
    n_q = length(alpha); % number of wind directions
else
    err = errordlg(['Unexpected number of columns in HURR_struct.Vsim (' num2str(size(Vsim,2)) '); 16 expected.'], err_dlgname);
    uiwait(err); MRI_struct.error = 1; return;
end
Z = 10;    % elevation (m) corresponding to wind speeds in Vsim
Z0 = 0.03; % roughness length (m) corresponding to wind speeds in Vsim
H = bldg_struct.d0(3); % eave height (units should be 'ft')
if strcmp(bldg_struct.length_units,'ft')
    H = H*.3048; % convert eave height from 'ft' to 'm'
else
    err = errordlg('Length units in building input file must be ''ft''', err_dlgname);
    uiwait(err); MRI_struct.error = 1; return;
end
if strcmp(bldg_struct.terrain, 'Open_Country')
    z0 = 0.03*ones(size(alpha)); %roughness lengths (m)
    DIF.mu = ones(size(alpha));
elseif strcmp(bldg_struct.terrain, 'Suburban')
    z0 = 0.3*ones(size(alpha)); % roughness lengths (m)
    DIF.mu = ones(size(alpha));
elseif strcmp(bldg_struct.terrain, 'Directional')
    z0 = interp1(0:45:360, [bldg_struct.z0 bldg_struct.z0(1)], alpha)*.3048; % roughness lengths (m)
    if length(DIF_struct)==1
        % DIFs from either 'Open_Country' or 'Suburban' terrain will be used to compute MRIs (no interpolation)
        DIF.mu = ones(size(alpha));
    elseif length(DIF_struct)==2
        % Interpolation between DIFs for 'Open_Country' and 'Suburban' terrain
        for i=1:length(DIF_struct)
            % Define weighting factors for interpolation between 'Open_Country' and 'Suburban' DIFs:
            if strcmp(DIF_struct(i).terrain, 'Open_Country')
                DIF(i).mu = (.3-z0)/(0.3-0.03);
                DIF(i).mu(find(z0>0.3)) = 0;
                DIF(i).mu(find(z0<0.03)) = 1;
            elseif strcmp(DIF_struct(i).terrain, 'Suburban')
                DIF(i).mu = (z0-0.03)/(0.3-0.03);
                DIF(i).mu(find(z0>0.3)) = 1;
                DIF(i).mu(find(z0<0.03)) = 0;
            else
                err = errordlg(['The terrain specified in the DIF file (''' DIF_struct(i).terrain ...
                    ''') is unrecognized; either ''Open_Country'' or ''Suburban'' expected.'],err_dlgname);
                uiwait(err); MRI_struct.error = 1; return;
            end
        end
    end
else
    err = errordlg(['The terrain specified in the building input file (''' bldg_struct(i).terrain ...
        ''') is unrecognized; either ''Open_Country'', ''Suburban'', or ''Directional'' expected.'],err_dlgname);
    uiwait(err); MRI_struct.error = 1; return;
end
lambda = (z0/Z0).^0.0706.*log(H./z0)/log(Z/Z0); % scale factors for each direction
V_H = Vsim*diag(lambda); % wind speeds scaled to eave height over appropriate terrain

h = size(V_H,1); % number of hurricanes
mu = HURR_struct.recurrence; % hurricane recurrence rate
MRI_rank = (1-exp(-mu*(1:h)/(h+1))).^-1; % MRIs associated with ranked responses

MRI_struct.milepost = HURR_struct.milepost;

if ~exist('MRI_list') || isempty(MRI_list)
    % Prompt user to select MRIs for which responses are desired:
    MRI_list = [50 100 250 500 750 1000];
    MRI_list = MRI_list(find(MRI_list<max(MRI_rank))); % make sure maximum return period is not exceeded 
    for i=1:length(MRI_list)
        MRI_cell{i} = num2str(MRI_list(i));
    end

    [MRI_sel,ok] = listdlg('ListString',MRI_cell, 'InitialValue', [1 4], ...
        'PromptString',{'Select Mean Recurrence','Intervals (MRIs) for which','responses are required:'} , ...
        'Name','MRI selection','ListSize',[160 140]);
    if ok
        MRI_list = MRI_list(MRI_sel);
    else
        err = errordlg('MRI selection cancelled by user.','Analysis cancelled');
        uiwait(err); MRI_struct.error = 1; return;
    end
end

MRI_struct.orientations = alpha_0;
MRI_struct.MRI_list = MRI_list;

for o=1:n_o % loop through list of building orientations
    
    theta_hat = alpha - alpha_0(o); % wind direction relative to reference axis of building
    theta_hat(find(theta_hat<0))   = theta_hat(find(theta_hat<0))+360;     % make sure theta_hat is >=0
    theta_hat(find(theta_hat>=360)) = theta_hat(find(theta_hat>=360))-360;   % make sure theta_hat is <360
    DIF_hat(o).theta_hat = theta_hat;
    
    % define upper and lower boundaries of sector for each wind direction:
    theta_hat_max = theta_hat+delta_alpha/2;
    theta_hat_min = theta_hat-delta_alpha/2;
    
    % RESAMPLE DIFs:
    % initialize empty arrays:
    if max_obs
        DIF_hat(o).X_max_obs = zeros(n_r,n_q,n_f);
    end
    if min_obs
        DIF_hat(o).X_min_obs = zeros(n_r,n_q,n_f);
    end
    if max_est
        DIF_hat(o).X_max_est = zeros(n_r,n_q,n_f);
    end
    if min_est
        DIF_hat(o).X_min_est = zeros(n_r,n_q,n_f);
    end
        
    for q=1:length(theta_hat) % loop through wind directions

        overlap = min(theta_wrap_max, theta_hat_max(q)*ones(size(theta_wrap_max)))...
                - max(theta_wrap_min, theta_hat_min(q)*ones(size(theta_wrap_min)));
        overlap(find(overlap<0))=0;
        overlap_fraction = overlap/delta_alpha;

        for f=1:n_f % loop through frames
            for i=1:length(DIF_struct)
                if max_obs
                    DIF_hat(o).X_max_obs(:,q,f) = DIF_hat(o).X_max_obs(:,q,f)...
                        + DIF(i).mu(q)*sum(DIF(i).X_max_obs_wrap(:,:,f)*diag(overlap_fraction),2);
                end
                if min_obs
                    DIF_hat(o).X_min_obs(:,q,f) = DIF_hat(o).X_min_obs(:,q,f)...
                        + DIF(i).mu(q)*sum(DIF(i).X_min_obs_wrap(:,:,f)*diag(overlap_fraction),2);
                end
                if max_est
                    DIF_hat(o).X_max_est(:,q,f) = DIF_hat(o).X_max_est(:,q,f)...
                        + DIF(i).mu(q)*sum(DIF(i).X_max_est_wrap(:,:,f)*diag(overlap_fraction),2);
                end
                if min_est
                    DIF_hat(o).X_min_est(:,q,f) = DIF_hat(o).X_min_est(:,q,f)...
                        + DIF(i).mu(q)*sum(DIF(i).X_min_est_wrap(:,:,f)*diag(overlap_fraction),2);
                end
            end
        end
    end

    for f=1:n_f % loop through frames
        for r=1:n_r % loop through responses
            if max_obs
                r_rank = sort(max(V_H.^2*diag(DIF_hat(o).X_max_obs(r,:,f)),[],2),1,'descend'); % rank ordered max responses in each hurricane
                MRI_struct.MRI_max_obs(:,f,o,r) = interp1(MRI_rank, r_rank, MRI_list ,'cubic'); 
            end
            if min_obs
                r_rank = sort(min(V_H.^2*diag(DIF_hat(o).X_min_obs(r,:,f)),[],2),1,'ascend'); % rank ordered min responses in each hurricane
                MRI_struct.MRI_min_obs(:,f,o,r) = interp1(MRI_rank, r_rank, MRI_list ,'cubic'); 
            end
            if max_est
                r_rank = sort(max(V_H.^2*diag(DIF_hat(o).X_max_est(r,:,f)),[],2),1,'descend'); % rank ordered max responses in each hurricane
                MRI_struct.MRI_max_est(:,f,o,r) = interp1(MRI_rank, r_rank, MRI_list ,'cubic'); 
            end
            if min_est
                r_rank = sort(min(V_H.^2*diag(DIF_hat(o).X_min_est(r,:,f)),[],2),1,'ascend'); % rank ordered min responses in each hurricane
                MRI_struct.MRI_min_est(:,f,o,r) = interp1(MRI_rank, r_rank, MRI_list ,'cubic'); 
            end
        end
    end
end

MRI_struct.DIF_hat = DIF_hat;

% Unknown orientation -- expected peaks for uniform probability of each orientation:
if max_obs
    MRI_struct.MRI_UKO_max_obs = squeeze(mean(MRI_struct.MRI_max_obs,3)); 
end
if min_obs
    MRI_struct.MRI_UKO_min_obs = squeeze(mean(MRI_struct.MRI_min_obs,3)); 
end
if max_est
    MRI_struct.MRI_UKO_max_est = squeeze(mean(MRI_struct.MRI_max_est,3)); 
end
if min_est
    MRI_struct.MRI_UKO_min_est = squeeze(mean(MRI_struct.MRI_min_est,3)); 
end