computemripeaks.m

低矮房屋风压系数、风荷载计算分析matlab程序

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