From eb4f6b424f65b681b4ba316df686a3b8774855fe Mon Sep 17 00:00:00 2001
From: LorenzGit001 <lorenz.baer@stud.uni-hannover.de>
Date: Wed, 27 May 2026 21:05:56 +0200
Subject: [PATCH] Refactor SimulateEarthquake for scaling and output

Updated SimulateEarthquake function to include acceleration mismatch output and scaling methods. Enhanced earthquake signal processing with new parameters for length scale, Young's modulus ratio, and material density ratio.
---
 MATLAB/Methods/SimulateEarthquake.m | 105 ++++++++++++++++++++++------
 1 file changed, 84 insertions(+), 21 deletions(-)

diff --git a/MATLAB/Methods/SimulateEarthquake.m b/MATLAB/Methods/SimulateEarthquake.m
index d627047..3f9bc90 100644
--- a/MATLAB/Methods/SimulateEarthquake.m
+++ b/MATLAB/Methods/SimulateEarthquake.m
@@ -1,35 +1,98 @@
-function [pos,t,name] = SimulateEarthquake()
+function [pos,t,name,accMissmatch] = SimulateEarthquake(~)
 %SIMULATEEARTHQUAKE To apply a in amplitude scaled signal
-%   This function should apply an earthquake signal to the table motion,
-%   the position is roughly approximated from an rathquake's accelerogram
-%   However this function does only scale the amplitude and doesn't concern
-%   the energy content.
+% the scaling method has to be changed accordingly (Cauchy/Froude), this might change the outcome significantly 
 
 name = 'EQ_ElCentro';
-%threshold set for maximum displacement achieved with ElCentro
-threshold_displacements = 30;
+                                                                        
 %File seperator for the current system
 fsep = filesep;
 %Location of El Centro EQ signal
 eq_record_el_Centro = fullfile([pwd, fsep, 'Methods', fsep, 'earthquake_records'], 'elcentrons.txt');
+
 %Load the El Centro time and accelerogram into an array
-M_elcentro = readmatrix(eq_record_el_Centro);
+M_elcentro = readmatrix(eq_record_el_Centro); 
+
+%loaded acceleration in g
+accOrigin = M_elcentro(:,2); 
+
+%converting acceleration from g in mm/s^2 
+acc = accOrigin.*(9810);
 
 %loaded time in seconds
-t = M_elcentro(:,1);
+t = M_elcentro(:,1); 
 
-%loaded acceleration in g
-acc = M_elcentro(:,2);
-%converting acceleration from g in mm/s^2
-acc = acc*9810;
-%speed in mm/s
-speed = cumtrapz(t, acc);
-%position in mm
-pos = cumtrapz(t, speed);
-a0 = t\pos;
-pos = pos-a0*t;
-%scale maximum displacement to be 10mm
-pos = (pos./max(pos)).*threshold_displacements;
+%% Scaling
+scalingMethod = 'Cauchy'; % could also be defined in MainCreateSimulation
+
+%input length scale 
+lengthScale = 0.53; 
+
+%input Youngs modulus ratio of specimen to prototype 
+EScale = 1; 
+
+%input material density ratio of specimen to prototype 
+rhoScale = 1; 
+
+%place holder for Cauchy scaling
+accMissmatch = 9.81; 
 
+
+switch scalingMethod
+    case 'Froude'                    
+
+        % Time scaling only 
+        t= t .* sqrt(lengthScale);
+        % Acceleration is *not* scaled 
+        
+    case 'Cauchy'                  
+        % Time scaling
+        t = t .* ( lengthScale * (1/sqrt(rhoScale)) * sqrt(EScale) );
+        % Acceleration scaling 
+        acc = acc .* ( EScale * (1/rhoScale) * lengthScale ); 
+        % Acceleration miss match (g distortion) 
+        accMissmatch = 9.81/( EScale * (1./rhoScale) * lengthScale); 
+        
+         if t == sqrt(lengthScale) 
+           disp('Froude-Cauchy condition true')
+         end 
+    otherwise
+        error('Unsupported scaling method. Use ''Froude'' or ''Cauchy''.');
+end
+
+%% integration 
+% cumaltive trapizoids
+%speed in mm/s and detrend 
+speed_int = cumtrapz(t, acc);
+speed = detrend(speed_int,1);
+
+%position in mm and detrend 
+pos_int= cumtrapz(t, speed);
+pos = detrend(pos_int,1); 
+
+% filter starting displacement 0 Linear interpol
+pos(1:5) = 0; 
+for i = 6:15
+    % Linear interpolation between zero and the actual value
+    interp_factor = (i - 5) / (15 - 5);
+    pos(i) = pos(i) * interp_factor/5;  % abitrary devison for softer entry 
 end
 
+                                                                             
+% OLD SAFTY CAP WAS 30mm CHECK CURRENT SHAKING TABLE
+%% safety restricted displacement - search for max. displacement in current Namazu framework
+ if max(abs(pos)) > 50
+     pos = NaN; 
+     disp('Scaling to big, maximal displacement should be under 50 [mm]')
+ end 
+
+%% Plot for reference 
+% plot(t,pos); 
+% box off
+% xlabel('t in $$ [s] $$')
+% ylabel('displacement $$ [mm] $$')
+% fname = 'Displacement input namazu';
+
+%% Output acceleration missmatch
+disp(['Required gravitational acceleration: ',num2str(accMissmatch),' in [g]' ...
+    ' the acceleration distorted by factor ',num2str(accMissmatch/9.81),]); 
+end