
%% File Names

datafile='DAT0000.HDR';

fil1file='mass_fil1';
ol_mass='mass1';
fil2file='pzt_fil1';
olfile='ol2';

calib=(8.6+10)/2*4*pi/1.064e-6;

%% Data read

[d1,time,date,prop]=ykdl(datafile, 1);
[d2,time,date,prop]=ykdl(datafile, 2);

l=max(find(d2<=-13));

if not(isempty(l))
d1=d1(l+1:end);
d2=d2(l+1:end);
end;

%% Spectrogram

nfft=2^17;   % FFT length
Fs=1e3;      % Sampling rate
window=hanning(nfft); Be=8/3; % Window Function
noverlap=nfft/2;
 
[yo,fr,to] = specg_data(d1,nfft,Fs,window,noverlap);
[avg,avgr]=spe_avg(yo,0.5);

[r,c]=size(yo);

%% Calibration factor

%[cf,fr]=c_factor(fr,filfile,olfile,calib);
[cf,fr]=c_factor2(fr,fil1file,ol_mass,fil2file,olfile,calib);

%% Calibrate to m/sqrt(Hz)

dis=sqrt(yo).*cf(:,ones(1,c)).*sqrt(Be);
disa=sqrt(avg).*cf.*sqrt(Be);
disar=sqrt(avgr).*cf.*sqrt(Be);

%% Plot noise spectrum

%d_dat=[fr,d0,fr,disp];
%save disp100719_mi1.dat d_dat -ascii

figure(1)
loglog(fr,dis)
grid on
hold on
loglog(fr,disa,fr,disar,'LineWidth',3); 
hold off
xlabel('Frequency  [Hz]')
ylabel('Disp. noise [m/Hz^{1/2}]')
axis([1e-2 1e2 1e-14 1e-4])

%% Plot sensitivity drift

f003=find(fr>=0.03, 1 );
f01=find(fr>=0.1, 1 );
f03=find(fr>=0.3, 1 );
ti=to/3600;
dones=ones(length(ti),1);

figure(2)
semilogy(ti,dis(f003,:),ti,dis(f01,:),ti,dis(f03,:));
grid on
hold on
semilogy(ti,disar(f003,:).*dones,ti,disar(f01,:).*dones,ti,disar(f03,:).*dones);
hold off
xlabel('Time  [Hour]')
ylabel('Disp. noise [m/Hz^{1/2}]')
title(['Avg. value  0.03Hz: ',num2str(disar(f003)),', 0.1Hz: ',...
       num2str(disar(f01)),', 0.3Hz: ',num2str(disar(f03)),' [m/Hz^{1/2}]'])

%% Plot spectrogram

figure(8)
imagesc(to/3600,fr(3:50),log10(dis(3:50,:)))
colorbar
xlabel('Time  [Hour]')
ylabel('Frequency  [Hz]')

%% Data save
fr=fr(:);

dat=[fr,disar];
save disa100904_1.dat dat -ascii

