%==================================================
% Noise Plotter
% Loads ascii darm & as_q data from the DTT
% exported text file and uses the loop gain model's
% State Space Object 'darmsys' to make the noise curve
%===================================================


S1darm
S1darm_08
[A,B,C,D]=linmod('S1darm_08');
darmsys=ss(A,B,C,D);

%load noise_020825.txt
%z=noise_020825;
%f = z(:,1);
%dc = z(:,2);
%asqc = z(:,3);

load 020906_2253_cal %load measured spectra asqc and dc, in counts/rtHz,
                     %and frequency f in Hertz.

w0 = 2 * pi * 0.76; %pendulum frequency
Q = 10;   
fc = 87.3;          %cavity pole

darmcal = 1.72e-9; %1.72nm/count, measured calibration for DARM_CTRL
asqcal = 2.34e-15;  %m/count, to meet DARM_CTRL at ugf at 250 Hz

darmtf = darmcal * tf(w0^2,[1 w0/Q w0^2]);
asqtf = asqcal * tf([1 2*pi*fc],[2*pi*fc]);
at = mybodesys(asqtf,f);
dt = mybodesys(darmtf,f);

%noise in m/rtHz
darmn = abs(dt) .* dc;    
asqn = abs(at) .* asqc;

Lmin2asqout=abs(mybodesys(darmsys(3,1),f));
Lmin2darmout=abs(mybodesys(darmsys(4,1),f));
DigFlt=abs(asq2darm*GW_K*mybodesys(DARMDF,f));

gwp1 = asqc./Lmin2asqout;
gwp2 = dc./Lmin2darmout;

calib=1./mybodesys(darmsys(3,1),f);
calib=calib/4e3;

%fid=fopen('ASQcalibration.txt','w');
%fprintf(fid,'%s\t','Freq(Hz)');
%fprintf(fid,'%s\t','Mag(strain/ct)');
%fprintf(fid,'%s\n','Phase(rad)');
%
%for i=1:length(f)
%  fprintf(fid,'%s\t',num2str(f(i)));
%  fprintf(fid,'%s\t',num2str(abs(calib(i))));
%  fprintf(fid,'%s\n',num2str(angle(calib(i))));
%end
%fclose(fid)

figure(1)
clf
%loglog(f,darmn,'b',f,asqn,'k',f,gwp1,'r',f,gwp2,'g')
loglog(f,darmn,'b',f,asqn,'k',f,gwp1,'r')
%loglog(f,darmn,f,gwp2)
grid
%axis tight
axis([9 5000 1e-18 1.1e-10])
set(gca,'YTick',10.^[-17:-10])
ylabel('Displacement (m/rtHz)')
title('Noise')
xlabel('Frequency (Hz)')
legend('measured DARM','measured AS\_Q',...
       'calibrated AS\_Q');
%'modelled noise from DARM\_CTRL')
shg

figure(2)
clf
%loglog(f,darmn,'b',f,asqn,'k',f,gwp1,'r',f,gwp2,'g')
loglog(f,darmn,'b',f,asqn,'k',f,gwp1,'r')
%loglog(f,darmn,f,gwp2)
grid
%axis tight
axis([50 2000 1e-18 1.1e-14])
set(gca,'YTick',10.^[-17:-10])
ylabel('Displacement (m/rtHz)')
title('Noise')
xlabel('Frequency (Hz)')
legend('measured DARM','measured AS\_Q',...
       'calibrated AS\_Q');
%'modelled noise from DARM\_CTRL')

shg

%save noise_0825 f gwp1 darmn asqn

figure(3)
%check the model has the right digital filters

loglog(f,dc./asqc,f,DigFlt,'.-',f,Lmin2darmout./Lmin2asqout)
legend('DARM\_CTRL/AS\_Q','asq2darm*DARMDF*GW\_K','darmsys(4,1)/darmsys(3,1)')
shg