% Noise budget script for bKAGRA Phase 1 Michelson % Author: Yutaro Enomoto, Yuta Michimura %% Add paths clear all close all findNbSVNroot; % find the root of Simulink NB addpath(genpath([NbSVNroot 'Common/Utils'])); addpath(genpath([NbSVNroot 'Dev/Utils/'])); %% Define some parameters, filters, and noises % see the Simulink file (noiseModel) for each parameter freq = logspace(-1,4,1000); noiseModel = 'Michelson'; IFO = 'K1'; site = 'kamioka'; % get online data (see http://klog.icrr.u-tokyo.ac.jp/osl/index.php?r=1357) OnlineData = GWData; %GPStime = 1208513418; % time to get online data (2018Apr23 19:10 JST) [~,GPStime]=system('tconvert Apr 30 2018 18:30 JST'); [~,GPStime]=system('tconvert May 06 2018 12:01 JST'); figdir = './results/'; % directory to save figures GPStime=str2num(GPStime); duration = 64; % data get time duration freqsamp = 16384; % sampling frequency freqmin = 0.1; % minimum frequency (= frequency reslution) Ndata = ceil(freqsamp/freqmin); ndata = 2^nextpow2(Ndata); duration2 = 64; % data get time duration freqsamp2 = 2048; % sampling frequency freqmin2 = 0.1; % minimum frequency (= frequency reslution) Ndata2 = ceil(freqsamp2/freqmin2); ndata2 = 2^nextpow2(Ndata2); figdir = ['./results/GPS',num2str(GPStime),'/']; % directory to save figures mkdir(figdir); %% get LiveParts parameters liveParts(noiseModel,GPStime,duration,freq); % online data not available anymore (as of 2017.2.23) % Manually input LiveParts parameters instead %PD_DOF_MTRX = [0,0,0,1,0,0]; %OUTPUT_MTRX = [1;-1;0]; % 1x3 or 3x1 depends on OS???? %% % Constants c = 299792458.; % speed of light in vacuum (m/s) hP = 6.62606957e-34; % Planck's Constant [J*s] lamb = 1064e-9; % laser wavelength (m) nu = c/lamb; % laser frequency (Hz) Larm = 3e3; % arm length [m] % Dictionaries Filt = containers.Map; % digital filter transfer functions [cnts/cnts] Susp = containers.Map; % suspension actuation transfer functions [m/N] Noise = containers.Map; % all sorts of noises PdResp = containers.Map; % PD rensponse [A/W] TransImp = containers.Map; % PD trans impedance [V/A] Elec = containers.Map; % electronics gain [V/V] CoilDriver = containers.Map; % coil driver V-I conversion [Ohm] Coil = containers.Map; % coil gain [N/V] % ADC and DAC ADC_V2C = 2^16/40; % ADC factor [cnts/V] DAC_C2V = 20/2^16; % DAC factor [V/cnts] % PD signal chain and noises for pd = {'REFL'} pd = char(pd); PdResp([pd,'_DC']) = 0.3; % A/W (Thorlabs PDA100A) TransImp([pd,'_DC']) = 0.75e3; % V/A (0dB setting; see https://www.thorlabs.co.jp/thorcat/13000/PDA100A-Manual.pdf) PdResp([pd,'_RF']) = 0.8; % A/W (PD is PerkinElmer C30642G (see JGW-D1201280-v2) / ) TransImp([pd,'_RF']) = 200/0.8; % V/A (see http://granite.phys.s.u-tokyo.ac.jp/internal/2012_m_ishigaki.pdf) noise = load('./Noises/20160402/LSC-REFL_PDA1_DC_OUT_DQ_dark.txt'); % dark noise measured when PR3 misaligned Noise([pd,'_DC_PD']) = interp1(noise(:,1), noise(:,2), freq, 'linear'); % PD dark noise in [counts/rtHz]; noise = load('./Noises/20180425/LSC-REFLAIR_A_RF17_Q_ERR_dark_adc.txt'); % dark noise measured when IFI viewport was covered, on 20180424 Noise([pd,'_RF_PD']) = interp1(noise(:,1), noise(:,2), freq, 'linear'); % RF PD dark noise in [counts/rtHz]; end DeModGain = 10.9; % Demodulation gain (see LIGO-F1100004-v4) WhiteningGain = 16; % http://klog.icrr.u-tokyo.ac.jp/osl/?r=4695, http://klog.icrr.u-tokyo.ac.jp/osl/?r=4831 noise = load('./Noises/adc_noise.txt'); %noise = load('./Noises/20180425/LSC-REFLAIR_A_RF17_Q_ERR_adc.txt'); %[data,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-ADCNOISE_GND_OUT_DQ'); % use online data %[noisep,noisef] = pwelch(data,hanning(ndata),[],ndata,freqsamp); % calculate power spectral density %noise = [noisef sqrt(noisep)]; %Noise('ADC') = interp1(noise(:,1), noise(:,2), freq, 'linear'); % ADC noise [cnts/rtHz] Noise('ADC') = 0; % since ADC noise is included in PD dark noise now % Intensity noise (see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1209) %% TO BE UPDATED FOR bKAGRA PHASE 1 % [data,time] = OnlineData.fetch(GPStime,duration,'K1:IMC-MCL_TRANS_OUT_DQ'); % use online data data = importdata('./Noises/GPS1144857600/K1:IMC-MCL_TRANS_OUT_DQ.txt'); [noisep,noisef] = pwelch(data,hanning(ndata),[],ndata,freqsamp); % calculate power spectral density noise = [noisef sqrt(noisep)]; Noise('RIN') = interp1(noise(:,1), noise(:,2), freq, 'linear')/0.4; % relative intensity noise [/rtHz] INT_NOISE_COUP = 0.0096; % intensity noise coupling to REFL DC PD [W] fitted to fit with measured spectrum % Optical lever control loops noise = load('./Noises/20180430_oplev_BS_PR2/20180430_oplev_BS_PR2.txt'); Noise('BSOPLEV_P') = interp1(noise(:,1), noise(:,2), freq, 'linear'); % BS pitch control noise in counts [counts/rtHz] Noise('BSOPLEV_Y') = interp1(noise(:,1), noise(:,3), freq, 'linear'); % BS yaw control noise in counts [counts/rtHz] Noise('PR2OPLEV_P') = interp1(noise(:,1), noise(:,4), freq, 'linear'); % PR2 pitch control noise in counts [counts/rtHz] Noise('PR2OPLEV_Y') = interp1(noise(:,1), noise(:,5), freq, 'linear'); % PR2 yaw control noise in counts [counts/rtHz] BSoplev_pitch_2_MICH = 1.1e-11; % coupling coefficient as measured on the 30th BSoplev_yaw_2_MICH = 9.3e-11; % coupling coefficient as PR2oplev_pitch_2_MICH = 3.6e-11; PR2oplev_yaw_2_MICH = 5e-12; % IFO Optical responses and noises %% TO BE UPDATED FOR bKAGRA PHASE 1 las = 3.3299; % Schnupp asymmetry [m] P0 = 250e-3; % input power to BS [W] %Ppd = 8e-3; % power at REFL PD [W] ~3000 counts Ppd = 56e-6; % power at REFL PD [W] ( http://klog.icrr.u-tokyo.ac.jp/osl/?r=4903 ) Freq2MICH = las/nu; % Laser frequency noise to MICH coupling [m/Hz] noise = load('./Noises/IMCnoise/freqnoise.txt'); Noise('LaserFreq') = interp1(noise(:,1), noise(:,2), freq, 'linear'); % Laser frequency noise [Hz/rtHz]; calculated with /kagranoisebudget/FSS/KAGRA_FSS.mdl model %MICHOpticalGain = -2*pi*nu/c*P0; % W/m (theoretical) MICHOpticalGainDC = 3.4e10/ADC_V2C/PdResp('REFL_DC')/TransImp('REFL_DC'); % measured value (see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1340); this is nominal value, adjusted later using calibration peak Noise('REFL_DC_SHOT') = sqrt(2*hP*nu*Ppd); % W/rtHz %MICHOpticalGainRF = 8500*4*pi/lamb/ADC_V2C/PdResp('REFL_RF')/TransImp('REFL_RF')/DeModGain; % measured value (http://klog.icrr.u-tokyo.ac.jp/osl/?r=4795) MICHOpticalGainRF = 8500*4*pi/lamb/ADC_V2C/PdResp('REFL_RF')/TransImp('REFL_RF')/DeModGain/WhiteningGain; % measured value (http://klog.icrr.u-tokyo.ac.jp/osl/?r=4795) Noise('REFL_RF_SHOT') = sqrt(2*hP*nu*Ppd); % W/rtHz THIS IS A FAKE (power at PD should be measured / Shot noise should be calculated in a case of dark fringe) % Suspension chain and noises %% TO BE UPDATED FOR bKAGRA PHASE 1 CoilDriver('HighPower') = 40*2; % Ohm; See https://dcc.ligo.org/LIGO-T0900232, https://dcc.ligo.org/LIGO-D0902747 Noise('HighPower') = 57e-9./freq+2.1e-9; % coil driver noise [V/rtHz]; See Figure 3 of https://dcc.ligo.org/LIGO-T080014 R1 = 750; R2 = 3900; C = 0.68e-6; CoilDriver('LowPower') = myzpk('zpk',[1/(2*pi*C*R1)],[1/(2*pi*C*(R1+R2))],R2*2); %Ohm; See https://dcc.ligo.org/LIGO-D070481, https://dcc.ligo.org/LIGO-T0900233 R2 = 700; CoilDriver('LowPowerMod') = myzpk('zpk',[1/(2*pi*C*R1)],[1/(2*pi*C*(R1+R2))],R2*2); %Ohm; modified version for Type-A MN and IM Noise('LowPower') = 170e-9./freq+6.4e-9; % coil driver noise [V/rtHz]; See https://dcc.ligo.org/LIGO-T0900233 (1/f noise at low frequency is assumed) % BS BSIMCoilN=1; BSTMCoilN=4; BSEUL2OSEM=0.25; % for EUL2OSEM matrix element (added by K. Komori on Apr 24, 2018 http://klog.icrr.u-tokyo.ac.jp/osl/?r=4811) Rcoil = 13; % resistance of coil CoilDriver('BS_IM') = CoilDriver('HighPower'); CoilDriver('BS_TM') = CoilDriver('HighPower'); Coil('BS_IM') = -BSIMCoilN*1.25/(CoilDriver('BS_IM')+Rcoil); % N/V; See http://gwdoc.icrr.u-tokyo.ac.jp/cgi-bin/DocDB/ShowDocument?docid=3239 Coil('BS_TM') = BSTMCoilN*0.0138/(CoilDriver('BS_TM')+Rcoil); % N/V; See http://gwdoc.icrr.u-tokyo.ac.jp/cgi-bin/DocDB/ShowDocument?docid=3239 Noise('BS_IM_Coil') = Noise('HighPower'); Noise('BS_TM_Coil') = Noise('HighPower'); Elec('BS_TM_Whitening') = 1; % no whitening? Elec('BS_TM_DeWhitening') = 1; % no dewhitening? Elec('BS_IM_Whitening') = 1; % no whitening? Elec('BS_IM_DeWhitening') = 1; % no dewhitening? % ETMY/ETMX ETMIMCoilN=2; ETMTMCoilN=4; for mir = {'ETMX','ETMY'} Rcoil = 13; % resistance of coil (TO BE CONFIRMED) mir = char(mir); CoilDriver('BS_IM') = CoilDriver('HighPower'); CoilDriver('BS_TM') = CoilDriver('HighPower'); Coil([mir,'_IM']) = ETMIMCoilN*0.015/(CoilDriver('BS_IM')+Rcoil); % N/V; See http://gwdoc.icrr.u-tokyo.ac.jp/cgi-bin/private/DocDB/ShowDocument?docid=5938 Coil([mir,'_TM']) = ETMTMCoilN*0.0015/(CoilDriver('BS_TM')+Rcoil); % N/V; See http://gwdoc.icrr.u-tokyo.ac.jp/cgi-bin/private/DocDB/ShowDocument?docid=5938 Noise([mir,'_IM_Coil']) = Noise('HighPower'); Noise([mir,'_TM_Coil']) = Noise('HighPower'); Elec([mir,'_TM_Whitening']) = 1; % no whitening? Elec([mir,'_TM_DeWhitening']) = 1; % no dewhitening? Elec([mir,'_IM_Whitening']) = 1; % no whitening? Elec([mir,'_IM_DeWhitening']) = 1; % no dewhitening? end % Filt('BS_ISCINF') = LvFilt_BS_ISCINF.gain; % TO BE REPLACED WITH LIVE PARTS % Filt('ETMX_ISCINF') = LvFilt_ETMX_ISCINF.gain; % TO BE REPLACED WITH LIVE PARTS % Filt('ETMY_ISCINF') = LvFilt_ETMY_ISCINF.gain; % TO BE REPLACED WITH LIVE PARTS Filt('BS_ISCINF') = 1; Filt('ETMX_ISCINF') = 1; Filt('ETMY_ISCINF') = 1; Filt('BS_IM_LOCK') = myzpk('zpk',[],[0],0.0016); % somehow live filter for this does not work % TM act to TM [m/N] see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1340 %% TO BE UPDATED FOR bKAGRA PHASE 1 Susp('BS_IM_TM') = -1*myzpk('zpk',[0.4487 20;0.6982 100;1.067 100;1.6 100],[0.4246 20;0.5495 100;0.7244 100;1.009 100;1.282 100;1.629 100],0.000744); % fitted TF Susp('BS_TM_TM') = -1*myzpk('zpk',[],[0.6918 1000],0.0028); % fitted TF Susp('ETMX_IM_TM') = myzpk('zpk',[0.7259 100;0.7428 100;0.9795 100;1.518 100;2.072 10;2.195 10],[0.6177 100;0.7343 100;0.7514 100;0.9572 100;1.416 100;1.626 100;2.12 100;2.52 100],0.0011); % fitted Tf Susp('ETMX_TM_TM') = myzpk('zpk',[0.7176 100;0.7428 100;0.9682 100;1.277 100;1.589 100;2.145 100;2.491 100],[0.6177 100;0.7343 100;0.7514 100;0.9572 100;1.416 100;1.626 100;2.12 100;2.52 100],0.0023); % fitted TF Susp('ETMY_IM_TM') = Susp('ETMX_IM_TM'); Susp('ETMY_TM_TM') = Susp('ETMX_TM_TM'); noise = load('./Noises/adc_noise.txt'); Noise('DAC') = interp1(noise(:,1), noise(:,2), freq, 'linear')*DAC_C2V/2; % DAC noise [V/rtHz]; noise level is OK, but fake (see http://gwwiki.icrr.u-tokyo.ac.jp/JGWwiki/CLIO/Tasks/DigitalControl/Caltech_setup?action=AttachFile&do=get&target=analog_system_investigation.pdf) % seismic motion to mirror motion [m/m] dataH = load('../Suspensions/suspensionTFs/BS_seismic_TF.dat'); dataV = load('../Suspensions/suspensionTFs/BS_seism_vertical.dat'); gg = interp1(dataH(:,1), dataH(:,2), freq, 'linear');gg(601:end)=gg(600)*freq(600)^10./freq(601:end).^10; % extrapolation ph = interp1(dataH(:,1), dataH(:,3), freq, 'linear');ph(601:end)=-180; % extrapolation Susp('BS_Seism_TM') = frd(gg.*exp(i*ph/180*pi),freq,'FrequencyUnit','Hz'); % Seismic noise to TM [m/m] gg = interp1(dataV(:,1), dataV(:,2), freq, 'linear');gg(601:end)=gg(600)*freq(600)^10./freq(601:end).^10; % extrapolation ph = interp1(dataV(:,1), dataV(:,3), freq, 'linear');ph(601:end)=-180; % extrapolation Susp('BS_Seism_Vert') = frd(gg.*exp(i*ph/180*pi),freq,'FrequencyUnit','Hz'); % Seismic noise to TM from vertical coupling [m/m] dataH = load('../Suspensions/suspensionTFs/TypeA_seismic_TF.dat'); dataV = load('../Suspensions/suspensionTFs/TypeA_seismic_vertical.dat'); ggH = interp1(dataH(:,1), dataH(:,2), freq, 'linear');ggH(601:end)=ggH(600)*freq(600)^10./freq(601:end).^10; % extrapolation phH = interp1(dataH(:,1), dataH(:,3), freq, 'linear');phH(601:end)=-180; % extrapolation ggV = interp1(dataV(:,1), dataV(:,2), freq, 'linear');ggV(601:end)=ggV(600)*freq(600)^10./freq(601:end).^10; % extrapolation phV = interp1(dataV(:,1), dataV(:,3), freq, 'linear');phV(601:end)=-180; % extrapolation for mir = {'ETMX','ETMY'} mir = char(mir); Susp([mir,'_Seism_TM']) = frd(ggH.*exp(i*phH/180*pi),freq,'FrequencyUnit','Hz'); % Seismic noise to TM [m/m] Susp([mir,'_Seism_Vert']) = frd(ggV.*exp(i*phV/180*pi),freq,'FrequencyUnit','Hz'); % Seismic noise to TM from vertical coupling [m/m] end VertCoup = 0.01; % vertical coupling % seismic vibration data from seismometer [data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-IY0_SEIS_NS_SENSINF_OUT_DQ'); % in unit of [m/s] %data = importdata('./Noises/GPS1144857600/K1:PEM-BS_SEIS_NS_SENSINF_OUT_DQ.txt'); [noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density noise = [noisef sqrt(noisep)./(2*pi*noisef)]; Noise('BS_SEIS') = interp1(noise(:,1), noise(:,2), freq, 'linear'); %[data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-EX_SEIS_WE_SENSINF_OUT_DQ'); % in unit of [m/s] [data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-EX1_SEIS_WE_SENSINF_OUT_DQ'); % in unit of [m/s] %data = importdata('./Noises/GPS1144857600/K1:PEM-BS_SEIS_WE_SENSINF_OUT_DQ.txt'); [noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density noise = [noisef sqrt(noisep)./(2*pi*noisef)]; Noise('ETMX_SEIS') = interp1(noise(:,1), noise(:,2), freq, 'linear'); [data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-EY1_SEIS_NS_SENSINF_OUT_DQ'); % in unit of [m/s] %data = importdata('./Noises/GPS1144857600/K1:PEM-BS_SEIS_NS_SENSINF_OUT_DQ.txt'); [noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density noise = [noisef sqrt(noisep)./(2*pi*noisef)]; Noise('ETMY_SEIS') = interp1(noise(:,1), noise(:,2), freq, 'linear'); % thermal noise noise = load('./Noises/spectrum_DRSE.txt'); % Latest estimated sensitivity Aug2017 Noise('ETMY_THERM') = interp1(noise(:,1), Larm*sqrt(noise(:,3).^2+noise(:,4).^2), freq, 'linear'); % 22 K thermal noise hsusp = load('./Noises/300K/hsusp.dat'); % 300K suspension thermal noise hmirror = load('./Noises/300K/hmirror.dat'); % 300K suspension thermal noise Noise('ETMX_THERM') = interp1(hsusp(:,1), Larm*sqrt(hsusp(:,2).^2+hmirror(:,2).^2), freq, 'linear'); % 300 K thermal noise % LSC filters and matrices FiltGain = 50; % filter gain Filt('LSC_MICH') = FiltGain*myzpk('zpk',[70;0.2],[2000;2000;0.02],10); % WILL BE REPLACED WITH LIVE PARTS % AA and AI filters % Anaglog AA/AI: 3rd-order Butterworth at 10kHz, notch at 65535Hz (LIGO-T070038, LIGO-D070081) % Digital AA/AI: /opt/rtcds/rtscore/release/src/fe/controller.c [z,p,k]=butter(3,2*pi*1e4,'low','s'); Elec('Digital_AA') = zpk(z,p,k); Elec('Analog_AA') = zpk(z,p,k); Elec('Digital_AI') = zpk(z,p,k); Elec('Analog_AI') = zpk(z,p,k); % Whitening and Dewhitening filters Elec('REFL_DC_Whitening') = 1; % no whitening for iKAGRA Elec('REFL_DC_DeWhitening') = 1; % no dewhitening for iKAGRA Elec('REFL_RF_Whitening') = myzpk('zpk',[1],[10],WhiteningGain); % added whitening on Apr 25, 2018 Elec('REFL_RF_DeWhitening') = myzpk('zpk',[10],[1],1); % added whitening on Apr 25, 2018 % Measured noise at K1:LSC-MICH1_OUT_DQ [data,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-MICH1_OUT_DQ'); % use online data [noisep,noisef] = pwelch(data,hanning(ndata),[],ndata,freqsamp); % calculate power spectral density noise = [noisef sqrt(noisep)]; MeasuredNoise = interp1(noise(:,1), noise(:,2), freq, 'linear'); % UGF servo UGFSERVO_GAIN = 1; % DOF Switches loopNames = {'total','BSTM','BSIM'}; SW = ones(1,length(loopNames)); % turn on the loops swtot = 1; % switches for the total loop swact = 2:3; % switches for each actuator loop %% Plot actuator transfer function % set(0, 'DefaultAxesFontSize',10); % figure(1001) % plotbode(freq,[1e6*DAC_C2V*Coil('BS_IM')*Susp('BS_IM_TM'),1e6*DAC_C2V*Coil('BS_TM')*Susp('BS_TM_TM')]); % % % From Shoda code % addpath('/kagra/Dropbox/Subsystems/VIS/Scripts/SuspensionControlModel/utility/'); % addpath('/kagra/Dropbox/Subsystems/VIS/Scripts/SuspensionControlModel/script/'); % sysc = VISss('BS', 'safe'); % [mag1,phs1]=bodesus(sysc, 'actLIM', 'LTM', freq); % [mag2,phs2]=bodesus(sysc, 'actLTM', 'LTM', freq); % % subplot(2,1,1) % loglog(freq,mag1,'--'); % loglog(freq,mag2,'--'); % legend({'BS IM','BS TM','BS IM Shoda','BS TM Shoda'}) % ylim([1e-10,1]) % set(gca,'yTick',logspace(-10,0,11)); % ylabel('Gain [um/count]') % set(gca,'XTick',logspace(log10(freq(1)),log10(freq(end)),log10(freq(end))-log10(freq(1))+1)); % subplot(2,1,2) % loglog(freq,phs1,':'); % loglog(freq,phs2,':'); % set(gca,'XTick',logspace(log10(freq(1)),log10(freq(end)),log10(freq(end))-log10(freq(1))+1)); % saveas(gcf,[figdir,'BSActuator.png']) %% Plot nominal openloop transfer function meas=load('./TransferFunctions/20180424/MICHOLTF20180424.txt'); % measured OLTF data meas(:,3)=angle(-exp(i*meas(:,3)/180*pi))/pi*180; OLTF=[]; % OLTF for total loop SW = zeros(1,length(loopNames)); SW(swact) = ones(1,length(swact)); [A,B,C,D]=linmod(noiseModel); systm=ss(A,B,C,D); OLTF=[OLTF,systm(swtot,swtot,:)]; % OLTF for each actuator loop SW = zeros(1,length(loopNames)); SW(swtot) = ones(1,length(swtot)); [A,B,C,D]=linmod(noiseModel); systm=ss(A,B,C,D); for kk=swact OLTF=[OLTF,systm(kk,kk,:)]; end set(0, 'DefaultAxesFontSize',10); figure(101) plotbode(freq,OLTF); subplot(2,1,1) loglog(meas(:,1),meas(:,2),'.'); loglog([freq(1),freq(end)],[1,1],'k') ylim([1e-6,1e6]); xlim([freq(1),freq(end)]); set(gca,'YTick',logspace(-6,6,7)); set(gca,'XTick',logspace(log10(freq(1)),log10(freq(end)),log10(freq(end))-log10(freq(1))+1)); title(sprintf('MICH OLTF (uncalibrated, optical gain: %10.2e counts/m)',abs(MICHOpticalGainRF)*ADC_V2C*PdResp('REFL_RF')*TransImp('REFL_RF')*DeModGain)); subplot(2,1,2) loglog(meas(:,1),meas(:,3),'.'); xlim([freq(1),freq(end)]); set(gca,'XTick',logspace(log10(freq(1)),log10(freq(end)),log10(freq(end))-log10(freq(1))+1)); legend(loopNames{swtot},loopNames{swact},'Measured 20180424 w/o IM loop'); text(1.2e4,1-180,['GPStime: ',num2str(GPStime)],'Rotation',90); saveas(gcf,[figdir,'MICHOLTF_uncalibrated.png']) SW = ones(1,length(loopNames)); % turn on the loops %% Compute noises and save cache % Compute noises [noises, sys] = nbFromSimulink(noiseModel, freq, 'dof', 'MICH'); % save cached output saveFunctionCache(); %% Plot noise budget set(0, 'DefaultAxesFontSize',10); disp('Plotting noises') nb = nbGroupNoises(noiseModel, noises, sys); nb.sortModel(); matlabNoisePlot(nb); figure(2) ylim([1e-20,1e-9]); xlim([1e-1,1e4]); set(gca,'YTick',logspace(-20,-9,12)); set(gca,'XTick',logspace(-1,4,6)); text(1.2e4,1e-20,['GPStime: ',num2str(GPStime)],'Rotation',90); saveas(gcf,[figdir,'MICHNoiseBudgetActuator.png']) figure(3) ylim([1e-20,1e-9]); xlim([1e-1,1e4]); set(gca,'YTick',logspace(-20,-9,12)); set(gca,'XTick',logspace(-1,4,6)); text(1.2e4,1e-20,['GPStime: ',num2str(GPStime)],'Rotation',90); saveas(gcf,[figdir,'MICHNoiseBudgetSeismic.png']) figure(4) ylim([1e-20,1e-9]); xlim([1e-1,1e4]); set(gca,'YTick',logspace(-20,-9,12)); set(gca,'XTick',logspace(-1,4,6)); text(1.2e4,1e-20,['GPStime: ',num2str(GPStime)],'Rotation',90); saveas(gcf,[figdir,'MICHNoiseBudgetSensor.png']) figure(1) ylim([1e-20,1e-9]); xlim([1e-1,1e4]); set(gca,'YTick',logspace(-20,-9,12)); set(gca,'XTick',logspace(-1,4,6)); ikagra = importdata('./Noises/MICHNoiseBudget_iKAGRA.txt'); loglog(ikagra(:,1),ikagra(:,2)); leg = legend(gca); legend({leg.String{:},'iKAGRA'},'Interpreter','None'); text(1.2e4,1e-20,['GPStime: ',num2str(GPStime)],'Rotation',90); saveas(gcf,[figdir,'MICHNoiseBudget.png']) %% Output NB data nbdata=[freq;nb.referenceNoises{1}.noiseData.asd]; header=['%GPStime:',num2str(GPStime),'% freq, sum']; for kk=1:length(nb.modelNoises) nbdata=[nbdata;nb.modelNoises{kk}.asd]; hd=strsplit(nb.modelNoises{kk}.name,'}{') try hdd=char(hd(2)); hdd=hdd(1:end-1); catch hdd=hd; end header=[header,', ',hdd] end dlmwrite([figdir,'MICHNoiseBudget.txt'],header,'delimiter',''); dlmwrite([figdir,'MICHNoiseBudget.txt'],nbdata','-append','delimiter',' ','precision','%.6e');