%% 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 g = 9.81; %% IMPORT SUSPENSION MODEL matfile='typeB1susmdl'; load([matfile,'.mat']); %% TUNING DAMPER % This part compensates the failure in converting the structural damping to % viscous damping. % REDUCE DAMPING ON RIM sys1.a(43,43)=sys1.a(43,43)/30; % RIM sys1.a(49,49)=sys1.a(49,49)/30; % RRM sys1.a(55,55)=sys1.a(55,55)/30; % RTM % INCREASE DAMPING ON YF0 sys1.a(15,15)=sys1.a(15,15)*3000; % YF0 % INCREASE DAMPING ON YRM, YTM sys1.a(51,51)=sys1.a(51,51)*30; % YRM sys1.a(57,57)=sys1.a(57,57)*300; % YTM %% FREQUENCY freq =logspace(-2,1,1001); %% IMPORT SERVO FILTERS addpath('servofilter'); % Add path to servo typeB1proto_no_control_150618; % NO CONTROL rmpath ('servofilter'); % Remove path to servo mdlfile='typeB1simctrl_150616'; % typeB1 ver.150616 % SETTING FILTERS dmpflt0=myzpk(0,[1e1,1e1],(1e1*pp)^2); servo_LF0=dmpflt0*3; 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); %% PLOT bodesusplotopt(sysc0,'actLF0','LVDT_LF0',freq,... 'calibration',gain_act_LF0*servo_LF0,... 'title','Open loop gain') %% DAMPING GAIN CHECK addpath('servofilter'); % Add path to servo typeB1proto_servo_lock_150618; % LOCK ACQUISITION MODE rmpath ('servofilter'); % Remove path to servo fltgain={0.1,0.3,1.0,3.0,10.,30,100}; nlp=length(fltgain); lgdname=cell(1,nlp); for i=1:nlp lgdname{i}=['gain=',num2str(fltgain{i})]; end % ONLY LVDT USED blend_LVDT=1; blend_GEO =0; % TF LIST SET mags1=cell(1,nlp); times1=cell(1,nlp); tml=1:0.1:100; mdlfile='typeB1simctrl_150616'; % typeB1 ver.150616 % LOOP START for i=1:nlp servo_LF0=dmpflt0*fltgain{i}; 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); [mags1{i},~]=bodesus(sysc,'accLGND','IFO_LTM',freq); mags1{i}=mags1{i}.*freq.*freq*pp*pp; nin1=strcmp(sysc.InputName,'actLF0'); nout1=strcmp(sysc.OutputName,'IFO_LTM'); times1{i}=impulse(sysc(nout1,nin1),tml); end %% PLOT GAIN mypsdplotopt(mags1,freq,... 'title','Amplitude of transfer functions with different gain',... 'legend',lgdname,... 'ylim',[1e-3,1e2]) export_fig('figure/typeB1proto_IPdamp_gain_150618.pdf') %% IMPULSE RESPONSE fig0=figure; times2=cell(1,nlp*2); for i=1:nlp; times2{2*i-1}=tml;times2{2*i}=times1{i}*(fltgain{i}+10)+i*2e-2; end plot(times2{:},'LineWidth',2) ylabel('Amplitude','FontSize',12,'FontWeight','bold','FontName','Times New Roman') xlabel('Time [s]','FontSize',12,'FontWeight','bold','FontName','Times New Roman') set(gca,'FontSize',12,'FontName','Times New Roman') legend(lgdname) title('Impulse Response',... 'FontSize',12,'FontWeight','bold','FontName','Times New Roman') set(fig0,'Position',[50, 50, 850, 650]); export_fig('figure/typeB1proto_IPdamp_impulse_150618.pdf') %% BLENDING % SETTING FILTERS addpath('servofilter'); % Add path to servo typeB1proto_servo_lock_150618; % LOCK ACQUISITION MODE rmpath ('servofilter'); % Remove path to servo dmpflt0=myzpk(0,[1e1,1e1],(1e1*pp)^2); servo_LF0=dmpflt0*3; georesp = zpk([-2.13+1i*5.19;-2.13-1i*5.19],[0;0],1); vel2disp = zpk([],0,1); n_blend = 7; % 7th order blending nbd = (n_blend+1)/2; cf_poly = zeros(1,n_blend+1); fblendlist={0.02,0.05,0.1,0.2,0.5,1.0}; nlp=length(fblendlist); lgdname=cell(1,nlp+1); for i=1:nlp lgdname{i}=['f(blend)=',num2str(fblendlist{i})]; end lgdname{nlp+1}='NO CONTROL'; mdlfile='typeB1simctrl_150616'; % typeB1 ver.150616 % NOISE 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); % [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_LF0 = noise_GEO/sqrt(3); % [m/rtHz] data_seism = importdata('../../noise/KamiokaSeismicHighNoise.dat'); % m/rtHz noise_seism = interp1(data_seism(:,1),data_seism(:,2),freq')';% m//rtHz % MAGNITUDE SET mags2L=cell(1,nlp); mags2G=cell(1,nlp); mags2S=cell(1,nlp+1); mags2T=cell(1,nlp+1); rms2T=cell(1,nlp+1); rms2Tv=cell(1,nlp+1); % LOOP START for i=1:nlp f_blend = fblendlist{i}; w_blend = f_blend*pp; for n=0:n_blend; cf_poly(n+1)=nchoosek(n_blend,n)*w_blend^(n); end blend_HP = tf([cf_poly(1:nbd),zeros(1,nbd)],cf_poly); blend_LP = tf(cf_poly(nbd+1:n_blend+1),cf_poly); blend_LVDT = blend_LP; blend_GEO = minreal(blend_HP*georesp*vel2disp); 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); [mags2L{i},~]=bodesus(sysc,'n_LVDT_LF0','IFO_LTM',freq); [mags2G{i},~]=bodesus(sysc,'n_GEO_LF0','IFO_LTM',freq); [mags2S{i},~]=bodesus(sysc,'accLGND','IFO_LTM',freq); mags2L{i}=mags2L{i}.*n_LVDT_LF0; mags2G{i}=mags2G{i}.*n_GEO_LF0; mags2S{i}=mags2S{i}.*freq.*freq*pp*pp.*noise_seism; mags2T{i}=sumpsd({mags2L{i},mags2G{i},mags2S{i}}); rms2T{i}=makerms(freq,mags2T{i}); rms2Tv{i}=makerms(freq,mags2T{i}.*freq*pp); end servo_LF0=servo_LF0*0; st =linmod(mdlfile); sysc =ss(st.a,st.b,st.c,st.d,'inputname',invl,'outputname',outvl); [mags2S{nlp+1},~]=bodesus(sysc,'accLGND','IFO_LTM',freq); mags2S{nlp+1}=mags2S{nlp+1}.*freq.*freq*pp*pp.*noise_seism; mags2T{nlp+1}=mags2S{nlp+1}; rms2T{nlp+1}=makerms(freq,mags2T{nlp+1}); rms2Tv{nlp+1}=makerms(freq,mags2T{nlp+1}.*freq*pp); %% TOTAL mypsdplotopt(mags2T,freq,... 'title','Total TM displacement with different blending frequency [m/rtHz]',... 'legend',lgdname,... 'ylim',[1e-12,1e-5]) export_fig('figure/typeB1proto_IPdamp_blend_total_150618.pdf') %% LVDT mypsdplotopt(mags2L,freq,... 'title','LVDT noise coupling [m/rtHz]',... 'legend',lgdname,... 'ylim',[1e-12,1e-5]) export_fig('figure/typeB1proto_IPdamp_blend_lvdt_150618.pdf') %% GEO mypsdplotopt(mags2G,freq,... 'title','Geophone noise coupling [m/rtHz]',... 'legend',lgdname,... 'ylim',[1e-12,1e-5]) export_fig('figure/typeB1proto_IPdamp_blend_geo_150618.pdf') %% SEISM mypsdplotopt(mags2S,freq,... 'title','Seismic noise coupling [m/rtHz]',... 'legend',lgdname,... 'ylim',[1e-12,1e-5]) export_fig('figure/typeB1proto_IPdamp_blend_seism_150618.pdf') %% RMS mypsdplotopt(rms2T,freq,... 'title','Total RMS [m]',... 'legend',lgdname,... 'ylim',[1e-12,1e-3]) export_fig('figure/typeB1proto_IPdamp_blend_rms_150618.pdf') %% RMS mypsdplotopt(rms2Tv,freq,... 'title','Total RMS [m/sec]',... 'legend',lgdname,... 'ylim',[1e-12,1e-3]) export_fig('figure/typeB1proto_IPdamp_blend_rmsV_150618.pdf') %% compare for i=1:nlp mypsdplotopt({mags2G{i},mags2L{i},mags2S{i},mags2T{i},rms2T{i}},freq,... 'title',['Blending at ',num2str(fblendlist{i}),' Hz'],... 'legend',{'GEO','LVDT','Seismic','total','RMS'},... 'ylim',[1e-12,1e-5],'ylabel','Magnitude [m/rtHz]/[m]') export_fig(['figure/typeB1proto_IPdamp_blend_',num2str(fblendlist{i}*1000),'mHz_150618.pdf']) end %% Blend f_blend = 0.5; w_blend = f_blend*pp; for n=0:n_blend; cf_poly(n+1)=nchoosek(n_blend,n)*w_blend^(n); end blend_HP = tf([cf_poly(1:nbd),zeros(1,nbd)],cf_poly); blend_LP = tf(cf_poly(nbd+1:n_blend+1),cf_poly); mybodeplot({blend_LP,blend_HP},freq) %% NOISE MODEL [maggr,~]=mybode(georesp*vel2disp,freq); mypsdplotopt({n_LVDT_LF0,n_GEO_LF0.*maggr,noise_seism},freq,... 'title','Noise model',... 'legend',{'LVDT','Geophone','Seismic'},... 'ylim',[1e-14,1e-4],'unit','[m/rtHz]') export_fig('figure/typeB1proto_IPdamp_noisemodel_150618.pdf') %% No ctrl mypsdplotopt({mags2T{nlp+1},rms2T{nlp+1},rms2Tv{nlp+1},noise_seism},freq,... 'title','Seismic noise coupling',... 'color',{'r-','m-','m-.','k-'},... 'legend',{'PSD','RMS','RMS velocity','Ground PSD'},... 'ylabel','Magnitude [m/rtHz][m]/[m/sec]',... 'ylim',[1e-12,1e-5]) export_fig('figure/typeB1proto_IPdamp_noctrl_150618.pdf')