diff --git a/gtm_design/FlightCockpit_GUI.mlapp b/gtm_design/FlightCockpit_GUI.mlapp
new file mode 100644
index 0000000..8bb80c3
Binary files /dev/null and b/gtm_design/FlightCockpit_GUI.mlapp differ
diff --git a/gtm_design/Oscillatory Spin Trajectory.JPG b/gtm_design/Oscillatory Spin Trajectory.JPG
new file mode 100644
index 0000000..70ace8f
Binary files /dev/null and b/gtm_design/Oscillatory Spin Trajectory.JPG differ
diff --git a/gtm_design/OscillatorySpin.asv b/gtm_design/OscillatorySpin.asv
new file mode 100644
index 0000000..799d9ad
--- /dev/null
+++ b/gtm_design/OscillatorySpin.asv
@@ -0,0 +1,237 @@
+%% --------------------- Oscillatory Trajectory Simulation --------------------------
+%
+% This MATLAB script models a targeted oscillatory spin condition for a GTM_T2 aircraft, 
+% simulates the trajectory from that condition, and plots various flight 
+% parameters over time.
+%
+% Simulation Details:
+% - The script begins by initializing the simulation with nominal design conditions for 
+%   a GTM_T2 aircraft.
+% - It then proceeds to set up conditions for an oscillatory spin as specified in
+%   init_spin_design function.
+% - The trajectory is simulated for a duration of 100 seconds.
+% - The resulting trajectory is plotted, showing the behavior of key flight parameters 
+%   such as alpha, beta, flight path angle, airspeed, angular rates, and Euler angles 
+%   over time.
+% - The trajectory is also visualized in a 3D plot with an orientated aircraft shape.
+%
+%
+%% About
+%
+% Author:     Omar Mourad
+% Email:      <mailto: st181901@stud.uni-stuttgart.de>
+% Created:    20.03.2024 
+
+
+%% Conditions
+
+% Start with some altitude,otherwise nominal init
+MWS = init_design('GTM_T2');
+loadmws(MWS,'gtm_design');
+printStates(MWS);
+
+% Proceed with the oscillatory spin conditions
+MWS = init_spin_design('GTM_T2');
+printStates(MWS);
+loadmws(MWS,'gtm_design');
+
+%% Simulate
+%Determine the Simulation time (in Secs)
+time_of_simulation = 10; 
+[t,x,y]=sim('gtm_design',[0 time_of_simulation]); 
+
+%% Convert from lat/lon to ft 
+Xeom=y(:,7:18);
+dist_lat=(Xeom(:,7)-Xeom(1,7)) * 180/pi*364100.79;
+dist_lon=(Xeom(:,8)-Xeom(1,8)) * 180/pi*291925.24;
+alt=Xeom(:,9);
+tplot=[0:.1:max(t)];
+
+%% Save the data to file
+save('simulationData.mat', 't', 'y', 'Xeom', 'dist_lat', 'dist_lon', 'alt');
+
+%% Define simple airplane shape
+scale=1;
+x1=scale*([0.0,-.5, -2.0, -3.0,-4.0, -3.25, -5.5, -6.0, -6.0]+3.0);
+y1=scale*[0.0, 0.5,  0.5, 4.25, 4.5,  0.5,   0.5,  1.5,  0.0];
+Vehicletop=[ [x1,fliplr(x1)]; [y1,fliplr(-y1)]; -.01*ones(1,2*length(x1))];
+Vehiclebot=[ [x1,fliplr(x1)]; [y1,fliplr(-y1)];  .01*ones(1,2*length(x1))];
+
+
+%%NEW PLOTS
+% Interpolate data to match the tplot timeline
+X_ani = interp1(t, Xeom, tplot);
+Y_ani = interp1(t, y(:,1:5), tplot);
+lat_ani = interp1(t, dist_lat, tplot);
+lon_ani = interp1(t, dist_lon, tplot);
+alt_ani = interp1(t, alt, tplot);
+
+%% Setup Dynamic Plots
+h = figure(1);
+set(h, 'Position', [20, 20, 1200, 800]);
+clf;
+
+% Alpha/Beta Plot Initialization
+%hAlphaBeta = subplot(4, 2, 1);
+axes('position',[.1 .75 .23 .13])
+pAlpha = plot(tplot(1), X_ani(1, 3), 'r', 'LineWidth', 1.5); % Alpha
+hold on;
+pBeta = plot(hAlphaBeta, tplot(1), X_ani(1, 4), 'b', 'LineWidth', 1.5); % Beta
+grid on;
+legend({'\alpha', '\beta'}, 'Location', 'SouthEast');
+xlabel('Time (sec)');
+ylabel('\alpha (deg), \beta (deg)');
+title('Alpha/Beta');
+ylim([-50 50]);
+
+% Flight Path and Airspeed Initialization
+%hFlightPath = subplot(4, 2, 2);
+axes('position',[.1 .55 .23 .13])
+hGamma = plot(tplot(1), Y_ani(1, 5), 'r', 'LineWidth', 1.5); % Gamma
+hold on;
+hEAS = plot(hFlightPath, tplot(1), Y_ani(1, 1), 'b', 'LineWidth', 1.5); % Air-Speed
+grid on;
+legend({'\gamma', 'EAS'}, 'Location', 'SouthEast');
+xlabel('Time (sec)');
+ylabel('\gamma (deg), Equivalent Airspeed (knots)');
+title('Flight Path Angle and Airspeed');
+
+% Angular Rates Plot Initialization
+hAngularRates = subplot(4, 2, 3);
+axes('position',[.1 .35 .23 .13])
+pP = plot(hAngularRates, tplot(1), X_ani(1, 4) * 180 / pi, 'r', 'LineWidth', 1.5); % p
+hold on;
+pQ = plot(hAngularRates, tplot(1), X_ani(1, 5) * 180 / pi, 'g', 'LineWidth', 1.5); % q
+pR = plot(hAngularRates, tplot(1), X_ani(1, 6) * 180 / pi, 'b', 'LineWidth', 1.5); % r
+legend({'p', 'q', 'r'}, 'Location', 'SouthEast');
+xlabel('Time (sec)');
+ylabel('Angular Rates (deg/sec)');
+title('Angular Rates');
+grid on;
+
+% Euler Angles Plot Initialization
+hEulerAngles = subplot(4, 2, 4);
+pRoll = plot(hEulerAngles, tplot(1), X_ani(1, 10) * 180 / pi, 'r', 'LineWidth', 1.5); % Roll (phi)
+hold on;
+pPitch = plot(hEulerAngles, tplot(1), X_ani(1, 11) * 180 / pi, 'g', 'LineWidth', 1.5); % Pitch (theta)
+pYaw = plot(hEulerAngles, tplot(1), X_ani(1, 12) * 180 / pi, 'b', 'LineWidth', 1.5); % Yaw (psi)
+legend({'Roll', 'Pitch', 'Yaw'}, 'Location', 'NorthEast');
+xlabel('Time (sec)');
+ylabel('\phi (deg), \theta (deg), \psi (deg)');
+title('Euler Angles');
+grid on;
+
+% 3D Trajectory Plot Initialization
+hTrajectory = subplot(4, 2, [6 7]);
+plot3(hTrajectory, dist_lat, dist_lon, alt, 'k', 'LineWidth', 1.5); % Entire Trajectory
+hold on;
+h3DAircraft = plot3(hTrajectory, lat_ani(1), lon_ani(1), alt_ani(1), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
+grid on;
+xlabel('Lat. Crossrange (ft)');
+ylabel('Long. Crossrange (ft)');
+zlabel('Altitude (ft)');
+title('3D Oscillatory Spin Trajectory');
+view([25, 10]);
+
+%% Start Dynamic Plot Updates
+for i = 1:length(tplot)
+    % Update Alpha/Beta plot
+    set(pAlpha, 'XData', tplot(1:i), 'YData', X_ani(1:i, 3));
+    set(pBeta, 'XData', tplot(1:i), 'YData', X_ani(1:i, 4));
+    
+    % Update Flight Path and Airspeed plot
+    set(hGamma, 'XData', tplot(1:i), 'YData', Y_ani(1:i, 5));
+    set(hEAS, 'XData', tplot(1:i), 'YData', Y_ani(1:i, 1));
+    
+    % Update Angular Rates plot
+    set(pP, 'XData', tplot(1:i), 'YData', X_ani(1:i, 4) * 180 / pi);
+    set(pQ, 'XData', tplot(1:i), 'YData', X_ani(1:i, 5) * 180 / pi);
+    set(pR, 'XData', tplot(1:i), 'YData', X_ani(1:i, 6) * 180 / pi);
+    
+    % Update Euler Angles plot
+    set(pRoll, 'XData', tplot(1:i), 'YData', X_ani(1:i, 10) * 180 / pi);
+    set(pPitch, 'XData', tplot(1:i), 'YData', X_ani(1:i, 11) * 180 / pi);
+    set(pYaw, 'XData', tplot(1:i), 'YData', X_ani(1:i, 12) * 180 / pi);
+    
+    % Update 3D Trajectory plot
+    set(h3DAircraft, 'XData', lat_ani(i), 'YData', lon_ani(i), 'ZData', alt_ani(i));
+    
+    % Plot Aircraft Shape at Current Position
+    % Top View
+    fill3(Vehicletop(1,:) + lat_ani(i), Vehicletop(2,:) + lon_ani(i), Vehicletop(3,:) + alt_ani(i), 'r', 'FaceAlpha', 0.5);
+    % Bottom View
+    fill3(Vehiclebot(1,:) + lat_ani(i), Vehiclebot(2,:) + lon_ani(i), Vehiclebot(3,:) + alt_ani(i), 'b', 'FaceAlpha', 0.5);
+
+    % Pause for animation effect
+    pause(0.1);
+end
+
+
+
+% %% ------------------------Plots---------------------------
+% h=figure(1);,set(h,'Position',[20,20,1200,800]);clf
+% 
+% % Alpha/Beta
+% axes('position',[.1 .75 .23 .13])
+% plot(t,[y(:,3),y(:,4)]); % y(:,3) ----> Alphas // y(:,4) ----> Betas
+% legend({'\alpha','\beta'},'Location','SouthEast');,grid on
+% xlabel('time (sec)'),
+% ylabel('\alpha (deg), \beta (deg)');
+% title('Alpha/Beta');
+% 
+% % Flight Path Angle and Airspeed
+% axes('position',[.1 .55 .23 .13])
+% [ax,h1,h2]=plotyy(t,y(:,5),t,y(:,1));grid on % y(:,5) ----> Gammas //
+% xlabel('time (sec)'),                        % %y(:,1) ----> Airspeeds
+% ylabel(ax(1),'Flight Path,  \gamma (deg)');
+% ylabel(ax(2),'Equivalent Airspeed (konts)')
+% legend([h1;h2],{'\gamma','eas'},'Location','SouthEast');
+% title('Flight Path Angle and Airspeed');
+% 
+% % Angular Rates
+% axes('position',[.1 .35 .23 .13])
+% plot(t,Xeom(:,4:6)*180/pi);grid on
+% legend({ 'p','q','r'},'Location','SouthEast');
+% xlabel('time (sec)'),ylabel('angular rates (deg/sec)')
+% title('Angular Rates');
+% 
+% % Euler Angles
+% axes('position',[.1 .15 .23 .13]);
+% plot(t,[y(:,16)*180/pi, y(:,17)*180/pi, y(:,18)*180/pi]);
+% % y(:,16) ---> Phi  // y(:,17) ---> Theta // y(:,18) ---> Psi
+% legend({'roll','pitch','yaw'},'Location','NorthEast');,grid on
+% xlabel('time (sec)');
+% ylabel('\phi (deg), \theta (deg), \psi (deg)');
+% title('Euler Angles');
+% 
+% % Trajectory: 3D plot with orientated vehicle
+% axes('position',[.45,.15,.5,.7])
+% plot3(dist_lat,dist_lon,alt);grid on, axset=axis; % Just get axis limits
+% %cnt=0;
+% % resample at equally spaced points for animation plot
+% tplot=[0:.1:max(t)];
+% X_ani=interp1(t,Xeom,tplot);
+% lat_ani=interp1(t,dist_lat,tplot);
+% lon_ani=interp1(t,dist_lon,tplot);
+% alt_ani=interp1(t,alt,tplot);
+% tic
+% for i=[1:length(tplot)]
+%   plot3(dist_lat,dist_lon,alt);grid on
+%   Offset=repmat([lat_ani(i);lon_ani(i);alt_ani(i)],1,size(Vehicletop,2));
+%   Ptmp=diag([1,1,-1])*transpose(euler321(X_ani(i,10:12)))*Vehicletop + Offset;
+%   patch(Ptmp(1,:),Ptmp(2,:),Ptmp(3,:),'g');
+%   Ptmp=diag([1,1,-1])*transpose(euler321(X_ani(i,10:12)))*Vehiclebot + Offset;
+%   patch(Ptmp(1,:),Ptmp(2,:),Ptmp(3,:),'c');
+%   view(25,10),axis(axset),hold off
+%   xlabel('Lat. Crossrange(ft)');
+%   ylabel('Long. Crossrange(ft)');
+%   zlabel('Altitude(ft)');
+%   title('Simulation of Oscillatory Spin Trajectory');
+%   pause(.1);
+% end
+% toc
+% if(exist('AutoRun','var'))
+%     pause(.2);
+%     orient portrait; print -dpng OscillatorySpin;
+% end
+% 
diff --git a/gtm_design/OscillatorySpin.m b/gtm_design/OscillatorySpin.m
new file mode 100644
index 0000000..6764002
--- /dev/null
+++ b/gtm_design/OscillatorySpin.m
@@ -0,0 +1,51 @@
+%% --------------------- Oscillatory Trajectory Simulation --------------------------
+%
+% This MATLAB script models a targeted oscillatory spin condition for a GTM_T2 aircraft, 
+% simulates the trajectory from that condition, and plots various flight 
+% parameters over time.
+%
+% Simulation Details:
+% - The script begins by initializing the simulation with nominal design conditions for 
+%   a GTM_T2 aircraft.
+% - It then proceeds to set up conditions for an oscillatory spin as specified in
+%   init_spin_design function.
+% - The trajectory is simulated for a duration of 100 seconds.
+% - The resulting trajectory is plotted, showing the behavior of key flight parameters 
+%   such as alpha, beta, flight path angle, airspeed, angular rates, and Euler angles 
+%   over time.
+% - The trajectory is also visualized in a 3D plot with an orientated aircraft shape.
+%
+%
+%% About
+%
+% Author:     Omar Mourad
+% Email:      <mailto: st181901@stud.uni-stuttgart.de>
+% Created:    20.03.2024 
+
+
+%% Conditions
+
+% Start with some altitude,otherwise nominal init
+MWS = init_design('GTM_T2');
+loadmws(MWS,'gtm_design');
+printStates(MWS);
+
+% Proceed with the oscillatory spin conditions
+MWS = init_spin_design('GTM_T2');
+printStates(MWS);
+loadmws(MWS,'gtm_design');
+
+%% Simulate
+%Determine the Simulation time (in Secs)
+time_of_simulation = 100; 
+[t,x,y]=sim('gtm_design',[0 time_of_simulation]); 
+
+%% Convert from lat/lon to ft 
+Xeom=y(:,7:18);
+dist_lat=(Xeom(:,7)-Xeom(1,7)) * 180/pi*364100.79;
+dist_lon=(Xeom(:,8)-Xeom(1,8)) * 180/pi*291925.24;
+alt=Xeom(:,9);
+tplot=0:.1:max(t);
+
+%% Save the data to file
+save('OscillatorySpinData.mat','t', 'tplot' , 'y', 'Xeom', 'dist_lat', 'dist_lon', 'alt');
\ No newline at end of file
diff --git a/gtm_design/OscillatorySpinData.mat b/gtm_design/OscillatorySpinData.mat
new file mode 100644
index 0000000..0303da7
Binary files /dev/null and b/gtm_design/OscillatorySpinData.mat differ
diff --git a/gtm_design/Steep Spiral Trajectory.JPG b/gtm_design/Steep Spiral Trajectory.JPG
new file mode 100644
index 0000000..942fd7a
Binary files /dev/null and b/gtm_design/Steep Spiral Trajectory.JPG differ
diff --git a/gtm_design/SteepSpiral.m b/gtm_design/SteepSpiral.m
new file mode 100644
index 0000000..a1b6b4d
--- /dev/null
+++ b/gtm_design/SteepSpiral.m
@@ -0,0 +1,51 @@
+%% --------------------- Steep Spiral Trajectory Simulation --------------------------
+%
+% This MATLAB script models a targeted steep spiral condition for a GTM_T2 aircraft, 
+% simulates the trajectory from that condition, and plots various flight 
+% parameters over time.
+%
+% Simulation Details:
+% - The script begins by initializing the simulation with nominal design conditions for 
+%   a GTM_T2 aircraft.
+% - It then proceeds to set up conditions for a steep spiral as specified in
+%   init_spiral_design function.
+% - The trajectory is simulated for a duration of 100 seconds.
+% - The resulting trajectory is plotted, showing the behavior of key flight parameters 
+%   such as alpha, beta, flight path angle, airspeed, angular rates, and Euler angles 
+%   over time.
+% - The trajectory is also visualized in a 3D plot with an orientated aircraft shape.
+%
+%
+%% About
+%
+% Author:     Omar Mourad
+% Email:      <mailto: st181901@stud.uni-stuttgart.de>
+% Created:    15.03.2024 
+
+
+%% Conditions
+
+% Start with some altitude,otherwise nominal init
+MWS = init_design('GTM_T2');
+loadmws(MWS,'gtm_design');
+printStates(MWS);
+
+% Proceed with the steep spiral conditions
+MWS = init_spiral_design('GTM_T2');
+printStates(MWS);
+loadmws(MWS,'gtm_design');
+
+%% Simulate
+%Determine the Simulation time (in Secs)
+time_of_simulation = 100;    
+[t,x,y]=sim('gtm_design',[0 time_of_simulation]); 
+intermediate = y;
+%% Convert from lat/lon to ft 
+Xeom=y(:,7:18);
+dist_lat=(Xeom(:,7)-Xeom(1,7)) * 180/pi*364100.79;
+dist_lon=(Xeom(:,8)-Xeom(1,8)) * 180/pi*291925.24;
+alt=Xeom(:,9);
+tplot=0:.1:max(t);
+
+%% Save the data to file
+save('SteepSpiralData.mat','t', 'tplot' , 'y', 'Xeom', 'dist_lat', 'dist_lon', 'alt');
diff --git a/gtm_design/SteepSpiralData.mat b/gtm_design/SteepSpiralData.mat
new file mode 100644
index 0000000..4a64d89
Binary files /dev/null and b/gtm_design/SteepSpiralData.mat differ
diff --git a/gtm_design/mfiles/CreateSimOut.m b/gtm_design/mfiles/CreateSimOut.m
index 84331bc..f77346a 100644
--- a/gtm_design/mfiles/CreateSimOut.m
+++ b/gtm_design/mfiles/CreateSimOut.m
@@ -1,62 +1,62 @@
-function [SimOut] = CreateSimOut(logsout)
-% Create structure from logged bus signal that is named SimOut.  All
-% signals in SimOut must be the same sample time and rate.  SimOut signal
-% must be inside a block on the top level of the diagram named
-% "Format_Outputs".  
-
-% $Id: CreateSimOut.m 2 2015-05-07 21:03:59Z cox $
-
-% Paste the commented code below into the Model StopFcn callback:
-
-% if SimIn.Switches.LogData
-%     if exist('logsout','var')
-%         clear SimOut
-%         SimOut = CreatSimOut(logsout);
-%     else
-%         Disp('No Log Data Variable')
-%     end
-% end
-% clear logsout
-
-% Eugene Heim, NASA Langley Research Center
-% Modified, david.e.cox NASA Langley Research Center
-
-
-if (isa(logsout,'Simulink.ModelDataLogs')==1),   % Old ModelDataLogs format for logsout
-  temp = logsout.whos('all');% Get all field names
-  index = find(strcmp('Timeseries',{temp(:).simulinkClass}));% Find signals
-  for ii = 1:length(index)
-       % Remove blockname hierarchy. Cut string before first occurance of SimOut 
-      tmpstr=temp(index(ii)).name;
-      VariableName = tmpstr([min(strfind(tmpstr,'SimOut')):end]);
-      % Grab TimeSeries data, remove singletons and make time vector first dimension, if necessary
-      eval(sprintf('data = squeeze(logsout.%s.Data);', temp(index(ii)).name));
-      eval(sprintf('timelen = length(logsout.%s.Time);',temp(index(ii)).name))
-      if ndims(data) > 2 && timelen>1 % for N-D matrices time is last dim, unless time is singleton.
-          eval([VariableName, ' = permute(data,[ndims(data),[1:ndims(data)-1]]);'])
-      else
-          eval([VariableName,' = data;']);
-      end
-   end
-elseif (isa(logsout,'Simulink.SimulationData.Dataset')==1),  % New Dataset format for logsout
-  tmp=logsout.getElement(1);  %Should only be one.
-  SimOut_TS=tmp.Values;       %Already a data struct
-  SimOut=ts2vec(SimOut_TS);   %Convert timeseries to simple vector/matrix, recursive function
-end
-
-
-function sout = ts2vec(sin);
-%function sout = ts2vec(sin);
-
-sout=sin;
-fn=fieldnames(sin);
-for i=1:length(fn),
-    if (isa(sin.(fn{i}),'timeseries')==1),
-        tstmp=sin.(fn{i});
-        sout.(fn{i})=tstmp.Data;
-    elseif (isa(sin.(fn{i}),'struct')==1),
-        sout.(fn{i})=ts2vec(sin.(fn{i})); % Recursive voodoo...
-    else
-        fprintf(1,'WARNING: Action for class %s at field %s  is undefined, ignoring\n',class(sin.(fn{i})),fn{i});
-    end
-end
+function [SimOut] = CreateSimOut(logsout)
+% Create structure from logged bus signal that is named SimOut.  All
+% signals in SimOut must be the same sample time and rate.  SimOut signal
+% must be inside a block on the top level of the diagram named
+% "Format_Outputs".  
+
+% $Id: CreateSimOut.m 2 2015-05-07 21:03:59Z cox $
+
+% Paste the commented code below into the Model StopFcn callback:
+
+% if SimIn.Switches.LogData
+%     if exist('logsout','var')
+%         clear SimOut
+%         SimOut = CreatSimOut(logsout);
+%     else
+%         Disp('No Log Data Variable')
+%     end
+% end
+% clear logsout
+
+% Eugene Heim, NASA Langley Research Center
+% Modified, david.e.cox NASA Langley Research Center
+
+
+if (isa(logsout,'Simulink.ModelDataLogs')==1),   % Old ModelDataLogs format for logsout
+  temp = logsout.whos('all');% Get all field names
+  index = find(strcmp('Timeseries',{temp(:).simulinkClass}));% Find signals
+  for ii = 1:length(index)
+       % Remove blockname hierarchy. Cut string before first occurance of SimOut 
+      tmpstr=temp(index(ii)).name;
+      VariableName = tmpstr([min(strfind(tmpstr,'SimOut')):end]);
+      % Grab TimeSeries data, remove singletons and make time vector first dimension, if necessary
+      eval(sprintf('data = squeeze(logsout.%s.Data);', temp(index(ii)).name));
+      eval(sprintf('timelen = length(logsout.%s.Time);',temp(index(ii)).name))
+      if ndims(data) > 2 && timelen>1 % for N-D matrices time is last dim, unless time is singleton.
+          eval([VariableName, ' = permute(data,[ndims(data),[1:ndims(data)-1]]);'])
+      else
+          eval([VariableName,' = data;']);
+      end
+   end
+elseif (isa(logsout,'Simulink.SimulationData.Dataset')==1),  % New Dataset format for logsout
+  tmp=logsout.getElement(1);  %Should only be one.
+  SimOut_TS=tmp.Values;       %Already a data struct
+  SimOut=ts2vec(SimOut_TS);   %Convert timeseries to simple vector/matrix, recursive function
+end
+
+
+function sout = ts2vec(sin);
+%function sout = ts2vec(sin);
+
+sout=sin;
+fn=fieldnames(sin);
+for i=1:length(fn),
+    if (isa(sin.(fn{i}),'timeseries')==1),
+        tstmp=sin.(fn{i});
+        sout.(fn{i})=tstmp.Data;
+    elseif (isa(sin.(fn{i}),'struct')==1),
+        sout.(fn{i})=ts2vec(sin.(fn{i})); % Recursive voodoo...
+    else
+        fprintf(1,'WARNING: Action for class %s at field %s  is undefined, ignoring\n',class(sin.(fn{i})),fn{i});
+    end
+end
diff --git a/gtm_design/mfiles/appendmws.m b/gtm_design/mfiles/appendmws.m
index b5e390c..6767079 100644
--- a/gtm_design/mfiles/appendmws.m
+++ b/gtm_design/mfiles/appendmws.m
@@ -1,35 +1,35 @@
-function appendmws(MWS,model)
-%function appendmws(MWS,model)
-%
-% Appends model workspace with variables 
-% given by the fields of the structure MWS.  
-%
-% ! Warning !  Will overwrite indentically named variables without warning!
-%
-% Inputs:
-%  MWS   - Model Workspace Structure, contains simulation parameters
-%  model - Name of model to load into, default is 'gtm_design'
-%
-
-% $Id: appendmws.m 4852 2013-08-06 22:12:54Z cox $
-% Author: Austin Murch (Austin.M.Murch@nasa.gov)
-% Adpated from code by Dave Cox (d.e.cox@larc.nasa.gov)
-
-
-% By default use the bdroot model
-if ( ~exist('model','var') || isempty(model) ),
-  model=bdroot;
-end
-
-mws=get_param(model,'modelworkspace');
-
-fn=fieldnames(MWS);
-for i=1:length(fn)
-    mws.assignin(fn{i},MWS.(fn{i}));
-end
-
-
-
-
-
-  
+function appendmws(MWS,model)
+%function appendmws(MWS,model)
+%
+% Appends model workspace with variables 
+% given by the fields of the structure MWS.  
+%
+% ! Warning !  Will overwrite indentically named variables without warning!
+%
+% Inputs:
+%  MWS   - Model Workspace Structure, contains simulation parameters
+%  model - Name of model to load into, default is 'gtm_design'
+%
+
+% $Id: appendmws.m 4852 2013-08-06 22:12:54Z cox $
+% Author: Austin Murch (Austin.M.Murch@nasa.gov)
+% Adpated from code by Dave Cox (d.e.cox@larc.nasa.gov)
+
+
+% By default use the bdroot model
+if ( ~exist('model','var') || isempty(model) ),
+  model=bdroot;
+end
+
+mws=get_param(model,'modelworkspace');
+
+fn=fieldnames(MWS);
+for i=1:length(fn)
+    mws.assignin(fn{i},MWS.(fn{i}));
+end
+
+
+
+
+
+  
diff --git a/gtm_design/mfiles/clearmws.m b/gtm_design/mfiles/clearmws.m
index 10f3103..cac532b 100644
--- a/gtm_design/mfiles/clearmws.m
+++ b/gtm_design/mfiles/clearmws.m
@@ -1,21 +1,21 @@
-function clearmws(model);
-%function clearmws(model);
-%
-% Clears the current model workspace.
-%  
-% Inputs:
-%    model   - name of model, defaults to 'gtm_design'
-%  
-
-% d.e.cox@larc.nasa.gov
-% $Id: clearmws.m 4852 2013-08-06 22:12:54Z cox $
-
-% By default use the bdroot model
-if ( ~exist('model','var') || isempty(model) ),
-  model=bdroot;
-end
-
-mws=get_param(model,'modelworkspace');
-mws.clear;
-
-
+function clearmws(model);
+%function clearmws(model);
+%
+% Clears the current model workspace.
+%  
+% Inputs:
+%    model   - name of model, defaults to 'gtm_design'
+%  
+
+% d.e.cox@larc.nasa.gov
+% $Id: clearmws.m 4852 2013-08-06 22:12:54Z cox $
+
+% By default use the bdroot model
+if ( ~exist('model','var') || isempty(model) ),
+  model=bdroot;
+end
+
+mws=get_param(model,'modelworkspace');
+mws.clear;
+
+
diff --git a/gtm_design/mfiles/euler321.m b/gtm_design/mfiles/euler321.m
index 49a8428..dd76727 100644
--- a/gtm_design/mfiles/euler321.m
+++ b/gtm_design/mfiles/euler321.m
@@ -1,23 +1,23 @@
-function DCM = euler321(p);
-%function DCM = euler321(p);
-%
-% Creates directional cosine matrix for euler 321 rotation 
-% sequence. 
-% 
-% Inputs:
-%   p    - phi/theta/psi vector, in radians
-% 
-% Outputs:
-%   DCM  - Directional cosine matrix
-% 	
-
-% $Id: euler321.m 4852 2013-08-06 22:12:54Z cox $
-
-p=p(:);
-cp=cos(p);
-sp=sin(p);
-
-DCM=[cp(2)*cp(3), cp(2)*sp(3), -sp(2);
-sp(1)*sp(2)*cp(3)-cp(1)*sp(3),sp(1)*sp(2)*sp(3)+cp(1)*cp(3), sp(1)*cp(2);
-cp(1)*sp(2)*cp(3)+sp(1)*sp(3),cp(1)*sp(2)*sp(3)-sp(1)*cp(3), cp(1)*cp(2)];
-
+function DCM = euler321(p);
+%function DCM = euler321(p);
+%
+% Creates directional cosine matrix for euler 321 rotation 
+% sequence. 
+% 
+% Inputs:
+%   p    - phi/theta/psi vector, in radians
+% 
+% Outputs:
+%   DCM  - Directional cosine matrix
+% 	
+
+% $Id: euler321.m 4852 2013-08-06 22:12:54Z cox $
+
+p=p(:);
+cp=cos(p);
+sp=sin(p);
+
+DCM=[cp(2)*cp(3), cp(2)*sp(3), -sp(2);
+sp(1)*sp(2)*cp(3)-cp(1)*sp(3),sp(1)*sp(2)*sp(3)+cp(1)*cp(3), sp(1)*cp(2);
+cp(1)*sp(2)*cp(3)+sp(1)*sp(3),cp(1)*sp(2)*sp(3)-sp(1)*cp(3), cp(1)*cp(2)];
+
diff --git a/gtm_design/mfiles/grabmws.m b/gtm_design/mfiles/grabmws.m
index f60f26f..441c2b0 100644
--- a/gtm_design/mfiles/grabmws.m
+++ b/gtm_design/mfiles/grabmws.m
@@ -1,28 +1,28 @@
-function MWS = grabmws(model);
-%function MWS = grabmws(model);
-%
-% Grabs current model workspace and returns as
-% fields in a strucure
-% 
-% Inputs:
-%    model  - name of model, defaults to 'gtm_design'
-% 
-% Outputs:
-%    MWS    - simlation parameters from the model workspace
-% 
-
-% d.e.cox@larc.nasa.gov
-% $Id: grabmws.m 4852 2013-08-06 22:12:54Z cox $
-
-% By default use the bdroot model
-if ( ~exist('model','var') || isempty(model) ),
-  model=bdroot;
-end
-
-mws=get_param(model,'modelworkspace');
-mws_vars=mws.whos;
-for i=[1:length(mws_vars)],
-  MWS.(mws_vars(i).name)=mws.evalin(mws_vars(i).name);
-end
-
-
+function MWS = grabmws(model);
+%function MWS = grabmws(model);
+%
+% Grabs current model workspace and returns as
+% fields in a strucure
+% 
+% Inputs:
+%    model  - name of model, defaults to 'gtm_design'
+% 
+% Outputs:
+%    MWS    - simlation parameters from the model workspace
+% 
+
+% d.e.cox@larc.nasa.gov
+% $Id: grabmws.m 4852 2013-08-06 22:12:54Z cox $
+
+% By default use the bdroot model
+if ( ~exist('model','var') || isempty(model) ),
+  model=bdroot;
+end
+
+mws=get_param(model,'modelworkspace');
+mws_vars=mws.whos;
+for i=[1:length(mws_vars)],
+  MWS.(mws_vars(i).name)=mws.evalin(mws_vars(i).name);
+end
+
+
diff --git a/gtm_design/mfiles/gtm_design.slxc b/gtm_design/mfiles/gtm_design.slxc
new file mode 100644
index 0000000..f46da14
Binary files /dev/null and b/gtm_design/mfiles/gtm_design.slxc differ
diff --git a/gtm_design/mfiles/init_design.m b/gtm_design/mfiles/init_design.m
index f1c03a6..e0fcea6 100644
--- a/gtm_design/mfiles/init_design.m
+++ b/gtm_design/mfiles/init_design.m
@@ -1,143 +1,143 @@
-function MWS = init_design(vehicle,noise_on,fuelburn_on)
-%function MWS = init_design(vehicle);
-%
-% Creates Model Workspace Structure for gtm_design simulation.
-%
-% Inputs:
-%   vehicle     - type of vehicle, currently only 'GTM_T2'  (default = 'GTM_T2')
-%   noise_on    - flag for turning sensor noise models on  (default = 0)
-%   fuelburn_on - flag for turning fuel burn model on  (default = 0)
-%
-% Outputs:
-%   MWS     - Simulation parameters, to be loaded into Model Workspace
-%             with loadmws(MWS) or appendmws(MWS).
-%
-
-% d.e.cox@larc.nasa.gov
-% $Id: init_design.m 4852 2013-08-06 22:12:54Z cox $
-
-if ~exist('vehicle','var') || isempty(vehicle)
-    vehicle='GTM_T2';
-end
-
-if strcmp(vehicle,'GTM'),
-    vehicle='GTM_T2';  % Backwards compatible
-end
-
-if ~exist('noise_on','var') || isempty(noise_on)
-    noise_on = 0;
-end
-
-if ~exist('fuelburn_on','var') || isempty(fuelburn_on)
-    fuelburn_on = 0;
-end
-
-% Initialize MWS structure
-MWS=[];
-
-% Winds
-MWS.WindDir = 270;  %deg,   true heading wind is coming FROM
-MWS.WindSpd = 0;    %kts,   wind speed
-MWS.WindShearOn = 1;
-
-% Sim timestep
-MWS.Timestep = 1/200; % sec
-
-% Initial Conditions
-ini.Altitude =     800; % ft
-ini.tas      =     75;      % knots
-ini.alpha    =     3;        % deg
-ini.gamma    =     1;        % deg
-ini.stab     =     0;
-% ini.lat      =     37.02814654;  %Smithfield, deg
-% ini.lon      =    -76.58696588;  %Smithfield, deg
-% Flying at Wallops
-ini.lat        =  37.827926944;
-ini.lon        = -75.494061666;
-MWS.StatesInp  = SetICs(ini);
-MWS.runway_alt = 12; %ft
-
-MWS.fuel_in_use  = fuelburn_on; % fuel burn on/off
-
-
-switch(upper(vehicle))
-    case 'GTM_T2',
-      %MWS.Aero=load('T2_restricted_aerodatabase');
-      MWS.Aero=load('T2_polynomial_aerodatabase');
-
-       % Load aircraft parameters and noise models
-       MWS = AC_baseparams_T2(MWS);
-       MWS = NoiseParams_T2(MWS,noise_on); % Second parameter is on/off
-
-       % Generate Damage Models
-       MWS = DamageModels_T2(MWS);
-    otherwise,
-       error('No Parameters defined for vehicles: %s',vehicle);
-end
-
-% Surface/Throttle Offsets for trim
-MWS.bias.aileron    =0;
-MWS.bias.speedbrake =0;
-MWS.bias.elevator   =0;
-MWS.bias.flaps      =0;
-MWS.bias.rudder     =0;
-MWS.bias.stabilizer =0;
-MWS.bias.throttle   =20;
-MWS.bias.geardown   =0;
-
-% Basic Table Options
-MWS.LinearizeModeOn=0;
-MWS.TrimModeOn=0;
-MWS.SurfaceDeadbandOn=0;
-MWS.TrimWithStab=0;
-MWS.symmetric_aero_on  = 1;
-MWS.DamageCase = 0;
-MWS.DamageOnsetTime=10; % In secs.
-
-% Engine on/off parameters
-MWS.LengON = 1;
-MWS.RengON = 1;
-
-%Engine Ram Drag Coefficient
-MWS.ram_drag_coef = 0.010;
-
-% Set turbulence model parameters
-MWS = init_turbulence(MWS);
-
-% The following increase asymetric response at stall and
-% have been used for pilot training
-MWS.stall_cl_asym_enabled         = 0;     % turn the Cl asymetry at stall on/off
-%MWS.stall_cl_asym_add_uncertainty = 0;    % add random uncertainties to Cl
-%MWS.stall_cl_asym_vary_sign       = 0;    % apply random +/-1 gain to Cl
-
-
-%--------------------SubFunctions--------------------
-
-function [StatesInp] = SetICs(ini)
-StatesInp= zeros(12,1);
-
-knots2fps = 1/0.592487;
-d2r=pi/180;
-
-ub = (knots2fps)*ini.tas*cos(ini.alpha*d2r);
-wb = (knots2fps)*ini.tas*sin(ini.alpha*d2r);
-
-theta    = ini.alpha + ini.gamma;
-
-StatesInp(1)  = ub;              %  1 - ub (ft/s)
-StatesInp(2)  = 0;		 %  2 - vb (ft/s)
-StatesInp(3)  = wb;		 %  3 - wb (ft/s)
-StatesInp(4)  = 0;		 %  4 - pb (rad/s)
-StatesInp(5)  = 0;		 %  5 - qb (rad/s)
-StatesInp(6)  = 0;		 %  6 - rb (rad/s)
-StatesInp(7)  = ini.lat*d2r;	 %  7 - lat (rad), +north,
-StatesInp(8)  = ini.lon*d2r;	 %  8 - lon (rad), +east,
-StatesInp(9)  = ini.Altitude;	 %  9 - h (ft)
-StatesInp(10) = 0;	             % 10 - phi (rad)
-StatesInp(11) = theta*d2r;         % 11 - theta (rad)
-StatesInp(12) = pi/2;	           % 12 - psi (rad)
-
-
-
-
-
+function MWS = init_design(vehicle,noise_on,fuelburn_on)
+%function MWS = init_design(vehicle);
+%
+% Creates Model Workspace Structure for gtm_design simulation.
+%
+% Inputs:
+%   vehicle     - type of vehicle, currently only 'GTM_T2'  (default = 'GTM_T2')
+%   noise_on    - flag for turning sensor noise models on  (default = 0)
+%   fuelburn_on - flag for turning fuel burn model on  (default = 0)
+%
+% Outputs:
+%   MWS     - Simulation parameters, to be loaded into Model Workspace
+%             with loadmws(MWS) or appendmws(MWS).
+%
+
+% d.e.cox@larc.nasa.gov
+% $Id: init_design.m 4852 2013-08-06 22:12:54Z cox $
+
+if ~exist('vehicle','var') || isempty(vehicle)
+    vehicle='GTM_T2';
+end
+
+if strcmp(vehicle,'GTM'),
+    vehicle='GTM_T2';  % Backwards compatible
+end
+
+if ~exist('noise_on','var') || isempty(noise_on)
+    noise_on = 0;
+end
+
+if ~exist('fuelburn_on','var') || isempty(fuelburn_on)
+    fuelburn_on = 0;
+end
+
+% Initialize MWS structure
+MWS=[];
+
+% Winds
+MWS.WindDir = 270;  %deg,   true heading wind is coming FROM
+MWS.WindSpd = 0;    %kts,   wind speed
+MWS.WindShearOn = 1;
+
+% Sim timestep
+MWS.Timestep = 1/200; % sec
+
+% Initial Conditions
+ini.Altitude =     800; % ft
+ini.tas      =     75;      % knots
+ini.alpha    =     3;        % deg
+ini.gamma    =     1;        % deg
+ini.stab     =     0;
+% ini.lat      =     37.02814654;  %Smithfield, deg
+% ini.lon      =    -76.58696588;  %Smithfield, deg
+% Flying at Wallops
+ini.lat        =  37.827926944;
+ini.lon        = -75.494061666;
+MWS.StatesInp  = SetICs(ini);
+MWS.runway_alt = 12; %ft
+
+MWS.fuel_in_use  = fuelburn_on; % fuel burn on/off
+
+
+switch(upper(vehicle))
+    case 'GTM_T2',
+      %MWS.Aero=load('T2_restricted_aerodatabase');
+      MWS.Aero=load('T2_polynomial_aerodatabase');
+
+       % Load aircraft parameters and noise models
+       MWS = AC_baseparams_T2(MWS);
+       MWS = NoiseParams_T2(MWS,noise_on); % Second parameter is on/off
+
+       % Generate Damage Models
+       MWS = DamageModels_T2(MWS);
+    otherwise,
+       error('No Parameters defined for vehicles: %s',vehicle);
+end
+
+% Surface/Throttle Offsets for trim
+MWS.bias.aileron    =0;
+MWS.bias.speedbrake =0;
+MWS.bias.elevator   =0;
+MWS.bias.flaps      =0;
+MWS.bias.rudder     =0;
+MWS.bias.stabilizer =0;
+MWS.bias.throttle   =20;
+MWS.bias.geardown   =0;
+
+% Basic Table Options
+MWS.LinearizeModeOn=0;
+MWS.TrimModeOn=0;
+MWS.SurfaceDeadbandOn=0;
+MWS.TrimWithStab=0;
+MWS.symmetric_aero_on  = 1;
+MWS.DamageCase = 0;
+MWS.DamageOnsetTime=10; % In secs.
+
+% Engine on/off parameters
+MWS.LengON = 1;
+MWS.RengON = 1;
+
+%Engine Ram Drag Coefficient
+MWS.ram_drag_coef = 0.010;
+
+% Set turbulence model parameters
+MWS = init_turbulence(MWS);
+
+% The following increase asymetric response at stall and
+% have been used for pilot training
+MWS.stall_cl_asym_enabled         = 0;     % turn the Cl asymetry at stall on/off
+%MWS.stall_cl_asym_add_uncertainty = 0;    % add random uncertainties to Cl
+%MWS.stall_cl_asym_vary_sign       = 0;    % apply random +/-1 gain to Cl
+
+
+%--------------------SubFunctions--------------------
+
+function [StatesInp] = SetICs(ini)
+StatesInp= zeros(12,1);
+
+knots2fps = 1/0.592487;
+d2r=pi/180;
+
+ub = (knots2fps)*ini.tas*cos(ini.alpha*d2r);
+wb = (knots2fps)*ini.tas*sin(ini.alpha*d2r);
+
+theta    = ini.alpha + ini.gamma;
+
+StatesInp(1)  = ub;              %  1 - ub (ft/s)
+StatesInp(2)  = 0;		 %  2 - vb (ft/s)
+StatesInp(3)  = wb;		 %  3 - wb (ft/s)
+StatesInp(4)  = 0;		 %  4 - pb (rad/s)
+StatesInp(5)  = 0;		 %  5 - qb (rad/s)
+StatesInp(6)  = 0;		 %  6 - rb (rad/s)
+StatesInp(7)  = ini.lat*d2r;	 %  7 - lat (rad), +north,
+StatesInp(8)  = ini.lon*d2r;	 %  8 - lon (rad), +east,
+StatesInp(9)  = ini.Altitude;	 %  9 - h (ft)
+StatesInp(10) = 0;	             % 10 - phi (rad)
+StatesInp(11) = theta*d2r;         % 11 - theta (rad)
+StatesInp(12) = pi/2;	           % 12 - psi (rad)
+
+
+
+
+
diff --git a/gtm_design/mfiles/init_spin_design.m b/gtm_design/mfiles/init_spin_design.m
new file mode 100644
index 0000000..1f15f6f
--- /dev/null
+++ b/gtm_design/mfiles/init_spin_design.m
@@ -0,0 +1,168 @@
+function MWS = init_spin_design(vehicle,noise_on,fuelburn_on)
+%
+% Initializes the Model Workspace Structure for simulating a targeted oscillatory 
+% spin condition.
+%
+% Inputs:
+%   vehicle      - Type of vehicle, currently only 'GTM_T2' (default = 'GTM_T2').
+%   noise_on     - Flag for turning sensor noise models on (default = 0).
+%   fuelburn_on  - Flag for turning the fuel burn model on (default = 0).
+%
+% Outputs:
+%   MWS          - Simulation parameters to be loaded into the Model Workspace 
+%                  with loadmws(MWS) or appendmws(MWS).
+%
+%
+% This function initializes the simulation environment for an oscillatory spin condition 
+% of the specified vehicle type. It sets up parameters such as wind direction and speed, 
+% simulation timestep, initial conditions, aircraft dynamics, noise models, damage models,
+% and other simulation options.
+%
+% The resulting MWS structure contains all the necessary parameters for simulating the 
+% targeted oscillatory spin condition.
+%
+%% About
+%
+% Author:     Omar Mourad
+% Email:      <mailto: st181901@stud.uni-stuttgart.de>
+% Created:    20.03.2024
+
+if ~exist('vehicle','var') || isempty(vehicle)
+    vehicle='GTM_T2';
+end
+
+if strcmp(vehicle,'GTM'),
+    vehicle='GTM_T2';  % Backwards compatible
+end
+
+if ~exist('noise_on','var') || isempty(noise_on)
+    noise_on = 0;
+end
+
+if ~exist('fuelburn_on','var') || isempty(fuelburn_on)
+    fuelburn_on = 0;
+end
+
+% Initialize MWS structure
+MWS=[];
+
+% Winds
+MWS.WindDir = 270;  %deg,   true heading wind is coming FROM
+MWS.WindSpd = 0;    %kts,   wind speed
+MWS.WindShearOn = 1;
+
+% Sim timestep
+MWS.Timestep = 1/200; % sec
+
+% Oscillatory Spin Initial Conditions
+ini.Altitude =     0;          % ft
+
+ini.tas      =     30.42;      % (meter/sec)
+
+ini.alpha    =     34.032;           % deg
+ini.beta     =     9.934;           % deg
+
+ini.p        =     -195.7;           % (deg/s)
+ini.q        =       15.56;           % (deg/s)
+ini.r        =      -137;           % (deg/s)
+
+ini.bank       =       5.069;            % deg
+ini.pitch      =       -54.72;            % deg
+ini.heading    =           0;            % deg
+
+
+ini.p = -250; 
+
+ini.stab     =     0;
+
+% Flying at Wallops
+ini.lat        =  37.827926944;            % deg
+ini.lon        = -75.494061666;            % deg
+MWS.StatesInp  = SetICs(ini);
+MWS.runway_alt = 12;                       %ft
+
+MWS.fuel_in_use  = fuelburn_on;            % fuel burn on/off
+
+
+switch(upper(vehicle))
+    case 'GTM_T2',
+      %MWS.Aero=load('T2_restricted_aerodatabase');
+      MWS.Aero=load('T2_polynomial_aerodatabase');
+
+       % Load aircraft parameters and noise models
+       MWS = AC_baseparams_T2(MWS);
+       MWS = NoiseParams_T2(MWS,noise_on); % Second parameter is on/off
+
+       % Generate Damage Models
+       MWS = DamageModels_T2(MWS);
+    otherwise,
+       error('No Parameters defined for vehicles: %s',vehicle);
+end
+
+% Surface/Throttle Offsets for trim
+MWS.bias.aileron    =0;
+MWS.bias.speedbrake =0;
+MWS.bias.elevator   =0;
+MWS.bias.flaps      =0;
+MWS.bias.rudder     =0;
+MWS.bias.stabilizer =0;
+MWS.bias.throttle   =20;
+MWS.bias.geardown   =0;
+
+% Basic Table Options
+MWS.LinearizeModeOn=0;
+MWS.TrimModeOn=0;
+MWS.SurfaceDeadbandOn=0;
+MWS.TrimWithStab=0;
+MWS.symmetric_aero_on  = 1;
+MWS.DamageCase = 0;
+MWS.DamageOnsetTime=10; % In secs.
+
+% Engine on/off parameters
+MWS.LengON = 1;
+MWS.RengON = 1;
+
+%Engine Ram Drag Coefficient
+MWS.ram_drag_coef = 0.010;
+
+% Set turbulence model parameters
+MWS = init_turbulence(MWS);
+
+% The following increase asymetric response at stall and
+% have been used for pilot training
+MWS.stall_cl_asym_enabled         = 0;     % turn the Cl asymetry at stall on/off
+%MWS.stall_cl_asym_add_uncertainty = 0;    % add random uncertainties to Cl
+%MWS.stall_cl_asym_vary_sign       = 0;    % apply random +/-1 gain to Cl
+
+
+%--------------------SubFunctions--------------------
+
+function [StatesInp] = SetICs(ini)
+StatesInp= zeros(12,1);
+
+metersPerSecond2fps = 3.28084;
+d2r=pi/180;
+
+ub = (metersPerSecond2fps)*ini.tas*cos(ini.alpha*d2r)*cos(ini.beta*d2r);
+vb = (metersPerSecond2fps)*ini.tas*sin(ini.beta*d2r);
+wb = (metersPerSecond2fps)*ini.tas*sin(ini.alpha*d2r)*cos(ini.beta*d2r);
+
+ini.gamma    = ini.pitch - ini.alpha;
+
+StatesInp(1)  = ub;      %  1 - ub (ft/s)
+StatesInp(2)  = vb;		 %  2 - vb (ft/s)
+StatesInp(3)  = wb;		 %  3 - wb (ft/s)
+StatesInp(4)  = ini.p*d2r;		 %  4 - pb (rad/s)
+StatesInp(5)  = ini.q*d2r;		 %  5 - qb (rad/s)
+StatesInp(6)  = ini.r*d2r;		 %  6 - rb (rad/s)
+StatesInp(7)  = ini.lat*d2r;	 %  7 - lat (rad), +north,
+StatesInp(8)  = ini.lon*d2r;	 %  8 - lon (rad), +east,
+StatesInp(9)  = ini.Altitude;	 %  9 - h (ft)
+StatesInp(10) = ini.bank*d2r;	      % 10 - phi (rad)
+StatesInp(11) = ini.pitch*d2r;        % 11 - theta (rad)
+StatesInp(12) = ini.heading*d2r;	  % 12 - psi (rad)
+
+
+
+
+
diff --git a/gtm_design/mfiles/init_spiral_design.m b/gtm_design/mfiles/init_spiral_design.m
new file mode 100644
index 0000000..3ffda0b
--- /dev/null
+++ b/gtm_design/mfiles/init_spiral_design.m
@@ -0,0 +1,167 @@
+function MWS = init_spiral_design(vehicle,noise_on,fuelburn_on)
+%
+% Initializes the Model Workspace Structure for simulating a targeted steep spiral condition.
+%
+% Inputs:
+%   vehicle      - Type of vehicle, currently only 'GTM_T2' (default = 'GTM_T2').
+%   noise_on     - Flag for turning sensor noise models on (default = 0).
+%   fuelburn_on  - Flag for turning the fuel burn model on (default = 0).
+%
+% Outputs:
+%   MWS          - Simulation parameters to be loaded into the Model Workspace 
+%                  with loadmws(MWS) or appendmws(MWS).
+%
+%
+% This function initializes the simulation environment for a steep spiral condition 
+% of the specified vehicle type. It sets up parameters such as wind direction and speed, 
+% simulation timestep, initial conditions, aircraft dynamics, noise models, damage models,
+% and other simulation options.
+%
+% The resulting MWS structure contains all the necessary parameters for simulating the 
+% targeted steep spiral condition.
+%
+%% About
+%
+% Author:     Omar Mourad
+% Email:      <mailto: st181901@stud.uni-stuttgart.de>
+% Created:    15.03.2024
+
+if ~exist('vehicle','var') || isempty(vehicle)
+    vehicle='GTM_T2';
+end
+
+if strcmp(vehicle,'GTM'),
+    vehicle='GTM_T2';  % Backwards compatible
+end
+
+if ~exist('noise_on','var') || isempty(noise_on)
+    noise_on = 0;
+end
+
+if ~exist('fuelburn_on','var') || isempty(fuelburn_on)
+    fuelburn_on = 0;
+end
+
+% Initialize MWS structure
+MWS=[];
+
+% Winds
+MWS.WindDir = 270;  %deg,   true heading wind is coming FROM
+MWS.WindSpd = 0;    %kts,   wind speed
+MWS.WindShearOn = 1;
+
+% Sim timestep
+MWS.Timestep = 1/200; % sec
+
+% Steep Spiral Initial Conditions
+ini.Altitude =     0;          % ft
+
+ini.tas      =     43.32;      % (meter/sec)
+
+ini.alpha    =     20.05;           % deg
+ini.beta     =     -7.55;           % deg
+
+ini.p        =     -122.9;           % (deg/s)
+ini.q        =       49.1;           % (deg/s)
+ini.r        =      -45.4;           % (deg/s)
+
+ini.bank       =       -47.3;            % deg
+ini.pitch      =       -61.5;            % deg
+ini.heading    =           0;            % deg
+
+
+ini.p = -250; 
+
+ini.stab     =     0;
+
+% Flying at Wallops
+ini.lat        =  37.827926944;            % deg
+ini.lon        = -75.494061666;            % deg
+MWS.StatesInp  = SetICs(ini);
+MWS.runway_alt = 12;                       %ft
+
+MWS.fuel_in_use  = fuelburn_on;            % fuel burn on/off
+
+
+switch(upper(vehicle))
+    case 'GTM_T2',
+      %MWS.Aero=load('T2_restricted_aerodatabase');
+      MWS.Aero=load('T2_polynomial_aerodatabase');
+
+       % Load aircraft parameters and noise models
+       MWS = AC_baseparams_T2(MWS);
+       MWS = NoiseParams_T2(MWS,noise_on); % Second parameter is on/off
+
+       % Generate Damage Models
+       MWS = DamageModels_T2(MWS);
+    otherwise,
+       error('No Parameters defined for vehicles: %s',vehicle);
+end
+
+% Surface/Throttle Offsets for trim
+MWS.bias.aileron    =0;
+MWS.bias.speedbrake =0;
+MWS.bias.elevator   =0;
+MWS.bias.flaps      =0;
+MWS.bias.rudder     =0;
+MWS.bias.stabilizer =0;
+MWS.bias.throttle   =20;
+MWS.bias.geardown   =0;
+
+% Basic Table Options
+MWS.LinearizeModeOn=0;
+MWS.TrimModeOn=0;
+MWS.SurfaceDeadbandOn=0;
+MWS.TrimWithStab=0;
+MWS.symmetric_aero_on  = 1;
+MWS.DamageCase = 0;
+MWS.DamageOnsetTime=10; % In secs.
+
+% Engine on/off parameters
+MWS.LengON = 1;
+MWS.RengON = 1;
+
+%Engine Ram Drag Coefficient
+MWS.ram_drag_coef = 0.010;
+
+% Set turbulence model parameters
+MWS = init_turbulence(MWS);
+
+% The following increase asymetric response at stall and
+% have been used for pilot training
+MWS.stall_cl_asym_enabled         = 0;     % turn the Cl asymetry at stall on/off
+%MWS.stall_cl_asym_add_uncertainty = 0;    % add random uncertainties to Cl
+%MWS.stall_cl_asym_vary_sign       = 0;    % apply random +/-1 gain to Cl
+
+
+%--------------------SubFunctions--------------------
+
+function [StatesInp] = SetICs(ini)
+StatesInp= zeros(12,1);
+
+metersPerSecond2fps = 3.28084;
+d2r=pi/180;
+
+ub = (metersPerSecond2fps)*ini.tas*cos(ini.alpha*d2r)*cos(ini.beta*d2r);
+vb = (metersPerSecond2fps)*ini.tas*sin(ini.beta*d2r);
+wb = (metersPerSecond2fps)*ini.tas*sin(ini.alpha*d2r)*cos(ini.beta*d2r);
+
+ini.gamma    = ini.pitch - ini.alpha;
+
+StatesInp(1)  = ub;      %  1 - ub (ft/s)
+StatesInp(2)  = vb;		 %  2 - vb (ft/s)
+StatesInp(3)  = wb;		 %  3 - wb (ft/s)
+StatesInp(4)  = ini.p*d2r;		 %  4 - pb (rad/s)
+StatesInp(5)  = ini.q*d2r;		 %  5 - qb (rad/s)
+StatesInp(6)  = ini.r*d2r;		 %  6 - rb (rad/s)
+StatesInp(7)  = ini.lat*d2r;	 %  7 - lat (rad), +north,
+StatesInp(8)  = ini.lon*d2r;	 %  8 - lon (rad), +east,
+StatesInp(9)  = ini.Altitude;	 %  9 - h (ft)
+StatesInp(10) = ini.bank*d2r;	      % 10 - phi (rad)
+StatesInp(11) = ini.pitch*d2r;        % 11 - theta (rad)
+StatesInp(12) = ini.heading*d2r;	  % 12 - psi (rad)
+
+
+
+
+
diff --git a/gtm_design/mfiles/linmodel.m b/gtm_design/mfiles/linmodel.m
index 5b2769a..6952d3f 100644
--- a/gtm_design/mfiles/linmodel.m
+++ b/gtm_design/mfiles/linmodel.m
@@ -1,217 +1,217 @@
-function [sys,londyn,latdyn] = linmodel(MWS,vabflag,use_all_inputs,Ts)
-%function [sys,londyn,latdyn] = linmodel(MWS,vabflag,use_all_inputs,Ts)
-%
-% Linearize gtm_design simulation at the current trim point.  The
-% linear model is continuous by default, discrete if the 'Ts' argument
-% is non-zero.
-%
-% Inputs:
-%   MWS    - simulation parameters, defaults to pre-loaded model workspace
-%   vabflag - flag to get linear models in terms of V, alpha, beta (0 default)
-%   use_all_inputs -for using full set of inputs rather than groups (0 default)
-%   Ts      - Timestep for generating a discrete linear model (0 default)
-%
-% Outputs:
-%   sys    - 6dof system with control surface inputs/state outputs
-%   londyn - Approximate (4th order) longitudinal dynamics
-%   latdyn - Approximate (4th order) lateral dynamics
-%
-
-% d.e.cox@larc.nasa.gov
-% $Id: linmodel.m 4852 2013-08-06 22:12:54Z cox $
-
-%% Parse Input Arguments
-% Load new model workspace if supplied
-if (exist('MWS','var') && ~isempty(MWS))
-    if isstruct(MWS)
-        loadmws(MWS,'gtm_design');
-    else
-        error('MWS input must be a data structure')
-    end
-end
-
-% Set optional flags
-if ~exist('vabflag','var') || isempty(vabflag),
-    vabflag = 0;
-end
-
-if ~exist('use_all_inputs','var') || isempty(use_all_inputs)
-    use_all_inputs = 0;
-end
-
-Args = [];
-if ~exist('Ts','var') || isempty(Ts) || Ts <= 0
-    Ts   = 0;
-    Args = 'IgnoreDiscreteStates';
-end
-
-%% Initializations and sim setup
-% Set Inline Params to off.
-% If not "off", trim results degrade (why??)
-dirtyflag =get_param('gtm_design','Dirty');
-inlineflag=get_param('gtm_design','InlineParams');
-set_param('gtm_design','InlineParams','off');
-
-% Grab state names
-[tmp1,tmp2,Statename]=gtm_design;
-
-% Get index for Equations of Motion states
-EOM =find(strcmp(Statename,['gtm_design/GTM_T2/EOM/Integrator'])==1);
-
-% Set LinearizationFlag
-appendmws(struct('LinearizeModeOn',1),'gtm_design');
-
-%% Run sim for 0.1 seconds.
-% This initializes some states that are not set by trimgtm, e.g.
-% memory blocks, filter states, etc.
-[t,Xtime,y]=sim('gtm_design',[0 0.1]);
-Xo=Xtime(end,:);
-
-%% Check to ensure non-accelerating set point
-tol=1e-3;
-delta_vel=Xtime(1,EOM(1:6))-Xtime(end,EOM(1:6));
-if max(abs(delta_vel))>tol,
-    fprintf(1,'\n  at t=%3.2f sec Xvel=[%5.2e %5.2e %5.2e %5.2e %5.2e %5.2e]',t(1),  Xtime(1,EOM(1:6)));
-    fprintf(1,'\n  at t=%3.2f sec Xvel=[%5.2e %5.2e %5.2e %5.2e %5.2e %5.2e]\n',t(end),Xtime(end,EOM(1:6)));
-    warning('Model does not appear to be at a stationary point, deltaV=%3.2f',max(abs(delta_vel)));
-end
-
-%% Extract the Linear Model using dlinmod.  For Ts > 0, the model is discrete.
-warning('off','Simulink:tools:dlinmodIgnoreDiscreteStates');
-[a,b,c,d]=dlinmod('gtm_design',Ts,Xo,Args);
-sys=ss(a,b,c,d,Ts);
-warning('on','Simulink:tools:dlinmodIgnoreDiscreteStates');
-
-% Define states to retain, remove others
-keep=EOM(1)+[0:11]';
-elim=setdiff([1:length(a)]',keep);
-sys=modred(sys,elim,'del');
-
-% Remove some of the outputs that are used for trim only
-% First 6 are from AUX, last 12 are from the EOM
-sys=sys([2:4,7:18],:);
-
-% Remove firsst four inputs (trim inputs)
-sys=sys(:,[5:23]);
-
-% Convert units on selected outputs
-sys.c(1,:) = sys.c(1,:)*1.6878; % convert tas to feet per second
-sys.c(2,:) = sys.c(2,:)*pi/180; % convert alpha to radians
-sys.c(3,:) = sys.c(3,:)*pi/180; % convert beta to radians
-
-% Set names, hardwired from ordering in block diagram
-
-Inames={'ElevLOB','ElevLIB','ElevROB','ElevRIB',              ...%  1  2  3  4
-        'AileronL',   'AileronR',                             ...%  5  6
-        'RudderUpper','RudderLower',                          ...%  7  8
-        'SpoilerLIB', 'SpoilerLOB', 'SpoilerRIB','SpoilerROB',...%  9 10 11 12
-        'FlapsLIB',   'FlapsLOB',   'FlapsRIB',  'FlapsROB',  ...% 13 14 15 16
-        'Stabilizer',                                         ...% 17
-        'L Throttle', 'R Throttle'}; 		                 % 18 19
-
-Snames={'u',   'v',     'w',    ... %  1  2  3
-        'p',   'q',     'r',    ... %  4  5  6
-        'Lat', 'Lon',   'Alt',  ... %  7  8  9
-        'phi', 'theta', 'psi'};     % 10 11 12
-
-Onames={'tas', 'alpha', 'beta', ... %  1  2  3
-        'u',   'v',     'w',    ... %  4  5  6
-        'p',   'q',     'r',    ... %  7  8  9
-        'Lat', 'Lon',   'Alt',  ... % 10 11 12
-        'phi', 'theta', 'psi'};     % 13 14 15
-
-set(sys,'Statename',Snames,'Outputname',Onames,'Inputname',Inames);
-
-% Create 4th order longitudinal/lateral models
-Xlon=[1 3 5 11]; % States to keep (u,w,q,theta)
-Xlat=[2 4 6 10]; % States to keep (v,p,r,phi)
-
-
-if use_all_inputs
-
-  % Keep the inputs using indices for the full set of control surface inputs
-  Ilon=[1 2 3 4 9 10 11 12 18 19]; % inputs to keep: elev, spoilers, thrott
-  Ilat=[5 6 7 8 9 10 11 12];       % inputs to keep: ail, rudder, spoilers
-
-else
-
-  % Keep the inputs using indices for a subset of control surface inputs
-
-  % Model inputs are: elevator, aileron, rudder, LSpoiler, RSpoiler, 
-  %                   Flap, Stab, LThrottle, RThrottle
-  % 
-  % Reduce/group the full set of control surfaces into the subsets:
-  %
-  sys.b(:,1) = sys.b(:,1)  + sys.b(:,2) + sys.b(:,3) + sys.b(:,4);    %1 elev
-  sys.b(:,2) = sys.b(:,5)  + sys.b(:,6);                              %2 ail
-  sys.b(:,3) = sys.b(:,7)  + sys.b(:,8);                              %3 rud
-  sys.b(:,4) = sys.b(:,9)  + sys.b(:,10);                             %4 Lspoil
-  sys.b(:,5) = sys.b(:,11) + sys.b(:,12);                             %5 Rspoil
-  sys.b(:,6) = sys.b(:,13) + sys.b(:,14) + sys.b(:,15) + sys.b(:,16); %6 flap
-  sys.b(:,7) = sys.b(:,17);                                           %7 stab
-  sys.b(:,8) = sys.b(:,18);                                           %8 Lthrot
-  sys.b(:,9) = sys.b(:,19);                                           %9 Rthrot
-
-  sys.d(:,1) = sys.d(:,1)  + sys.d(:,2) + sys.d(:,3) + sys.d(:,4);    %1 elev  
-  sys.d(:,2) = sys.d(:,5)  + sys.d(:,6);			      %2 ail   
-  sys.d(:,3) = sys.d(:,7)  + sys.d(:,8);			      %3 rud   
-  sys.d(:,4) = sys.d(:,9)  + sys.d(:,10);			      %4 Lspoil
-  sys.d(:,5) = sys.d(:,11) + sys.d(:,12);			      %5 Rspoil
-  sys.d(:,6) = sys.d(:,13) + sys.d(:,14) + sys.d(:,15) + sys.d(:,16); %6 flap  
-  sys.d(:,7) = sys.d(:,17);					      %7 stab  
-  sys.d(:,8) = sys.d(:,18);					      %8 Lthrot
-  sys.d(:,9) = sys.d(:,19);					      %9 Rthrot
-
-  Ilon=[1 4 5 8 9];    % inputs to keep: elev, spoilers, throttle
-  Ilat=[2 3 4 5];      % inputs to keep: aileron, rudder, spoilers
-
-  sys(:,10:19) = []; % Remove the unused inputs
-  
-  Inames={'Elevator',  'Aileron','Rudder',... % 1 2 3
-          'L Spoiler', 'R Spoiler',       ... % 4 5
-          'Flaps',     'Stabilizer',      ... % 6 7
-          'L Throttle','R Throttle'};         % 8-9
-
-  set(sys,'Inputname',Inames);
-
-end
-
-
-if vabflag
-    Ylon = [1 2 8 14]; % outputs to keep (tas,alpha,q,theta)
-    Ylat = [3 7 9 13]; % outputs to keep (beta,p,r,phi)
-    londyn=modred(sys(Ylon,Ilon),setdiff(1:12,Xlon),'del');
-    londyn.a = londyn.c*londyn.a*inv(londyn.c);
-    londyn.b = londyn.c*londyn.b;
-    londyn.c = eye(4);
-    set(londyn,'Statename',Onames([1 2 8 14]))
-    
-    latdyn=modred(sys(Ylat,Ilat),setdiff(1:12,Xlat),'del');
-    latdyn.a = latdyn.c*latdyn.a*inv(latdyn.c);
-    latdyn.b = latdyn.c*latdyn.b;
-    latdyn.c = eye(4);
-    set(latdyn,'Statename',Onames([3 7 9 13]))
-    
-    sys = sys([1:3, 7:15],:);
-    sys.a = sys.c*sys.a*inv(sys.c);
-    sys.b = sys.c*sys.b;
-    sys.c = eye(12);
-    set(sys,'Statename',Onames([1:3, 7:15]))    
-    
-else
-    Ylon = [4 6 8 14]; % outputs to keep (u,w,q,theta)
-    Ylat = [5 7 9 13]; % outputs to keep (v,p,r,phi)  
-    londyn=modred(sys(Ylon,Ilon),setdiff(1:12,Xlon),'del');
-    latdyn=modred(sys(Ylat,Ilat),setdiff(1:12,Xlat),'del');
-    sys = sys(4:15,:);
-end
-
-%% Return to initial parameters
-set_param('gtm_design','InlineParams',inlineflag);
-set_param('gtm_design','Dirty',dirtyflag);
-
-% reset Linearization Flag
-appendmws(struct('LinearizeModeOn',0),'gtm_design');
-
-
-
+function [sys,londyn,latdyn] = linmodel(MWS,vabflag,use_all_inputs,Ts)
+%function [sys,londyn,latdyn] = linmodel(MWS,vabflag,use_all_inputs,Ts)
+%
+% Linearize gtm_design simulation at the current trim point.  The
+% linear model is continuous by default, discrete if the 'Ts' argument
+% is non-zero.
+%
+% Inputs:
+%   MWS    - simulation parameters, defaults to pre-loaded model workspace
+%   vabflag - flag to get linear models in terms of V, alpha, beta (0 default)
+%   use_all_inputs -for using full set of inputs rather than groups (0 default)
+%   Ts      - Timestep for generating a discrete linear model (0 default)
+%
+% Outputs:
+%   sys    - 6dof system with control surface inputs/state outputs
+%   londyn - Approximate (4th order) longitudinal dynamics
+%   latdyn - Approximate (4th order) lateral dynamics
+%
+
+% d.e.cox@larc.nasa.gov
+% $Id: linmodel.m 4852 2013-08-06 22:12:54Z cox $
+
+%% Parse Input Arguments
+% Load new model workspace if supplied
+if (exist('MWS','var') && ~isempty(MWS))
+    if isstruct(MWS)
+        loadmws(MWS,'gtm_design');
+    else
+        error('MWS input must be a data structure')
+    end
+end
+
+% Set optional flags
+if ~exist('vabflag','var') || isempty(vabflag),
+    vabflag = 0;
+end
+
+if ~exist('use_all_inputs','var') || isempty(use_all_inputs)
+    use_all_inputs = 0;
+end
+
+Args = [];
+if ~exist('Ts','var') || isempty(Ts) || Ts <= 0
+    Ts   = 0;
+    Args = 'IgnoreDiscreteStates';
+end
+
+%% Initializations and sim setup
+% Set Inline Params to off.
+% If not "off", trim results degrade (why??)
+dirtyflag =get_param('gtm_design','Dirty');
+inlineflag=get_param('gtm_design','InlineParams');
+set_param('gtm_design','InlineParams','off');
+
+% Grab state names
+[tmp1,tmp2,Statename]=gtm_design;
+
+% Get index for Equations of Motion states
+EOM =find(strcmp(Statename,['gtm_design/GTM_T2/EOM/Integrator'])==1);
+
+% Set LinearizationFlag
+appendmws(struct('LinearizeModeOn',1),'gtm_design');
+
+%% Run sim for 0.1 seconds.
+% This initializes some states that are not set by trimgtm, e.g.
+% memory blocks, filter states, etc.
+[t,Xtime,y]=sim('gtm_design',[0 0.1]);
+Xo=Xtime(end,:);
+
+%% Check to ensure non-accelerating set point
+tol=1e-3;
+delta_vel=Xtime(1,EOM(1:6))-Xtime(end,EOM(1:6));
+if max(abs(delta_vel))>tol,
+    fprintf(1,'\n  at t=%3.2f sec Xvel=[%5.2e %5.2e %5.2e %5.2e %5.2e %5.2e]',t(1),  Xtime(1,EOM(1:6)));
+    fprintf(1,'\n  at t=%3.2f sec Xvel=[%5.2e %5.2e %5.2e %5.2e %5.2e %5.2e]\n',t(end),Xtime(end,EOM(1:6)));
+    warning('Model does not appear to be at a stationary point, deltaV=%3.2f',max(abs(delta_vel)));
+end
+
+%% Extract the Linear Model using dlinmod.  For Ts > 0, the model is discrete.
+warning('off','Simulink:tools:dlinmodIgnoreDiscreteStates');
+[a,b,c,d]=dlinmod('gtm_design',Ts,Xo,Args);
+sys=ss(a,b,c,d,Ts);
+warning('on','Simulink:tools:dlinmodIgnoreDiscreteStates');
+
+% Define states to retain, remove others
+keep=EOM(1)+[0:11]';
+elim=setdiff([1:length(a)]',keep);
+sys=modred(sys,elim,'del');
+
+% Remove some of the outputs that are used for trim only
+% First 6 are from AUX, last 12 are from the EOM
+sys=sys([2:4,7:18],:);
+
+% Remove firsst four inputs (trim inputs)
+sys=sys(:,[5:23]);
+
+% Convert units on selected outputs
+sys.c(1,:) = sys.c(1,:)*1.6878; % convert tas to feet per second
+sys.c(2,:) = sys.c(2,:)*pi/180; % convert alpha to radians
+sys.c(3,:) = sys.c(3,:)*pi/180; % convert beta to radians
+
+% Set names, hardwired from ordering in block diagram
+
+Inames={'ElevLOB','ElevLIB','ElevROB','ElevRIB',              ...%  1  2  3  4
+        'AileronL',   'AileronR',                             ...%  5  6
+        'RudderUpper','RudderLower',                          ...%  7  8
+        'SpoilerLIB', 'SpoilerLOB', 'SpoilerRIB','SpoilerROB',...%  9 10 11 12
+        'FlapsLIB',   'FlapsLOB',   'FlapsRIB',  'FlapsROB',  ...% 13 14 15 16
+        'Stabilizer',                                         ...% 17
+        'L Throttle', 'R Throttle'}; 		                 % 18 19
+
+Snames={'u',   'v',     'w',    ... %  1  2  3
+        'p',   'q',     'r',    ... %  4  5  6
+        'Lat', 'Lon',   'Alt',  ... %  7  8  9
+        'phi', 'theta', 'psi'};     % 10 11 12
+
+Onames={'tas', 'alpha', 'beta', ... %  1  2  3
+        'u',   'v',     'w',    ... %  4  5  6
+        'p',   'q',     'r',    ... %  7  8  9
+        'Lat', 'Lon',   'Alt',  ... % 10 11 12
+        'phi', 'theta', 'psi'};     % 13 14 15
+
+set(sys,'Statename',Snames,'Outputname',Onames,'Inputname',Inames);
+
+% Create 4th order longitudinal/lateral models
+Xlon=[1 3 5 11]; % States to keep (u,w,q,theta)
+Xlat=[2 4 6 10]; % States to keep (v,p,r,phi)
+
+
+if use_all_inputs
+
+  % Keep the inputs using indices for the full set of control surface inputs
+  Ilon=[1 2 3 4 9 10 11 12 18 19]; % inputs to keep: elev, spoilers, thrott
+  Ilat=[5 6 7 8 9 10 11 12];       % inputs to keep: ail, rudder, spoilers
+
+else
+
+  % Keep the inputs using indices for a subset of control surface inputs
+
+  % Model inputs are: elevator, aileron, rudder, LSpoiler, RSpoiler, 
+  %                   Flap, Stab, LThrottle, RThrottle
+  % 
+  % Reduce/group the full set of control surfaces into the subsets:
+  %
+  sys.b(:,1) = sys.b(:,1)  + sys.b(:,2) + sys.b(:,3) + sys.b(:,4);    %1 elev
+  sys.b(:,2) = sys.b(:,5)  + sys.b(:,6);                              %2 ail
+  sys.b(:,3) = sys.b(:,7)  + sys.b(:,8);                              %3 rud
+  sys.b(:,4) = sys.b(:,9)  + sys.b(:,10);                             %4 Lspoil
+  sys.b(:,5) = sys.b(:,11) + sys.b(:,12);                             %5 Rspoil
+  sys.b(:,6) = sys.b(:,13) + sys.b(:,14) + sys.b(:,15) + sys.b(:,16); %6 flap
+  sys.b(:,7) = sys.b(:,17);                                           %7 stab
+  sys.b(:,8) = sys.b(:,18);                                           %8 Lthrot
+  sys.b(:,9) = sys.b(:,19);                                           %9 Rthrot
+
+  sys.d(:,1) = sys.d(:,1)  + sys.d(:,2) + sys.d(:,3) + sys.d(:,4);    %1 elev  
+  sys.d(:,2) = sys.d(:,5)  + sys.d(:,6);			      %2 ail   
+  sys.d(:,3) = sys.d(:,7)  + sys.d(:,8);			      %3 rud   
+  sys.d(:,4) = sys.d(:,9)  + sys.d(:,10);			      %4 Lspoil
+  sys.d(:,5) = sys.d(:,11) + sys.d(:,12);			      %5 Rspoil
+  sys.d(:,6) = sys.d(:,13) + sys.d(:,14) + sys.d(:,15) + sys.d(:,16); %6 flap  
+  sys.d(:,7) = sys.d(:,17);					      %7 stab  
+  sys.d(:,8) = sys.d(:,18);					      %8 Lthrot
+  sys.d(:,9) = sys.d(:,19);					      %9 Rthrot
+
+  Ilon=[1 4 5 8 9];    % inputs to keep: elev, spoilers, throttle
+  Ilat=[2 3 4 5];      % inputs to keep: aileron, rudder, spoilers
+
+  sys(:,10:19) = []; % Remove the unused inputs
+  
+  Inames={'Elevator',  'Aileron','Rudder',... % 1 2 3
+          'L Spoiler', 'R Spoiler',       ... % 4 5
+          'Flaps',     'Stabilizer',      ... % 6 7
+          'L Throttle','R Throttle'};         % 8-9
+
+  set(sys,'Inputname',Inames);
+
+end
+
+
+if vabflag
+    Ylon = [1 2 8 14]; % outputs to keep (tas,alpha,q,theta)
+    Ylat = [3 7 9 13]; % outputs to keep (beta,p,r,phi)
+    londyn=modred(sys(Ylon,Ilon),setdiff(1:12,Xlon),'del');
+    londyn.a = londyn.c*londyn.a*inv(londyn.c);
+    londyn.b = londyn.c*londyn.b;
+    londyn.c = eye(4);
+    set(londyn,'Statename',Onames([1 2 8 14]))
+    
+    latdyn=modred(sys(Ylat,Ilat),setdiff(1:12,Xlat),'del');
+    latdyn.a = latdyn.c*latdyn.a*inv(latdyn.c);
+    latdyn.b = latdyn.c*latdyn.b;
+    latdyn.c = eye(4);
+    set(latdyn,'Statename',Onames([3 7 9 13]))
+    
+    sys = sys([1:3, 7:15],:);
+    sys.a = sys.c*sys.a*inv(sys.c);
+    sys.b = sys.c*sys.b;
+    sys.c = eye(12);
+    set(sys,'Statename',Onames([1:3, 7:15]))    
+    
+else
+    Ylon = [4 6 8 14]; % outputs to keep (u,w,q,theta)
+    Ylat = [5 7 9 13]; % outputs to keep (v,p,r,phi)  
+    londyn=modred(sys(Ylon,Ilon),setdiff(1:12,Xlon),'del');
+    latdyn=modred(sys(Ylat,Ilat),setdiff(1:12,Xlat),'del');
+    sys = sys(4:15,:);
+end
+
+%% Return to initial parameters
+set_param('gtm_design','InlineParams',inlineflag);
+set_param('gtm_design','Dirty',dirtyflag);
+
+% reset Linearization Flag
+appendmws(struct('LinearizeModeOn',0),'gtm_design');
+
+
+
diff --git a/gtm_design/mfiles/loadmws.m b/gtm_design/mfiles/loadmws.m
index ef8dd71..290006c 100644
--- a/gtm_design/mfiles/loadmws.m
+++ b/gtm_design/mfiles/loadmws.m
@@ -1,36 +1,36 @@
-function loadmws(MWS,model);
-%function loadmws(MWS,model);
-%
-% Clears model workspace and replaces with variables 
-% given by the fields of the structure MWS
-%
-% Inputs:
-%  MWS   - Model Workspace Structure, contains simulation parameters
-%  model - Name of model to load into, default is 'gtm_design'
-%
-
-% d.e.cox@larc.nasa.gov
-% $Id: loadmws.m 4852 2013-08-06 22:12:54Z cox $
-
-% By default use the bdroot model
-if ( ~exist('model','var') || isempty(model) ),
-  model=bdroot;
-end
-
-mws=get_param(model,'modelworkspace');
-
-if ( ~exist('MWS') || isempty(MWS) ),
-  mws.clear;
-  return
-else
-  mws.clear;
-  fn=fieldnames(MWS);
-  for i=[1:length(fn)],
-    mws.assignin(fn{i},MWS.(fn{i}));
-  end
-end
-
-
-
-
-  
+function loadmws(MWS,model);
+%function loadmws(MWS,model);
+%
+% Clears model workspace and replaces with variables 
+% given by the fields of the structure MWS
+%
+% Inputs:
+%  MWS   - Model Workspace Structure, contains simulation parameters
+%  model - Name of model to load into, default is 'gtm_design'
+%
+
+% d.e.cox@larc.nasa.gov
+% $Id: loadmws.m 4852 2013-08-06 22:12:54Z cox $
+
+% By default use the bdroot model
+if ( ~exist('model','var') || isempty(model) ),
+  model=bdroot;
+end
+
+mws=get_param(model,'modelworkspace');
+
+if ( ~exist('MWS') || isempty(MWS) ),
+  mws.clear;
+  return
+else
+  mws.clear;
+  fn=fieldnames(MWS);
+  for i=[1:length(fn)],
+    mws.assignin(fn{i},MWS.(fn{i}));
+  end
+end
+
+
+
+
+  
diff --git a/gtm_design/mfiles/printStates.m b/gtm_design/mfiles/printStates.m
new file mode 100644
index 0000000..7e1be80
--- /dev/null
+++ b/gtm_design/mfiles/printStates.m
@@ -0,0 +1,66 @@
+function printStates(MWS)
+%% --------------------- Print Aircraft States Information --------------------------------
+%
+% This function, printStates, is designed to display the state information 
+% of an aircraft stored in a Model Workspace Structure (MWS). It prints 
+% various states such as velocity components, angular velocities, geographic 
+% coordinates, altitude, and Euler angles in a human-readable format.
+%
+% The function takes an input parameter MWS, which is expected to be a 
+% structure containing the aircraft state information. It then prints out 
+% each state variable along with its corresponding value.
+%
+% The printed states include:
+%   1. Linear velocity components (ub, vb, wb) in meters per second (m/s).
+%   2. Angular velocity components (pb, qb, rb) in degrees per second (deg/s).
+%   3. Geographic coordinates (lat, lon) in degrees.
+%   4. Altitude (h) in meters.
+%   5. Euler angles (phi, theta, psi) in degrees.
+%   6. Derived states (alpha, beta) in degrees, calculated from linear velocity components.
+%
+% This function is useful for debugging and monitoring the state of an aircraft 
+% during simulation or analysis.
+%
+% Syntax:
+%   printStates(MWS)
+%
+% Inputs:
+%   - MWS: Model Workspace Structure containing aircraft state information.
+%
+% Example:
+%   % Define an example Model Workspace Structure
+%   MWS.StatesInp = [100, 50, 20, 0.1, 0.2, 0.3, 0.5, 1.0, 5000, 0.2, 0.1, 0.3];
+%   
+%   % Print the state information
+%   printStates(MWS)
+%
+% Output:
+%   The function prints the state information to the MATLAB command window.
+%
+%% Author: Omar Mouard
+% Email: <mailto: st181901@stud.uni-stuttgart.de>
+% Date: 15.03.2024
+    
+    fps2metersPerSecond=0.3048;
+    
+    fprintf('Aircraft States Information:\n');
+    fprintf('  1. Linear Velocities:\n');
+    fprintf('     - u (m/s): %.2f\n', MWS.StatesInp(1)*fps2metersPerSecond);
+    fprintf('     - v (m/s): %.2f\n', MWS.StatesInp(2)*fps2metersPerSecond);
+    fprintf('     - w (m/s): %.2f\n', MWS.StatesInp(3)*fps2metersPerSecond);
+    fprintf('  2. Angular Velocities:\n');
+    fprintf('     - p (deg/s): %.2f\n', rad2deg(MWS.StatesInp(4)));
+    fprintf('     - q (deg/s): %.2f\n', rad2deg(MWS.StatesInp(5)));
+    fprintf('     - r (deg/s): %.2f\n', rad2deg(MWS.StatesInp(6)));
+    fprintf('  3. Geographic Coordinates:\n');
+    fprintf('     - latitude (deg): %.6f\n', rad2deg(MWS.StatesInp(7)));
+    fprintf('     - longitude (deg): %.6f\n', rad2deg(MWS.StatesInp(8)));
+    fprintf('  4. Altitude (m): %.2f\n', MWS.StatesInp(9));
+    fprintf('  5. Euler Angles:\n');
+    fprintf('     - roll (phi) (deg): %.2f\n', rad2deg(MWS.StatesInp(10)));
+    fprintf('     - pitch (theta) (deg): %.2f\n', rad2deg(MWS.StatesInp(11)));
+    fprintf('     - yaw (psi) (deg): %.2f\n', rad2deg(MWS.StatesInp(12)));
+    fprintf('  6. Derived States:\n');
+    fprintf('     - angle of attack (alpha) (deg): %.2f\n', atan2(MWS.StatesInp(3), MWS.StatesInp(1)) * 180/pi);
+    fprintf('     - sideslip angle (beta) (deg): %.2f\n', atan2(MWS.StatesInp(2), sqrt((MWS.StatesInp(1))^2 + (MWS.StatesInp(3))^2)) * 180/pi);
+end
diff --git a/gtm_design/mfiles/seteomic.m b/gtm_design/mfiles/seteomic.m
index 303b880..97749e9 100644
--- a/gtm_design/mfiles/seteomic.m
+++ b/gtm_design/mfiles/seteomic.m
@@ -1,107 +1,107 @@
-function MWS_out = seteomic(MWS,varargin);
-%function MWS_out = seteomic(MWS,statename,value,statename,value,...)
-%
-% This function allows the initial condtions of the gtm_design simulation's 
-% equation of motion (MWS.StatesInp) to be specified by named values. 
-%
-% Specifically statename is one of:
-%     tas   - total airspeed (knots)
-%     alpha - angle of attack (deg)
-%     beta  - sideslip (deg)
-%     p     - body rate (deg/sec)
-%     q     - body rate (deg/sec)
-%     r     - body rate (deg/sec)
-%     lat   - latitude (deg)
-%     lon   - longitude (deg)
-%     alt   - altitude (ft)
-%     phi   - Euler angle (deg)
-%     theta - Euler angle (deg)
-%     psi   - Euler angle (deg)
-%
-% Values not explictly specified remain at the value they had in MWS
-%
-% Examples: 
-%     MWS=seteomic(init_design(),'alt',3000);
-%     MWS=seteomic(MWS,'tas',90,'beta',3,'p',15);  
-%
-
-% d.e.cox@larc.nasa.gov
-% $Id: seteomic.m 4852 2013-08-06 22:12:54Z cox $
-
-% Constants
-fps2knots = 0.592487;
-knots2fps = 1/fps2knots;
-d2r=pi/180;
-r2d=180/pi;
-
-% default to incoming parameters
-MWS_out=MWS;
-
-% put named IC values into a structure
-ini=struct(varargin{:});
-
-
-% Only set these if one of tas/alpha/beta is being specified
-fn=fieldnames(ini);
-if ~isempty(strmatch('tas',fn)) || ~isempty(strmatch('alpha',fn))|| ~isempty(strmatch('beta',fn)),
-  V=MWS.StatesInp(1:3);
-  Vt=sqrt(sum(V.^2));
-  % For unspecified but coupled parameters, preserve incoming condition.
-  if ~isfield(ini,'alpha'), ini.alpha=atan2(V(3),V(1))*r2d; end
-  if ~isfield(ini,'beta');  ini.beta =atan2(V(2)*cos(ini.alpha),V(1))*r2d;  end
-  if ~isfield(ini,'tas'),   ini.tas=Vt/knots2fps; end
-end
-
-% Error Check input, kick out on unknown names
-fn=fieldnames(ini); 
-knownames={ 'alpha','beta','tas','p','q','r','alt','lat','lon','phi','theta','psi'};
-for i=[1:length(fn)],
-  if isempty(strmatch(fn{i},knownames,'exact')),
-    error(sprintf('Unknown name: %s\n',fn{i}));
-  end
-end
-
-for i=[1:length(fn)],
-  switch(fn{i})
-    case {'alpha','beta','tas'},  % coupled params, this sets IC multiple times,... inefficient but harmless
-      ub = (knots2fps)*ini.tas*cos(ini.alpha*d2r)*cos(ini.beta*d2r);
-      vb = (knots2fps)*ini.tas*sin(ini.beta*d2r);
-      wb = (knots2fps)*ini.tas*sin(ini.alpha*d2r)*cos(ini.beta*d2r);
-      MWS_out.StatesInp(1)  = ub;        %  1 - ub (ft/s)                       
-      MWS_out.StatesInp(2)  = vb;	  %  2 - vb (ft/s)                  
-      MWS_out.StatesInp(3)  = wb;	  %  3 - wb (ft/s)                  
-      
-    case 'p'
-      MWS_out.StatesInp(4)=ini.p*d2r;
-    case 'q'
-      MWS_out.StatesInp(5)=ini.q*d2r;
-    case 'r'
-      MWS_out.StatesInp(6)=ini.r*d2r;
-      
-    case 'lat'
-      MWS_out.StatesInp(7)=ini.lat*d2r;
-    case 'lon'
-      MWS_out.StatesInp(8)=ini.lon*d2r;
-    case 'alt'
-      MWS_out.StatesInp(9)=ini.alt;
-      
-    case 'phi',
-      MWS_out.StatesInp(10)=ini.phi*d2r;
-    case 'theta',
-      MWS_out.StatesInp(11)=ini.theta*d2r;
-    case 'psi',
-      MWS_out.StatesInp(12)=ini.psi*d2r;
-
-  end
-end
-
-
-
-
-
-
-
-
-
-
-
+function MWS_out = seteomic(MWS,varargin);
+%function MWS_out = seteomic(MWS,statename,value,statename,value,...)
+%
+% This function allows the initial condtions of the gtm_design simulation's 
+% equation of motion (MWS.StatesInp) to be specified by named values. 
+%
+% Specifically statename is one of:
+%     tas   - total airspeed (knots)
+%     alpha - angle of attack (deg)
+%     beta  - sideslip (deg)
+%     p     - body rate (deg/sec)
+%     q     - body rate (deg/sec)
+%     r     - body rate (deg/sec)
+%     lat   - latitude (deg)
+%     lon   - longitude (deg)
+%     alt   - altitude (ft)
+%     phi   - Euler angle (deg)
+%     theta - Euler angle (deg)
+%     psi   - Euler angle (deg)
+%
+% Values not explictly specified remain at the value they had in MWS
+%
+% Examples: 
+%     MWS=seteomic(init_design(),'alt',3000);
+%     MWS=seteomic(MWS,'tas',90,'beta',3,'p',15);  
+%
+
+% d.e.cox@larc.nasa.gov
+% $Id: seteomic.m 4852 2013-08-06 22:12:54Z cox $
+
+% Constants
+fps2knots = 0.592487;
+knots2fps = 1/fps2knots;
+d2r=pi/180;
+r2d=180/pi;
+
+% default to incoming parameters
+MWS_out=MWS;
+
+% put named IC values into a structure
+ini=struct(varargin{:});
+
+
+% Only set these if one of tas/alpha/beta is being specified
+fn=fieldnames(ini);
+if ~isempty(strmatch('tas',fn)) || ~isempty(strmatch('alpha',fn))|| ~isempty(strmatch('beta',fn)),
+  V=MWS.StatesInp(1:3);
+  Vt=sqrt(sum(V.^2));
+  % For unspecified but coupled parameters, preserve incoming condition.
+  if ~isfield(ini,'alpha'), ini.alpha=atan2(V(3),V(1))*r2d; end
+  if ~isfield(ini,'beta');  ini.beta =atan2(V(2)*cos(ini.alpha),V(1))*r2d;  end
+  if ~isfield(ini,'tas'),   ini.tas=Vt/knots2fps; end
+end
+
+% Error Check input, kick out on unknown names
+fn=fieldnames(ini); 
+knownames={ 'alpha','beta','tas','p','q','r','alt','lat','lon','phi','theta','psi'};
+for i=[1:length(fn)],
+  if isempty(strmatch(fn{i},knownames,'exact')),
+    error(sprintf('Unknown name: %s\n',fn{i}));
+  end
+end
+
+for i=[1:length(fn)],
+  switch(fn{i})
+    case {'alpha','beta','tas'},  % coupled params, this sets IC multiple times,... inefficient but harmless
+      ub = (knots2fps)*ini.tas*cos(ini.alpha*d2r)*cos(ini.beta*d2r);
+      vb = (knots2fps)*ini.tas*sin(ini.beta*d2r);
+      wb = (knots2fps)*ini.tas*sin(ini.alpha*d2r)*cos(ini.beta*d2r);
+      MWS_out.StatesInp(1)  = ub;        %  1 - ub (ft/s)                       
+      MWS_out.StatesInp(2)  = vb;	  %  2 - vb (ft/s)                  
+      MWS_out.StatesInp(3)  = wb;	  %  3 - wb (ft/s)                  
+      
+    case 'p'
+      MWS_out.StatesInp(4)=ini.p*d2r;
+    case 'q'
+      MWS_out.StatesInp(5)=ini.q*d2r;
+    case 'r'
+      MWS_out.StatesInp(6)=ini.r*d2r;
+      
+    case 'lat'
+      MWS_out.StatesInp(7)=ini.lat*d2r;
+    case 'lon'
+      MWS_out.StatesInp(8)=ini.lon*d2r;
+    case 'alt'
+      MWS_out.StatesInp(9)=ini.alt;
+      
+    case 'phi',
+      MWS_out.StatesInp(10)=ini.phi*d2r;
+    case 'theta',
+      MWS_out.StatesInp(11)=ini.theta*d2r;
+    case 'psi',
+      MWS_out.StatesInp(12)=ini.psi*d2r;
+
+  end
+end
+
+
+
+
+
+
+
+
+
+
+
diff --git a/gtm_design/simulationData.mat b/gtm_design/simulationData.mat
new file mode 100644
index 0000000..4381fdd
Binary files /dev/null and b/gtm_design/simulationData.mat differ