diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000000000000000000000000000000000000..602d83c24fb3ac8feb821b397d650650183d792c
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,28 @@
+*.fig binary
+*.mat binary
+*.mdl binary
+*.mdlp binary
+*.mexa64 binary
+*.mexw64 binary
+*.mexmaci64 binary
+*.mlapp binary
+*.mldatx binary
+*.mlproj binary
+*.mlx binary
+*.p binary
+*.sfx binary
+*.sldd binary
+*.slreqx binary
+*.slmx binary
+*.sltx binary
+*.slxc binary
+*.slx binary merge=mlAutoMerge
+*.slxp binary
+
+## Other common binary file types
+*.docx binary
+*.exe binary
+*.jpg binary
+*.pdf binary
+*.png binary
+*.xlsx binary
diff --git a/README.md b/README.md
index c32ebf9e38c3900e03241ca8e93c1aaf346be0d8..30b56e743f4bd770cf2b51eeca032e8be72df3be 100644
--- a/README.md
+++ b/README.md
@@ -1,33 +1,3 @@
# ACME
-Accelerometer Models Extended: is a MATLAB/Simulink toolbox designed to model and simulate accelerometers. You can use it as a workbench for your accelerometer model alone, or integrate with XHPS for in-flight simulation.
-This toolbox is intended to simulate a variety of accelerometer types, with different geometries, sensor technology, actuation schemes, and so on. Hence, it is modular and parametric.
-
-## Installation
-
-Download all files from the repository and add them to your Matlab path.
-All blocks are in the library (`core/libAcme.slx`), to use it, drag to your project and then use the appropriate setter function to initialize its parameters.
-The tutorials are in the `simulations` folder.
-
-## Requirements
-Matlab 2021a
-
-XHPS 2022 (Optional)
-
-## Features
-The blocks implemented so far will allow you to model a 1DoF GRACE-like ACC, with capacitive sensing and electrostatic actuation. Furthermore:
-* parameters accessible from script;
-* noise generation from ASD expression;
-* compatibility with XHPS for in-flight simulation;
-
-Please follow the tutorials to check its features.
-Let me know of any bugs, questions and suggestions.
-
-## Proposed Features
-
-* noise source models;
-* a more complete dynamics system including
- + actuators with multiple electrodes;
- + rotational DOFs;
- + linearized DOF couplings;
-* novel optical sensors;
+Accelerometer Models Extended: package with the framework to simulate accelerometers for space-based missions, from GRACE-like to after LISA.
\ No newline at end of file
diff --git a/core/AcmeSim.m b/core/AcmeSim.m
new file mode 100644
index 0000000000000000000000000000000000000000..d0c9b3f1ede49b3af8af70bc5ceef8165f867b8e
--- /dev/null
+++ b/core/AcmeSim.m
@@ -0,0 +1,59 @@
+classdef AcmeSim
+ %ACMESIM Summary of this class goes here
+ % Detailed explanation goes here
+
+ properties
+ mdl
+ stopTime = 1e4;
+ nSteps = 1e4;
+ freqs
+ tStep
+ fSample
+ t
+ enbw
+ eps0 = 8.86e-12; %F/m, vacuum permittivity
+ kb = 1.3806e-23; %boltzmann
+ T = 300; % K
+ lpsdopts
+% lpsd
+ end
+
+ methods
+ function obj = AcmeSim(mdl,varargin)
+ %ACMESIM Construct an instance of this class
+ % Detailed explanation goes here
+ obj.mdl = mdl;
+ if length(varargin) == 2
+ obj.stopTime = varargin{1};
+ obj.nSteps = varargin{2};
+ end
+ obj.tStep = obj.stopTime/obj.nSteps;
+ obj.fSample = 1/obj.tStep;
+ obj.freqs = logspace(-log10(obj.stopTime),log10(obj.fSample/2),100);
+ obj.enbw = 1/obj.stopTime;
+ obj.t = (0:obj.tStep:(obj.stopTime-obj.tStep))';
+
+ obj.lpsdopts = {'fs',obj.fSample,...
+ 'Kdes', 100,...
+ 'olap', 30,...
+ 'Jdes', 1000,...
+ 'Lmin', 1,...
+ 'order', 0,...
+ 'scale','ASD',...
+ 'LTPDAmode',false,...
+ 'win','hanning',...
+ };
+ end
+
+ function initModel(obj)
+ load_system(obj.mdl);
+ end
+
+ function TM = addTestMass(obj,mass,sideLength,gap,geometricFactor,block,varargin)
+ TM = TestMass(mass,sideLength,gap,geometricFactor);
+ TM = TM.setDoubleIntegrator(obj.mdl,block,varargin);
+
+ end
+ end
+end
+
diff --git a/core/component/TIA.slx b/core/component/TIA.slx
index 4e52a08a4a10b9803d1929772f30054c8de4b263..55e759f63e94c7298c765ef45aca20eea4554b52 100644
Binary files a/core/component/TIA.slx and b/core/component/TIA.slx differ
diff --git a/core/component/capBridge.slx b/core/component/capBridge.slx
index c95da3abc3f78e3d12682c9569af3f59b7f5869a..150da5b76f81e580208b1fdb3b639fdc9c3345a5 100644
Binary files a/core/component/capBridge.slx and b/core/component/capBridge.slx differ
diff --git a/core/component/electrode.slx b/core/component/electrode.slx
index 8c31a0345d5f1296cd67d99cd81889e034e05ecb..43deb2eca6605e5960de656c18fcee8b4e2ca39f 100644
Binary files a/core/component/electrode.slx and b/core/component/electrode.slx differ
diff --git a/core/component/planeCapacitor.slx b/core/component/planeCapacitor.slx
index 97f96224a5ca4ddb18d40eb0fb762668f180c5d6..c61f399ccf37c91fab3a6212d5a85340d6f609bc 100644
Binary files a/core/component/planeCapacitor.slx and b/core/component/planeCapacitor.slx differ
diff --git a/core/dynamics/dblIntTransFS.slx b/core/dynamics/dblIntTransFS.slx
index 62b77aea70790f32dc517079734b033d8d70c896..2e050c7244f3fd1b2c1ff98cb5c1728991884fe7 100644
Binary files a/core/dynamics/dblIntTransFS.slx and b/core/dynamics/dblIntTransFS.slx differ
diff --git a/core/dynamics/dblIntTransVS.slx b/core/dynamics/dblIntTransVS.slx
index 82fbc661ad07d9d4e6a58223a32b5c0b8468517a..eb855ed99ca2c3eea0afbc146eb9c916c6012d65 100644
Binary files a/core/dynamics/dblIntTransVS.slx and b/core/dynamics/dblIntTransVS.slx differ
diff --git a/core/dynamics/setDoubleIntegrator.m b/core/dynamics/setDoubleIntegrator.m
index 80ecf30d49999db2f4a6b55384f0a55ae776b722..160b57806368e10746493ca02fbf4eaffcaaaf15 100644
--- a/core/dynamics/setDoubleIntegrator.m
+++ b/core/dynamics/setDoubleIntegrator.m
@@ -1,7 +1,4 @@
-function blockHandle = setDoubleIntegrator(model, block, varargin)
-% ========================================================================
-% DEPRECATED, PLZ USE SEPARATED LINEAR AND ROTATION VERSIONS
-% ========================================================================
+function blockHandle = setDblIntTrans(model, block, varargin)
%SETDOUBLEINTEGRATOR - Sets the initial parameters of a double interator block.
%Sets the initial parameters for one instance of the double integrator SISO,
%linear or rotational, block from Acme library.
@@ -11,13 +8,11 @@ function blockHandle = setDoubleIntegrator(model, block, varargin)
% Inputs:
% model - Name of the Simulink model, string
% block - Name of the target instance of the Double Integrator Block, string
-%
-% varargin - optional
-% m - mass, kg
+% varargin - m - mass, kg
% initial position
% - dq0, initial velocity
%
-% See also: setDoubleIntegrator_6DOF
+% See also:
%
% Author: Arthur Reis
% email: arthur.reis@aei.mpg.de
@@ -26,31 +21,28 @@ function blockHandle = setDoubleIntegrator(model, block, varargin)
blockHandle = Simulink.findBlocks(strcat(model),...
'Name',block);
-%% process inputs
defaultq0 = 0;
defaultdq0 = 0;
-%defaultmode = 'linear';
-%expectedmodes = {'linear', 'rotational'};
+defaultdt = 1;
+%defaultSolver = 'VariableStep';
-% parameter parsing, required parms will throw error w/ defaults;
p = inputParser;
-%addrequired
-%add
-
-addRequired(p, 'm' , @(x) x>0);
-addOptional(p, 'q0' , defaultq0 , @(x) isnumeric(x));
-addOptional(p, 'dq0' , defaultdq0, @(x) isnumeric(x));
-addParameter(p, 'mode', defaultmode);
+addRequired(p, 'm', @(x) (x > 0) && isnumeric(x) && isscalar(x));
+addRequired(p, 'gap', @(x) (x > 0) && isnumeric(x) && isscalar(x));
+addOptional(p, 'q0', defaultq0, @(x) isnumeric(x) && isscalar(x));
+addOptional(p, 'dq0', defaultdq0, @(x) isnumeric(x) && isscalar(x));
+%addParameter(p, 'solver', defaultsolver)
+%addOptional(p, 'dt', defaultdt, @(x) isnumeric(x) && isscalar(x));
parse(p, varargin{:});
+
+set_param(blockHandle,'m', num2str(p.Results.m));
+set_param(blockHandle,'gap', num2str(p.Results.gap));
+set_param(blockHandle,'q0', num2str(p.Results.q0));
+set_param(blockHandle,'dq0',num2str(p.Results.dq0));
+%set_param(blockHandle,'dt_solver',num2str(p.Results.dt));
+%if p.Results.solver
+
-switch p.Results.mode
- case 'linear'
- set_param(blockHandle,'m', num2str(p.Results.m));
- set_param(blockHandle,'q0', num2str(p.Results.q0));
- set_param(blockHandle,'dq0',num2str(p.Results.dq0));
- case 'rotational'
- ...
-end
end
diff --git a/core/libAcme.slx b/core/libAcme.slx
index f2cef80e53e25ac68cd6233f8ce60a11c4302944..f813d1140bd5f8d75a889447dc020e95f8899b67 100644
Binary files a/core/libAcme.slx and b/core/libAcme.slx differ
diff --git a/core/noises/gasBrownianNoise.m b/core/noises/gasBrownianNoise.m
new file mode 100644
index 0000000000000000000000000000000000000000..57f76aeb33087b28b6c490d53ece468b108b1dda
--- /dev/null
+++ b/core/noises/gasBrownianNoise.m
@@ -0,0 +1,69 @@
+function varargout = gasBrownianNoise(varargin)
+%THERMALNOISE uses the 4 kb Beta model for white, thermal noise
+% Considering a body inserted on a infinite large cavity filled with Knudsen
+% gas, the ASD of forces acting on the body due to collisions with the gas can
+% be calculated with:
+%
+% S^1/2 = sqrt(4*kb*T*beta) (white noise)
+%
+% Also includes a semi-empirical correction for finite cavity.
+%
+% Syntax:
+% S_1/2 = gasBrownianNoise(vp, T , s, d[, m0])
+% or
+% [Sx_1/2 Sy_1/2 Sz_1/2] =
+% gasBrownianNoise(vp, T , [sx sy sz], [dx dy dz], m0)
+%
+% Input:
+% p - pressure, Pa
+% T - temperature, K
+% s - length of side, m
+% or [sx, sy, sz] - length of sides, m
+% d - gap size, m
+% or [dx, dy, dz] - gap sizes, m
+% m0 - mass of molecule
+%
+% Output:
+% S^1/2 - thermal noise ASD, constant over freqs, N/sqrtHz
+% or [S_x^1/2, S_y^1/2, S_z^1/2] - thermal noise ASD in the 3 axis
+%
+% See also:
+% Cavalleri, et al. Phys Let A (2010)
+% Cavalleri, et al. PRL (2009)
+%
+% Author: Arthur Reis
+% email: arthur.reis@aei.mpg.de
+% Last revision: 03.Nov.2022;
+
+kb = 1.38e-23; %J K-1
+m0_default = 3e-26; %water molecule weight in kg;
+
+par = inputParser;
+par.addRequired('p');
+par.addRequired('T');
+par.addRequired('s');
+par.addRequired('d');
+par.addOptional('m0',m0_default);
+parse(par,varargin{:});
+p = par.Results.p;
+T = par.Results.T;
+s = par.Results.s;
+d = par.Results.d;
+m0 = par.Results.m0;
+
+parse(par,varargin{:});
+
+if (length(s) == 1) && (length(d) == 1)
+ thermal = sqrt(4*kb*T*gasDampingCoeff(p,T,s,d,m0));
+ varargout = {thermal};
+else
+ ax = s(2)*s(3);
+ ay = s(1)*s(3);
+ az = s(1)*s(2);
+ thermalx = sqrt(4*kb*T*gasDampingCoeff(p,T,s(1),d(1),m0,ax,ay,az));
+ thermaly = sqrt(4*kb*T*gasDampingCoeff(p,T,s(2),d(2),m0,ay,ax,az));
+ thermalz = sqrt(4*kb*T*gasDampingCoeff(p,T,s(3),d(3),m0,az,ax,ay));
+ varargout = {thermalx, thermaly, thermalz};
+end
+end
+
diff --git a/core/noises/gasDampingCoeff.m b/core/noises/gasDampingCoeff.m
new file mode 100644
index 0000000000000000000000000000000000000000..633d76825e60ab538ba8bd41ea439e8f61f29949
--- /dev/null
+++ b/core/noises/gasDampingCoeff.m
@@ -0,0 +1,75 @@
+function beta_coeff = gasDampingCoeff(varargin)
+%GASDAMPINGCOEFF calculates the coefficient for 4 kb Beta model
+% Considering a body inserted on a infinite large cavity filled with Knudsen
+% gas, the ASD of forces acting on the body due to collisions with the gas can
+% be calculated with:
+%
+% S^1/2 = sqrt(4*kb*T*beta)
+%
+% Also includes a semi-empirical correction for finite cavity.
+%
+% Syntax:
+% beta_coeff = gasDampingCoeff(p,T,s,d [, m0, A, A2, A3])
+%
+% Input:
+% p - pressure, Pa
+% T - temperature, K
+% s - length of side, m
+% d - gap size, d
+% m0 - kg, mass of molecule, default water
+% A - area of face (normal), m2, optional, if not provided will assume A = s2;
+% A2 - area of face (perp), m2, optional, if not provided assume A2 = A3 = A;
+% A3 - area of face (perp)
+%
+% Output:
+% beta_coeff - gas damping coefficient for that axis
+%
+% See also:
+% Cavalleri, et al. Phys Let A (2010)
+% Cavalleri, et al. PRL (2009)
+%
+% Author: Arthur Reis
+% email: arthur.reis@aei.mpg.de
+% Last revision: 03.Nov.2022;
+
+
+kb = 1.38e-23; %J K-1
+m0_default = 3e-26; %water molecule weight in kg;
+
+par = inputParser;
+par.addRequired('p');
+par.addRequired('T');
+par.addRequired('s');
+par.addRequired('d');
+par.addOptional('m0',m0_default);
+par.addOptional('A',-1);
+par.addOptional('A2',-1);
+par.addOptional('A3',-1);
+
+parse(par,varargin{:});
+
+p = par.Results.p;
+T = par.Results.T;
+s = par.Results.s;
+d = par.Results.d;
+m0 = par.Results.m0;
+A = par.Results.A;
+A2 = par.Results.A2;
+A3 = par.Results.A3;
+
+if A==-1
+ % areas not provided, assume square
+ A = s^2;
+end
+
+
+if (A2 == -1) && (A3 == -1)
+ % one area provided, assume cube
+ A2 = A;
+ A3 = A;
+end
+beta_inf = p*sqrt(2*m0/(pi*kb*T))*((2+pi/2)*A + A2 + A3);
+beta_coeff = beta_inf/(log(s/d)*(d/s)^2); % correction for finite chamber
+end
+
+
diff --git a/core/noises/outgassing.m b/core/noises/outgassing.m
new file mode 100644
index 0000000000000000000000000000000000000000..7da907d7aa4b89528adbfa294a9e6ffdafb75ea8
--- /dev/null
+++ b/core/noises/outgassing.m
@@ -0,0 +1,7 @@
+function [outputArg1,outputArg2] = outgassing(inputArg1,inputArg2)
+%OUTGASSING Summary of this function goes here
+% Detailed explanation goes here
+outputArg1 = inputArg1;
+outputArg2 = inputArg2;
+end
+
diff --git a/core/noises/radiometricNoise.m b/core/noises/radiometricNoise.m
new file mode 100644
index 0000000000000000000000000000000000000000..8339b9c8bd3fddfb7012861dc1d1e0eb741622d0
--- /dev/null
+++ b/core/noises/radiometricNoise.m
@@ -0,0 +1,21 @@
+function dF = radiometricNoise(A, P, T0, dT, varargin)
+%RADIOMETRICNOISE Summary of this function goes here
+% Detailed explanation goes here
+% INPUT
+% A - area of the faces, m2;
+% P - pressure
+% T0 - temperature of the system
+% dT - difference between average
+% OUTPUT
+
+% Reference: Carbone 2007
+
+parser = inputParser;
+parser.addOptional('kappaR', 1); % correction factor
+parser.parse;
+kappaR = parser.Results.kappaR;
+
+dF = kappaR*A*P*dT/(4*T0);
+
+end
+
diff --git a/core/noises/thermalRadiationNoise.m b/core/noises/thermalRadiationNoise.m
new file mode 100644
index 0000000000000000000000000000000000000000..67e0a85cf4226659ac0ce8ceda1a562521d1cc7f
--- /dev/null
+++ b/core/noises/thermalRadiationNoise.m
@@ -0,0 +1,22 @@
+function dF = thermalRadiationNoise(A, P, T0, dT, varargin)
+%RADIATION Summary of this function goes here
+% Detailed explanation goes here
+% INPUT
+% A - area of the faces, m2;
+% P - pressure
+% T0 - temperature of the system
+% dT - difference between average
+% OUTPUT
+
+% Reference: Carbone 2007
+sig = 5.67e-8;
+c = 299792000;
+parser = inputParser;
+parser.addOptional('kappaRP', 1); % correction factor
+parser.parse;
+kappaRP = parser.Results.kappaR;
+
+dF = kappaRP*A*8*sig/(3*c)*T0^3*dT;
+
+end
+
diff --git a/core/sensors/diffTransformer.slx b/core/sensors/diffTransformer.slx
index 0e93622c1954d149a468e702bed7c4608b40aa6b..a899e43e806db5ccb66cca5bd12ea6f1b31640f2 100644
Binary files a/core/sensors/diffTransformer.slx and b/core/sensors/diffTransformer.slx differ
diff --git a/core/sensors/lockIn.slx b/core/sensors/lockIn.slx
index 36884818b130d9520df787d7a4109d67d01a94ab..a7f9fac19abaf20f3d77bd5be6278ee05d901cc3 100644
Binary files a/core/sensors/lockIn.slx and b/core/sensors/lockIn.slx differ
diff --git a/simulations/Electrostatic/Copy_of_LISA_EA.slx b/simulations/Electrostatic/Copy_of_LISA_EA.slx
new file mode 100644
index 0000000000000000000000000000000000000000..751e59f932cdd8d26704c6557669c5e754f47469
Binary files /dev/null and b/simulations/Electrostatic/Copy_of_LISA_EA.slx differ
diff --git a/simulations/Electrostatic/Copy_of_LISA_EA_script.m b/simulations/Electrostatic/Copy_of_LISA_EA_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..33c8eddf7e6f2342608119d3836cd83b8bd7ade9
--- /dev/null
+++ b/simulations/Electrostatic/Copy_of_LISA_EA_script.m
@@ -0,0 +1,102 @@
+%% Electrostatic ACC Model
+% Following the description in Davor Mance's 2012 thesis. References are from
+% that, unless otherwise mentioned. The capacitance sensor model
+% Author: Arthur Reis
+% email: arthur.reis@aei.mpg.de
+% last review: 19.10.2022
+clear;
+
+mdl = 'LISA_EA';
+load_system(mdl);
+
+%% Simulation parameters
+% frequency and time
+stopTime = 1e4;
+nSteps = 1e6;
+tStep = stopTime/nSteps;
+fSample = 1/(stopTime/nSteps);
+t = (0:tStep:(stopTime-tStep))';
+enbw = 1/stopTime;
+freq = logspace(-4,0,100);
+% Physical constants
+eps0 = 8.86e-12; %F/m, vacuum permittivity
+kb = 1.3806e-23; %boltzmann
+T = 300; % K
+
+
+%% Accelerometer Parameters (Ch.2)
+% TM
+l_TM = 0.046; %m,
+mass = 2; %kg
+gap0 = 0.004; % m
+area = l_TM^2; %m2
+
+dblint = setDblIntTrans(mdl, 'dblInt', mass, gap0);
+
+% Electrodes
+% There are four electrodes, 2 on the faces on the cage, at potential -V and +V,
+% and 2 faces of the TM, at same potential
+% setElectrode(model, block, area, pos0, V_DC, [V_AC], [f], [solver])
+Vd = 5; %V detection rms voltage
+Vp = 1; %V polarizing voltage
+fCar = 1e5;
+el1 = setElectrode(mdl,'electrode1',area,-(l_TM/2+gap0),0 ,0 ,fCar); %Cage left
+el2 = setElectrode(mdl,'electrode2',area,-l_TM/2 ,Vp,Vd,fCar); %TM left
+el3 = setElectrode(mdl,'electrode3',area,+l_TM/2 ,Vp,Vd,fCar); %TM right
+el4 = setElectrode(mdl,'electrode4',area,+(l_TM/2+gap0),0 ,0 ,fCar); %Cage right
+
+% Plane capacitors, each formed by two electrodes
+% setPlaneCapacitor(model, block, eps0)
+cap1 = setPlaneCapacitor(mdl, 'planeCapacitor1', eps0);
+cap2 = setPlaneCapacitor(mdl, 'planeCapacitor2', eps0);
+capSensorGain = 1.8e13;
+capSensorNoiseASD = @(f) 5.0124e-19 + 1.852e-22./f;
+capSensorNoiseTS = [t timeseriesFromASD(capSensorNoiseASD, enbw, fSample, nSteps).Data];
+% Capacitive sensing
+nominal_cap = 1.15e-12;% F
+Cmax = 0.12e-12; %Tab. B-1
+fTypical = 1e-2; %Hz
+
+% Control
+Kp = 1;
+Ki = 0;
+Kd = 0;
+driverGain = 1;
+
+%% Simulating
+mdlOutput = sim(mdl);
+
+%% Analysing the spectrum
+
+[axx,fxx] = flpsd('ts', mdlOutput.simout.Data, 'fs', fSample, 'Kdes', 100, ...
+ 'Lmin', 1, 'Jdes', 1000, 'win', 'hanning', 'olap', 30, ...
+ 'order', 0, 'scale', 'ASD', 'LTPDAmode', false);
+
+% Comparing the results
+% Capacitance sensor
+io(1) = linio(strcat(mdl,'/Minus1'),1,'input');
+io(2) = linio(strcat(mdl,'/simplifiedCapSensor'),1,'output');
+readout2capTF = getTFresponse(mdl, io, fxx); % TF between readout and acc signal
+figure(1);
+loglog(fxx,axx./readout2capTF,fxx,capSensorNoiseASD(fxx),'linewidth',2);
+legend('Simulated measurement','Simulated white noise input');
+xlabel('Frequency [Hz]');
+ylabel('Capacitance [F]');
+
+% Accelerometer readout
+io(1) = linio(strcat(mdl,'/Gain'),1,'input');
+io(2) = linio(strcat(mdl,'/simplifiedCapSensor'),1,'output');
+readout2accTF = getTFresponse(mdl, io, fxx); % TF between readout and acc signal
+figure(2);
+loglog(fxx,axx./readout2accTF);
+legend('Simulated measurement','Simulated white noise input');
+xlabel('Frequency [Hz]');
+ylabel('Acceleration [m/s2]');
+
+%%
+dof = 'x';
+[noises, sys] = nbFromSimulink(mdl, freq, 'dof', dof);
+nb = nbGroupNoises(mdl, noises, sys);
+matlabNoisePlot(nb);
+%loglog(fxx, axx./readout2accTF)
+
diff --git a/simulations/LISA_EA/GRS.slx b/simulations/LISA_EA/GRS.slx
new file mode 100644
index 0000000000000000000000000000000000000000..66f8b62f6fdd6fb35bb9f53dbffb4409f08ca962
Binary files /dev/null and b/simulations/LISA_EA/GRS.slx differ
diff --git a/simulations/LISA_EA/GRS_script.m b/simulations/LISA_EA/GRS_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..6807af582785a16655a36a4b77dc93d95ffbfb3d
--- /dev/null
+++ b/simulations/LISA_EA/GRS_script.m
@@ -0,0 +1,123 @@
+%% Electrostatic ACC Model
+% Following the description in Davor Mance's 2012 thesis. References are from
+% that, unless otherwise mentioned. The capacitance sensor model
+% Author: Arthur Reis
+% email: arthur.reis@aei.mpg.de
+% last review: 19.10.2022
+clear;
+
+mdl = 'LISA_EA';
+load_system(mdl);
+
+%% Simulation parameters
+% frequency and time
+stopTime = 1e4;
+nSteps = 1e4;
+tStep = stopTime/nSteps;
+fSample = 1/(stopTime/nSteps);
+t = (0:tStep:(stopTime-tStep))';
+enbw = 1/stopTime;
+freq = logspace(-4,0,100);
+% Physical constants
+eps0 = 8.86e-12; %F/m, vacuum permittivity
+kb = 1.3806e-23; %boltzmann
+T = 300; % K
+
+
+%% Accelerometer Parameters (Ch.2)
+% TM
+l_TM = 0.030; %m,
+mass = .54; %kg
+gap0 = 0.001; % m
+area = l_TM^2; %m2
+p = 1e-5; %Pa
+dblint = setDblIntTrans(mdl, 'dblInt', mass, gap0);
+
+% Electrodes
+% There are four electrodes, 2 on the faces on the cage, at potential -V and +V,
+% and 2 faces of the TM, at same potential
+% setElectrode(model, block, area, pos0, V_DC, [V_AC], [f], [solver])
+% Vd = 5; %V detection rms voltage
+% Vp = 1; %V polarizing voltage
+% fCar = 1e5;
+% Vxbias = 16;
+% el1 = setElectrode(mdl,'electrode1',area,-(l_TM/2+gap0),0 ,0 ,fCar); %Cage left
+% el2 = setElectrode(mdl,'electrode2',area,-l_TM/2 ,Vp,Vd,fCar); %TM left
+% el3 = setElectrode(mdl,'electrode3',area,+l_TM/2 ,Vp,Vd,fCar); %TM right
+% el4 = setElectrode(mdl,'electrode4',area,+(l_TM/2+gap0),0 ,0 ,fCar); %Cage right
+%
+Vd = 5; %V detection rms voltage "The TM is modulated at f0 = 100 kHz using the
+% injection electrodes, which induce an oscillating voltage in the TM, VTM ,
+% of about one tenth of the injected voltage" "VTM = 0.5 V"
+V_TM = 0.5;
+Vp = 0; %V polarizing voltage #TODO
+Vx0 = 16; %V rm
+fCar = 1e5;
+el1 = setElectrode(mdl,'electrode1',area,-(l_TM/2+gap0), 0, 0 ,fCar); %Cage left
+el2 = setElectrode(mdl,'electrode2',area,-l_TM/2 , 0, V_TM,fCar); %TM left
+el3 = setElectrode(mdl,'electrode3',area,+l_TM/2 , 0, V_TM,fCar); %TM right
+el4 = setElectrode(mdl,'electrode4',area,+(l_TM/2+gap0), 0, 0 ,fCar); %Cage right
+
+% Plane capacitors, each formed by two electrodes
+% setPlaneCapacitor(model, block, eps0)
+cap1 = setPlaneCapacitor(mdl, 'planeCapacitor1', eps0);
+cap2 = setPlaneCapacitor(mdl, 'planeCapacitor2', eps0);
+%capSensorGain = 1.8e13;
+Cf = 10e-12; %F
+capSensorGain = 10*4*V_TM/Cf; %IA, BPF, cap sensor
+%capSensorNoiseASD = @(f) 5.0124e-19 + 1.852e-22./f;
+%capSensorNoiseASD = @(f) 5.0124e-19*f.^0;
+capSensorNoiseASD = @(f) 5.0124e-21*f.^0;
+capSensorNoiseTS = [t timeseriesFromASD(capSensorNoiseASD, enbw, fSample, nSteps).Data];
+% Capacitive sensing
+nominal_cap = 1.15e-12;% F
+Cmax = 0.12e-12; %Tab. B-1
+fTypical = 1e-2; %Hz
+% Control
+Kp = 1;
+Ki = 0.0;
+Kd = 10;
+driverGain =-1;
+
+gasThermalNoiseASD = @(f) gasBrownianNoise(p,T,l_TM,gap0)/mass + 1.6e-15./sqrt(f);
+gasThermalNoiseTS = [t timeseriesFromASD(gasThermalNoiseASD, enbw, fSample, nSteps).Data];
+%%
+%bbb = [all_non_grav_forces_A_x(1:50000); all_non_grav_forces_A_x(50000:-1:1)];
+%aaa = [t bbb];
+%% Simulating
+mdlOutput = sim(mdl);
+
+%% Analysing the spectrum
+
+[axx,fxx] = flpsd('ts', mdlOutput.simout.Data, 'fs', fSample, 'Kdes', 100, ...
+ 'Lmin', 1, 'Jdes', 1000, 'win', 'hanning', 'olap', 30, ...
+ 'order', 0, 'scale', 'ASD', 'LTPDAmode', false);
+
+% Comparing the results
+% Capacitance sensor
+io(1) = linio(strcat(mdl,'/Minus1'),1,'input');
+io(2) = linio(strcat(mdl,'/simplifiedCapSensor'),1,'output');
+readout2capTF = getTFresponse(mdl, io, fxx); % TF between readout and acc signal
+% figure(1);
+% loglog(fxx,axx./readout2capTF,fxx,capSensorNoiseASD(fxx),'linewidth',2);
+% legend('Simulated measurement','Simulated white noise input');
+% xlabel('Frequency [Hz]');
+% ylabel('Capacitance [F]');
+%
+% % Accelerometer readout
+io(1) = linio(strcat(mdl,'/Gain'),1,'input');
+io(2) = linio(strcat(mdl,'/simplifiedCapSensor'),1,'output');
+readout2accTF = getTFresponse(mdl, io, fxx); % TF between readout and acc signal
+% figure(2);
+% loglog(fxx,axx./readout2accTF);
+% legend('Simulated measurement','Simulated cap noise input');
+% xlabel('Frequency [Hz]');
+% ylabel('Acceleration [m/s2]');
+
+%%
+dof = 'x';
+[noises, sys] = nbFromSimulink(mdl, freq, 'dof', dof);
+nb = nbGroupNoises(mdl, noises, sys);
+matlabNoisePlot(nb);
+loglog(fxx, axx./readout2accTF)
+
diff --git a/simulations/LISA_EA/capSensor_complete.slx b/simulations/LISA_EA/capSensor_complete.slx
new file mode 100644
index 0000000000000000000000000000000000000000..7c5f95350042260ac8fa2d39b6c53e7a1242c843
Binary files /dev/null and b/simulations/LISA_EA/capSensor_complete.slx differ
diff --git a/simulations/LISA_EA/capSensor_complete_script.m b/simulations/LISA_EA/capSensor_complete_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..67fc0bc6771127feb5b1dbcdf25053684042f785
--- /dev/null
+++ b/simulations/LISA_EA/capSensor_complete_script.m
@@ -0,0 +1,136 @@
+%% ACME - Capacitive Sensor Simulator
+% All references are from Mance's 2012 thesis, unless otherwise noted.
+% This model is calculated in the frequency range 1e-4 - 1Hz. Effects of the
+% 100kHz injection signal are accounted for but the signal is not fully
+% simulated.
+% Arthur Reis, 29.11.2022
+clear;
+
+mdl = 'capSensor_complete';
+load_system(mdl);
+%% Simulation parameters
+% frequency and time
+stopTime = 1e4;
+nSteps = 1e6;
+tStep = stopTime/nSteps;
+fSample = 1/(stopTime/nSteps);
+t = (0:tStep:(stopTime-tStep))';
+enbw = 1/stopTime;
+freq = logspace(-4,0,100);
+% Physical constants
+eps0 = 8.86e-12; %F/m, vacuum permittivity
+kb = 1.3806e-23; %boltzmann
+T = 300; % K
+
+%% Accelerometer Parameters
+% TM, Ch.2
+TM_side = 0.046; %m,
+mass = 2; %kg
+nominal_gap = 0.004; % m
+area = TM_side^2; %m2
+
+% Capacitive sensing
+nominal_cap = 1.15e-12;% F
+Cmax = 0.12e-12; %Tab. B-1
+fTypical = 1e-2; %Hz
+
+%% Voltage Reference
+Vd = 4.8; % Vpp, Ch. 2.6
+fDetection = 1e5; %Hz
+capAttenuation = 0.123; %Eq 2.60, from ratio between injection
+% and all other capacitances
+VRefUNoiseASD = @(f) (125e-9.*f.^0) + 1.5e-7.*f.^-.5; %V/sqrt(Hz)
+%Ch. 2.5, LT1021 datasheet
+VRefUNoiseTS = [t timeseriesFromASD(VRefUNoiseASD, enbw, fSample, nSteps).Data];
+% DAC Noise - AD660
+DACNoiseASD = @(f) (0.125e-6.*f.^0) + 0.800e-6.*f.^-.5; %V/sqrt(Hz)
+% AD660 20V span, int ref
+DACNoiseTS = [t timeseriesFromASD(DACNoiseASD, enbw, fSample, nSteps).Data];
+
+%% Capacitance Bridge
+% Tab. B-1
+L = 4.2e-3; %H
+Cp = 320e-12; %F
+C0 = 1.15e-12; %Mance tab B1 ???
+Q = 200; %Mance tab B1
+
+[brBlk, Zbr, f0br] = setCapBridge(mdl, 'capBridge', L, Cp, C0, Q);
+% Bridge Impedance Thermal Noise - Mance 2012 eq 2.24
+brRThNoise = sqrt(4*kb*T*real(Zbr(f0br)));
+brRThNoiseASD = @(f) brRThNoise.*f.^0;
+brRThNoiseTS = [t timeseriesFromASD(brRThNoiseASD, enbw, fSample, nSteps).Data];
+
+%% Transimpedance Amplifier
+% Tab. B-1
+tiaCfb = 3.3e-12; %F
+tiaRfb = 10e6; %ohm
+%brR = real(brZ(fDetection));
+[tiaBlk, tiaZfb,tiaGs] = setTIA(mdl, 'TIAmp', Zbr(f0br), tiaCfb, tiaRfb);
+%tiaR = real(tiaZfb(brf0));
+
+% TIA Feedback Impedance Thermal Noise - Mance 2012
+tiaRfbThNoise = 2*sqrt(4*kb*T*real(tiaZfb(f0br))); %V,sqrt(hz) eq 2.94
+tiaRfbThNoiseASD = @(f) tiaRfbThNoise.*f.^0;
+tiaRfbThNoiseTS = [t timeseriesFromASD(tiaRfbThNoiseASD, enbw, fSample, nSteps).Data];
+
+% TIA Current Noise - eq 2.95
+tiaINoise = 2*1.6e-15*real(tiaZfb(f0br)); % tab B1
+tiaINoiseASD = @(f) tiaINoise.*f.^0;
+tiaINoiseTS = [t timeseriesFromASD(tiaINoiseASD, enbw, fSample, nSteps).Data];
+
+% TIA Voltage Noise - eq 2.95
+NGs = 1 + 2*abs(tiaZfb(f0br)/Zbr(f0br));
+tiaUNoise = 4.5e-9*NGs*sqrt(2); % tab B1
+tiaUNoiseASD = @(f) tiaUNoise.*f.^0;
+tiaUNoiseTS = [t timeseriesFromASD(tiaUNoiseASD, enbw, fSample, nSteps).Data];
+
+
+%% AC amplifier
+acGain = 45.4; % Ch. 2.8.5
+% amplifier voltage noise
+acUNoise = 5.2e-9; %V/sqrt(Hz) tab 2.15
+acUNoiseASD = @(f) acUNoise.*f.^0;
+acUNoiseTS = [t timeseriesFromASD(acUNoiseASD, enbw, fSample, nSteps).Data];
+
+%% Demodulator
+% DC amp voltage noise
+dcampgain = 1.2; %calibrate so max +- 2.5V
+DCampUNoiseASD = @(f) (0.03e-6.*f.^0) ; %V/sqrt(Hz) AD8629, mance
+%DCampUNoiseASD = @(f) (14e-9.*f.^0) + 6.5e-9.*f.^-.7; %V/sqrt(Hz) OP97 (arbitrary)
+DCampUNoiseTS = [t timeseriesFromASD(DCampUNoiseASD, enbw, fSample, nSteps).Data];
+
+%% ADC
+% quantization noise
+ADCqNoiseASD = @(f) (2.4e-6.*f.^0) ; %V/sqrt(Hz) from eq 2.104
+ADCqNoiseTS = [t timeseriesFromASD(ADCqNoiseASD, enbw, fSample, nSteps).Data];
+
+%% Noise Budget
+% this uses the SimulinkNB toolbox.
+% it doesn't run the simulation. Instead, it finds the transfer functions of the
+% linearized model and plots the noise budged
+[noises, sys] = nbFromSimulink(mdl, freq, 'dof', 'x');
+nb = nbGroupNoises(mdl, noises, sys);
+matlabNoisePlot(nb);
+
+%% Time domain simulation
+mdlOutput = sim(mdl);
+% lpsd
+Kdes = 100;
+olap = 30;
+Jdes = 1000;
+Lmin = 1;
+order = 0;
+
+[axx,fxx] = flpsd( 'ts', mdlOutput.simout.Data, 'fs', fSample, 'Kdes', Kdes, 'Lmin', Lmin, ...
+ 'Jdes', Jdes, 'win', 'hanning', 'olap', olap, 'order', order,...
+ 'scale', 'ASD', 'LTPDAmode',false);
+
+%%
+% finding the TF between displacement readout and input perturbation
+io(1) = linio(strcat(mdl,'/C1'),1,'input');
+io(2) = linio(strcat(mdl,'/Sum13'),1,'output');
+readout2capTF = getTFresponse(mdl, io, fxx); % TF between readout and acc signal
+
+loglog(fxx, axx./readout2capTF,'k','LineWidth',1);
+xlabel('Frequency [Hz]');
+ylabel('Measured capacitance [F]');
\ No newline at end of file
diff --git a/simulations/LISA_EA/capSensor_simplified.slx b/simulations/LISA_EA/capSensor_simplified.slx
new file mode 100644
index 0000000000000000000000000000000000000000..20f8439da5a8a59d40cbe687715761ff078a1abb
Binary files /dev/null and b/simulations/LISA_EA/capSensor_simplified.slx differ
diff --git a/simulations/LISA_EA/capSensor_simplified_script.m b/simulations/LISA_EA/capSensor_simplified_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..44d99828077020a850ab71a6353ba1f74139cb72
--- /dev/null
+++ b/simulations/LISA_EA/capSensor_simplified_script.m
@@ -0,0 +1,117 @@
+%% ACME - Simplified Capacitive Sensor Simulator
+% All references are from Mance's 2012 thesis, unless otherwise noted.
+% This model is calculated in the frequency range 1e-4 - 1Hz. Effects of the
+% 100kHz injection signal are accounted for but the signal is not fully
+% simulated.
+% Only the dominant noise sources from the complete model are included
+% Arthur Reis, 29.11.2022
+clear;
+
+mdl = 'capSensor_simplified';
+load_system(mdl);
+%% Simulation parameters
+% frequency and time
+stopTime = 1e4;
+nSteps = 1e6;
+tStep = stopTime/nSteps;
+fSample = 1/(stopTime/nSteps);
+t = (0:tStep:(stopTime-tStep))';
+enbw = 1/stopTime;
+freq = logspace(-4,0,100);
+% Physical constants
+eps0 = 8.86e-12; %F/m, vacuum permittivity
+kb = 1.3806e-23; %boltzmann
+T = 300; % K
+
+%% Accelerometer Parameters
+% TM, Ch.2
+TM_side = 0.046; %m,
+mass = 2; %kg
+nominal_gap = 0.004; % m
+area = TM_side^2; %m2
+
+% Capacitive sensing
+nominal_cap = 1.15e-12;% F
+Cmax = 0.12e-12; %Tab. B-1
+fTypical = 1e-2; %Hz
+
+%% Voltage Reference
+Vd = 4.8; % Vpp, Ch. 2.6
+fDetection = 1e5; %Hz
+capAttenuation = 0.123; %Eq 2.60, from ratio between injection
+% and all other capacitances
+% DAC Noise - AD660
+DACNoiseASD = @(f) (0.125e-6.*f.^0) + 0.800e-6.*f.^-.5; %V/sqrt(Hz)
+% AD660 20V span, int ref
+DACNoiseTS = [t timeseriesFromASD(DACNoiseASD, enbw, fSample, nSteps).Data];
+
+%% Capacitance Bridge
+% Tab. B-1
+L = 4.2e-3; %H
+Cp = 320e-12; %F
+C0 = 1.15e-12; %Mance tab B1 ???
+Q = 200; %Mance tab B1
+
+[brBlk, Zbr, f0br] = setCapBridge(mdl, 'capBridge', L, Cp, C0, Q);
+% Bridge Impedance Thermal Noise - Mance 2012 eq 2.24
+brRThNoise = sqrt(4*kb*T*real(Zbr(f0br)));
+brRThNoiseASD = @(f) brRThNoise.*f.^0;
+brRThNoiseTS = [t timeseriesFromASD(brRThNoiseASD, enbw, fSample, nSteps).Data];
+
+%% Amplifier
+% arbitrary gain for the simplified model, adjust if needed
+% output should be max around +-2.5V
+ampgain = 100;
+%% Noise Budget
+% this uses the SimulinkNB toolbox.
+% it doesn't run the simulation. Instead, it finds the transfer functions of the
+% linearized model and plots the noise budged
+[noises, sys] = nbFromSimulink(mdl, freq, 'dof', 'x');
+nb = nbGroupNoises(mdl, noises, sys);
+
+%% Time domain simulation
+mdlOutput = sim(mdl);
+% lpsd options
+Kdes = 100;
+olap = 30;
+Jdes = 1000;
+Lmin = 1;
+order = 0;
+% lpsd
+[axx,fxx] = flpsd( 'ts', mdlOutput.simout.Data, 'fs', fSample, 'Kdes', Kdes, 'Lmin', Lmin, ...
+ 'Jdes', Jdes, 'win', 'hanning', 'olap', olap, 'order', order,...
+ 'scale', 'ASD', 'LTPDAmode',false);
+%% Finding the transfer functions
+% finding the transfer function between recorded signal and the signal we wish
+% to measure
+io(1) = linio(strcat(mdl,'/C1'),1,'input');
+io(2) = linio(strcat(mdl,'/Amp'),1,'output');
+readout2capTF = getTFresponse(mdl, io, fxx); % TF between readout and Cap1
+sensorTFexp = fit(fxx,readout2capTF, 'poly2')
+
+%% Expression for the total noise ASD
+% assuming ASD = A/f + B, where B is the high freq limit
+x = nb.sumNoise.f';
+y = nb.sumNoise.asd';
+sensorASDexp = fit(x,y, fittype(str2func(['@(a, x) a./x + ' num2str(y(end))])))
+
+%from time domain simulation
+xsim = fxx;
+ysim = axx./readout2capTF;
+sensorASDexpsim = fit(xsim,ysim, fittype(...
+ str2func(['@(a, x) a./x + ' num2str(ysim(end))])))
+
+%% Plots
+matlabNoisePlot(nb);
+xlabel('Frequency [Hz]');
+ylabel('Measured capacitance [F]');
+loglog(freq, sensorASDexp(freq),'b.', 'displayName', 'Fit from lin model');
+loglog(freq, sensorASDexpsim(freq),'k.','displayName', 'Fit from sim');
+loglog(fxx, axx./readout2capTF,'k-','LineWidth',1,'DisplayName', ...
+ 'ASD from time domain simulation');
+
+%% Finding the transfer functions
+% lpsd
+[axxs,fxxs] = flpsd( 'ts', mdlOutput.simout1.Data, 'fs', fSample, 'Kdes', Kdes, 'Lmin', Lmin, ...
+ 'Jdes', Jdes, 'win', 'hanning', 'olap', olap, 'order', order,...
+ 'scale', 'ASD', 'LTPDAmode',false);
\ No newline at end of file
diff --git a/simulations/SGRS/capnoise2acctf.m b/simulations/SGRS/capnoise2acctf.m
new file mode 100644
index 0000000000000000000000000000000000000000..9489354b7ade79bf3cb3ad9d202a4b187556dc55
--- /dev/null
+++ b/simulations/SGRS/capnoise2acctf.m
@@ -0,0 +1,36 @@
+G = 40
+Vtm = .5;
+m = 0.54;
+A= (0.021)^2;
+eps0 = 8.89e-12;
+Cf = 10e-12;
+d = 0.001;
+C0 = 3.9e-12
+
+ntofV = tf(51e-9*[-m*Cf*d^3 0 2*eps0*A*Cf*Vtm^2]./m,[2*eps0*A*d*Vtm]); %acc noise from the 51nV/sqrtHz contribution
+freqs = logspace(-5,0,100);
+ntofVfreq = frd(ntofV, freqs, 'unit', 'Hz');
+ntofVfreq = squeeze(abs(ntofVfreq.ResponseData));
+loglog(freqs,ntofVfreq,'.k','DisplayName','charge amplifier voltage noise')
+thermalasd_fromdavila = @(f) 4e-13 *sqrt(1 + 700e-6./f + (300e-6./f).^2);
+hold on
+loglog(freqs, thermalasd_fromdavila(freqs),'--','DisplayName','Thermal+electric Noise from Davila Alvarez')
+legend()
+xlabel('Frequency [Hz]');
+xlim([1e-4 1e-1]);
+ylim([1e-15 1e-8]);
+ylabel('acceleration ASD [m/s2 sqrt(hz)]');
+title('contribution of capacitive noise')
+
+2*G*Vtm*C0/(Cf*d)
+
+(2*eps0*A*Vtm)/(Cf*d^2)
+
+%%
+x = tf(1.3e-10*[1 0 -1e-6/m],[1]); %acc noise from the 51nV/sqrtHz contribution
+freqs = logspace(-5,0,100);
+xtoacc = frd(x, freqs, 'unit', 'Hz');
+xtoacc = squeeze(abs(xtoacc.ResponseData));
+loglog(freqs,xtoacc,'DisplayName','charge amplifier voltage noise')
+
+k = 2*eps0*A*Vtm^2/d^3
\ No newline at end of file
diff --git a/simulations/SGRS/getSGRSCapSensorASD.m b/simulations/SGRS/getSGRSCapSensorASD.m
new file mode 100644
index 0000000000000000000000000000000000000000..69841c0e8110cc76440068ed3ed5eb526ce412aa
--- /dev/null
+++ b/simulations/SGRS/getSGRSCapSensorASD.m
@@ -0,0 +1,45 @@
+function sensorASDexp = getSGRSCapSensorASD(acme, TM, C0, Cf, IAGain,chargeAmpNoiseLevel)
+%GETSGRSCAPSENSORASD Summary of this function goes here
+% Detailed explanation goes here
+
+
+
+ C0% F, capacitance between tm and elec. when dx = 0
+%fTypical = 1e-2; %Hz
+Cf = 10e-12; %F feedback capacitor
+%% Voltage Reference
+Vd = 5; % Vpp
+fDetection = 1e5; %Hz
+capAttenuation = 0.1;
+
+
+%% Amplifier
+
+chargeAmpNoiseASD = @(f) chargeAmpNoiseLevel.*f.^0; %V/sqrtHz
+assignin('base','chargeAmpNoiseASD',chargeAmpNoiseASD);
+assignin('base','IA_BPF_gain',IAGain);
+assignin('base','nominal_cap',C0);
+assignin('base','Cf',Cf);
+assignin('base','fDetection',fDetection);
+assignin('base','capAttenuation',capAttenuation);
+assignin('base','Vd',Vd);
+
+
+
+%% Noise Budget
+% this uses the SimulinkNB toolbox.
+% it doesn't run the simulation. Instead, it finds the transfer functions of the
+% linearized model and plots the noise budged
+load_system('sgrs_capSensor')
+[noises, sys] = nbFromSimulink('sgrs_capSensor', acme.freqs, 'dof', 'x');
+nb = nbGroupNoises('sgrs_capSensor', noises, sys);
+
+
+%% Expression for the total noise ASD
+% assuming ASD = A/f + B, where B is the high freq limit
+x = nb.sumNoise.f';
+y = nb.sumNoise.asd';
+sensorASDexp = fit(x,y, fittype(str2func(['@(a, x) a./x + ' num2str(y(end))])));
+clear('chargeAmpNoiseASD','IA_BPF_gain','nominal_cap');
+end
+
diff --git a/simulations/SGRS/sgrs_acc.slx b/simulations/SGRS/sgrs_acc.slx
new file mode 100644
index 0000000000000000000000000000000000000000..b8cac3d00643435a828973431f288cc4adb7880c
Binary files /dev/null and b/simulations/SGRS/sgrs_acc.slx differ
diff --git a/simulations/SGRS/sgrs_acc_parametrization.m b/simulations/SGRS/sgrs_acc_parametrization.m
new file mode 100644
index 0000000000000000000000000000000000000000..a9893209c549f4ac0408a8b353befc6cfa0091d8
--- /dev/null
+++ b/simulations/SGRS/sgrs_acc_parametrization.m
@@ -0,0 +1,101 @@
+% Electrostatic ACC Model
+% Following the Simplified Gravitational Reference Sensor model
+% The S-GRS is a scaled-down version of the LPF GRS
+% from: arXiv:2107.08545v1
+% Author: Arthur Reis
+% email: arthur.reis@aei.mpg.de
+% last review: 29.11.2022
+clear;
+acme = AcmeSim('sgrs_acc');
+acme.initModel;
+
+% If is non-drag-free:
+driverGain = -20;
+V_TM = 0.5;
+Vp = 0; %V polarizing
+Vx0 = 16; %V rms
+
+% If is drag-free:
+% driverGain = 0;
+% V_TM = 0.5;
+% Vp = 0; %V polarizing
+% Vx0 = 0; %V rms
+
+tmgap = 1e-3*[1, 2, 4];
+TM = acme.addTestMass(0.54,30e-3,tmgap(1),1,'dblInt',0.07);
+
+fCar = 1e5;
+
+% Plane capacitors, each formed by two electrodes
+cap1 = setPlaneCapacitor(acme.mdl, 'planeCapacitor1', acme.eps0);
+cap2 = setPlaneCapacitor(acme.mdl, 'planeCapacitor2', acme.eps0);
+%CapNoiseASD = @(f) 1.854e-22 ./f + 5.0124e-19; %F/sqrt(hz)
+CapNoiseASD = @(f) 51e-9.*f.^0; %F/sqrt(hz)
+CapNoiseTS = [acme.t timeseriesFromASD(CapNoiseASD, acme.enbw, acme.fSample, acme.nSteps).Data];
+Cf = 10e-12; %F
+%Cf = 10e-12; %F
+%getSGRSCapSensorASD(acme, TM,1e-12,Cf,40,51e-9)
+ActuationNoiseASD = @(f) 20e-6.*f.^0; %V/sqrt(hz)
+ActuationNoiseTS = [acme.t timeseriesFromASD(ActuationNoiseASD, acme.enbw, acme.fSample, acme.nSteps).Data];
+io(2) = linio(strcat(acme.mdl,'/BPF'),1,'output');
+io(1) = linio(strcat(acme.mdl,'/Gain1'),1,'input');
+dof = 'x';
+% Control
+Kp = 10;
+Ki = 0.0;
+Kd = 10;
+mdlOutput = [];
+vqs = [];
+rd = {};
+as = {};
+fs = {};
+nb = [];
+%%
+for i=1:3
+
+ TM = acme.addTestMass(0.54,30e-3,tmgap(i),1,'dblInt',0.07);
+ el1 = setElectrode(acme.mdl,'electrode1',TM.area.x,-(TM.l.x/2+TM.gap.x), 0, 0 ,fCar); %Cage left
+ el2 = setElectrode(acme.mdl,'electrode2',TM.area.x,- TM.l.x/2 , 0, V_TM,fCar); %TM left
+ el3 = setElectrode(acme.mdl,'electrode3',TM.area.x,+ TM.l.x/2 , 0, V_TM,fCar); %TM right
+ el4 = setElectrode(acme.mdl,'electrode4',TM.area.x,+(TM.l.x/2+TM.gap.x), 0, 0 ,fCar); %Cage right
+
+
+ p = 1e-5; %Pa
+ gasThermalNoiseASD = @(f) gasBrownianNoise(p,acme.T,TM.l.x,TM.gap.x).*f.^0 + TM.mass*1.6e-15./sqrt(f);
+ gasThermalNoiseTS = [acme.t timeseriesFromASD(gasThermalNoiseASD, acme.enbw, acme.fSample, acme.nSteps).Data];
+ mdlOutput = [mdlOutput sim(acme.mdl)]; % might not work for drag free
+
+
+ [noises, sys] = nbFromSimulink(acme.mdl, acme.freqs, 'dof', dof);
+ nb = [nb nbGroupNoises(acme.mdl, noises, sys)];
+ vq = interp1(mdlOutput(i).simout.Time,mdlOutput(i).simout.Data,acme.t);
+ vqs = [vqs interp1(mdlOutput(i).simout.Time,mdlOutput(i).simout.Data,acme.t)];
+
+ % ts = [ts mdlOutput(i).simout.Data];%
+
+ [axx,fxx] = flpsd('ts', vq, acme.lpsdopts{:});
+ rd(end+1) = {getTFresponse(acme.mdl, io, fxx)}; %
+ as(end+1) = {axx};
+ fs(end+1) = {fxx};
+end
+
+%%
+
+ loglog(fs{1},as{1}./rd{1},'b', ...
+ fs{2},as{2}./rd{2},'r',...
+ fs{3},as{3}./rd{3},'g',...
+ nb(1).f,nb(1).asd,'--b',...
+ nb(2).f,nb(2).asd,'--r',...
+ nb(3).f,nb(3).asd,'--g')
+
+ legend('1mm gap (time domain simulation)', ...
+ '2mm ', ...
+ '4mm ', ...
+ '1mm gap (freq domain model linearization)', ...
+ '2mm',...
+ '4mm')
+ title('Acc asd gap parametrization')
+ ylabel('m/s² sqrt(hz)')
+ xlabel('frequency hz')
+ xlim([1e-4 1e-1]);
+ ylim([1e-15 1e-8]);
\ No newline at end of file
diff --git a/simulations/SGRS/sgrs_acc_script.m b/simulations/SGRS/sgrs_acc_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..df99e20bc66ff38e73fe1860b71c87e8290f0cbd
--- /dev/null
+++ b/simulations/SGRS/sgrs_acc_script.m
@@ -0,0 +1,76 @@
+% Electrostatic ACC Model
+% Following the Simplified Gravitational Reference Sensor model
+% The S-GRS is a scaled-down version of the LPF GRS
+% from: arXiv:2107.08545v1
+% Author: Arthur Reis
+% email: arthur.reis@aei.mpg.de
+% last review: 29.11.2022
+clear;
+acme = AcmeSim('sgrs_acc');
+acme.initModel;
+
+%% Accelerometer Parameters
+TM = acme.addTestMass(0.54,30e-3,1e-3,1,'dblInt',0.07);
+
+
+% If is non-drag-free:
+% driverGain = -20;
+% V_TM = 0.5;
+% Vp = 0; %V polarizing
+% Vx0 = 16; %V rms
+
+% If is drag-free:
+driverGain = 0;
+V_TM = 0.5;
+Vp = 0; %V polarizing
+Vx0 = 0; %V rms
+
+
+fCar = 1e5;
+el1 = setElectrode(acme.mdl,'electrode1',TM.area.x,-(TM.l.x/2+TM.gap.x), 0, 0 ,fCar); %Cage left
+el2 = setElectrode(acme.mdl,'electrode2',TM.area.x,- TM.l.x/2 , 0, V_TM,fCar); %TM left
+el3 = setElectrode(acme.mdl,'electrode3',TM.area.x,+ TM.l.x/2 , 0, V_TM,fCar); %TM right
+el4 = setElectrode(acme.mdl,'electrode4',TM.area.x,+(TM.l.x/2+TM.gap.x), 0, 0 ,fCar); %Cage right
+
+% Plane capacitors, each formed by two electrodes
+cap1 = setPlaneCapacitor(acme.mdl, 'planeCapacitor1', acme.eps0);
+cap2 = setPlaneCapacitor(acme.mdl, 'planeCapacitor2', acme.eps0);
+
+%CapNoiseASD = @(f) 1.854e-22 ./f + 5.0124e-19; %F/sqrt(hz)
+CapNoiseASD = @(f) 51e-9.*f.^0; %F/sqrt(hz)
+CapNoiseTS = [acme.t timeseriesFromASD(CapNoiseASD, acme.enbw, acme.fSample, acme.nSteps).Data];
+
+Cf = 10e-12; %F
+getSGRSCapSensorASD(acme, TM,1e-12,Cf,40,51e-9)
+ActuationNoiseASD = @(f) 20e-6.*f.^0; %V/sqrt(hz)
+ActuationNoiseTS = [acme.t timeseriesFromASD(ActuationNoiseASD, acme.enbw, acme.fSample, acme.nSteps).Data];
+
+p = 1e-5; %Pa
+gasThermalNoiseASD = @(f) gasBrownianNoise(p,acme.T,TM.l.x,TM.gap.x).*f.^0 + TM.mass*1.6e-15./sqrt(f);
+gasThermalNoiseTS = [acme.t timeseriesFromASD(gasThermalNoiseASD, acme.enbw, acme.fSample, acme.nSteps).Data];
+
+% Control
+Kp = 10;
+Ki = 0.0;
+Kd = 10;
+
+
+%%
+mdlOutput = sim(acme.mdl); % might not work for drag free
+%%
+vq = interp1(mdlOutput.simout.Time,mdlOutput.simout.Data,acme.t);
+ts = mdlOutput.simout.Data;%(10000:end);
+
+[axx,fxx] = flpsd('ts', vq, acme.lpsdopts{:});
+io(2) = linio(strcat(acme.mdl,'/BPF'),1,'output');
+io(1) = linio(strcat(acme.mdl,'/Gain1'),1,'input');
+readout2accTF = getTFresponse(acme.mdl, io, fxx); %
+
+%%
+dof = 'x';
+[noises, sys] = nbFromSimulink(acme.mdl, acme.freqs, 'dof', dof);
+nb = nbGroupNoises(acme.mdl, noises, sys);
+matlabNoisePlot(nb);
+loglog(fxx,axx./readout2accTF)
+xlim([1e-4 1e-1]);
+ylim([1e-15 1e-8]);
diff --git a/simulations/SGRS/sgrs_capSensor.slx b/simulations/SGRS/sgrs_capSensor.slx
new file mode 100644
index 0000000000000000000000000000000000000000..2b1f6d987869cd46aef649c438f3df9ba1df9927
Binary files /dev/null and b/simulations/SGRS/sgrs_capSensor.slx differ
diff --git a/simulations/SGRS/sgrs_capSensor_script.m b/simulations/SGRS/sgrs_capSensor_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..774ca2f778b953776193e2392c33581c17136898
--- /dev/null
+++ b/simulations/SGRS/sgrs_capSensor_script.m
@@ -0,0 +1,80 @@
+%% ACME - Simplified Capacitive Sensor Simulator
+% All references are from Alvarez's 2022 paper, unless otherwise noted.
+% This model is calculated in the frequency range 1e-4 - 1Hz. Effects of the
+% 100kHz injection signal are accounted for but the signal is not fully
+% simulated.
+% Only the dominant noise sources from the complete model are included
+% Arthur Reis, 28.11.2022
+clear;
+%%
+
+%mdl = 'sgrs_capSensor';
+%load_system(mdl);
+%% Simulation parameters
+% frequency and time
+acme = AcmeSim('sgrs_capSensor');
+
+%% Accelerometer Parameters
+% TM
+TM = acme.addTestMass(0.54,0.03,0.001,0.7,'');
+
+% Capacitive sensing
+nominal_cap = 3.9e-12;% F, capacitance between tm and elec. when dx = 0
+fTypical = 1e-2; %Hz
+Cf = 10e-12; %F feedback capacitor
+%% Voltage Reference
+Vd = 5; % Vpp
+fDetection = 1e5; %Hz
+capAttenuation = 0.1;
+
+
+%% Amplifier
+IA_BPF_gain = 40;
+chargeAmpNoiseASD = @(f) 51e-9.*f.^0; %V/sqrtHz
+chargeAmpNoiseTS = [acme.t timeseriesFromASD(chargeAmpNoiseASD, ...
+ acme.enbw, acme.fSample, acme.nSteps).Data];
+%% Noise Budget
+% this uses the SimulinkNB toolbox.
+% it doesn't run the simulation. Instead, it finds the transfer functions of the
+% linearized model and plots the noise budged
+[noises, sys] = nbFromSimulink(acme.mdl, acme.freqs, 'dof', 'x');
+nb = nbGroupNoises(acme.mdl, noises, sys);
+
+%% Time domain simulation
+mdlOutput = sim(acme.mdl);
+
+%% Finding the transfer functions
+% finding the transfer function between recorded signal and the signal we wish
+% to measure
+[axx,fxx] = flpsd( 'ts', mdlOutput.simout.Data,acme.lpsdopts{:});
+io(1) = linio(strcat(acme.mdl,'/Gain1'),1,'input');
+io(2) = linio(strcat(acme.mdl,'/IA_BPF'),1,'output');
+readout2capTF = getTFresponse(acme.mdl, io, fxx); % TF between readout and Cap1
+sensorTFexp = fit(fxx,readout2capTF, 'poly2')
+
+%% Expression for the total noise ASD
+% assuming ASD = A/f + B, where B is the high freq limit
+x = nb.sumNoise.f';
+y = nb.sumNoise.asd';
+sensorASDexp = fit(x,y, fittype(str2func(['@(a, x) a./x + ' num2str(y(end))])))
+
+%from time domain simulation
+xsim = fxx;
+ysim = axx./readout2capTF;
+sensorASDexpsim = fit(xsim,ysim, fittype(...
+ str2func(['@(a, x) a./x + ' num2str(ysim(end))])))
+
+%% Plots
+matlabNoisePlot(nb);
+xlabel('Frequency [Hz]');
+ylabel('Measured capacitance [F]');
+loglog(acme.freqs, sensorASDexp(acme.freqs),'b.', 'displayName', 'Fit from lin model');
+loglog(acme.freqs, sensorASDexpsim(acme.freqs),'k.','displayName', 'Fit from sim');
+loglog(fxx, axx./readout2capTF,'k-','LineWidth',1,'DisplayName', ...
+ 'ASD from time domain simulation');
+
+%% Finding the transfer functions
+% lpsd
+[axxs,fxxs] = flpsd( 'ts', mdlOutput.simout1.Data, 'fs', fSample, 'Kdes', Kdes, 'Lmin', Lmin, ...
+ 'Jdes', Jdes, 'win', 'hanning', 'olap', olap, 'order', order,...
+ 'scale', 'ASD', 'LTPDAmode',false);
\ No newline at end of file
diff --git a/utils/SimulinkNbLite/SimulinkNb/NbLibrary.mdl b/utils/SimulinkNbLite/SimulinkNb/NbLibrary.mdl
index 81f61a08ba9d3258b9cd1f90a423cbc38a620ed0..8d0b48a487cdb09b9cedd7e61376fce311f01eb0 100644
Binary files a/utils/SimulinkNbLite/SimulinkNb/NbLibrary.mdl and b/utils/SimulinkNbLite/SimulinkNb/NbLibrary.mdl differ
diff --git a/utils/SimulinkNbLite/SimulinkNb/example/OLvsCL/OLvsCL.mdl b/utils/SimulinkNbLite/SimulinkNb/example/OLvsCL/OLvsCL.mdl
index adc6b822102ec078da4da3307a51b371db1bd22a..707609b8f10b837ee07ed2e572fbcb66c77e109f 100644
Binary files a/utils/SimulinkNbLite/SimulinkNb/example/OLvsCL/OLvsCL.mdl and b/utils/SimulinkNbLite/SimulinkNb/example/OLvsCL/OLvsCL.mdl differ
diff --git a/utils/SimulinkNbLite/SimulinkNb/nbFromSimulink.m b/utils/SimulinkNbLite/SimulinkNb/nbFromSimulink.m
index ef12f2ed0132ef9690bc7f8872059ae2af948bea..be832eda21db1d8147d73183ee5288432daebea0 100644
--- a/utils/SimulinkNbLite/SimulinkNb/nbFromSimulink.m
+++ b/utils/SimulinkNbLite/SimulinkNb/nbFromSimulink.m
@@ -142,10 +142,12 @@ end
io = [ioSink ioCal ioSource];
%% Perform the linearization using FlexTf functions
-
+warning('off')
% Don't abbreviate I/O block names (it's faster that way)
linOpt = linoptions('UseFullBlockNameLabels', 'on');
+
[sys, flexTfs] = linFlexTf(mdl, io, linOpt);
+
% Attempt to improve numerical accuracy with prescale
minPosFreq = min(freq(freq>0));
maxPosFreq = max(freq(freq>0));
@@ -160,7 +162,7 @@ sys = frd(sys, freq, 'Units', 'Hz');
% (UseFullBlockNameLabels appends signal names to the block names)
sys.InputName = [{nbNoiseSink nbNoiseCal} nbNoiseSources'];
sys.OutputName = nbNoiseSink;
-
+warning('on')
%% Apply noise/calibration TFs to each NbNoiseSource's spectrum
cal = 1/sys(2);
diff --git a/utils/convert2OlderMatlab.m b/utils/convert2OlderMatlab.m
index a8728985e594e3dcbabfb114e4d9bfc691a48631..1f652f8da95cdc580ec3d6565a7535511af90c21 100644
--- a/utils/convert2OlderMatlab.m
+++ b/utils/convert2OlderMatlab.m
@@ -1,2 +1,25 @@
proj = currentProject;
-Simulink.exportToVersion(proj,'Acme_R2020a','R2020a');
\ No newline at end of file
+Simulink.exportToVersion(proj,'Acme_R2020a','R2020a');
+
+%unzip
+
+%remove files
+
+%copy to other folder
+
+%save simulations "manually"
+
+ % make temp folder with simulations files
+ % proj.addFile
+ % copy xhps toolbox to temp folder
+ % add xhps to project (how to do that programmatically??)
+ % exportToVersion on temp folder
+ % unzip
+ % copy simulation files
+ % exclude xhps files (%git .ignore?)
+ % system('git .....')
+
+%git commit
+
+% add starter script
+ % add xhps path/start up
\ No newline at end of file