function [finalRecords, finalDirectios, finalScaleFactors,...
          F_x, F1_y1, F1_y1a, F1_y2, F1_y3, F1_yn, F2_y1, F3_y1, F3_y2,...  
          F4_x1, F4_y1, F5_y1, F6_y1, F7_y1, errCode] = Select_Ground_Motions(...  
    idRef, Sa_1, Sa_2, sVs30, ...  
    nGM, T1, T1dej, ecSoilCode, isScaled, maxScale, weight_mean, weight_sd, nLoop, penalty,...
    seedValue, M_bar, R_bar, eps_bar, nX, useVar) 
%compileReplaceIncludes 
errCode = 0; 
%finalRecords_tmp=0; 
finalDirectios = 0; 
coder.varsize('recID'); 
%coder.varsize('notAllowed'); 
coder.extrinsic('clock'); 
recID = zeros(round(nGM),1); 
% 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. The target means and variances can be 
% modified by the user. The ground motion selection algorithm used here is 
% based on a single arbitrary ground motion component (for two-dimensional 
% structural models), instead of the geometric mean of two ground motion 
% components (for three-dimensional structural models), as described in 
% Jayaram et al. (2010). This is used to select ground motions for
% two-dimensional structural models, where a single ground motion component
% is required as an input for every time history analysis.
%
% Nirmal Jayaram, Ting Lin, Jack W. Baker
% Department of Civil and Environmental Engineering
% Stanford University
% Last Updated: 7 June 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

S=coder.load('rec_selection_meta_data','Periods'); % Information from the NGA database
% Each horizontal component of a recording from the same station in the NGA
% database is treated separately as an individual ground motion.
SaKnown = [Sa_1; Sa_2];
%%%
RSN = [(1:size(Sa_1,1))';(1:size(Sa_1,1))'];
%Filename = [Filename_1; Filename_2]; 
soil_Vs30 = [sVs30; sVs30];
% mechanism = [mechanism; mechanism];
%magnitude = [magnitude_; magnitude_];
% lowest_usable_freq = [lowest_usable_freq; lowest_usable_freq];
% distance_jb = [distance_jb; distance_jb];
% distance_hyp = [distance_hyp; distance_hyp];
% distance_epi = [distance_epi; distance_epi];
% distance_closest = [distance_closest_; distance_closest_];
% distance_campbell = [distance_campbell; distance_campbell];
%%%

perKnown   = S.Periods;
%nGM        = 10;
%T1         = 2.63;
%isScaled   = 1;
%maxScale   = 4;
weights    = [weight_mean weight_sd];
%nLoop      = 2;
%penalty    = 0;
notAllowed = [];
%checkCorr  = 1;
%seedValue  = 1;
%outputFile = 'Output_File.dat';
% NOTE: MORE user input required in the next cell

% if ecSoilCode==1
%     S1=0; % is there shallow alluvium
%     S2=0; % is there deep alluvium
%     coef=4; % 1-fault dist,raw coef; 2-fault dist,smooth coef 3-epicentral dist,raw coef; 4-epicentral dist,smooth coef
% elseif ecSoilCode==2
%     S1=0 % is there shallow alluvium
%     S2=1 % is there deep alluvium
%     coef=4 % 1-fault dist,raw coef; 2-fault dist,smooth coef 3-epicentral dist,raw coef; 4-epicentral dist,smooth coef  
% elseif ecSoilCode==3
%     S1=1; % is there shallow alluvium
%     S2=0; % is there deep alluvium
%     coef=4; % 1-fault dist,raw coef; 2-fault dist,smooth coef 3-epicentral dist,raw coef; 4-epicentral dist,smooth coef  
% else
%     S1=0;% is there shallow alluvium
%     S2=1; % is there deep alluvium
%     coef=4; % 1-fault dist,raw coef; 2-fault dist,smooth coef 3-epicentral dist,raw coef; 4-epicentral dist,smooth coef  
% end

% Limit records without acc. data.
% List of missing records: Usability testing of application programming
% interface for SAP2000, appendix B; Žiga Šebenik 2013
%notAllowed = [notAllowed; 260; 930; 1171; 1298; 2136];
%notAllowed = [notAllowed; 29; 168; 177; 267; 361; 365; 378; 473; 474; 475; 500; 599; 660; 805; 1635; 1753; 1774;];
%notAllowed = [notAllowed; 43; 61; 78; 98; 104; 105; 106; 256; 275; 361; 390; 425; 436; 437; 438; 439; 442; 443; 444; 446; 496; 517; 530; 532; 651; 719; 720; 722; 723; 724; 725; 726; 727; 803; 822; 973; 1048; 1122; 1132; 1134; 1299; 1591; 1628; 1629; 1694; 1696; 1701; 1720; 1732; 3549; 3550; 3551];

% Limiting the records to be considered using the `notAllowed' variable

% 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

%M_bar     = 7;
%R_bar     = 10;
%eps_bar   = 2;
%Vs30      = 400;
%Ztor      = 0;
%delta     = 90;
%lambda    = 180;
%Zvs       = 2;
%arb       = 1; % Arbitrary component
PerTgt_tmp    = linspace(0,4,nX);
%showPlots = 0;
%useVar    = 1;

% Modify perTgt to include T1
if ~any(PerTgt_tmp == T1)
    PerTgt_tmp2 = zeros(1,size(PerTgt_tmp,2)+1);
    PerTgt_tmp2 = [PerTgt_tmp(PerTgt_tmp<T1) T1 PerTgt_tmp(PerTgt_tmp>T1)];
else
    PerTgt_tmp2 = PerTgt_tmp;
end
% Modify perTgt to include T1dej
if ~any(PerTgt_tmp2 == T1dej)
    PerTgt = zeros(1,size(PerTgt_tmp2,2)+1);
    PerTgt = [PerTgt_tmp2(PerTgt_tmp2<T1dej) T1dej PerTgt_tmp2(PerTgt_tmp2>T1dej)];
else
    PerTgt = PerTgt_tmp2;
end

% Find position of T1 in PerTgt
posT1=find(PerTgt==T1);

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

% sa = zeros(1,length(PerTgt));
% sigma = zeros(1,length(PerTgt));
% for i = 1:length(PerTgt)
%     [temp1, temp2] = SP_1996(PerTgt(i), M_bar, R_bar, S1, S2, coef);
%     sa(1,i)=temp1(1);
%     sigma(1,i)=temp2(1);
% end

sa    = zeros(1,length(PerTgt));
sigma = zeros(1,length(PerTgt));
for i = 1:length(PerTgt)
    %[saSP(1,i) sigmaSP(1,i)] = SP_1996(PerTgt(i), M_mean, R_mean, S1, S2, coef);
    [koefAtonacijsko]=SP1996hor(2,PerTgt(i),ecSoilCode);
    [temp1,temp2]=AtonacijskoTip1(M_bar,R_bar,koefAtonacijsko);
    sa(1,i)=temp1(1);
    sigma(1,i)=temp2(1);    
end
sa=exp(sa);


% User inputs end here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Use the CMS method to estimate target means and covariances

% (Log) Response Spectrum Mean: meanReq
rho = zeros(1,length(PerTgt));
for i = 1:length(PerTgt)
    rho(i) = baker_jayaram_correlation(PerTgt(i), T1);
end
meanReq = log(sa) + sigma.*eps_bar.*rho;
lnSa1 = meanReq(PerTgt == T1);

% (Log) Response Spectrum Covariance: covReq
covReq = zeros(length(PerTgt));
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);
        varT = sigma(rec)^2;
        
        sigma11 = [var1 baker_jayaram_correlation(Ti, Tj)*sqrt(var1*var2);baker_jayaram_correlation(Ti, Tj)*sqrt(var1*var2) var2];
        sigma22 = varT;
        sigma12 = [baker_jayaram_correlation(Ti, T1)*sqrt(var1*varT);baker_jayaram_correlation(T1, Tj)*sqrt(var2*varT)];
        
        sigmaCond = sigma11 - sigma12*inv(sigma22)*(sigma12)';
        
        covReq(i,j) = sigmaCond(1,2);
        
        if useVar == 0
            covReq(i,j) = 0.0;
        end
        
    end
end
%% 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
    rng(seedValue,'v4');
else
    rng(sum(100*clock),'v4');
end
%compileReplaceClock
devTotalSim = zeros(nTrials,1);
tempf=zeros(nGM,size(PerTgt,2),nTrials);

for j=1:nTrials
    for i=1:nGM
        tempf(i,:,j)=real(exp(mvnrnd_(meanReq,covReq)));
    end
    devMeanSim = mean(log(tempf(:,:,j))) - meanReq;
    % Round to 10 digits (to get the same results as in Matlab) and replace
    % NaNs with 0; zsebenik
    devSkewSim = skew(round(log(tempf(:,:,j))*1e10)/1e10, 1);
    devSkewSim(isnan(devSkewSim)) = 0;
    devSigSim = std(log(tempf(:,:,j))) - sqrt(diag(covReq))';
    devTotalSim(j) = weights(1) * sum(devMeanSim.^2,2) + weights(2) * sum(devSigSim.^2,2)+ 0.1 * (weights(1)+weights(2)) * sum(devSkewSim.^2,2);
end
[tmp, recUse] = min(abs(devTotalSim));
gm = tempf(:,:,recUse);

% if showPlots == 1
F_x = PerTgt;
F1_y1 = exp(meanReq);
F1_y1a = meanReq;
F1_y2 = exp(meanReq + 1.96*sqrt(diag(covReq))');
F1_yn = gm;
F1_y3 = exp(meanReq - 1.96*sqrt(diag(covReq))');
F2_y1 = exp(mean(log(gm)));
F3_y1 = sqrt(diag(covReq))';
F3_y2 = std(log(gm));
%     
%     % Plot simulated response spectra
%     figure
%     loglog(PerTgt, exp(meanReq), '-r', 'linewidth', 3)
%     hold on
%     loglog(PerTgt, exp(meanReq + 1.96*sqrt(diag(covReq))'), '--r', 'linewidth', 3)
%     for i=1:nGM
%         loglog(PerTgt,gm,'k');
%     end
%     loglog(PerTgt, exp(meanReq - 1.96*sqrt(diag(covReq))'), '--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))
%     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))')
%     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))))
    errCode=3390;   %'NaNs found in input response spectra'   
    finalRecords=0;
    finalScaleFactors=0;
    F_x=0;
    F1_y1=0;
    F1_y1a = 0;
    F1_y2=0;
    F1_y3=0;
    F1_yn=0;
    F2_y1=0;
    F3_y1=0;
    F3_y2=0;
    F4_x1=perKnown;
    F4_y1=0;
    F5_y1=0;
    F6_y1=0;
    F7_y1=0;
    return;
end

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

%% Find best matches to the simulated spectra from ground-motion database
err = zeros(nBig,1);
sampleSmall = zeros(nGM,size(PerTgt,2));
finalScaleFac = ones(nGM,1);
for i = 1:nGM
    err = zeros(nBig,1);
    
    scaleFac = ones(nBig,1);
    for j=1:nBig
        if (isScaled == 1)
            
            if all(exp(sampleBig(j,posT1)) == 0)
                scaleFac(j) = -1;
                err(j) = 1000000;
            else
                scaleFac(j) = exp(lnSa1)/exp(sampleBig(j,PerTgt == T1));
                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
        errCode=3389;    %d'Warning: Possible problem with simulated spectrum. No good matches found'
    end
    if (isScaled == 1)
        finalScaleFac(i) = scaleFac(recID(i));
    else
        finalScaleFac(i) = 1;
    end
    sampleSmall(i,:) = log(exp(sampleBig(recID(i),:))*scaleFac(recID(i)));
    
end

if (all(recID==1) && all(err==1000000))
    errCode=3391;   %'Error: no match spectra found'
    finalRecords=0;
    finalScaleFactors=0;
    F_x=0;
    F1_y1=0;
    F1_y1a = 0;
    F1_y2=0;
    F1_y3=0;
    F1_yn=0;
    F2_y1=0;
    F3_y1=0;
    F3_y2=0;
    F4_x1=perKnown;
    F4_y1=0;
    F5_y1=0;
    F6_y1=0;
    F7_y1=0;
    return;
end
%% Greedy subset modification procedure

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

minID=nan;
sdCovReq = sqrt(diag(covReq))';
emrsdc = exp(meanReq+3*sdCovReq);
for k=1:nLoop % Number of passes
    
    for i=1:nGM % Selects nSelect ground motions
        
        %compileReplacePercent
        %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
        scaleFac=zeros(nBig,1);
        for j=1:nBig
            if isScaled == 1
                if all(exp(sampleBig(j,posT1)) == 0)
                    scaleFac(j) = 1000000;
                else
                    scaleFac(j) = exp(lnSa1)/exp(sampleBig(j,PerTgt == T1));
                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;
            devSkew = skew(sampleSmall,1);
            devSig = std(sampleSmall) - sqrt(diag(covReq))';
            %devTotal = weights(1) * sum(devMean.^2) + weights(2) * sum(devSig.^2)+ weights(3) * sum(devSkew.^2);
            devTotal = weights(1) * sum(devMean.^2,2) + weights(2) * sum(devSig.^2,2);
            
            % Penalize bad spectra (set penalty to zero if this is not required)
            if penalty ~= 0
                for m=1:size(sampleSmall,1)
                    devTotal = devTotal + sum(abs(exp(sampleSmall(m,:))>emrsdc)) * penalty;
                end
            end
%             % 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+3*sqrt(diag(covReq))')),2) * 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 (all(devTotal < minDev) && ~any(recID == j))
                minID = j;
                minDev = devTotal(1);
            end
            sampleSmall = sampleSmall(1:end-1,:);
        end
        
        if ~isnan(minID)
            % 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
end

%display('100% done');
%compileReplace100Percent

if (size(recID,1)==0)
    errCode=3391;   %'Error: no match spectra found'
    finalRecords=0;
    finalScaleFactors=0;
    F_x=0;
    F1_y1=0;
    F1_y1a=0;
    F1_y2=0;
    F1_y3=0;
    F1_yn=0;
    F2_y1=0;
    F3_y1=0;
    F3_y2=0;
    F4_x1=perKnown;
    F4_y1=0;
    F5_y1=0;
    F6_y1=0;
    F7_y1=0;
    return;
end
% Output information
finalRecords_tmp = RSN(recID);
finalRecords_tmp = idRef(finalRecords_tmp);

finalDirections_tmp = ones(size(finalRecords_tmp,1),1);
for i=1:size(finalRecords_tmp,1)
    if recID(i) > size(Sa_1,1)
       finalDirections_tmp(i) = 2;
    end 
end

finalRecords = finalRecords_tmp;
finalDirectios = finalDirections_tmp;
finalScaleFactors = finalScaleFac;

%% Spectra Plots

F4_x1 = perKnown;
F4_y1 = SaKnown(recID,:).*repmat(finalScaleFactors,1,size(SaKnown,2));
F5_y1 = exp(sampleBig(recID,:)).*repmat(finalScaleFactors,1,length(PerTgt));
F6_y1 = exp(mean(sampleSmall));
F7_y1=std(sampleSmall);
% 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    : Target mean for the (log) response spectrum
%     % covReq     : Target covariance for the (log) response spectrum
%     %
%     
%     
%     % Plot at all periods
%     figure
%     loglog(PerTgt, exp(meanReq), 'b', 'linewidth', 3)
%     hold on
%     loglog(PerTgt, exp(meanReq + 1.96*sqrt(diag(covReq))'), '--b', 'linewidth', 3)
%     loglog(perKnown,SaKnown(finalRecords,:).*repmat(finalScaleFactors,1,size(SaKnown,2)),'k');
%     loglog(PerTgt, exp(meanReq - 1.96*sqrt(diag(covReq))'), '--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), 'b', 'linewidth', 3)
%     hold on
%     loglog(PerTgt, exp(meanReq + 1.96*sqrt(diag(covReq))'), '--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 - 1.96*sqrt(diag(covReq))'), '--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),'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))','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


%% Output data to file (best viewed with textpad)

% fin = fopen(outputFile,'w');
% fprintf(fin,'%s \t %s \t %s \t %s \n','Record Number','Scale Factor','File Name','URL'); % Title of relevant outputs
% for i = 1 : length(finalRecords)
%     rec = finalRecords(i);
%     url = ['http://peer.berkeley.edu/nga_files/ath/' Filename{rec}(1:end-3) 'AT2']; % URL of selected record
%     fprintf(fin,'%d \t %6.2f \t %s \t %s \n',i,finalScaleFactors(i),Filename{rec},url); % Print relevant outputs
% end
% 
% fclose(fin);
