
%% File Names

datafile='DAT0040.HDR';
filfile='fil1';
olfile='ol2135';

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

%% Data read

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

%%

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

[r,c]=size(yo);

%%

[cf,fr]=c_factor(fr,filfile,olfile,calib);

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

%%

%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])

%%

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}]'])

%%

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

%%
dat=[fr',disar];
save disa100715.dat dat -ascii

