v 0.2.5 Fixed actuation voltage calculation, added Satellite class

core/AcmeSim.m

@@ -198,6 +198,10 @@ classdef AcmeSim
             end
                 
         end
+
+        function openModel(obj)
+            open(obj.ModelName);
+        end
     end
 end
 

core/dynamics/Satellite.m

+classdef Satellite
+%SATELLITE Summary of this class goes here
+%   Detailed explanation goes here
+
+properties
+    mass (1,1) double = 470;
+    altitude (1,1) double = 450;
+    Cd (1,1) double = 2.25;
+    area (1,1) double = 0.9;
+end
+
+properties(Dependent)
+    vel_orb
+end
+
+methods
+    function obj = Satellite()
+        %SATELLITE Construct an instance of this class
+        %   Detailed explanation goes here
+
+    end
+    
+    function v = get.vel_orb(obj)
+        v = orbitalvelocity(obj.altitude);
+    end
+
+
+    function acc = typicalDragAcc(obj)
+        %UNTITLED Summary of this function goes here
+        %   Detailed explanation goes here
+        rho = Satellite.harrispriester(obj.altitude);
+        acc = 0.5*rho*obj.Cd*obj.area*obj.vel_orb^2 /obj.mass;
+    end
+end
+
+methods (Static)
+    function rho = harrispriester(h)
+        v = [100 497400.0 497400.0 
+        120 24900.0 24900.0
+        130 8377.0 8710.0 
+        140 3899.0 4059.0 
+        150 2122.0 2215.0 
+        160 1263.0 1344.0 
+        180 528.3 601.0 
+        190 361.7 429.7 
+        170 800.8 875.8
+        200 255.7 316.2 
+        210 183.9 239.6 
+        220 134.1 185.3
+        230 99.49 145.5 
+        240 74.88 115.7 
+        250 57.09 93.08 
+        260 44.03 75.55 
+        270 34.30 61.82 
+        280 26.97 50.95 
+        290 21.39 42.26 
+        300 17.08 35.26
+        320 10.99 25.11
+        340 7.214 18.19
+        360 4.824 13.37 
+        380 3.274 9.955
+        400 2.249 7.492 
+        420 1.558 5.684
+        440 1.091 4.355
+        460 0.7701 3.362
+        480 0.5474 2.612
+        500 0.3916 2.042
+        520 0.2819 1.605
+        540 0.2042 1.267
+        560 0.1488 1.005
+        580 0.1092 0.7997
+        600 0.08070 0.6390
+        620 0.06012 0.5123
+        640 0.04519 0.4121
+        660 0.03430 0.3325
+        680 0.02632 0.2691
+        700 0.02043 0.2185
+        720 0.01607 0.1779
+        740 0.01281 0.1452
+        760 0.01036 0.1190
+        780 0.008496 0.09776
+        800 0.007069 0.08059
+        840 0.004680 0.05741
+        880 0.003200 0.04210
+        920 0.002210 0.03130
+        960 0.001560 0.02360
+        1000 0.001150 0.01810];
+        rho = interp1(v(:,1),v(:,3),h)*1e-12; %rho in kg/m³
+     end
+    end
+end
+

core/noises/CPDnoise.m

-function fp = CPDnoise(varvoltrms, gap)
+function fp = CPDnoise(varvoltrms, gap0, varargin)
 %CPDNOISE Summary of this function goes here
 %   Detailed explanation goes here
 eps0 = 8.8e-12;
 
-kpeak = 1.3/gap;
+if length(varargin) == 1
+    l = varargin{1};
+    k0 = 1/l;
+else
+    k0 = 1.3/gap0;
+end
 
-kmax = 1.2*kpeak;
-kmin = 0.8*kpeak;
+kmax = 1.2*k0;
+kmin = 0.8*k0;
 
-func = @(k) k.^3./(sinh(gap.*k)).^2;
+func = @(k) k.^3./(sinh(gap0.*k)).^2;
 
 intr = integral(func,kmin,kmax);
 

core/noises/actuationRequirement.m

@@ -10,7 +10,8 @@ function [acc, volt_required] = actuationRequirement(f_ng, a_ng, TM, cap)
 
 %psd = a_ng.^2;
 %loglog(fxx_input_forces, psd)
-acc = sqrt(trapz(f_ng,a_ng.^2)); %mean? acceleration m/s2
+
+acc = sqrt(trapz(f_ng(10:60),a_ng(10:60).^2)); %mean? acceleration m/s2
 forces_ng = acc*TM.mass;
 volt_required = cap.voltage_from_force(forces_ng);
 end

core/noises/actuationRequirementSimplified.m

+function [acc, volt_required] = actuationRequirementSimplified(sat, TM, cap)
+%ACTUATIONREQUIREMENT Calculates required actuation voltage to keep test mass
+% stable
+%   Detailed explanation goes here
+% Input:
+%    [f_ng, a_ng] - ASD (frequencies and amplitudes) of the non gravitational
+%    forces that impact the satellite
+% Output:
+%
+
+%psd = a_ng.^2;
+%loglog(fxx_input_forces, psd)
+arguments
+sat Satellite
+TM TestMass
+cap Capacitor
+end
+acc = sat.typicalDragAcc;
+volt_required = cap.voltage_from_force(acc*TM.mass);
+end
+

simulations/SGRS/Copy_of_sgrs_acc_script_tests.m

+% 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: 03.02.2023
+clear;
+acme = AcmeSim('sgrs_acc_tests',1e4,1e4);
+%acme.openModel
+%% Accelerometer Parameters 
+%TM & EH Geometry
+side = 0.03;
+gap = 0.001;
+elSides = [side, side];
+density = 20099;
+
+% % % If is drag-free:
+% driverGain = 0;
+% V_TM = 0.5;
+% Vp = 0; %V polarizing 
+
+% If is non-drag-free:1000
+driverGain = -1;
+V_TM = 0.5;
+%Vp = 0; %V polarizing 
+
+fCar = 1e5;
+
+sat = Satellite;
+%sat.altitude = 450; %edit here if needed. These are the defaults for GFO
+% sat.mass = 470
+% sat.Cd = 2.25
+% sat.area = 0.9
+
+TM   = acme.addTestMass('dblInt',side,gap,'density',density);
+el1  = acme.addElectrode('electrode1', elSides, -(side/2 + gap), 0,    0, fCar);
+el2  = acme.addElectrode('electrode2', elSides, -side/2        , 0, V_TM, fCar);
+el3  = acme.addElectrode('electrode3', elSides,  side/2        , 0, V_TM, fCar);
+el4  = acme.addElectrode('electrode4', elSides, +(side/2 + gap), 0,    0, fCar);
+cap1 = acme.addCapacitor('planeCapacitor1', 'ElectrodePair', [el1, el2]);
+cap2 = acme.addCapacitor('planeCapacitor2', 'ElectrodePair', [el3, el4]);
+
+%CapNoiseASD = @(f) 1.854e-22 ./f + 5.0124e-19; %F/sqrt(hz)
+CapNoiseASD = @(f) 51e-9.*f.^0; %F/sqrt(hz)
+CapNoiseTS  = acme.noiseTS(CapNoiseASD);
+
+CPDnoiseASD = @(f) cap1.force_from_CPD(10e-3)*f.^0; %worst case
+CPDnoiseTS = acme.noiseTS(CPDnoiseASD);
+Cf = 10e-12; %F
+
+[a, v] = actuationRequirementSimplified(sat,TM,cap1);
+ActuationNoiseASD = @(f) v * 1e-6*f.^0; %V/sqrt(hz)
+ActuationNoiseTS  = acme.noiseTS(ActuationNoiseASD);
+Vx0 = 2*v;
+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.noiseTS(gasThermalNoiseASD);
+
+% Control
+Kp = 10;
+Ki = 0.0;
+Kd = 10;
+
+%%
+mdlOutput = sim(acme.ModelName);     % might not work for drag free
+%%
+
+vq = interp1(mdlOutput.simout.Time,mdlOutput.simout.Data,acme.Time);
+ts = mdlOutput.simout.Data;%(10000:end);
+%vq = interp1(out.simout.Time,mdlOutput.simout.Data,acme.Time);
+%ts = out.simout.Data;%(10000:end);
+[axx,fxx] = flpsd('ts', vq, acme.flpsdOpts{:});
+io(2) = linio(strcat(acme.ModelName,'/BPF'),1,'output');
+io(1) = linio(strcat(acme.ModelName,'/Gain1'),1,'input');
+readout2accTF = getTFresponse(acme.ModelName, io, fxx); % 
+
+%%
+dof = 'x';  
+[noises, sys] = nbFromSimulink(acme.ModelName, acme.Frequencies, 'dof', dof);  
+nb = nbGroupNoises(acme.ModelName, noises, sys);  
+matlabNoisePlot(nb);  
+loglog(fxx,axx./readout2accTF)
+xlim([1e-4 1e-1]);
+ylim([1e-15 1e-8]);
+
+%%
+% [agt,fgt] = flpsd('ts', gasThermalNoiseTS(:,2), acme.flpsdOpts{:});
+% io(2) = linio(strcat(acme.ModelName,'/Sum6'),1,'output');
+% gast2accTF = getTFresponse(acme.ModelName, io, fgt); % 
+% 
+[acp,fcp] = flpsd('ts', CapNoiseTS(:,2), acme.flpsdOpts{:});
+io(2) = linio(strcat(acme.ModelName,'/Sum'),1,'output');
+cap2accTF = getTFresponse(acme.ModelName, io, fcp); % 
+
+[aac,fac] = flpsd('ts', ActuationNoiseTS(:,2), acme.flpsdOpts{:});
+io(2) = linio(strcat(acme.ModelName,'/Sum7'),1,'output');
+act2accTF = getTFresponse(acme.ModelName, io, fgt); % 
+
+hold on
+%%
+% loglog(fgt,agt./gast2accTF,'DisplayName','Gas Thermal');
+matlabNoisePlot(nb);  
+loglog(fxx,axx./readout2accTF,'DisplayName','Measurements')
+loglog(fcp,acp./cap2accTF,'DisplayName','Cap');
+loglog(fac,aac./act2accTF,'DisplayName','Act');
+