program2.m

关于高层建筑强风条件下风压系数计算

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Program2.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Calculates the mean-recurrence intervals (MRI) of wind effects for each member

global  textA axesA_text textB axesB_text input_coords pos_on pos_off fig
% 1. Building Information
global  Nfloors H_bldg  Fdofs  Ndofs   RNmodes  AnL_A  AnL_B  anal
global  T  T_str  D  D_str  flnMode_shape  evectors  omega  h
%% 2. Wind Tunnel Information
global  Ntaps  flnTaps  flnCp  Vm  freq  ms  Npoints
% 3. Program Inputs
global  flnMem  flDif  flnProps  mem_list  props  WS  WS_str  WD  WD_str
global  MRI  MRI_str  mri  IMdir_A  IMdir_B  intmeth  mem_out
% 4. Loads and Reactions
global  flnMass  DL  SDL  LL
global  flnDLr  flnSDLr  flnLLr  frames_DL  frames_SDL  frames_LL
global  DLf  SDLf  LLf  WLf
% 5. Wind Effect Selection
global  WLE  Hmp  Hfile  Vth  rf  tstep  Npointsout
global  Qs  flHfile
% 6. Calculations
global  saveBij  flnSaveBij  flnSavemaxBij  flnSaveMRI
% 7. Output
global  AnL_A  AnL_B  IMdir_A  IMdir_B  Hmp  Vth  DLf  SDLf  LLf  WLf
global  mem_out  mem_list
% 8. Save/Load Data
global  fnSAVE  flSAVE  flnSAVE  fnLOAD  flLOAD  flnLOAD

%%%%%%%%%%%%%%%%%%%%%%%% Member list : mem_list, Page 3, [1 x any] %%%%%%%%%%%%%%%%%%%%%%%%%%%

w = waitbar(0,'Running Program2.m   Please wait...'); c=0;fighan=w;

%%%%%%%%%%%%%%%%%%%%%%%% Member list : mem_list, Page 3, [1 x any] %%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop through each member; load the saved minimum, maximum, observed, estimated peak wind
% effect matrix; map each wind speed in the hurricane database to a peak wind effect from
% each matrix; calculate the mean recurrence interval
for p=1:length(mem_list)
    mem = mem_list(p);
    load ([flnSavemaxBij '_' num2str(mem)]);
    for axpos_index = 1:3
        eval(['max_Bij_obs = Bij_obs_max_' num2str(mem) '(:,:,axpos_index);']);
        eval(['min_Bij_obs = Bij_obs_min_' num2str(mem) '(:,:,axpos_index);']);
        eval(['max_Bij_est = Bij_est_max_' num2str(mem) '(:,:,axpos_index);']);
        eval(['min_Bij_est = Bij_est_min_' num2str(mem) '(:,:,axpos_index);']);
        % This function reads in the NIST hurricane wind speed matrix
        [V_lookup,dir_lookup,nhurr,recurrence] = HurrData;
        for a=1:nhurr
            for b=1:size(V_lookup,2)
                if V_lookup(a,b) > Vth             % Impose minimum threshold of wind speeds, Vth
                    V_kept(a,b)=V_lookup(a,b);
                    if IMdir_A == 1                % Interpolation Method A: between WS and WD
                        if length(WD)==1
                            Bij_Omax_lookup(a,b) = interp1(WS,max_Bij_obs,V_kept(a,b));
                            Bij_Omin_lookup(a,b) = interp1(WS,min_Bij_obs,V_kept(a,b));
                            Bij_Emax_lookup(a,b) = interp1(WS,max_Bij_est,V_kept(a,b));
                            Bij_Emin_lookup(a,b) = interp1(WS,min_Bij_est,V_kept(a,b));
                        else
                            Bij_Omax_lookup(a,b) = interp2(WS,WD,max_Bij_obs,V_kept(a,b),dir_lookup(b));
                            Bij_Omin_lookup(a,b) = interp2(WS,WD,min_Bij_obs,V_kept(a,b),dir_lookup(b));
                            Bij_Emax_lookup(a,b) = interp2(WS,WD,max_Bij_est,V_kept(a,b),dir_lookup(b));
                            Bij_Emin_lookup(a,b) = interp2(WS,WD,min_Bij_est,V_kept(a,b),dir_lookup(b));
                        end
                    elseif IMdir_B == 1           % Interpolation Method B: between WS, larger WD
                            for u=1:size(max_Bij_obs,1)
                                Bij_Omax_V(u)        = interp1(WS,max_Bij_obs(u,:),V_kept(a,b));
                                Bij_Omin_V(u)        = interp1(WS,min_Bij_obs(u,:),V_kept(a,b));
                                Bij_Emax_V(u)        = interp1(WS,max_Bij_est(u,:),V_kept(a,b));
                                Bij_Emin_V(u)        = interp1(WS,min_Bij_est(u,:),V_kept(a,b));
                            end
                            D_index                  = 1:size(max_Bij_obs,1);   % 1:36 to match table directions
                            D_V_index = interp1(WD,D_index,dir_lookup(b));
                            Bij_Omax_lookup(a,b) = max(Bij_Omax_V(floor(D_V_index)),Bij_Omax_V(ceil(D_V_index)));
                            Bij_Omin_lookup(a,b) = min(Bij_Omin_V(floor(D_V_index)),Bij_Omin_V(ceil(D_V_index)));
                            Bij_Emax_lookup(a,b) = max(Bij_Emax_V(floor(D_V_index)),Bij_Emax_V(ceil(D_V_index)));
                            Bij_Emin_lookup(a,b) = min(Bij_Emin_V(floor(D_V_index)),Bij_Emin_V(ceil(D_V_index)));
                    end     % loop on method 1 or 2
                    if Bij_Emin_lookup(a,b) < 0
                        Bij_Emin_lookup(a,b) = 0;
                    end
                else
                    V_kept(a,b)=0;
                    Bij_Omax_lookup(a,b)     = 0;
                    Bij_Omin_lookup(a,b)     = 0;
                    Bij_Emax_lookup(a,b)     = 0;
                    Bij_Emin_lookup(a,b)     = 0;
                end               % loop V_lookup > Vth

            end                 % loop on 'b'
            % Max Bij for each event (999 in total)
            Bij_Obs_max_event(a)         = max(Bij_Omax_lookup(a,:));
            Bij_Obs_min_event(a)         = min(Bij_Omin_lookup(a,:));
            Bij_Est_max_event(a)         = max(Bij_Emax_lookup(a,:));
            Bij_Est_min_event(a)         = min(Bij_Emin_lookup(a,:));
        end                   % loop on 'a'

        % Rank order the vector to prepare for MRI calculation
        sorted_Bij_Obs_max_mat(:,axpos_index)      = -sort( -Bij_Obs_max_event );
        sorted_Bij_Obs_min_mat(:,axpos_index)      =  sort(  Bij_Obs_min_event );
        sorted_Bij_Est_max_mat(:,axpos_index)      = -sort( -Bij_Est_max_event );
        sorted_Bij_Est_min_mat(:,axpos_index)      =  sort(  Bij_Est_min_event );
        
        for i = 1:length(MRI)
            rank                = 999 / recurrence / MRI(i);
            Mem_Bij_Obs_max(i)  = interp1(sorted_Bij_Obs_max_mat(:,axpos_index), rank);
            Mem_Bij_Obs_min(i)  = interp1(sorted_Bij_Obs_min_mat(:,axpos_index), rank);
            Mem_Bij_Est_max(i)  = interp1(sorted_Bij_Est_max_mat(:,axpos_index), rank);
            Mem_Bij_Est_min(i)  = interp1(sorted_Bij_Est_min_mat(:,axpos_index), rank);
        end

        Mem_Bij_Obs_max_mat(:,axpos_index) = Mem_Bij_Obs_max';
        Mem_Bij_Obs_min_mat(:,axpos_index) = Mem_Bij_Obs_min';
        Mem_Bij_Est_max_mat(:,axpos_index) = Mem_Bij_Est_max';
        Mem_Bij_Est_min_mat(:,axpos_index) = Mem_Bij_Est_min';
    end  %loop on member axial position index (axpos_index)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SAVE MRI FOR EACH MEMBER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Save the sorted vectors and all associated MRI's for each member using 'flnSaveMRI'
        K3 = findstr(flnSaveMRI,'.mat');
        if ~isempty(K3)
            flnSaveMRI = flnSaveMRI(1:K3-1);
        end
        save ([flnSaveMRI '_' num2str(mem) ],...
            'sorted_Bij_Obs_max_mat','sorted_Bij_Obs_min_mat','sorted_Bij_Est_max_mat','sorted_Bij_Est_min_mat',...
            'Mem_Bij_Obs_max_mat','Mem_Bij_Obs_min_mat','Mem_Bij_Est_max_mat','Mem_Bij_Est_min_mat',...
            'MRI', 'recurrence');
    c=c+1;
    counter=c/( length(mem_list)  );
    waitbar(counter,fighan);
end     % loop on member
close(fighan)