% This file is the Master to plot everything!
% Curt Haselton
% 9-26-05

clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the data that is the same for all plots
    % Define the anchor period that all the analyses were run at (i.e. what T1
    % was used for defining Sa in the model; this is the same T1 that is used
    % for the "Sa" value used in the folder naming for all output)
        originalTOneUsedWhenAnalysesWereRun = 1.0; % seconds
        originalDampRatForSaTOneUsedWhenAnalysesWereRun = 0.05;
    
    % Input the damping ratio to use for the new Sa values
        currentDampRatToUseForSaCol = 0.05;
        
    % Define the period list used for the sigmaLn calculations (trying to
    % find the period with the lowest sigmaLn(Sa,col) value.
        periodLISTForCalcs = [0.10:0.10:5.00];
        
    % Define the periods to use for computing Sa,col statictics.  For each
    % period, we compute the collapse CDFs (lognormal fit) for each record 
    % set and save a file for the period. 
        periodsForCollapseCDFProcessingLIST = [0.5, 1.0, 1.5, 2.0, 2.4, 3.0, 4.0];
    
    % Define the periods to use for plotting the collapse CDFs (empirical
    % now) for each record set.
        periodsForCollapseCDFPlottingLIST = [0.5, 1.0, 1.5, 2.0, 2.4, 3.0, 4.0];
    
    % Input the information needed for the collapse risk calculations
    % (probability at the 2% in 50 year event and the MAF calcs.
        % 2% in 50 year GM level (this needs to have values corresponding
        % with the periods in the
        % "periodThatRecordSetsWereSelectedForLIST").  I put -1.0 for the
        % position of the non-epsilon set.  See hand notes on 9-23-05 for EERI paper with Baker)
            % Based on the Code spectrum, but this does not agree with the
            % non-epsilon set spectrum
                %twoPercentInFiftySaForEachPeriodLIST = [-1.0, 1.50, 0.90, 0.60, 0.45, 0.38, 0.36, 0.30, 0.23];    % With error from before - WRONG - DO NOT USE
                %twoPercentInFiftySaForEachPeriodLIST = [-1.0, 1.50, 0.90, 0.60, 0.45, 0.38, 0.30, 0.23];    % With error corrected
            % Based on assuming that the mean of the non-epsilon set
            % spectrum is the UHS
                %twoPercentInFiftySaForEachPeriodLIST = [-1.0, 1.496, 0.90, 0.629, 0.386, 0.291, 0.203, 0.139];  
                factorOnMAF = 0.75;
                twoPercentInFiftySaForEachPeriodLIST = factorOnMAF * [-1.0, 1.496, 0.90, 0.629, 0.386, 0.291, 0.203, 0.139];  % Make it 1.2 as large to test something! 

    
        % Input the exponent used for the hazard curve
            k = 3.0;  % This is general value from FEMA 450 (I think that is the pub. #).  I 
                                        %   may want to use 4-5 instead (to be consistent with the LA Bulk Mail Site)
        
    % Define the lists of EQ record sets and the T1 values that each set was
    % selected for.  Note that the first set is part of ATC-63 Set A and was
    % selected without considering epsilon (hence the -1.0 in the periodLIST).
        recordSetsLIST = {  [11011,	11012,	11021,	11022,	11031,	11032,	11041,	11042,	11051,	11052,	11061,	11062,	11071,	11072,	11081,	11082,	11091,	11092,	11101,	11102,	11111,	11112,	11121,	11122,	11131,	11132],
                        [1202	1203	1204	1212	1213	1215	1220	1222	1226	1230	1231	1232	1238	1239	1243	1244	1245	1251	1252	1253],
                        [1209	1212	1213	1214	1215	1218	1221	1222	1224	1225	1226	1227	1228	1230	1231	1232	1239	1240	1244	1251],
                        [1201	1206	1210	1211	1212	1215	1218	1219	1221	1222	1230	1231	1232	1235	1236	1238	1239	1249	1250	1251],
                        [1205	1208	1212	1215	1218	1219	1222	1230	1236	1237	1238	1239	1247	1249	1251	1252	1257	1258	1259	1261],
                        [1201	1207	1212	1215	1218	1219	1222	1223	1224	1233	1234	1236	1239	1242	1249	1250	1251	1252	1257	1258],
                        [1201	1212	1215	1217	1218	1223	1224	1225	1226	1230	1233	1234	1249	1258	1259	1260	1262	1263	1264	1265],
                        [1205	1212	1215	1216	1217	1218	1224	1226	1229	1233	1234	1241	1246	1248	1249	1251	1252	1254	1255	1256] };

        periodThatRecordSetsWereSelectedForLIST = [-1.0, 0.5, 1.0, 1.5, 2.0, 2.4, 3.0, 4.0];  % associated with recordSetsLIST
        markerTypeForPlotsLIST_indivPoints = {'k-', 'bo-', 'bv--', 'bx-', 'b+--', 'b*-', 'bs--', 'bd-'};
        markerTypeForPlotsLIST_lognormalFit = {'k-', 'k-', 'k-', 'k-', 'k-', 'k-', 'k-', 'k-'};
        legendEntriesLIST = {'Record selection without epsilon', 'Record selection at 0.5 sec.', 'Record selection at 1.0 sec.', 'Record selection at 1.5 sec.', 'Record selection at 2.0 sec.', 'Record selection at 2.4 sec.', 'Record selection at 3.0 sec.', 'Record selection at 4.0 sec.'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over SDOF results that we want to look at.  Send all of the above
% defined variables to each procedure (so we won't get lots or strange
% variables in this file).
%     outputFolderNameLIST = {'(SDOF_EQandCycDisp_v.4)_(TOne_1.00,Duct_4.00,PostK_-0.10)'};
%     currentFundPeriodOfSDOFLIST = [1.00];

    outputFolderNameLIST = {'(SDOF_EQandCycDisp_v.4)_(TOne_1.00,Duct_2.00,PostK_-0.10)',...
        '(SDOF_EQandCycDisp_v.4)_(TOne_1.00,Duct_2.00,PostK_-0.30)', '(SDOF_EQandCycDisp_v.4)_(TOne_1.00,Duct_4.00,PostK_-0.10)',...
        '(SDOF_EQandCycDisp_v.4)_(TOne_1.00,Duct_4.00,PostK_-0.30)', '(SDOF_EQandCycDisp_v.4)_(TOne_1.00,Duct_6.00,PostK_-0.10)',...
        '(SDOF_EQandCycDisp_v.4)_(TOne_1.00,Duct_6.00,PostK_-0.30)', '(SDOF_EQandCycDisp_v.4)_(TOne_2.00,Duct_2.00,PostK_-0.10)',...
        '(SDOF_EQandCycDisp_v.4)_(TOne_2.00,Duct_2.00,PostK_-0.30)', '(SDOF_EQandCycDisp_v.4)_(TOne_2.00,Duct_4.00,PostK_-0.10)',...
        '(SDOF_EQandCycDisp_v.4)_(TOne_2.00,Duct_4.00,PostK_-0.30)', '(SDOF_EQandCycDisp_v.4)_(TOne_2.00,Duct_6.00,PostK_-0.10)',...
        '(SDOF_EQandCycDisp_v.4)_(TOne_2.00,Duct_6.00,PostK_-0.30)'};
    currentFundPeriodOfSDOFLIST = [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00];

% Loop over all 
for outputFolderIndex = 1:length(outputFolderNameLIST)
    outputFolderName = outputFolderNameLIST{outputFolderIndex};
    currentFundPeriodOfSDOF = currentFundPeriodOfSDOFLIST(outputFolderIndex);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Process the collapse CDFs for all record sets using a single period for the
        % IM (definition of period is hard-coded in this file).
            isPlotThisOneNow = 0;
            if(isPlotThisOneNow == 1)
                for currentPeriodIndex = 1:length(periodsForCollapseCDFProcessingLIST)
                    currentPeriodToUseForSaCol = periodsForCollapseCDFProcessingLIST(currentPeriodIndex);
                    Driver_ProcessColCDFsForAllRecordSetsAtAPeriod(currentPeriodToUseForSaCol, outputFolderName, currentFundPeriodOfSDOF, originalTOneUsedWhenAnalysesWereRun, originalDampRatForSaTOneUsedWhenAnalysesWereRun, currentDampRatToUseForSaCol, periodLISTForCalcs, recordSetsLIST, periodThatRecordSetsWereSelectedForLIST, markerTypeForPlotsLIST_indivPoints, markerTypeForPlotsLIST_lognormalFit, legendEntriesLIST);
                end
            end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Plot the collapse CDFs for all record sets using a single period for the
        % IM.
        % NOTICE: Right now this reprocesses the data to make the plot.
        % This is probably okay just to leave it this way
            isPlotThisOneNow = 1;
            if(isPlotThisOneNow == 1)
                for currentPeriodIndex = 1:length(periodsForCollapseCDFPlottingLIST)
                    currentPeriodToUseForSaCol = periodsForCollapseCDFPlottingLIST(currentPeriodIndex);
                    figure(1);
                    hold on;
                    Driver_PlotColCDFsForAllRecordSetsAtAPeriod(currentPeriodToUseForSaCol, outputFolderName, currentFundPeriodOfSDOF, originalTOneUsedWhenAnalysesWereRun, originalDampRatForSaTOneUsedWhenAnalysesWereRun, currentDampRatToUseForSaCol, periodLISTForCalcs, recordSetsLIST, periodThatRecordSetsWereSelectedForLIST, markerTypeForPlotsLIST_indivPoints, markerTypeForPlotsLIST_lognormalFit, legendEntriesLIST);
                    close(1);
                end
                
            end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Compute the P[col | 2% in 50] and the MAF of collapse for each
        % ground motion selection option (i.e. with the period for the IM
        % and the period that the GM set was selected for being the same)
        % and for the non-epsilon set for all periods.
            isPlotThisOneNow = 0;
            if(isPlotThisOneNow == 1)
                % Note: Start at an idex of two two skip the -1 at the
                % start of the list.
                for currentPeriodIndex = 2:length(periodThatRecordSetsWereSelectedForLIST)
                    currentPeriodForColRiskCalcs = periodThatRecordSetsWereSelectedForLIST(currentPeriodIndex);
                    twoPercentInFiftySa = twoPercentInFiftySaForEachPeriodLIST(currentPeriodIndex);
                    
                    % Pass this to the function to compute the P[col | 2%
                    % in 50] and the MAF of collapse for this period; for
                    % all GM sets.
                    [probOfCollapseAtTwoPercentInFifty_AllGMSetsForCurrentPeriod, MAFofCollapse_AllGMSetsForCurrentPeriod] = ComputeAndSaveCollapseRiskForAPeriodAndAllGMSets(currentPeriodForColRiskCalcs, twoPercentInFiftySa, k, outputFolderName, currentFundPeriodOfSDOF, originalTOneUsedWhenAnalysesWereRun, originalDampRatForSaTOneUsedWhenAnalysesWereRun, currentDampRatToUseForSaCol, periodLISTForCalcs, recordSetsLIST, periodThatRecordSetsWereSelectedForLIST, markerTypeForPlotsLIST_indivPoints, markerTypeForPlotsLIST_lognormalFit, legendEntriesLIST);
                
                end
            end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Retreieve and output the P[col | 2% in 50] and the MAF of
        % collapse for each period (outputing the value for both the
        % non-epsilon and epsilon set at that period).
        % 
        % If the option is set, also make plots!
        %
        %
            isPlotThisOneNow = 0;
            isMakePlots = 1;
            if(isPlotThisOneNow == 1)
                % Initialize a few things for the plots
                probOfCollapse_nonEpsilonSet_LIST = zeros(1, (length(periodThatRecordSetsWereSelectedForLIST)));
                mafOfCollapse_nonEpsilonSet_LIST = zeros(1, (length(periodThatRecordSetsWereSelectedForLIST)));
                probOfCollapse_epsilonSet_LIST = zeros(1, (length(periodThatRecordSetsWereSelectedForLIST)));
                mafOfCollapse_epsilonSet_LIST = zeros(1, (length(periodThatRecordSetsWereSelectedForLIST)));
            
                % Note: Start at an index of two two skip the -1 at the
                % start of the list.
                for currentPeriodIndex = 2:length(periodThatRecordSetsWereSelectedForLIST)
                    currentPeriodForColRiskCalcs = periodThatRecordSetsWereSelectedForLIST(currentPeriodIndex);

                    % Call the function to print the results to the
                    % screen.  Note that the function also returns the
                    % results, and now I am using it to make plots.
                    [probOfCollapse_nonEpsilonSet, mafOfCollapse_nonEpsilonSet, probOfCollapse_epsilonSet, mafOfCollapse_epsilonSet] = PrintCollapseRiskForAPeriodForEpsilonAndNonEpsilonGMSets(currentPeriodIndex, currentPeriodForColRiskCalcs, outputFolderName);

                    % Save the results into a vector
                    probOfCollapse_nonEpsilonSet_LIST(currentPeriodIndex) = probOfCollapse_nonEpsilonSet;
                    mafOfCollapse_nonEpsilonSet_LIST(currentPeriodIndex) = mafOfCollapse_nonEpsilonSet;
                    probOfCollapse_epsilonSet_LIST(currentPeriodIndex) = probOfCollapse_epsilonSet;
                    mafOfCollapse_epsilonSet_LIST(currentPeriodIndex) = mafOfCollapse_epsilonSet;
                end
                
                % If the option is set to do so, make some plots
                if(isMakePlots == 1)
                    figure(1);
                    hold on;
                    PlotProbOfCollapseForMultGMSets(probOfCollapse_nonEpsilonSet_LIST, probOfCollapse_epsilonSet_LIST, periodThatRecordSetsWereSelectedForLIST, outputFolderName);
                    close(1);
                    figure(1);
                    hold on;
                    PlotMAFOfCollapseForMultGMSets(mafOfCollapse_nonEpsilonSet_LIST, mafOfCollapse_epsilonSet_LIST, periodThatRecordSetsWereSelectedForLIST, k, outputFolderName);
                    close(1);
                end
            end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Plot the sigmaLn(Sa,col) for each record set over many periods.
        % This takes a while to do.  ****** It now both makes the plots and saves
        % the results, so I may want to seperate these two pieces later.
            isPlotThisOneNow = 0;
            if(isPlotThisOneNow == 1)
                figure(1);
                hold on;
                Driver_PlotSigmaLnSaColVSPeriodForMultRecordSets(outputFolderName, currentFundPeriodOfSDOF, originalTOneUsedWhenAnalysesWereRun, originalDampRatForSaTOneUsedWhenAnalysesWereRun, currentDampRatToUseForSaCol, periodLISTForCalcs, recordSetsLIST, periodThatRecordSetsWereSelectedForLIST, markerTypeForPlotsLIST_indivPoints, markerTypeForPlotsLIST_lognormalFit, legendEntriesLIST);
                close(1);
            end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Compute and plot the ratio between the median collapse capacity
        % and the 2% in 50 Sa level for the epsilon and non-epsilon sets
        % for each period.  Make a plot that shows this collapse median
        % collapse capacity for the non-epsilon set for all periods and for
        % the epsilon sets for each period (for the epsilon sets, the
        % period for selection and the period for the IM are consistent)
            isPlotThisOneNow = 0;
            if(isPlotThisOneNow == 1)
                % Call the procedure
                figure(1);
                hold on;
                PlotCollapseMarginForManyRecordSets(outputFolderName, periodThatRecordSetsWereSelectedForLIST, twoPercentInFiftySaForEachPeriodLIST);
                close(1);
            end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Compute and plot the sigmaLn (Sa,col) for the epsilon and non-epsilon sets
        % for each period.  Note that this was done for 
            isPlotThisOneNow = 0;
            if(isPlotThisOneNow == 1)
                % Call the procedure
                figure(1);
                hold on;
                PlotSigmaLnForManyRecordSets(outputFolderName, periodThatRecordSetsWereSelectedForLIST, twoPercentInFiftySaForEachPeriodLIST);
                close(1);
            end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        
        
        
        
        
        
        
        
        
        
end; % end of loop for all output folders
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the exp(mean(Ln(Sa))) spectra for all of the GM sets against the Code 1/T spectrum
        isPlotThisOneNow = 0;
        if(isPlotThisOneNow == 1)
            % Define plot options
            dampingRatioToPlot = 0.05;
            
            % Plot Code 1/T spectrum - see hnad notes on 9-23-05 for
            % numbers
            codePeriods = [0.00, 0.25, 0.50, 1.00, 1.50, 2.00, 2.40, 3.00, 4.00];
            codeTwoPercentInFiftyValues = factorOnMAF * [1.50, 1.50, 1.50, 0.90, 0.60, 0.45, 0.38, 0.30, 0.23];
            lineTypeForCodeLine = 'r-';
            lineSizeForCodeLine = 3;
            plot(codePeriods, codeTwoPercentInFiftyValues, lineTypeForCodeLine, 'LineWidth', lineSizeForCodeLine);
            
            figure(1);
            hold on;
            
            % Loop and plot the spectra for all the GM sets
            for recordSetIndex = 1:length(recordSetsLIST)
                recordSet_current = recordSetsLIST{recordSetIndex};

                % Check if we are plotting the epsilon-independent set and
                % make adjustments if this is the case.
                if(recordSetIndex == 1)
                    % We are plotting the record set selected without
                    % considering epsilon
                    lineTypeForMeanLnLine = 'b--';
                    lineSizeForMeanLnLine = 3;
                    
                    periodThatRecordSetWasSelectedFor_current = periodThatRecordSetsWereSelectedForLIST(3);     % Anchor it at 1.0 seconds for simplicity
                    twoPercentInFiftySaForPeriod_current = twoPercentInFiftySaForEachPeriodLIST(3);             % Hazard at 1.0 seconds
                
                    PlotMultSaSpectrumFunc_ScaledToSaLevel_Mean(recordSet_current, twoPercentInFiftySaForPeriod_current, periodThatRecordSetWasSelectedFor_current, dampingRatioToPlot, lineTypeForMeanLnLine, lineSizeForMeanLnLine)
                    
                else
                    % We are plotting one of the sets selected based on
                    % epsilon at a period
                    lineTypeForMeanLnLine = 'k-';
                    lineSizeForMeanLnLine = 1;
                    
                    periodThatRecordSetWasSelectedFor_current = periodThatRecordSetsWereSelectedForLIST(recordSetIndex);
                    twoPercentInFiftySaForPeriod_current = twoPercentInFiftySaForEachPeriodLIST(recordSetIndex);
                    
                    PlotMultSaSpectrumFunc_ScaledToSaLevel_Mean(recordSet_current, twoPercentInFiftySaForPeriod_current, periodThatRecordSetWasSelectedFor_current, dampingRatioToPlot, lineTypeForMeanLnLine, lineSizeForMeanLnLine)

                end
            end

            % Add legend
            legh = legend('2% in 50 year Code spectrum', 'Mean spectrum for non-epsilon set (anchored at 1.0 sec.)', 'Mean spectra for seven sets selected based on epsilon');
                
            % Save the figures in the output folder
            % Go to the output folder
            currentFolder_temp = [pwd];
            cd ..;
            cd ..;
            cd Output;
            cd(outputFolderName);
   
            % Save the plots in this folder
            % Save the plot as a .fig file
                plotName = sprintf('Spectra_AllRecSetsAndCode.fig');
                hgsave(plotName);
            % Export the plot as a .emf file (Matlab book page 455)
                exportName = sprintf('Spectra_AllRecSetsAndCode.emf');
                print('-dmeta', exportName);
    
            % Go back to original folder
            cd(currentFolder_temp);  
            
            hold off
            close(1);
            
        end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%












    
    




