%% ABOUT THIS FILE % ---------------------------------------------------------------------- % Type-B SAS for KAGRA % Coded by Y. Fujii on July 21th, 2018. % ---------------------------------------------------------------------- %% PRELIMINARY clear all; % Clear workspace close all; % Close plot windows addpath('../../utility'); % Add path to utilities % addpath('/usr/local/matlab/altmany-export_fig-76bd7fa'); g = 9.81; %% IMPORT SUSPENSION MODEL matfile='SR3_180502'; load([matfile,'mdl.mat']); %% TUNING DAMPER % This part compensates the failure in converting the structural damping to % viscous damping. % REDUCE DAMPING ON IM sys1.a(37,37)=sys1.a(37,37)/3.e3; % RIR sys1.a(43,43)=sys1.a(43,43)/3.e3; % RIM %sys1.a(40,40)=sys1.a(40,40)/30; % RRM % INCREASE DAMPING ON F0 sys1.a(13,13)=sys1.a(13,13)/3.e3; % LF0 sys1.a(14,14)=sys1.a(14,14)/3.e3; % TF0 sys1.a(15,15)=sys1.a(15,15)*300; % YF0 % INCREASE DAMPING ON MD %sys1.a(16,16)=sys1.a(16,16)/5000; % LMD %sys1.a(17,17)=sys1.a(17,17)/5000; % TMD % INCREASE DAMPING ON RM, TM %sys1.a(42,42)=sys1.a(42,42)*30; % YRM %sys1.a(48,48)=sys1.a(48,48)*300; % YTM % INCREASE DAMPING ON VTM %sys1.a(45,45)=sys1.a(45,45)*100; % VTM %% IMPORT SERVO FILTERS addpath('servofilter'); % Add path to servo DefaultParameters; % NO CONTROL MODE rmpath ('servofilter'); % Remove path to servo %% IMPORT SIMULINK MODEL mdlfile='typeBsimctrl'; % typeB_bKAGRA ver. 20180502 st =linmod(mdlfile); invl =strrep(st.InputName, [mdlfile,'/'],''); outvl =strrep(st.OutputName,[mdlfile,'/'],''); sysc0 =ss(st.a,st.b,st.c,st.d,'inputname',invl,'outputname',outvl); %% FREQUENCY freq1=logspace(-2,3,1001); freq=logspace(-2,2,1001); %% IMPORT NOISE ---------------------------------------------------------------------- data_OSEM = importdata('../../noise/OSEMnoiseworstproto_disp.dat'); % um/rtHz noise_OSEM = interp1(data_OSEM(:,1),data_OSEM(:,2)*1e-6,freq')'; %m/rtHz n_OSEM_LTM = noise_OSEM/sqrt(4); % [m/rtHz] n_OSEM_PTM = noise_OSEM/sqrt(2)/117e-3; % [rad/rtHz] n_OSEM_YTM = noise_OSEM/sqrt(2)/117e-3; % [rad/rtHz] n_OSEM_LIM = noise_OSEM/sqrt(1); % [m/rtHz] n_OSEM_TIM = noise_OSEM/sqrt(2); % [m/rtHz] n_OSEM_VIM = noise_OSEM/sqrt(3); % [m/rtHz] n_OSEM_RIM = noise_OSEM/sqrt(3/2)/96e-3;% [rad/rtHz] n_OSEM_PIM = noise_OSEM/sqrt(3/2)/96e-3;% [rad/rtHz] n_OSEM_YIM = noise_OSEM/sqrt(2)/87e-3; % [rad/rtHz] data_LVDT = importdata('../../noise/LVDTnoiseADC_disp.dat'); % um/rtHz noise_LVDT = interp1(data_LVDT(:,1),data_LVDT(:,2)*1e-6,freq')'; % m/rtHz n_LVDT_LIP = noise_LVDT/sqrt(3/2); % [m/rtHz] n_LVDT_TIP = noise_LVDT/sqrt(3/2); % [m/rtHz] n_LVDT_YIP = noise_LVDT/sqrt(3)/600e-3; % [rad/rtHz] n_LVDT_GASF0 = noise_LVDT*2; % [m/rtHz] n_LVDT_GASF1 = noise_LVDT*2; % [m/rtHz] n_LVDT_GASBF = noise_LVDT*2; % [m/rtHz] data_GEO = importdata('../../noise/GEOnoiseproto_vel.dat'); % um/sec/rtHz noise_GEO = interp1(data_GEO(:,1),data_GEO(:,2)*1e-6,freq')';% m/sec/rtHz n_GEO_LIP = noise_GEO/sqrt(3); % [m/rtHz] n_GEO_TIP = noise_GEO/sqrt(2); % [m/rtHz] n_GEO_YIP = noise_GEO/sqrt(3)/594e-3; % [rad/rtHz] %data_F7 = importdata('../../noise/F7LVDTprotonoise.dat'); % m/rtHz %noise_F7 = interp1(data_F7(:,1),data_F7(:,2),freq')'; % m/rtHz %n_LVDT_LBF = noise_F7/sqrt(2); %n_LVDT_TBF = noise_F7/sqrt(2); %n_LVDT_YBF = noise_F7/sqrt(2)/440e-3; %n_LVDT_RBF = noise_F7/sqrt(3/2)/440e-3; %n_LVDT_PBF = noise_F7/sqrt(3/2)/440e-3; %n_LVDT_VBF = noise_F7/sqrt(3); n_Oplev = zeros(size(freq))+1.e-7; %[rad/rtHz] (from Michimura's slide JGW-T1202403-v1) data_seism = importdata('../../noise/KamiokaSeismicHighNoise.dat'); % m/rtHz noise_seism = interp1(data_seism(:,1),data_seism(:,2),freq')';% m//rtHz %% NOISE MODEL ---------------------------------------------------------------------- % mypsdplotopt({noise_OSEM,noise_seism},freq,... % 'title','Noise Model',... % 'legend',{'OSEM','seismic'},... % 'color',{'r-','k-'},... % 'ylim',[1e-13,1e-5],'ylabel','Magnitude [m/rtHz] ') %% MODEL ---------------------------------------------------------------------------- addpath('servofilter'); % Add path to servo servo_damp2; % Damping MODE rmpath ('servofilter'); % Remove path to servo %% IMPORT SIMULINK MODEL ------------------------------------------------------------- mdlfile='typeBsimctrl'; st =linmod(mdlfile); invl =strrep(st.InputName, [mdlfile,'/'],''); outvl =strrep(st.OutputName,[mdlfile,'/'],''); sysc =ss(st.a,st.b,st.c,st.d,'inputname',invl,'outputname',outvl); %% CALC INITIAL ---------------------------------------------------------------------- [mag0Vg,~]=bodesus(sysc0,'accVGND','IFO_LTM',freq); [mag0S,~]=bodesus(sysc0,'accLGND','IFO_LTM',freq); mag0Vg=mag0Vg.*freq.*freq*pp*pp.*noise_seism; mag0S = mag0S.*freq.*freq*pp*pp.*noise_seism; mag0T=sumpsd({mag0Vg,mag0S}); % Mirror Displacement w/o Ctrl rms0T=makerms(freq,mag0T); % Mirror Displacement RMS w/o Ctrl rms0Tv=makerms(freq,mag0T.*freq*pp); % Mirror Velocity RMS w/o Ctrl %% CALC SEISM for LENGTH ------------------------------------------------------------- [magVg,~]=bodesus(sysc,'accVGND','IFO_LTM',freq); [magsS,~]=bodesus(sysc,'accLGND','IFO_LTM',freq); magVg=magVg.*freq.*freq*pp*pp.*noise_seism; magsS=magsS.*freq.*freq*pp*pp.*noise_seism; magSeisT=sumpsd({magsS,magVg}); % Seismic noise of mirror rmsSeisT=makerms(freq,magSeisT); % Seismic noise RMS of mirror rmsSeisTv=makerms(freq,magSeisT.*freq*pp); % Seismic noise velocity RMS of mirror %%% --- Noise to LENGTH --- %%% %% CALC GAS for LENGTH --------------------------------------------------------------- [magg0,~]=bodesus(sysc,'n_LVDT_GASF0','IFO_LTM',freq); [magg1,~]=bodesus(sysc,'n_LVDT_GASF1','IFO_LTM',freq); [magg2,~]=bodesus(sysc,'n_LVDT_GASBF','IFO_LTM',freq); magg0=magg0.*n_LVDT_GASF0; magg1=magg1.*n_LVDT_GASF1; magg2=magg2.*n_LVDT_GASBF; magGasTot=sumpsd({magg0, magg1, magg2}); %GAS LVDT noise to IFO_LTM %% CALC TOP TRANS -------------------------------------------------------------------- [magsL,~]=bodesus(sysc,'n_LVDT_LIP','IFO_LTM',freq); [magsG,~]=bodesus(sysc,'n_GEO_LIP','IFO_LTM',freq); magsL=magsL.*n_LVDT_LIP; magsG=magsG.*n_GEO_LIP; magsLG=sumpsd({magsL, magsG}); magsT=sumpsd({magsL, magsG, magsS, magVg}); rmsT=makerms(freq,magsT); %F0 LVDT&GEO noise to IFO_LTM rmsTv=makerms(freq,magsT.*freq*pp); %% CALC BF-LVDT TRANS for LENGTH ----------------------------------------------------------------- % [magfL,~]=bodesus(sysc,'n_LVDT_LBF','IFO_LTM',freq); % [magfT,~]=bodesus(sysc,'n_LVDT_TBF','IFO_LTM',freq); % [magfY,~]=bodesus(sysc,'n_LVDT_YBF','IFO_LTM',freq); % [magfV,~]=bodesus(sysc,'n_LVDT_VBF','IFO_LTM',freq); % [magfP,~]=bodesus(sysc,'n_LVDT_PBF','IFO_LTM',freq); % [magfR,~]=bodesus(sysc,'n_LVDT_RBF','IFO_LTM',freq); % magfL=magfL.*n_LVDT_LBF; % magfT=magfT.*n_LVDT_TBF; % magfY=magfY.*n_LVDT_YBF; % magfV=magfV.*n_LVDT_VBF; % magfP=magfP.*n_LVDT_PBF; % magfR=magfR.*n_LVDT_RBF; % magLvdtTot=sumpsd({magfL,magfT,magfY,magfV,magfP,magfR}); %F2 LVDT noise %% CALC ABOUT OSEM for LENGTH --------------------------------------------------------- [magoL,~]=bodesus(sysc,'n_OSEM_LIM','IFO_LTM',freq); [magoP,~]=bodesus(sysc,'n_OSEM_PIM','IFO_LTM',freq); [magoY,~]=bodesus(sysc,'n_OSEM_YIM','IFO_LTM',freq); [magoT,~]=bodesus(sysc,'n_OSEM_TIM','IFO_LTM',freq); [magoR,~]=bodesus(sysc,'n_OSEM_RIM','IFO_LTM',freq); [magoV,~]=bodesus(sysc,'n_OSEM_VIM','IFO_LTM',freq); magoL=magoL.*n_OSEM_LIM; magoT=magoT.*n_OSEM_TIM; magoV=magoV.*n_OSEM_VIM; magoR=magoR.*n_OSEM_RIM; magoP=magoP.*n_OSEM_PIM; magoY=magoY.*n_OSEM_YIM; magOsemTot=sumpsd({magoL,magoP,magoY,magoT,magoR,magoV}); %IM-OSEM noise to IFO_LTM % ---------------------------------------------------------------------- % [magtL,~]=bodesus(sysc,'n_OSEM_LTM','IFO_LTM',freq); % [magtP,~]=bodesus(sysc,'n_OSEM_PTM','IFO_LTM',freq); % [magtY,~]=bodesus(sysc,'n_OSEM_YTM','IFO_LTM',freq); % magtL=magtL.*n_OSEM_LTM; % magtP=magtP.*n_OSEM_PTM; % magtY=magtY.*n_OSEM_YTM; % magOsemTot=sumpsd({magoL,magoP,magoY,magoT,magoR,magoV,magtL,magtP,magtY}); %TM-OSEM noise to LTM %% CALC ABOUT OPLEV for LENGTH ---------------------------------------------------------------------- [magolL,~]=bodesus(sysc,'n_OpLev_LTM','IFO_LTM',freq); [magolP,~]=bodesus(sysc,'n_OpLev_PTM','IFO_LTM',freq); [magolY,~]=bodesus(sysc,'n_OpLev_YTM','IFO_LTM',freq); magolL=magolL.*n_Oplev; magolP=magolP.*n_Oplev; magolY=magolY.*n_Oplev; magOlTot=sumpsd({magolL, magolP, magolY}); %Oplev noise to LTM %% ALL Noise to LENGTH ---------------------------------------------------------------------- magnTot=sumpsd({magGasTot, magsT, magOsemTot, magOlTot}); %ALL Noise rmsnTot=makerms(freq,magnTot); rmsnTotv=makerms(freq,magnTot.*freq*pp); %% ALL ---------------------------------------------------------------------- %magAllT=sumpsd({magSeisT, magnTot}); magAllT=sumpsd({magnTot}); rmsAllT=makerms(freq,magAllT); rmsAllTv=makerms(freq,magAllT.*freq*pp); %%% --- Requirement --- %%% %reqPR=freq; %for i=1:length(freq) % if freq(i)<8 % reqPR(i)=0; % else % reqPR(i)=5e-15*(freq(i)/10)^(-2)+5e-16; % end %end % dataBRSE=importdata('../../SuspensionRequirement/BRSE/DisplacementNoiseRequirementCopy5Hz.dat'); % dataBRSE=[dataBRSE(:,1),dataBRSE(:,2),dataBRSE(:,3),dataBRSE(:,4),dataBRSE(:,5),dataBRSE(:,6),... % dataBRSE(:,7),dataBRSE(:,8),dataBRSE(:,9)]; % reqSR_BRSE = interp1(dataBRSE(:,1),dataBRSE(:,9),freq')';% m/rtHz dataDRSE=importdata('../../SuspensionRequirement/DRSE/DisplacementNoiseRequirementCopy5Hz.dat'); dataDRSE=[dataDRSE(:,1),dataDRSE(:,2),dataDRSE(:,3),dataDRSE(:,4),dataDRSE(:,5),dataDRSE(:,6),... dataDRSE(:,7),dataDRSE(:,8),dataDRSE(:,9)]; reqSR_DRSE = interp1(dataDRSE(:,1), dataDRSE(:,9)/2.,freq')';% m/rtHz %%% --- Preperation for PLOT --- %%% output_dir = 'figures/SensorNoiseCoupling'; %%% --- PLOT --- %%% mypsdplotopt1({magOsemTot, magGasTot, magsT, magOlTot, magAllT, mag0T, rmsAllT, rms0T, reqSR_DRSE},freq,... 'title','Control noise coupling during damping phase',... 'legend',{'OSEM noise','GAS control noise','F0 control noise','Oplev control noise',... 'With control total','No control','total with control (RMS)','No control (RMS)','requirement (DRSE)'},... 'color',{'m-','g-','c-','y-','r-','b-','r--','b--','k-','k--'},... 'ylabel','Magnitude [m/rtHz] or [m]') figname_disp = sprintf('%s/disp_damp.pdf', output_dir); export_fig(figname_disp) %%% --- CALC INITIAL for Pitch --- %%% [mag0pS,~]=bodesus(sysc0,'accLGND','PTM',freq); mag0pS=mag0pS.*freq.*freq*pp*pp.*noise_seism; % Pitch Angular Fluctuation w/o Ctrl rms0pT=makerms(freq,mag0pS); % Pitch RMS w/o Ctrl rms0pTv=makerms(freq,mag0pS.*freq*pp); %% CALC SEISM for Pitch ---------------------------------------------------------------------- [magpS,~]=bodesus(sysc,'accLGND','PTM',freq); magpS=magpS.*freq.*freq*pp*pp.*noise_seism; % %% CALC LVDT for Pitch ---------------------------------------------------------------------- % [magpfL,~]=bodesus(sysc,'n_LVDT_LBF','PTM',freq); % [magpfT,~]=bodesus(sysc,'n_LVDT_TBF','PTM',freq); % [magpfY,~]=bodesus(sysc,'n_LVDT_YBF','PTM',freq); % [magpfV,~]=bodesus(sysc,'n_LVDT_VBF','PTM',freq); % [magpfP,~]=bodesus(sysc,'n_LVDT_PBF','PTM',freq); % [magpfR,~]=bodesus(sysc,'n_LVDT_RBF','PTM',freq); % magpfL=magpfL.*n_LVDT_LBF; % magpfT=magpfT.*n_LVDT_TBF; % magpfY=magpfY.*n_LVDT_YBF; % magpfV=magpfV.*n_LVDT_VBF; % magpfP=magpfP.*n_LVDT_PBF; % magpfR=magpfR.*n_LVDT_RBF; % magpLvdtTot=sumpsd({magpfL, magpfT, magpfY, magpfV, magpfR, magpfP}); %% CALC TOP PITCH ---------------------------------------------------------------------- [magpL,~]=bodesus(sysc,'n_LVDT_LIP','PTM',freq); [magpG,~]=bodesus(sysc,'n_GEO_LIP','PTM',freq); magpL=magpL.*n_LVDT_LIP; magpG=magpG.*n_GEO_LIP; magpT=sumpsd({magpL, magpG, magpS}); rmspT=makerms(freq, magpT); rmspTv=makerms(freq, magpT.*freq*pp); %F0 LVDT&GEO noise to PITCH %% CALC OSEM for Pitch ---------------------------------------------------------------------- [magpoL,~]=bodesus(sysc,'n_OSEM_LIM','PTM',freq); [magpoP,~]=bodesus(sysc,'n_OSEM_PIM','PTM',freq); [magpoY,~]=bodesus(sysc,'n_OSEM_YIM','PTM',freq); [magpoT,~]=bodesus(sysc,'n_OSEM_TIM','PTM',freq); [magpoR,~]=bodesus(sysc,'n_OSEM_RIM','PTM',freq); [magpoV,~]=bodesus(sysc,'n_OSEM_VIM','PTM',freq); magpoL=magpoL.*n_OSEM_LIM; magpoP=magpoP.*n_OSEM_PIM; magpoY=magpoY.*n_OSEM_YIM; magpoT=magpoT.*n_OSEM_TIM; magpoR=magpoR.*n_OSEM_RIM; magpoV=magpoV.*n_OSEM_VIM; magpOsemTot=sumpsd({magpoL,magpoP,magpoY,magpoT,magpoR,magpoV}); % OSEM noise to PITCH % ---------------------------------------------------------------------- % [magptL,~]=bodesus(sysc,'n_OSEM_LTM','PTM',freq); % [magptP,~]=bodesus(sysc,'n_OSEM_PTM','PTM',freq); % [magptY,~]=bodesus(sysc,'n_OSEM_YTM','PTM',freq); % magptL=magptL.*n_OSEM_LTM; % magptP=magptP.*n_OSEM_PTM; % magptY=magptY.*n_OSEM_YTM; % magpOsemTot=sumpsd({magpoL,magpoP,magpoY,magpoT,magpoR,magpoV,magptL,magptP,magptY}); % ---------------------------------------------------------------------- %% CALC Oplev for Pitch ---------------------------------------------------------------------- [magpolL,~]=bodesus(sysc,'n_OpLev_LTM','PTM',freq); [magpolP,~]=bodesus(sysc,'n_OpLev_PTM','PTM',freq); [magpolY,~]=bodesus(sysc,'n_OpLev_YTM','PTM',freq); magpolL=magpolL.*n_Oplev; magpolP=magpolP.*n_Oplev; magpolY=magpolY.*n_Oplev; magpOlTot=sumpsd({magpolL,magpolP,magpolY}); % OpLev noise to PITCH %%% --- ALL --- %%% magpAllT=sumpsd({magpT, magpOsemTot, magpOlTot}); rmspAllT =makerms(freq, magpAllT); rmspAllTv=makerms(freq, magpAllT.*freq*pp); %%% --- PLOT PITCH --- %%% mypsdplotopt({mag0pS, magpT, magpOsemTot, magpOlTot, magpAllT, rms0pT, rmspAllT},freq,... 'title','Angular Fluctuation of the Mirror',... 'legend',{'No control','F0 control noise',... 'OSEM noise','Oplev noise','total with control','no control (RMS)',... 'with control (RMS)'},... 'color',{'b-','k-','g-','y-','r-','c--','m--'},... 'ylim',[1e-12,1e-4],'ylabel','Magnitude [rad/rtHz] or [rad]') figname_pitch = sprintf('%s/pitch_damp.pdf', output_dir); export_fig(figname_pitch) %%% --- PLOT Velocity --- %%% magAllTv = magAllT.*freq*pp; mag0Tv = mag0T.*freq*pp; mypsdplotopt({mag0Tv, magOsemTot.*freq*pp, magGasTot.*freq*pp, magsT.*freq*pp, magOlTot.*freq*pp, magAllTv, rms0Tv, rmsAllTv}, freq,... 'title','Velocity of the Mirror',... 'legend',{'No control','OSEM noise', 'GAS noise', 'F0 control noise','Oplev noise',... 'With control total',... 'No control RMS','With control RMS'},... 'color',{'b-','k-','k--','g-','y-','r-','c--','m--'},... 'ylim',[1e-12,1e-4],'ylabel','Magnitude [m/s/rtHz] or [m/s]') figname_vel = sprintf('%s/vel_damp.pdf', output_dir); export_fig(figname_vel) %%% --- Plot 1/e damping time with damping controls ON and OFF --- %%% vecpole0=pole(sysc0); vecpole0=vecpole0(real(vecpole0)<0); vecpole0=vecpole0(imag(vecpole0)>1e-6); vecpole0i=imag(vecpole0); vecpole0r=real(vecpole0); [vecpole0i,numpole0i]=sort(vecpole0i); vecpole0r=vecpole0r(numpole0i); vecpole0=[vecpole0i,-vecpole0r]; vecpolec=pole(sysc); %vecpolec=vecpolec(real(vecpolec)<0); %vecpolec=vecpolec(imag(vecpolec)>1e-6); vecpolec=vecpolec(real(vecpolec)<0); vecpolec=vecpolec(imag(vecpolec)>1e-10); vecpoleci=imag(vecpolec); vecpolecr=real(vecpolec); [vecpoleci,numpoleci]=sort(vecpoleci); vecpolecr=vecpolecr(numpoleci); vecpolec=[vecpoleci,-vecpolecr]; %% FREQUENCY ---------------------------------------------------------------------- freq1 =logspace(-2,3,1001); t_noCtrl = [vecpole0(:,1)/2./pi,1./vecpole0(:,2)]; t_Ctrl = [vecpolec(:,1)/2./pi,1./vecpolec(:,2)]; fig=figure; grid on loglog(vecpole0(:,1)/2./pi, 1./vecpole0(:,2),'b.',... vecpolec(:,1)/2./pi, 1./vecpolec(:,2),'r.',... freq1, 1000/pi./freq1,'--y',... freq1, 10/pi./freq1,'--g',... freq1, 60,'--m',... 'MarkerSize',10) grid on title('Quality factor with damping controls ON and OFF','FontSize',12, ... 'FontWeight','bold','FontName','Times New Roman',... 'interpreter','none') ylabel('1/e damping time [sec]','FontSize',12,'FontWeight','bold','FontName','Times New Roman') xlabel('Frequency [Hz]','FontSize',12,'FontWeight','bold','FontName','Times New Roman') set(gca,'FontSize',12,'FontName','Times New Roman') legend('Control OFF','Control ON','Q = 1000','Q = 10','1 min.') ylim([1e-1 5e5]) xlim([0.01,1000]) set(fig,'Color','white') figname_Q = sprintf('%s/decay_freq_damp.pdf', output_dir); export_fig(figname_Q) % ---- save point 2018.7.21 --- %