% Noise budget script for iKAGRA Michelson % plot for iKAGRA Paper % Author: 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/'])); figdir = './results/'; % directory to save figures %% Define some parameters, filters, and noises % see the Simulink file (noiseModel) for each parameter freq = logspace(-1,4,1000); noiseModel = 'Michelson_foriKAGRAPaper'; IFO = 'K1'; site = 'kamioka'; % get online data (see http://klog.icrr.u-tokyo.ac.jp/osl/index.php?r=1357) % OnlineData = GWData; % online data not available anymore (as of 2017.2.23) GPStime = 1144857600; % time to get online data 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); % 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]'; % 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 = (2991.6+2988.3)/2; % 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] 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/20160407/LSC-REFL_PDA1_RF17_Q_dark.txt'); % dark noise measured when MI misaligned 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) noise = load('./Noises/adc_noise.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) % [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 % IFO Optical responses and noises las = 3.3299; % Schnupp asymmetry [m] P0 = 250e-3; % input power to BS [W] Ppd = 8e-3; % power at REFL PD [W] ~3000 counts Freq2MICH = las/nu; % Laser frequency noise to MICH coupling [m/Hz] noise = load('./Noises/IMCOOL.dat'); 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 = 5.3e10/ADC_V2C/PdResp('REFL_RF')/TransImp('REFL_RF')/DeModGain; % measured value (see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1390) 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 for mir = {'BS','ETMX','ETMY'} mir = char(mir); Coil([mir,'_TM']) = 1e-4; % N/V see email from K. Yamamoto on Mar 22 (to be confirmed) Noise([mir,'_TM_Coil']) = 57e-9./freq+2.1e-9; % high power coil driver noise [V/rtHz]; See Figure 3 of https://dcc.ligo.org/LIGO-T080014 Elec([mir,'_TM_Whitening']) = 1; % no whitening for iKAGRA Elec([mir,'_TM_DeWhitening']) = 1; % no dewhitening for iKAGRA 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; % TM act to TM [m/N] see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1340 Susp('BS_TM_TM') = myzpk('zpk',[],[0.87 3.5],1e-2); % THIS IS FAKE Susp('ETMX_TM_TM') = myzpk('zpk',[],[0.92 5.1],5.1e-6/1e-4); Susp('ETMY_TM_TM') = myzpk('zpk',[],[0.94 4.6],3.3e-6/1e-4); 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) %Seis_TypeC = 8e-9*freq.^-2./(1+(2./freq).^-8); %Seis_TypeB = 1.5e-11*freq.^-9; %Seis_TypeBp=2.3e-10*freq.^-7; %Seis_TypeBpfixed=2e-10*freq.^-5; % seismic motion to mirror motion [m/m] Susp('BS_SEIS_TM') = myzpk('zpk',[],[0.93 5.0; 2.8 3.0],1); % rough caculation Susp('ETMX_SEIS_TM') = myzpk('zpk',[],[0.93 5.0; 2.8 3.0],1); Susp('ETMY_SEIS_TM') = myzpk('zpk',[],[0.93 5.0; 2.8 3.0],1); % seismic vibration data from seismometer % [data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-BS_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-BS_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-EY_SEIS_NS_SENSINF_OUT_DQ'); % in unit of [m/s] % [data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-BS_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'); % coupling from mirror tilt(Pitch/Yaw) fluctuation BeamMissAlignP = 1e-2; % [m] beamspot offset from the center of mirror (fitted to fit measured spectrum by using peak height) % [data,time] = OnlineData.fetch(GPStime,duration2,'K1:VIS-BS_TM_OPLEV_PIT_OUT_DQ'); % mirror tilt fluctuation (OpLev signal; [urad]) data = importdata('./Noises/GPS1144857600/K1_VIS-BS_TM_OPLEV_PIT_OUT_DQ.txt'); [noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density noise = [noisef sqrt(noisep)*1e-6]; % [rad/rtHz] Noise('BS_PIT') = interp1(noise(:,1), noise(:,2), freq, 'linear'); BeamMissAlignY = 1e-2; % [m] beamspot offset from the center of mirror (fitted to fit measured spectrum by using peak height) % [data,time] = OnlineData.fetch(GPStime,duration2,'K1:VIS-BS_TM_OPLEV_YAW_OUT_DQ'); % mirror tilt fluctuation (OpLev signal; [urad]) data = importdata('./Noises/GPS1144857600/K1_VIS-BS_TM_OPLEV_YAW_OUT_DQ.txt'); [noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density noise = [noisef sqrt(noisep)*1e-6]; % [rad/rtHz] Noise('BS_YAW') = interp1(noise(:,1), noise(:,2), freq, 'linear'); % LSC filters and matrices FiltGain = 50; % filter gain Filt('LSC_MICH') = FiltGain*myzpk('zpk',[70;0.2],[2000;2000;0.02],10); % TO BE REPLACED WITH LIVE PARTS % Filt('LSC_MICH') = LvFilt_LSC_MICH.gain * LvFilt_LSC_MICH.fm1 * LvFilt_LSC_MICH.fm2 * (-1); % factor -1 is for adjusting +/- (I don't know where this factor comes from. Please find it.) % 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') = 1; % no whitening for iKAGRA Elec('REFL_RF_DeWhitening') = 1; % no dewhitening for iKAGRA % Measured noise at K1:LSC_MICH_OUT_DQ %noise = load('./Noises/20160407/LSC-MICH_OUT_DQ.txt'); % [data,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-MICH_OUT_DQ'); % use online data data = importdata('./Noises/GPS1144857600/K1_LSC-MICH_OUT_DQ.txt'); [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 % [data,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-UGF_SERVO_OUT_DQ'); data = importdata('./Noises/GPS1144857600/K1_LSC-UGF_SERVO_OUT_DQ.txt'); UGFSERVO_GAIN = 1/mean(data); % Acoustic noise data = importdata('./Noises/20160430aco/excitation_test_spectrum.txt'); ff = data(:,1); acc = data(:,3);%m/s^2 data = importdata('./Noises/20160430aco/excitation_test_tf2MICH_OUT.txt'); ff2 = data(:,1); tf = data(:,2); tf = tf; Noise('BS_ACOUSTIC') = interp1(ff,tf.*acc,freq,'linear'); % calibration fcal = 80; % calibration line frequency UGFAdjustFactor = 0.8; %OLTCALIB_MTRX(1); % get value from Live parts %[dataAr,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-CAL_A1_REAL_OUT_DQ'); % demodulated calibration line after injection (real part) %[dataAi,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-CAL_A1_IMAG_OUT_DQ'); % demodulated calibration line after injection (imaginary part) %[dataBr,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-CAL_B1_REAL_OUT_DQ'); % demodulated calibration line before injection (real part) %[dataBi,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-CAL_B1_IMAG_OUT_DQ'); % demodulated calibration line before injection (imaginary part) dataAr = importdata('./Noises/GPS1144857600/K1_LSC-CAL_A1_REAL_OUT_DQ.txt'); dataAi = importdata('./Noises/GPS1144857600/K1_LSC-CAL_A1_IMAG_OUT_DQ.txt'); dataBr = importdata('./Noises/GPS1144857600/K1_LSC-CAL_B1_REAL_OUT_DQ.txt'); dataBi = importdata('./Noises/GPS1144857600/K1_LSC-CAL_B1_IMAG_OUT_DQ.txt'); FBAfterInj = sqrt(mean(dataAr)^2+mean(dataAi)^2); % use LSC-LKIN_CAL1_A_REAL and IMAG FBBeforeInj = sqrt(mean(dataBr)^2+mean(dataBi)^2); % use LSC-LKIN_CAL1_B_REAL and IMAG OLTF_cal = FBBeforeInj/FBAfterInj; % OpenLoop gain / calibration factor % DOF Switches MICH_SW = 1; %% Plot nominal openloop transfer function with measured one % set(0, 'DefaultAxesFontSize',10); % MICH_SW = 0; % turn off the loop to measure OLTF % [A,B,C,D]=linmod(noiseModel); % systm=ss(A,B,C,D); % OLTF=systm(1,1,:); % meas=load('./TransferFunctions/20160409/MICHOLTF-UGFservo.txt'); % meas(:,3)=angle(-exp(i*meas(:,3)/180*pi))/pi*180; % figure(100) % plotbode(freq,OLTF); % subplot(2,1,1) % loglog(meas(:,1),meas(:,2),'.'); % legend({'NB model','Measured 20160409'}) % 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)); % saveas(gcf,[figdir,'MICHOLTF.png']) % % MICH_SW = 1; % turn on the loop %% Plot actual openloop transfer function using calibration line meas=load('./TransferFunctions/20160409/MICHOLTF-UGFservo.txt'); % measured OLTF data meas(:,3)=angle(-exp(i*meas(:,3)/180*pi))/pi*180; MICH_SW = 0; % turn off the loop to measure OLTF [A,B,C,D]=linmod(noiseModel); systm=ss(A,B,C,D); OLTF=systm(1,1,:); [OLTF_uncal,~]=bode(OLTF,2*pi*fcal); % theoretical openloop transfer function gain at fcal MICHOpticalGainRF = MICHOpticalGainRF * OLTF_cal/(OLTF_uncal/UGFSERVO_GAIN); % calibrate optical gain using calibration line MICH_SW = 0; % turn off the loop to measure OLTF [A,B,C,D]=linmod(noiseModel); systm=ss(A,B,C,D); OLTF=systm(1,1,:); set(0, 'DefaultAxesFontSize',10); figure(101) plotbode(freq,OLTF); subplot(2,1,1) loglog(meas(:,1),meas(:,2),'.'); legend({'NB model','Measured 20160409'}) loglog([fcal],[1/UGFAdjustFactor],'k+') % calibration point 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 (calibrated, 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)); saveas(gcf,[figdir,'MICHOLTF_calibrated.png']) MICH_SW = 1; % turn on the loop %% Compute noises and save cache % Compute noises [noises, sys] = nbFromSimulink(noiseModel, freq, 'dof', 'MICH'); % save cached outputs saveFunctionCache(); %% Plot noise budget set(0, 'DefaultAxesFontSize',10); disp('Plotting noises') nb = nbGroupNoises(noiseModel, noises, sys); nb.sortModel(); matlabNoisePlot(nb); % figure(2) % ylim([1e-17,1e-5]); % xlim([1e-1,1e4]); % set(gca,'YTick',logspace(-17,-5,13)); % set(gca,'XTick',logspace(-1,4,6)); % saveas(gcf,[figdir,'MICHNoiseBudgetActuator_iKAGRAPaper.png']) % % figure(4) % ylim([1e-17,1e-5]); % xlim([1e-1,1e4]); % set(gca,'YTick',logspace(-17,-5,13)); % set(gca,'XTick',logspace(-1,4,6)); % saveas(gcf,[figdir,'MICHNoiseBudgetSensor_iKAGRAPaper.png']) bkagraseis=2*2.*sqrt((4.1e-17.*freq.^(-6.8)).^2+(1.9e-16.*freq.^(-8)).^2+(1.8e-18.*freq.^(-5.7)).^2)/Larm; figure(1) loglog(freq,bkagraseis,'c:','LineWidth',2); ylim([1e-18,1e-9]); xlim([1e-1,4e3]); set(gca,'YTick',logspace(-18,-9,10)); set(gca,'XTick',[1e-1,1e0,1e1,1e2,1e3]); title('') leg = legend(gca); legend({leg.String{1:end},'bKAGRA seismic noise'}); set(legend,'Box','off') saveas(gcf,[figdir,'MICHNoiseBudget_iKAGRAPaper.png']) saveas(gcf,[figdir,'MICHNoiseBudget_iKAGRAPaper.eps'],'epsc') saveas(gcf,[figdir,'MICHNoiseBudget_iKAGRAPaper.fig']) nbdata=[freq;nb.referenceNoises{1}.noiseData.asd]; for kk=1:length(nb.modelNoises) nbdata=[nbdata;nb.modelNoises{kk}.asd]; end dlmwrite([figdir,'MICHNoiseBudget.txt'],nbdata','delimiter','\t','precision',6);