Updated ResponseSpectrum script

etc/matlab/ResponseSpectrum.m

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-% Ground motion model (GMM) response spectra plots
+% Ground motion model (GMM) response spectra subplots
 %
 % Author: Peter Powers (pmpowers@usgs.gov)
+%         Demi Girot (dgirot@usgs.gov)
 %
-% Updated: 01/19/2023
+% Updated: 05/16/2023
 %
 % This script demonstrates how to call web services for the ground motion
 % models implemented in USGS National Sesmic Hazard Models and plot the
@@ -19,14 +20,15 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-%% Plot a figure of GMM reponse spectra and epistemic uncertainty
+%% Plot a figure of subplots for GMM response spectra, sigmas, and 
+% epistemic uncertainty
 
 clearvars
 
 % The root url for the response spectra web service
 urlbase = "https://earthquake.usgs.gov/ws/nshmp/data/gmm/spectra?";
 
-% If you are running nshmp-haz locally, use this URL:
+% If you are running nshmp-ws locally, use this URL:
 % urlbase = "http://localhost:8080/gmm/spectra?";
 
 % Struct of ground motion model parameters
@@ -48,9 +50,15 @@ plot_handles = [];
 epi = false;
 
 gmms_2023 = ["ASK_14" "BSSA_14" "CB_14" "CY_14"];
+
+subplot(2,1,1)
 means = readSpectra(urlbase, gmms_2023, in).response.means.data;
 plotGmms(means, in, epi);
 
+subplot(2,1,2)
+sigmas = readSpectra(urlbase, gmms_2023, in).response.sigmas.data;
+plotSigmas(sigmas, in, epi);
+
 %% Shared Functions
 
 function plotGmms(gmm_data, in, epi)
@@ -58,7 +66,7 @@ function plotGmms(gmm_data, in, epi)
     plot_handles = [];
     legend_labels = {};
 
-    % loop means_new response array
+    % loop means response array
     for j = 1 : length(gmm_data) 
     
         gmm_xs = gmm_data(j).data.xs;
@@ -111,7 +119,69 @@ function plotGmms(gmm_data, in, epi)
     
     legend( ...
         plot_handles,legend_labels, ...
-        'Location','southwest', ...
+        'Location','best', ...
+        'Interpreter', 'none');
+end
+
+function plotSigmas(gmm_data, in, epi)
+
+    plot_handles = [];
+    legend_labels = {};
+
+    % loop sigmas response array
+    for m = 1 : length(gmm_data) 
+    
+        sig_xs = gmm_data(m).data.xs;
+        sig_ys = gmm_data(m).data.ys;
+        sig_tree = gmm_data(m).tree;
+
+        % if 0.01 absent use PGA (0.001 is used for PGA)
+        if sig_xs(1) == 0.001 && sig_xs(2) ~= 0.01
+            sig_xs(1) = 0.01;
+        end
+
+        % Plot the total spectrum
+        ph = loglog(sig_xs, sig_ys, ...
+            'LineWidth',3, ...
+            'LineStyle','-');
+        
+        color = get(ph,'Color');
+        hold on;
+        
+        plot_handles = [plot_handles ph];
+        legend_labels{end+1} = gmm_data(m).label;
+        
+        % Plot epistemic spectra, if present
+        if epi && ~isempty(sig_tree)
+            for n = 1 : length(sig_tree)
+                id = sig_tree(n).id;
+                epi_ys = sig_tree(n).values;
+                
+                phe = loglog(sig_xs, epi_ys, ...
+                    'LineWidth',1, ...
+                    'LineStyle',':', ...
+                    'Color',color); 
+
+                plot_handles = [plot_handles phe];
+                legend_labels{end+1} = [gmm_data(m).label ' ' id];
+            end
+        end
+    end
+    
+    grid on
+
+    xlabel('Spectral Period (s)','FontSize',16)
+    ylabel('Standard Deviation','FontSize',16)
+    title(['M',num2str(in.Mw),', ',num2str(in.rRup),'km'],'FontSize',20)
+    xlim([0.01 10]);
+
+    set(gca,'FontSize',20);
+    set(gca,'XTick',[0.01 0.1 1 10]);
+    set(gca,'XTickLabel',{'0.01','0.1','1','10'});
+    
+    legend( ...
+        plot_handles,legend_labels, ...
+        'Location','best', ...
         'Interpreter', 'none');
 end
 
@@ -153,4 +223,4 @@ end
 
 response = webread(url);
 
-end
+end
\ No newline at end of file