Skip to content
Snippets Groups Projects
Commit f0b44dc1 authored by Arthur Reis's avatar Arthur Reis
Browse files

Added SGRS model with thermal and cap noise

parent 7c2108d8
No related branches found
No related tags found
No related merge requests found
Showing
with 403 additions and 58 deletions
*.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
# 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
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
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
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{:});
switch p.Results.mode
case 'linear'
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));
case 'rotational'
...
end
%set_param(blockHandle,'dt_solver',num2str(p.Results.dt));
%if p.Results.solver
end
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
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
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
function [outputArg1,outputArg2] = outgassing(inputArg1,inputArg2)
%OUTGASSING Summary of this function goes here
% Detailed explanation goes here
outputArg1 = inputArg1;
outputArg2 = inputArg2;
end
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
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
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
%% 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)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment