close all
clear
clc
set(0,'DefaultLineLineWidth',2)
set(0,'DefaultAxesFontSize',16)
set(0,'DefaultAxesXGrid','on')
set(0,'DefaultAxesYGrid','on')
set(0,'DefaultAxesZGrid','on')
%%
verbosePlotFlag = true;
printFigs = true;
saveFilters = false;
susType = 'HXDS';
svnDir.sus = '/ligo/svncommon/SusSVN/sus/trunk/';
svnDir.sei = '/ligo/svncommon/SeiSVN/seismic/';
svnDir.isc = '/ligo/svncommon/IscCVS/iscmodeling/';
toolsDir.sus = [svnDir.sus 'Common/MatlabTools/'];
toolsDir.sei = [svnDir.sei 'Common/MatlabTools/'];
bscDir.isc = [svnDir.isc 'LentickleAligo/SeismicIsolation/'];
saveDir = [svnDir.sus susType '/Common/FilterDesign/MatFiles/'];
plotDir = [svnDir.sus susType '/Common/FilterDesign/Results/'];
% LHO aLOG 61255
author = ', R. Kumar';
designDate = '2022-01-10';
dof = {'L';'T';'V';'R';'P';'Y'};
sampleFreq = 16384;
sampleTime = 1/sampleFreq;
digType = 'tustin';
freq = logspace(-1,2,1001);
buildType = [lower(susType) 'opt_doublep'];
susTag = 'H1:SUS-ZM6';
optic = lower(susTag(regexp(susTag,'\-')+1:end));
EPICSgains = [-1.5 0 0 0 -1.0 -1.0];
% Because the OSEM sensors come into the damping path pre-calibrated into
% um, so that damping filter units are [(DAC ct) / um]. But,
% createdampingfilterstruct is calibrating the filters as though they were
% in [(DAC ct) / (ADC ct)]. So, we need to invert that sensor chain
% calibration gain within the damping loops by include a dimensionful gain
% of [(ADC ct) / um]. This gain should nominally in the "from_um" bank.
% <---------------- Sensor Chain --------------->
% OSEM SatAmp
% Sensitivity TransImp. ADC
% [uA / mm] [mm / um] [A/uA] [V/A] [ADC ct/V]
from_um = (76.29/0.7)*(1e3/1e6)*(1/1e6)*(240e3)*(2^16/40);
gains = EPICSgains; % HXDS
useEllip50 = true;
figFileTag = ['dampingfilters_' susType '_' upper(optic) '_' designDate];
filterFileName = [saveDir figFileTag '.mat'];
%%
addpath(toolsDir.sei,toolsDir.sus,bscDir.isc)
%%
[trans, rot] = seisHAM(freq);
hamData.(optic).asd(:,1) = trans; % X
hamData.(optic).asd(:,2) = trans; % Y
hamData.(optic).asd(:,3) = rot; % RZ
hamData.(optic).asd(:,4) = trans; % Z
hamData.(optic).asd(:,5) = rot; % RX
hamData.(optic).asd(:,6) = rot; % RY
hamData.freq = freq;
%% Model Filters
rolloff.model = create_cdf_filter(0,pair(100,55),1/(2*pi),freq,sampleTime,digType);
% function [out1, out2] = mycheby1(freq,order,ripple)
[cheby5.proto,cheby5.zeros_Hz,cheby5.poles_Hz,cheby5.gain] = mycheby1(5,2,2);
cheby5.model = create_cdf_filter(cheby5.zeros_Hz,cheby5.poles_Hz,1,freq,sampleTime,digType);
% Gains pulled for ZM6 from FM6 foton file H1SUSSQZOUT.txt
%% LONG
gain_L = 182.15;
% Took the ratio between the measurement from LHO:62585 and model without
% the scale.
% Gain at | Meas. | Model | Ratio (Meas/Model)
% 0.5Hz | 0.15381 | 0.18357 | 0.83788
% 10Hz | 0.02176 | 0.02246 | 0.96883
% The average between the two ratios: 0.90336
measscale_L = 0.90336;
% Total Filter
% (DAC ct / ADC ct) = (um/ADC_ct) ( ----------- DAC ct / um ----------------- ) (meas/model gain)
totalFilter.L.c = (1/from_um) * gain_L * rolloff.model.c * cheby5.model.c * measscale_L;
%% TRANS
gain_T = 0;
totalFilter.T.c = (1/from_um) * gain_T * rolloff.model.c * cheby5.model.c;
%% VERT
gain_V = 0;
totalFilter.V.c = (1/from_um) * gain_V * rolloff.model.c * cheby5.model.c;
%% ROLL
gain_R = 0;
totalFilter.R.c = (1/from_um) * gain_R * rolloff.model.c * cheby5.model.c;
%% PITCH
gain_P = 0.1028;
% Took the ratio between the measurement from LHO:62585 and model without
% the scale.
% Gain at | Meas. | Model | Ratio (Meas/Model)
% 0.5Hz | 0.02123 | 0.02424 | 0.87583
% 10Hz | 0.00513 | 0.00563 | 0.91119
% The average between the two ratios: 0.89351
measscale_P = 0.89351;
totalFilter.P.c = (1/from_um) * gain_P * rolloff.model.c * cheby5.model.c * measscale_P;
%% YAW
gain_Y = 0.1028;
% Took the ratio between the measurement from LHO:62585 and model without
% the scale.
% Gain at | Meas. | Model | Ratio (Meas/Model)
% 0.5Hz | 0.01823 | 0.02049 | 0.8897
% 10Hz | 0.00743 | 0.00715 | 1.0392
% The average between the two ratios: 0.96445
measscale_Y = 0.96445;
totalFilter.Y.c = (1/from_um) * gain_Y * rolloff.model.c * cheby5.model.c * measscale_Y;
%%
% Total Filter *MUST* be in units of (DAC ct / ADC ct)
[totalFilter,calibFilter,calibration] = createdampingfilterstruct(susType,totalFilter,gains,freq,sampleTime,digType,designDate,verbosePlotFlag);
%%
[opticDisp, eulerDisp, st2isiDisp, plant, oltf, stability, cltf, supp, susModel] = plothxdsdampingcontroldesign(freq,buildType,susTag,dof,calibFilter,hamData,verbosePlotFlag);
%% DARM Expectations
%[darm, srcl, SRCL2DARM] = plotsrclvsdarmreqs_hxts(freq,susType,opticDisp,verbosePlotFlag);
%%
% close([1:46]) % Filters Only
% close([1 8:46]) % Seismic Input Motion Only
% close([1:7 9:14 16:21 23:28]) % Loop Design Only
% close([1:7 8:11 13:16 17:24 25:32 33:40 41:44 46:49 50:53 55:58 100 110 410]) % Global TFs Only
% close([1:7 8 10:16 17 19:24 25 27:32 33 35:40 41 43:49 50 52:58 100 110 410]) % Impulse Response Only
% close([1:7 8:9 11:16 17:18 20:24 25:26 28:32 33:34 36:40 41:42 44:49 50:51 53:58 100 110 410]) % Seismic Noise Only
% close([1:7 8:10 12:16 17:19 21:24 25:27 29:32 33:35 37:40 41:43 45:49 50:52 54:58 100 110 410]) % M1 Sensor Noise Only
% close([1:7 8:15 17:23 25:31 33:39 41:48 50:57 100 110 410]) % Total Noise Only
% close([1:7 10 13:15 17:24 25:32 33:40 42:44 46:48 50:58 100 110 410]) % L Design Only
% close([1:7 8:16 19 21:23 25:32 35 37:39 41:49 50:58 100 110 410]) % T Design Only
% close([1:7 8:16 17:24 27:31 33:40 41:49 50:58 100 110]) % V Design Only
% close([1:7 8:16 19 21:23 25:32 35 37:39 41:49 50:58 100 110 410]) % R Design Only
% close([1:7 10 13:15 17:24 25:32 33:40 42:44 46:48 50:58 100 110 410]) % P Design Only
% close([1:7 8:16 17:24 25:32 33:40 41:49 52:53 55:57 100 110 410]) % Y Design Only
%%
%%
if printFigs
figNames = {[plotDir figFileTag '_normalized.pdf'];...
[plotDir figFileTag '_calibrated.pdf'];...
[plotDir figFileTag '_loopdesign_L.pdf'];...
[plotDir figFileTag '_impulseresp_L.pdf'];...
[plotDir figFileTag '_seismicnoise_L.pdf'];...
[plotDir figFileTag '_sensornoise_L.pdf'];...
[plotDir figFileTag '_allstages2optic_L.pdf'];...
[plotDir figFileTag '_m2actuatornoise_L.pdf'];...
[plotDir figFileTag '_totalbudget_L.pdf'];...
[plotDir figFileTag '_loopdesign_P.pdf'];...
[plotDir figFileTag '_impulseresp_P.pdf'];...
[plotDir figFileTag '_seismicnoise_P.pdf'];...
[plotDir figFileTag '_sensornoise_P.pdf'];...
[plotDir figFileTag '_allstages2optic_P.pdf'];...
[plotDir figFileTag '_m2actuatornoise_P.pdf'];...
[plotDir figFileTag '_totalbudget_P.pdf'];...
[plotDir figFileTag '_loopdesign_Y.pdf'];...
[plotDir figFileTag '_impulseresp_Y.pdf'];...
[plotDir figFileTag '_seismicnoise_Y.pdf'];...
[plotDir figFileTag '_sensornoise_Y.pdf'];...
[plotDir figFileTag '_allstages2optic_Y.pdf'];...
[plotDir figFileTag '_m2actuatornoise_Y.pdf'];...
[plotDir figFileTag '_totalbudget_Y.pdf'];...
[plotDir figFileTag '_seismicinput_L.pdf'];...
[plotDir figFileTag '_seismicinput_T.pdf'];...
[plotDir figFileTag '_seismicinput_V.pdf'];...
[plotDir figFileTag '_seismicinput_R.pdf'];...
[plotDir figFileTag '_seismicinput_P.pdf'];...
[plotDir figFileTag '_seismicinput_Y.pdf'];...
[plotDir figFileTag '_allsensornoiseinput.pdf']};
mergeCommand = [svnDir.sus 'Common/Misc/pdfmerge ' plotDir figFileTag '.pdf '];
figCount = 1;
for iFig = [100 110 8:14 15:21 22:28 2:7 1]
figure(iFig)
FillPage('w')
IDfig(author)
saveas(gcf,figNames{figCount})
mergeCommand = [mergeCommand ' ' figNames{figCount}];
figCount = figCount+1;
end
system(mergeCommand);
end
if saveFilters
README = {['Designed by ' author ', on ' designDate];...
['using ${SusSVN}/sus/trunk/' susType '/Common/FilterDesign/design_damping_' susType '_' datestr(designDate,'yyyymmdd') '.m'];...
'gains are EPICs value gains in the M1 filter bank, rows are the 6 DOFs';...
'For totalFilter.DOF, .c = continuous zpk model, .ss = continuous state space model, .d = discrete zpk model, .fd = freqresp of .d';...
'calibFilter = calibration * gains * totalFilter';...
'totalFilter is what goes into foton, calibFilter is used to wrap around models'};
save(filterFileName,'boost*','rolloff*','ellip*','totalFilter','calibFilter','gains','calibration','freq','README');
save([filterFileName(1:end-4) '_model.mat'],'freq','calibFilter','opticDisp','eulerDisp','st2isiDisp','stability','plant','oltf','cltf','supp','susModel');
end