% 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 = [1.0];
    
    % Define the periods to use for plotting the collapse CDFs (empirical
    % now) for each record set.
        periodsForCollapseCDFPlottingLIST = [1.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 = {  [15011	15021	15031	15041	15051	15061	15071	15081	15091	15101	15012	15022	15032	15042	15052	15062	15072	15082	15092	15102],  
                            [15111	15121	15131	15141	15151	15161	15171	15181	15191	15201	15112	15122	15132	15142	15152	15162	15172	15182	15192	15202]};
                    
        periodThatRecordSetsWereSelectedForLIST = [-1.0, 0.5, 1.0, 1.5, 2.0, 2.4, 3.0, 4.0];  % associated with recordSetsLIST
        markerTypeForPlotsLIST_indivPoints = {'bo-', 'bs--'};
        markerTypeForPlotsLIST_lognormalFit = {'k-', 'k-'};
        legendEntriesLIST = {'Chi-Chi with 15cm/s < PGV < 16cm/s', 'Chi-Chi with 51cm/s < PGV < 66cm/s'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.5)_(TOne_1.00,Duct_6.00,PostK_-0.10)'};
    currentFundPeriodOfSDOFLIST = [1.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%












    
    




