%% ABOUT THIS FILE % ---------------------------------------------------------------------- % Type-B1 prototype for KAGRA % Coded by T. Sekiguchi on 2015/06/16 % ---------------------------------------------------------------------- %% PRELIMINARY clear all; % Clear workspace close all; % Close plot windows addpath('../../utility'); % Add path to utilities addpath('../../utility/altmany-export_fig-524e556'); g = 9.81; %% IMPORT SUSPENSION MODEL matfile='typeBpsus_BRmdl'; load([matfile,'.mat']); %% TUNING DAMPER % This part compensates the failure in converting the structural damping to % viscous damping. % REDUCE DAMPING ON RIM sys1.a(28,28)=sys1.a(28,28)/30; % RIM sys1.a(34,34)=sys1.a(34,34)/30; % RRM sys1.a(40,40)=sys1.a(40,40)/30; % RTM % INCREASE DAMPING ON YF0 %sys1.a(15,15)=sys1.a(15,15)*3000; % YF0 % INCREASE DAMPING ON YRM, YTM sys1.a(36,36)=sys1.a(36,36)*30; % YRM sys1.a(42,42)=sys1.a(42,42)*300; % YTM %% ZERO CONTROL CASE addpath('servofilter'); % Add path to servo typeBptest_no_control_BR; % NO CONTROL rmpath ('servofilter'); % Remove path to servo mdlfile='typeBpsimctrl_BR'; % typeBp ver.150717 % SETTING FILTERS 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 freq =logspace(-2,1.99,1001); %% NOISE % 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] n_OSEM_LF2 = noise_OSEM/sqrt(2); n_OSEM_TF2 = noise_OSEM/sqrt(2); n_OSEM_YF2 = noise_OSEM/sqrt(2)/300e-3; 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_LF0 = noise_LVDT/sqrt(3/2); % [m/rtHz] %n_LVDT_TF0 = noise_LVDT/sqrt(3/2); % [m/rtHz] %n_LVDT_YF0 = noise_LVDT/sqrt(3)/600e-3; % [rad/rtHz] %n_LVDT_VF0 = noise_LVDT*2; % [m/rtHz] n_LVDT_VF1 = noise_LVDT*2; % [m/rtHz] n_LVDT_VF2 = noise_LVDT*2; % [m/rtHz] n_Oplev = zeros(size(freq))+1.e-7; %[rad/rtHz] (from Michimura's slide JGW-T1202403-v1) % 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_LF0 = noise_GEO/sqrt(3/2); % [m/rtHz] %n_GEO_TF0 = noise_GEO/sqrt(3/2); % [m/rtHz] %n_GEO_YF0 = noise_GEO/sqrt(3)/594e-3; % [m/rtHz] 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 typeBptest_servo_lock_BR; % DAMPING MODE rmpath ('servofilter'); % Remove path to servo mdlfile='typeBpsimctrl_BR'; % typeB1 ver.150616 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 %% CALC GAS for LENGTH [magg1,~]=bodesus(sysc,'n_LVDT_VF1','IFO_LTM',freq); [magg2,~]=bodesus(sysc,'n_LVDT_VF2','IFO_LTM',freq); magg1=magg1.*n_LVDT_VF1; magg2=magg2.*n_LVDT_VF2; magGasTot=sumpsd({magg1,magg2}); %GAS LVDT noise %% CALC ABOUT OSEM for LENGTH [magfL,~]=bodesus(sysc,'n_OSEM_LF2','IFO_LTM',freq); [magfT,~]=bodesus(sysc,'n_OSEM_TF2','IFO_LTM',freq); [magfY,~]=bodesus(sysc,'n_OSEM_YF2','IFO_LTM',freq); magfL=magfL.*n_OSEM_LF2; magfT=magfT.*n_OSEM_TF2; magfY=magfY.*n_OSEM_YF2; [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; [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({magfL,magfT,magfY,magoL,magoP,magoY,magoT,magoR,magoV,magtL,magtP,magtY}); %OSEM noise %% CALC ABOUT OPLEV for LENGTH [magolP,~]=bodesus(sysc,'n_OpLev_PTM','IFO_LTM',freq); [magolY,~]=bodesus(sysc,'n_OpLev_YTM','IFO_LTM',freq); magolP=magolP.*n_Oplev; %Assume Oplev noise level is same as OSEM magolY=magolY.*n_Oplev; magOlTot=sumpsd({magolP, magolY}); %% ALL Noise magnTot=sumpsd({magGasTot, magOsemTot, magOlTot}); %ALL Noise %magnTot=sumpsd({magGasTot, magOsemTot}); %ALL Noise rmsnTot=makerms(freq,magnTot); rmsnTotv=makerms(freq,magnTot.*freq*pp); %% ALL magAllT=sumpsd({magSeisT,magnTot}); % magaT=magoTot; 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 %% PLOT % mypsdplotopt({magoTot,maggT,magsLG,magaT,mag0T,reqBS},freq,... % 'title','Control noise coupling during observation phase',... % 'legend',{'OSEM noise','GAS control noise','top control noise'... % 'With control total','No control',... % 'BS requirement'},... % 'color',{'m-','c-','g-','r-','b-','k--'},... % 'ylim',[1e-18,1e-5],'ylabel','Magnitude [m/rtHz] or [m/sec]') mypsdplotopt({magOsemTot,magGasTot,magOlTot,magAllT,mag0T,reqPR},freq,... 'title','Control noise coupling during observation phase',... 'legend',{'OSEM noise','GAS control noise','Oplev control noise',... 'With control total','No control', 'PR requirement'},... 'color',{'m-','g-','y-','r-','b-','k--'},... 'ylabel','Magnitude [m/rtHz] or [m/sec]'); filename = input('disp filename? > '); export_fig(strcat('./figure/',filename)); %% 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 OSEM for Pitch [magpfL,~]=bodesus(sysc,'n_OSEM_LF2','PTM',freq); [magpfT,~]=bodesus(sysc,'n_OSEM_TF2','PTM',freq); [magpfY,~]=bodesus(sysc,'n_OSEM_YF2','PTM',freq); magpfL=magpfL.*n_OSEM_LF2; magpfT=magpfT.*n_OSEM_TF2; magpfY=magpfY.*n_OSEM_YF2; [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; [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({magpfL, magpfT, magpfY, magpoL,magpoP,magpoY,magpoT,magpoR,magpoV,magptL,magptP,magptY}); %% CALC Oplev for Pitch [magpolP,~]=bodesus(sysc,'n_OpLev_PTM','PTM',freq); [magpolY,~]=bodesus(sysc,'n_OpLev_YTM','PTM',freq); magpolP=magpolP.*n_Oplev; %Assume that Oplev noise level is the same as OSEM magpolY=magpolY.*n_Oplev; magpOlTot=sumpsd({magpolP,magpolY}); %% ALL magpAllT=sumpsd({magpS,magpOsemTot,magpOlTot}); rmspAllT =makerms(freq, magpAllT); rmspAllTv=makerms(freq, magpAllT.*freq*pp); %% PLOT PITCH mypsdplotopt({mag0pS,magpOsemTot,magpOlTot,magpAllT,rms0pT,rmspAllT},freq,... 'title','Angular Fluctuation of the Mirror',... 'legend',{'No control',... 'OSEM noise','Oplev noise','total with control','no control RMS',... 'with control RMS'},... 'color',{'b-','g-','y-','r-','c--','m--'},... 'ylim',[1e-12,1e-4],'ylabel','Magnitude [rad/rtHz] or [rad]'); filename = input('pitch filename? > '); export_fig(strcat('./figure/',filename)); %% PLOT Velocity magAllTv = magAllT.*freq*pp; mag0Tv = mag0T.*freq*pp; %rms0Tv = makerms(freq, mag0Tv); mypsdplotopt({mag0Tv,magOsemTot.*freq*pp,magGasTot.*freq*pp,magOlTot.*freq*pp,magAllTv,rms0Tv,rmsAllTv},freq,... 'title','Velocity of the Mirror',... 'legend',{'No control','OSEM noise', 'GAS noise','Oplev noise',... 'With control total',... 'No control RMS','With control RMS'},... 'color',{'b-','k-','g-','y-','r-','c--','m--'},... 'ylim',[1e-12,1e-4],'ylabel','Magnitude [m/s/rtHz] or [m/s]'); filename = input('vel filename ? >'); export_fig(strcat('./figure/',filename)); %% % [mag0pS,~]=bodesus(sysc0,'accLGND','OpLev_PTM',freq); % mag0pS=mag0pS.*freq.*freq*pp*pp.*noise_seism; % rms0pT=makerms(freq,mag0pS); % mypsdplotopt({mag0pS,rms0pT},freq,... % 'title','Noise coupling from OSEMs',... % 'legend',{... % 'Suspended mirror w/o Ctrl','Suspended mirror RMS'},... % 'color',{'r-','m--'},... % 'ylim',[1e-12,1e-4],'ylabel','Angular Fluctuation [rad/rtHz]') %% Force on F2 [magLf2,~]=bodesus(sysc,'accLGND','act_LF2_mon',freq); magLf2 = magLf2.*freq.*freq.*pp.*pp.*noise_seism; rmsLf2 = makerms(freq, magLf2); mypsdplotopt({magLf2, rmsLf2},freq, 'title', 'Feedback Force on BF',... 'legend',{'Feedback force', 'RMS'}, 'color',{'b','c'},... 'ylim', [1e-10,1e-2],'ylabel','Force [N/rtHz] or [N]'); filename=input('FB force filename ? >'); export_fig(strcat('./figure/',filename));