% By: Mohsen Kohrangi
% PhD candidate, UME School, Pavia, Italy
% 11/08/2015

% This code is an extension of the code for conditional mean spctrum based record
% selection based on AvgSA. The original code was generated by:

% Nirmal Jayaram, Ting Lin, Jack W. Baker
% Department of Civil and Environmental Engineering
% Stanford University

% here we modifed the code to perform the record selection of CS(AvgSA)
% based on the the following manuscript:

% Kohrangi, M., Bazzurro, P., Vamvatsikos, D., Spillatura, A., 2017. 
% Conditional Spectrum based ground motion record selection using average 
% spectral acceleration. Earthquake Engineering and Structural Dynamics. 
% DOI: 10.1002/eqe.2873.

% In the code, the lines that are modified are shown in the script. The code is
% based on the selection of an arbitrary component of the ground motion using the NGA west1
% record database. The code is based on the exact method which incorporates
% all the causal earthquakes in computation of the target spectrum.

% Additional input with respect to the origianl code to be used in this function:
% nGM: number of records to be selected
% T1_r: the vector of the periods to be included in AvgSA
% scenarios: avector of nx4 where:
            % n is the total number of bins considerred in disaggregation
            % first column, i-th row: Magnitude_bar of i-th bin
            % second column, i-th row: closest-rupture-to-site-distance of
            % i-th bin
            % third column, i-th row: epsilon_bar of i-th bin
            % fourth column, i-th row: the contribution of the bin in total
            % hazard
% outputFile: the name of mat output file
% HZ_SA: The value of the conditioning AvgSA

% Outputs of the function:
% SF: Final scale factors
% Fname: output mat file name
% meanReq: a vector of target spectrum with the mean+/-std in three columns


% Comments from the original code
%--------------------------------------------------------------------------
% This code is used to select conditional (structure- and site-
% specific) ground motions. The target means and covariances are
% obtained corresponding to a pre-defined target scenario earthquake, and
% are obtained based on the CMS method
%
% Nirmal Jayaram, Ting Lin, Jack W. Baker
% Department of Civil and Environmental Engineering
% Stanford University
% Last Updated: 27 March 2010
%
% Referenced manuscripts:
%
% N. Jayaram, T. Lin and and Baker, J. W. (2010). A computationally
% efficient ground-motion selection algorithm for matching a target
% response spectrum mean and variance, Earthquake Spectra, (in press).
%
% N. Jayaram and Baker, J. W. (2010). Ground-motion selection for PEER
% Transportation Systems Research Program, 7th CUEE and 5th ICEE Joint
% Conference, Tokyo, Japan.
%
% J. W. Baker and Jayaram, N. (2008). Correlation of spectral acceleration
% values from NGA ground motion models, Earthquake Spectra, 24 (1), 299-317
%
% Baker, J.W. (2010). The Conditional Mean Spectrum: A tool for ground
% motion selection, ASCE Journal of Structural Engineering, (in press).
%
%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% OUTPUT VARIABLES
% finalRecords      : Record numbers of selected records
% finalScaleFactors : Scale factors
%
% The final cell in this m file shows how to plot the selected spectra
% using this information.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Load workspace containing ground-motion information. Here, the NGA
% database is used. Documentation of the NGA database workspace
% 'rec_selection_meta_data.mat' can be found at 'WorkspaceDocumentation.m'.
% For an alternate database, the minimum information to be provided
% includes the pseudo-acceleration spectra of the available, ground
% motions, periods at which the spectra are defined, and other information
% required to compute means and variances using ground-motion models.
% This cell can be modified by the user if desired.
%
% Variable definitions
% saKnown   : (N*P matrix)
%             This is a matrix of Sa values at different periods (P) for
%             available ground-motion time histories (N).
% perKnown  : The set of P periods.
% nGM       : Number of ground motions to be selected
% T1        : Period at which spectra should be scaled and matched.
%             Usually, the structure's fundamental period.
% isScaled  : Should spectra be scaled before matching (1 -YES, 0-NO).
% maxScale  : Maximum allowable scale factor.
% weights   : [Weight for error in mean, Weight for error in standard
%             deviation] e.g., [1.0,1.0] - equal weight for both errors.
% nLoop     : This is a meta-variable of the algorithm. The greedy
%             improvement technique does several passes over the available
%             data set and tries to improve the selected records. This
%             variable denotes the number of passes. Recommended value: 2.
% penalty   : If a penalty needs to be applied to avoid selecting spectra
%             that have spectral acceleration values beyond 3 sigma at any
%             of the periods, set a value here. Use 0 otherwise.
% notAllowed: List of record numbers that should not be considered for
%             selection. Use [] to consider all available records. This can
%             be used to prevent certain magnitude-distance-Vs30 records
%             from being selected, if desired. (see example below)
% checkCorr : If 1, this runs a code that compares the correlation
%             structure of the selected ground motions to the correlations
%             published by Baker and Jayaram (2008).
% seedValue : For repeatability. For a particular seedValue not equal to
%             zero, the code will output the same set of ground motions.
%             The set will change when the seedValue changes. If set to
%             zero, the code randomizes the algorithm and different sets of
%             ground motions (satisfying the target mean and variance) are
%             generated each time.
% outputFile: File name of the output file
%
% If a database other than the NGA database is used, also define the
% following variables:
%
% magnitude        : Magnitude of all the records
% distance_closest : Closest distance for all the records
% soil_Vs30        : Soil Vs30 values corresponding to all the records
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User inputs begin here
function [SF, Fname, PerTgt, meanReq]=Select_GM_AvgSA(nGM,T1_r,scenarios, outputFile, HZ_SA)

load rec_selection_meta_data % Information from the NGA database
SaKnown    = sqrt(Sa_1.*Sa_2);
Filename = [Filename_1]; 
soil_Vs30 = [soil_Vs30];
% mechanism = [mechanism; mechanism];
magnitude = [magnitude];
perKnown   = Periods;

isScaled   = 1;
maxScale   = 4;
weights    = [1.0 2.0];
nLoop      = 2;
penalty    = 0;
% load('existing_records.mat')
% m=[url_ext(:,1);url_ext(:,2)];
notAllowed = [];
checkCorr  = 0;
seedValue  = 1;

% Limiting the records to be considered using the `notAllowed' variable
% Vs30_min=520;
% Vs30_max=720;
% Only hard rock
% recInvalid = find(soil_Vs30<Vs30_min | soil_Vs30>Vs30_max);
% notAllowed = [notAllowed; recInvalid;notAllowed; recInvalid];

% Limits on magnitude, distance
% recInvalid = find(magnitude<5 | magnitude>7 | distance_closest>200);
% notAllowed = [notAllowed; recInvalid];

% Using a different ground-motion database:

% MORE user input in the next cell
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Determination of target mean and covariances

% The Campbell and Bozorgnia (2008) ground-motion model is used in this
% code. The input variables defined below are the inputs required for this
% model. The user can change the ground-motion model as long as any
% additional information that may be required by the new model is provided.

% Please refer to Baker (2010) for details on the conditional mean spectrum
% method which is used for obtaining the target means and covariances in
% this example. Alternately, the details are summarized in Jayaram et al.
% (2010).

% The code provides the user an option to not match the target variance.
% This is done by setting the target variance to zero so that each selected
% response spectrum matches the target mean response spectrum.

% Variable definitions

% M_bar     : Magnitude of the target scenario earthquake
% R_bar     : Distance corresponding to the target scenario earthquake
% eps_bar   : Epsilon value for the target scenario
% Vs30      : Average shear wave velocity in the top 30m of the soil, used
%             to model local site effects (m/s)
% Ztor      : Depth to the top of coseismic rupture (km)
% delta     : Average dip of the rupture place (degree)
% lambda    : Rake angle (degree)
% Zvs       : Depth to the 2.5 km/s shear-wave velocity horizon (km)
% arb       : 1 for arbitrary component sigma
%             0 for average component sigma
% PerTgt    : Periods at which the target spectrum needs to be computed
% showPlots : 1 to display plots, 0 otherwise
% useVar    : 0 to set target variance to zero, 1 to compute target
%             variance using ground-motion model

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User inputs begin here
Vs30      = 620;
Ztor      = 0;
delta     = 60;
lambda    = 180;
Zvs       = 1; %2
arb       = 0;
PerTgt    = logspace(log10(0.1),log10(4.5),50);
showPlots = 1;
useVar    = 1;
Fault_Type=2;

% Modify perTgt to include T1
for ii=1:length(T1_r)
if ~any(PerTgt == T1_r(ii))
    PerTgt = [PerTgt(PerTgt<T1_r(ii)) T1_r(ii) PerTgt(PerTgt>T1_r(ii))]; % Modified     
end
end

% Determine target spectra using ground-motion model (replace ground-motion
% model code if desired)

sa = zeros(1,length(PerTgt));
sigma = zeros(1,length(PerTgt));

% covReq = zeros(length(PerTgt),length(PerTgt),length(scenarios(:,1))); 
covReq = zeros(length(PerTgt),length(PerTgt));
rho = zeros(1,length(PerTgt));

for kk=1:length(scenarios(:,1))
        
    M_bar=scenarios(kk,1);
    eps_bar=scenarios(kk,3);
    R_bar=scenarios(kk,2); Rjb=R_bar;
    
    for i = 1:length(PerTgt)
        %[sa(1,i) sigma(1,i)] = CB_2008_nga (M_bar, PerTgt(i), Rrup, Rjb, Ztor, delta, lambda, Vs30, Zvs, arb);
        [sa(1,i) sigma(1,i)] =BA_2008_nga (M_bar, PerTgt(i), Rjb, Fault_Type, Vs30);
    end
    
    % User inputs end here
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Use the CMS method to estimate target means and covariances
    % (Log) Response Spectrum Mean: meanReq

    for i = 1:length(PerTgt)
        rho(i) = rho_AvgSA_SA(PerTgt(i), T1_r, M_bar, Rjb, Fault_Type, Vs30); %Modified
    end
    
    [lnSa1, Avg_sig] = Sa_avg (M_bar, T1_r, Rjb, Fault_Type, Vs30) ;

    eps_bar=(log(HZ_SA)-log(lnSa1))/Avg_sig; %backcalculate epsilon to ensure matching with AvgSA
    meanReq(kk,:) = log(sa) + sigma.*eps_bar.*rho; % Modified

    if rem(kk/length(scenarios(:,1))*100,5) == 0 %**
        disp ( strcat (num2str(kk/length(scenarios(:,1))*100) , '% of scenarios computed') ) %**
    end %**
    
    % (Log) Response Spectrum Covariance: covReq, it remains constant with
    % common GMPEs. if u want to consider that variation you need to
    % uncomment the lines with the symbols #* and comment the %**
    
    for i=1:length(PerTgt)
        for j=1:length(PerTgt)

            % Periods
            Ti = PerTgt(i);
            Tj = PerTgt(j);

            % Means and variances
            rec1 = find(PerTgt == Ti);
            rec2 = find(PerTgt == Tj);
            var1 = sigma(rec1)^2;
            var2 = sigma(rec2)^2;

            %rec = find(PerTgt == T1); %Modified
            varT = Avg_sig^2;

            sigma11 = [var1 baker_jayaram_correlation(Ti, Tj)*sqrt(var1*var2);baker_jayaram_correlation(Ti, Tj)*sqrt(var1*var2) var2]; %Modified
            sigma22 = varT;%Modified
            sigma12 = [rho(i)*sqrt(var1*varT);rho(j)*sqrt(var2*varT)]; %Modified

            sigmaCond = sigma11 - sigma12*inv(sigma22)*(sigma12)';

            covReq(i,j) = sigmaCond(1,2); %**
            
            if useVar == 0
                covReq(i,j) = 0.0; %**
            end

        end
    end
cov_mat_all(:,:,kk)=covReq;
end%**
%     if rem(kk/length(scenarios(:,1))*100,5) == 0  % #*
%     disp ( strcat (num2str(kk/length(scenarios(:,1))*100) , '% of scenarios computed') ) % #*
%     end % #*
%     
% end % #*

% length(meanReq(:,1))
% length(meanReq(1,:))
meanReq_fin=sum(meanReq.*repmat(scenarios(:,4),1,length(meanReq(1,:))),1);
cov_final=cov_mat_all(:,:,1); %all 2D matrices are the same for each kk scenario, since sigma is only T dependent
for ii=1:length(PerTgt)
    for kk=1:length(scenarios(:,1))
        
        cov_elms(ii,kk)= (cov_mat_all(ii,ii,kk)+(meanReq(kk,ii)-meanReq_fin(ii))^2) * scenarios(kk,4); %scenarios(kk,4)= contribution of the k-th scenario

    end
end
cov_diag=sum(cov_elms,2);
cov_final(eye(length(PerTgt))==1)=cov_diag;

id_avg=ismember(PerTgt,T1_r); %looking for the Period id of the range where we have computed the Sa_avg
lnSa1_fin=mean( meanReq_fin(1,id_avg));
% lnSa1_fin=1/length*sum(Sa(Ti); %mean di "meanReq_fin" per indici di T1_r poi sostituisci lnSa1 con lnSa1_fin
covReq_fin=cov_final;
%% Simulate response spectra using Monte Carlo Simulation

% 20 sets of response spectra are simulated and the best set (in terms of
% matching means, variances and skewness is chosen as the seed). The user
% can also optionally rerun this segment multiple times before deciding to
% proceed with the rest of the algorithm. It is to be noted, however, that
% the greedy improvement technique significantly improves the match between
% the means and the variances subsequently.

nTrials = 20;
% Setting initial seed for simulation
if seedValue ~= 0
    randn('seed',seedValue);
else
    randn('seed',sum(100*clock));
end
devTotalSim = zeros(nTrials,1);
for j=1:nTrials
    gmCell{j} = zeros(nGM,length(meanReq_fin));
    for i=1:nGM
        gmCell{j}(i,:) = exp(mvnrnd(meanReq_fin,covReq_fin));
    end
    devMeanSim = mean(log(gmCell{j})) - meanReq_fin;
    devSkewSim = skewness(log(gmCell{j}),1);
    devSigSim = std(log(gmCell{j})) - sqrt(diag(covReq_fin))';
    devTotalSim(j) = weights(1) * sum(devMeanSim.^2) + weights(2) * sum(devSigSim.^2)+ 0.1 * (weights(1)+weights(2)) * sum(devSkewSim.^2);
end
[tmp recUse] = min(abs(devTotalSim));
gm = gmCell{recUse};

if showPlots == 1
    
    % Plot simulated response spectra
    figure
    loglog(PerTgt, exp(meanReq_fin), '-r', 'linewidth', 3)
    hold on
    loglog(PerTgt, exp(meanReq_fin + 1.96*sqrt(diag(covReq_fin))'), '--r', 'linewidth', 3)
    for i=1:nGM
        loglog(PerTgt,gm,'k');
    end
    loglog(PerTgt, exp(meanReq_fin - 1.96*sqrt(diag(covReq_fin))'), '--r', 'linewidth', 3)
    axis([min(PerTgt) max(PerTgt) 1e-2 5])
    xlabel('T (s)')
    ylabel('S_a (g)')
    legend('Median response spectrum','2.5 and 97.5 percentile response spectra','Response spectra of simulated ground motions')
    title('Response spectra of simulated ground motions')
    
    % Plot target and simulated means
    figure
    loglog(PerTgt,exp(meanReq_fin))
    hold on
    loglog(PerTgt,exp(mean(log(gm))),'--')
    axis([min(PerTgt) max(PerTgt) 1e-2 5])
    xlabel('T (s)')
    ylabel('Median S_a (g)')
    legend('exp(Target mean lnS_a)','exp(Mean of simulated lnS_a)')
    title('Target and sample exponential logarithmic means (i.e., medians)')
    
    % Plot target and simulated standard deviations
    figure
    semilogx(PerTgt,sqrt(diag(covReq_fin))')
    hold on
    semilogx(PerTgt,std(log(gm)),'--')
    axis([min(PerTgt) max(PerTgt) 0 1])
    xlabel('T (s)')
    ylabel('Standard deviation of lnS_a')
    legend('Target standard deviation of lnS_a','Standard deviation of simulated lnS_a')
    title('Target and sample logarithmic standard deviations')
    
end
%% Arrange the available spectra in a usable format and check for invalid
% input
% Match periods (known periods and periods for error computations)
recPer = zeros(length(PerTgt),1);
for i=1:length(PerTgt)
    [tmp recPer(i)] = min(abs(perKnown - PerTgt(i)));
end

% Check for invalid input
sampleBig = SaKnown(:,recPer);
if (any(any(isnan(sampleBig))))
    error ('NaNs found in input response spectra')
end

% Processing available spectra
sampleBig = log(sampleBig);
nBig = size(sampleBig,1);

%% Find best matches to the simulated spectra from ground-motion database

recID = zeros(nGM,1);
sampleSmall = [];
finalScaleFac = ones(nGM,1);
for i = 1:nGM
    err = zeros(nBig,1);
    
    scaleFac = ones(nBig,1);
    for j=1:nBig
        
        ind=ismember(PerTgt, T1_r);
        rec_Avg=exp(sum(sampleBig(j,ind))/(length(T1_r)));
        
        if (isScaled == 1)
            
            %if exp(sampleBig(j,PerTgt == T1)) == 0 %Modified
            if   rec_Avg==0 % we could do it by ismember function
                scaleFac(j) = -1;
                err(j) = 1000000;
            else
                scaleFac(j) = exp(lnSa1_fin)/rec_Avg;% Modified
                if (scaleFac(j) > maxScale | soil_Vs30(j)==-1 | any(notAllowed==j))
                    err(j) = 1000000;
                else
                    err(j) = sum((log(exp(sampleBig(j,:))*scaleFac(j)) - log(gm(i,:))).^2); 
                end
            end
        else
            if (soil_Vs30(j)==-1 | any(notAllowed==j))
                err(j) = 1000000;
            else
                err(j) = sum((sampleBig(j,:) - log(gm(i,:))).^2);
                if err(j) == inf
                    err(j) = 1000000;
                end
            end
        end
        if (any(recID == j))
            err(j) = 1000000;
        end
    end
    [tmp recID(i)] = min(err);
    if tmp >= 1000000
        display('Warning: Possible problem with simulated spectrum. No good matches found');
        display(recID(i));
    end
    if (isScaled == 1)
        finalScaleFac(i) = scaleFac(recID(i));
    else
        finalScaleFac(i) = 1;
    end
    sampleSmall = [sampleSmall;log(exp(sampleBig(recID(i),:))*scaleFac(recID(i)))];
    
end
%% Greedy subset modification procedure

display('Please wait...This algorithm takes a few minutes depending on the number of records to be selected');

for k=1:nLoop % Number of passes
    
    for i=1:nGM % Selects nSelect ground motions
        
        display([num2str(round(((k-1)*nGM + i-1)/(nLoop*nGM)*100)) '% done']);
        
        minDev = 100000;
        
        sampleSmall(i,:) = [];
        recID(i,:) = [];
        
        % Try to add a new spectra to the subset list
        for j=1:nBig
           % rec_Avg2=exp(sum(sampleBig(j,PerTgt >= min(T1_r) & PerTgt <= max(T1_r)))/(length(T1_r)));
        ind2=ismember(PerTgt, T1_r);
        rec_Avg2=exp(sum(sampleBig(j,ind2))/(length(T1_r)));
         
           if isScaled == 1
                
                if rec_Avg2 == 0 % Modified
                    
                    scaleFac(j) = 1000000;
                else
                    scaleFac(j) = exp(lnSa1_fin)/rec_Avg2; %Modified
                end
                sampleSmall = [sampleSmall;sampleBig(j,:)+log(scaleFac(j))];
            else
                sampleSmall = [sampleSmall;sampleBig(j,:)];
                scaleFac(j) = 1;
            end
            % Compute deviations from target
            devMean = mean(sampleSmall) - meanReq_fin;
            devSkew = skewness(sampleSmall,1);
            devSig = std(sampleSmall) - sqrt(diag(covReq_fin))';
            %devTotal = weights(1) * sum(devMean.^2) + weights(2) * sum(devSig.^2)+ weights(3) * sum(devSkew.^2);
            devTotal = weights(1) * sum(devMean.^2) + weights(2) * sum(devSig.^2);
            
            % Penalize bad spectra (set penalty to zero if this is not required)
            for m=1:size(sampleSmall,1)
                devTotal = devTotal + sum(abs(exp(sampleSmall(m,:))>exp(meanReq_fin+3*sqrt(diag(covReq_fin))'))) * penalty;
            end
            
            if (scaleFac(j) > maxScale | soil_Vs30(j)==-1 | any(notAllowed==j))
                devTotal = devTotal + 1000000;
            end
            
            % Should cause improvement and record should not be repeated
            if (devTotal < minDev && ~any(recID == j))
                minID = j;
                minDev = devTotal;
            end
            sampleSmall = sampleSmall(1:end-1,:);
        end
        
        % Add new element in the right slot
        if isScaled == 1
            finalScaleFac(i) = scaleFac(minID);
        else
            finalScaleFac(i) = 1;
        end
        sampleSmall = [sampleSmall(1:i-1,:);sampleBig(minID,:)+log(scaleFac(minID));sampleSmall(i:end,:)];
        recID = [recID(1:i-1);minID;recID(i:end)];
    end
end

display('100% done');

% Output information
finalRecords = recID;
finalScaleFactors = finalScaleFac;

%% Spectra Plots

if (showPlots)
    
    % Variables used here
    
    % SaKnown    : As before, it contains the response spectra of all the
    %              available ground motions (N*P matrix) - N ground motions,
    %              P periods
    % sampleBig  : Same as SaKnown, but is only defined at PerTgt, the
    %              periods at which the target response spectrum properties
    %              are computed
    % sampleSmall: The response spectra of the selected ground motions,
    %              defined at PerTgt
    % meanReq_fin    : Target mean for the (log) response spectrum
    % covReq_fin     : Target covariance for the (log) response spectrum
    %
    
    
    % Plot at all periods
    figure
    loglog(PerTgt, exp(meanReq_fin), 'b', 'linewidth', 3)
    hold on
    loglog(PerTgt, exp(meanReq_fin + 1.96*sqrt(diag(covReq_fin))'), '--b', 'linewidth', 3)
    loglog(perKnown,SaKnown(finalRecords,:).*repmat(finalScaleFactors,1,size(SaKnown,2)),'k');
    loglog(PerTgt, exp(meanReq_fin - 1.96*sqrt(diag(covReq_fin))'), '--b', 'linewidth', 3)
    axis([min(PerTgt) max(PerTgt) 1e-2 5])
    xlabel('T (s)');
    ylabel('S_a (g)');
    legend('Median response spectrum','2.5 and 97.5 percentile response spectra','Response spectra of selected ground motions');
    title ('Response spectra of selected ground motions');
    
    % Plot spectra only at periods where error is minimized
    figure
    loglog(PerTgt, exp(meanReq_fin), 'b', 'linewidth', 3)
    hold on
    loglog(PerTgt, exp(meanReq_fin + 1.96*sqrt(diag(covReq_fin))'), '--b', 'linewidth', 3)
    loglog(PerTgt,exp(sampleBig(finalRecords,:)).*repmat(finalScaleFactors,1,length(PerTgt)),'color',[0.5 0.5 0.5],'linewidth',1)
    loglog(PerTgt, exp(meanReq_fin - 1.96*sqrt(diag(covReq_fin))'), '--b', 'linewidth', 3)
    axis([min(PerTgt) max(PerTgt) 1e-2 5])
    xlabel('T (s)');
    ylabel('S_a (g)');
    legend('Median response spectrum','2.5 and 97.5 percentile response spectra','Response spectra of selected ground motions');
    title ('Response spectra of selected ground motions at periods where error is minimized');
    
    % Sample and target means
    figure
    loglog(PerTgt,exp(meanReq_fin),'k','linewidth',1)
    hold on
    loglog(PerTgt,exp(mean(sampleSmall)),'b--','linewidth',1)
    axis([min(PerTgt) max(PerTgt) 1e-2 5])
    xlabel('T (s)')
    ylabel('Median S_a (g)')
    legend('exp(Target mean lnS_a)','exp(Mean of selected lnS_a)')
    title('Target and sample exponential logarithmic means (i.e., medians)')
    
    % Sample and target standard deviations
    figure
    semilogx(PerTgt,sqrt(diag(covReq_fin))','k','linewidth',1)
    hold on
    semilogx(PerTgt,std(sampleSmall),'b--','linewidth',1)
    axis([min(PerTgt) max(PerTgt) 0 1])
    xlabel('T (s)')
    ylabel('Standard deviation of lnS_a')
    legend('Target standard deviation of lnS_a','Standard deviation of selected lnS_a')
    title('Target and sample logarithmic standard deviations')
    
end

if (checkCorr)
    conditionalCovariance
end
%% Output data to file (best viewed with textpad)

fin = fopen(outputFile,'w');
fprintf(fin,'%s \t %s \t %s \t %s \t %s \t %s \t %s \n','Record Number','NGA Record Sequence Number','Scale Factor', 'Filename');
for i = 1 : length(finalRecords)
    rec = finalRecords(i);
    fprintf(fin,'%d \t %d \t %6.2f \t %s \t %s \t %s \t %s \n',i,rec,finalScaleFactors(i),Filename{rec});
    Fname{i} =  strcat(Filename{rec}(1:end-3), 'at2');
    SF(i) = finalScaleFactors (i);
end
fclose(fin);

meanrReq_minus = exp(meanReq_fin - 1.96*sqrt(diag(covReq_fin))');
meanReq_plus     = exp(meanReq_fin + 1.96*sqrt(diag(covReq_fin))');

meanReq= exp(meanReq_fin');
meanReq(:,2) = meanrReq_minus;
meanReq(:,3) = meanReq_plus;
end