function interp_DIF_struct = interpDIF( DIF_struct_cell )

err_dlgname = 'Error in DIF interpolation';

% Load saved sensitivity factors for building dimensions:
if exist('tmp_s_f.mat','file'), load 'tmp_s_f'; end
% If saved options are not found, set to default options:
if ~exist('s_f','var')
    err = errordlg('Sensitivity factors for interpolation not defined.', err_dlgname);
end
s_f = s_f;

n_m = length( DIF_struct_cell ); % number of models to interpolate between
if n_m<2
    err = errordlg(['Results from at least two building models are required for interpolation:' ...
        'DIF_struct_array must have length > 1'], err_dlgname);
    uiwait(err); interp_DIF_struct.error = 1; return;
end    

% Assemble list of fields common to all structures in DIF_struct_cell:
fields = fieldnames(DIF_struct_cell{1});
for i=2:n_m
    fields_i = fieldnames(DIF_struct_cell{i});
    fields = intersect(fields,fields_i);
end
% Create a structure array containing only the common fields:
for i=1:n_m
    for j=1:length(fields)
        DIF_struct_array(i).(fields{j}) = DIF_struct_cell{i}.(fields{j});
    end
end

% error checking:
% make sure DIF_struct_array has length > 1


% make sure building dimensions, terrain, response names, and frame
% locations coincide for all models:

force_units = DIF_struct_array(1).force_units;
length_units = DIF_struct_array(1).length_units;
ws_units = DIF_struct_array(1).ws_units;
terrain = DIF_struct_array(1).terrain;
theta = DIF_struct_array(1).theta;
frame_coords = DIF_struct_array(1).frame_coords;
attach_pts = DIF_struct_array(1).attach_pts;
d0 = DIF_struct_array(1).d0;
resp_names = DIF_struct_array(1).resp_names;
resp_units = DIF_struct_array(1).resp_units;
n_r = length(resp_names);

tol = 10^-3; % tolerance for coincidence of dimensions, frame locations, scale factors
for i = 1:n_m
    if ~strcmp(force_units,DIF_struct_array(i).force_units)
        err = errordlg('The force units must coincide for all models used in interpolation', err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif ~strcmp(length_units,DIF_struct_array(i).length_units)
        err = errordlg('The length units must coincide for all models used in interpolation', err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif ~strcmp(ws_units,DIF_struct_array(i).ws_units)
        err = errordlg('The wind speed units must coincide for all models used in interpolation', err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif ~strcmp(terrain,DIF_struct_array(i).terrain)
        err = errordlg('The terrain conditions must coincide for all models used in interpolation', err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif length(theta)~=length(DIF_struct_array(i).theta) || any(theta~=DIF_struct_array(i).theta)
        err = errordlg( 'The wind directions must coincide for all models used in interpolation', err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif any(size(frame_coords)~=size(DIF_struct_array(i).frame_coords)) || max(max(abs(frame_coords-DIF_struct_array(i).frame_coords)))>tol
        err = errordlg( 'The frame coordinates must coincide for all models used in interpolation', err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif any(size(attach_pts)~=size(DIF_struct_array(i).attach_pts)) || max(max(abs(attach_pts-DIF_struct_array(i).attach_pts)))>tol
        err = errordlg( 'The attachment locations must coincide for all models used in interpolation', err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif max(max(abs(d0-DIF_struct_array(i).d0)))>tol
        err = errordlg( 'The building dimensions must coincide for all models used in interpolation',err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif max(max(abs(d0-DIF_struct_array(i).d)))<tol
        err = errordlg( 'For interpolation, the model dimensions must not exactly coincide with the building dimensions',err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    elseif size(DIF_struct_array(i).d,1)>1
        err = errordlg( 'The DIF results used in interpolation must not have been obtained by interpolation themselves.',err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    else
        d0_mx(i,:) = DIF_struct_array(i).d0;
        d(i,:) = DIF_struct_array(i).d;
    end
end
interp_DIF_struct.force_units = force_units;
interp_DIF_struct.length_units = length_units;
interp_DIF_struct.ws_units = ws_units;
interp_DIF_struct.terrain = terrain;
interp_DIF_struct.theta = theta;
interp_DIF_struct.frame_coords = frame_coords;
interp_DIF_struct.attach_pts = attach_pts;
interp_DIF_struct.d0 = d0;
interp_DIF_struct.d = d;

% Form matrix of normalized deviations of model dimensions from building dimensions:
delta(:,1) = abs(d(:,1)./d0_mx(:,1)-1);                         % width: W_i/W_0 -1
delta(:,2) = abs((d(:,2)./d(:,1))./(d0_mx(:,2)./d0_mx(:,1))-1); % aspect ratio: (L_i/W_i)/(L_0/W_0)-1
delta(:,3) = abs((d(:,3)./d(:,1))./(d0_mx(:,3)./d0_mx(:,1))-1); % height ratio: (H_i/W_i)/(H_0/W_0)-1
delta(:,4) = abs((d(:,4)./d(:,1))./(d0_mx(:,4)./d0_mx(:,1))-1); % roof slope: (R_i/W_i)/(R_0/W_0)-1

norm_s_f = s_f*4/sum(s_f); % normalized sensitivity factors
S = diag(norm_s_f); % diagonal matrix of sensitivity factors

% scaled vector deviations of model dimensions from building dimensions:
delta_s = delta*S;
% norm of scaled deviation:
dev = diag(delta_s*delta_s').^.5;

gam = sum(dev.^-1)^-1; % constant factor in numerator of weighted average
wf = gam./dev; % weighting factors for each model in weighted average

interp_DIF_struct.wf = wf;

% Identify which quantities are available to interpolate:
% initialize boolean values for which quantities to interpolate:
intrp.X_max_obs = 1;
intrp.X_min_obs = 1;
intrp.X_max_est = 1;
intrp.X_min_est = 1;
intrp.load_max = 1;
intrp.load_min = 1;

for i = 1:n_m
    if length(DIF_struct_array(i).resp_names) ~= n_r
        err = errordlg('The number of response quantities must coincide for all models used in interpolation', err_dlgname);
        uiwait(err); interp_DIF_struct.error = 1; return;
    end
    for j = 1:n_r
        if ~strcmp(resp_names(j), DIF_struct_array(i).resp_names(j))
            err = errordlg('The response names must coincide for all models used in interpolation', err_dlgname);
            uiwait(err); interp_DIF_struct.error = 1; return;
        elseif ~strcmp(resp_units(j), DIF_struct_array(i).resp_units(j))
            err = errordlg('The response units must coincide for all models used in interpolation', err_dlgname);
            uiwait(err); interp_DIF_struct.error = 1; return;
        end
    end
    
    if ~isfield(DIF_struct_array, 'X_max_obs') || isempty(DIF_struct_array(i).X_max_obs)
        intrp.X_max_obs = 0;
    elseif intrp.X_max_obs==1
        if exist('X_max_obs','var')
            if any(size(X_max_obs)~=size(DIF_struct_array(i).X_max_obs))
                err = errordlg('Mismatch in dimensions of the field ''X_max_obs''.', err_dlgname  );
                uiwait(err); interp_DIF_struct.error =1; return;
            end
        else
            X_max_obs = zeros(size(DIF_struct_array(i).X_max_obs));
            X_max_obs_bnd = [];
        end
    end

    if ~isfield(DIF_struct_array, 'X_min_obs') || isempty(DIF_struct_array(i).X_min_obs)
        intrp.X_min_obs = 0;
    elseif intrp.X_min_obs==1
        if exist('X_min_obs','var')
            if any(size(X_min_obs)~=size(DIF_struct_array(i).X_min_obs))
                err = errordlg('Mismatch in dimensions of  the field ''X_min_obs''.', err_dlgname  );
                uiwait(err); interp_DIF_struct.error =1; return;
            end
        else
            X_min_obs = zeros(size(DIF_struct_array(i).X_min_obs));
            X_min_obs_bnd = [];
        end
    end

    if ~isfield(DIF_struct_array, 'X_max_est') || isempty(DIF_struct_array(i).X_max_est)
        intrp.X_max_est = 0;
    elseif intrp.X_max_est==1
        if exist('X_max_est','var')
            if any(size(X_max_est)~=size(DIF_struct_array(i).X_max_est))
                err = errordlg('Mismatch in dimensions of  the field ''X_max_est''.', err_dlgname  );
                uiwait(err); interp_DIF_struct.error =1; return;
            end
        else
            X_max_est = zeros(size(DIF_struct_array(i).X_max_est));
            X_max_est_bnd = [];
        end
    end

    if ~isfield(DIF_struct_array, 'X_min_est') || isempty(DIF_struct_array(i).X_min_est)
        intrp.X_min_est = 0;
    elseif intrp.X_min_est==1
        if exist('X_min_est','var')
            if any(size(X_min_est)~=size(DIF_struct_array(i).X_min_est))
                err = errordlg('Mismatch in dimensions of  the field ''X_min_est''.', err_dlgname  );
                uiwait(err); interp_DIF_struct.error =1; return;
            end
        else
            X_min_est = zeros(size(DIF_struct_array(i).X_min_est));
            X_min_est_bnd = [];
        end
    end
    
    if ~isfield(DIF_struct_array, 'load_max') || isempty(DIF_struct_array(i).load_max)
        intrp.load_max = 0;
    elseif intrp.load_max==1
        if exist('load_max','var')
            if any(size(load_max)~=size(DIF_struct_array(i).load_max))
                err = errordlg('Mismatch in dimensions of  the field ''load_max''.', err_dlgname  );
                uiwait(err); interp_DIF_struct.error =1; return;
            end
        else
            load_max = zeros(size(DIF_struct_array(i).load_max));
            load_max_bnd = [];
        end
    end

    if ~isfield(DIF_struct_array, 'load_min') || isempty(DIF_struct_array(i).load_min)
        intrp.load_min = 0;
    elseif intrp.load_min==1
        if exist('load_min','var')
            if any(size(load_min)~=size(DIF_struct_array(i).load_min))
                err = errordlg('Mismatch in dimensions of  the field ''load_min''.', err_dlgname  );
                uiwait(err); interp_DIF_struct.error =1; return;
            end
        else
            load_min = zeros(size(DIF_struct_array(i).load_min));
            load_min_bnd = [];
        end
    end
    
end    
interp_DIF_struct.resp_names = DIF_struct_array(1).resp_names;
interp_DIF_struct.resp_units = DIF_struct_array(1).resp_units;

% Perform the interpolation:

for i = 1:n_m
    if intrp.X_max_obs ==1
        X_max_obs = X_max_obs + DIF_struct_array(i).X_max_obs*wf(i);
        X_max_obs_bnd = cat(2,X_max_obs_bnd,DIF_struct_array(i).X_max_obs);
    end
    if intrp.X_min_obs ==1
        X_min_obs = X_min_obs + DIF_struct_array(i).X_min_obs*wf(i);
        X_min_obs_bnd = cat(2,X_min_obs_bnd,DIF_struct_array(i).X_min_obs);
    end
    if intrp.X_max_est ==1
        X_max_est = X_max_est + DIF_struct_array(i).X_max_est*wf(i);
        X_max_est_bnd = cat(2,X_max_est_bnd,DIF_struct_array(i).X_max_est);
    end
    if intrp.X_min_est ==1
        X_min_est = X_min_est + DIF_struct_array(i).X_min_est*wf(i);
        X_min_est_bnd = cat(2,X_min_est_bnd,DIF_struct_array(i).X_min_est);
    end
    if intrp.load_max ==1
        load_max = load_max + DIF_struct_array(i).load_max*wf(i);
        load_max_bnd = cat(3,load_max_bnd,DIF_struct_array(i).load_max);
    end
    if intrp.load_min ==1
        load_min = load_min + DIF_struct_array(i).load_min*wf(i);
        load_min_bnd = cat(3,load_min_bnd,DIF_struct_array(i).load_min);
    end
end

% Assemble the results in a data structure:
if intrp.X_max_obs ==1
    interp_DIF_struct.X_max_obs = X_max_obs;
    interp_DIF_struct.X_max_obs_bnd = X_max_obs_bnd;
end
if intrp.X_min_obs ==1
    interp_DIF_struct.X_min_obs = X_min_obs;
    interp_DIF_struct.X_min_obs_bnd = X_min_obs_bnd;
end
if intrp.X_max_est ==1
    interp_DIF_struct.X_max_est = X_max_est;
    interp_DIF_struct.X_max_est_bnd = X_max_est_bnd;
end
if intrp.X_min_est ==1
    interp_DIF_struct.X_min_est = X_min_est;
    interp_DIF_struct.X_min_est_bnd = X_min_est_bnd;
end
if intrp.load_max ==1
    interp_DIF_struct.load_max = load_max;
    interp_DIF_struct.load_max_bnd = load_max_bnd;
end
if intrp.load_min ==1
    interp_DIF_struct.load_min = load_min;
    interp_DIF_struct.load_min_bnd = load_min_bnd;
end