% I.Fiori, 8 April 2020
% This script: reads data  from a specified acoustic transfer function 
% file (written by combine_TF.m)
% then it gets frame data of Hrec and MICs at the given gps
% time, computes the acoustic projection for these data and plots and 
% save it.
%------------------------------

clear all
close all

set(0,'DefaultAxesFontSize',10);
set(0,'DefaultTextFontSize',10);
run /users/swinkels/deploy/MatlabVirgoTools/trunk/startup.m
format longg

%---- read AdV+ design curves
[fO4,hO4]=textread('/data/procdata/detchar/env/injections/acoustic/results/O3/AdV+_design_curves/O4baseline.txt');
[fO5,hO5]=textread('/data/procdata/detchar/env/injections/acoustic/results/O3/AdV+_design_curves/O5low.txt');

%--- directory where plot is saved
OUTdirectory='/data/procdata/detchar/env/injections/acoustic/results/O3/example_projection/';

%--- specify date when projection is to be computed
year=2020, month=3; day = 21; hour = 21; min = 30; sec = 0;
utc = datenum(year,month,day,hour,min,sec); % N = datenum(Y,MO,D,H,MI,S)
gps_current=utc2gps(utc)
datestring_current = [num2str(day),'/',num2str(month),'/',num2str(year)];

%--- specify which TF file to use
%
IDstring = 'October_6_2019';
directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/NEB_MIC/';
filename=[directory,'ENV_NEB_MIC_',char(IDstring),'_combined_TF.txt'];
% %
IDstring = 'August_19_2019';
directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/WEB_MIC/';
filename=[directory,'ENV_WEB_MIC_',char(IDstring),'_combined_TF.txt'];
% %
% IDstring = 'March_26_2020';
% directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/CEB_MIC/';
% filename=[directory,'ENV_CEB_MIC_',char(IDstring),'_combined_TF.txt'];
% %
% IDstring = 'no_SQZ_March_28_2020';
% directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/CEB_MIC/';
% filename=[directory,'ENV_CEB_MIC_',char(IDstring),'_combined_TF.txt'];
% %
% IDstring = 'March_28_2020';
% directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/SQZ_MIC/';
% filename=[directory,'ENV_SQZ_MIC_',char(IDstring),'_combined_TF.txt'];
% %
% IDstring = 'March_28_2020';
% directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/EDB_MIC/';
% filename=[directory,'ENV_EDB_MIC_',char(IDstring),'_combined_TF.txt'];
% %
% IDstring = 'no_SQZ_March_28_2020';
% directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/SQZ_MIC/';
% filename=[directory,'ENV_SQZ_MIC_',char(IDstring),'_combined_TF.txt'];
% %
% IDstring = 'no_SQZ_March_28_2020';
% directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/EDB_MIC/';
% filename=[directory,'ENV_EDB_MIC_',char(IDstring),'_combined_TF.txt'];
% %
% IDstring = 'April_26_2019';
% directory = '/data/procdata/detchar/env/injections/acoustic/results/O3/CEB_MIC/';
% filename=[directory,'ENV_CEB_MIC_',char(IDstring),'_combined_TF.txt'];
%

%---- reads parameters from 3rd line
[where,ch1,ch2,gps,dur,fs,T,THR_Hrec,THR_sensor,sk_hrec,sk_sensor]=textread(filename,'%s %s %s %f %f %f %f %f %f %f %f',1,'delimiter',',','commentstyle','matlab','headerlines',2);
[yy,mm,dd]=datevec(gps2utc(gps));
datestring = [num2str(dd),'/',num2str(mm),'/',num2str(yy)]
%---- reads values starting from line 7th
[f,hrec,mic,flagUL,proj,TF]=textread(filename,'','delimiter',',','commentstyle','matlab','headerlines',6);

%---- reads Hrec and MIC data 
dec = 4;
fs = 5000;
%T = 10; % use same T parameter as read from TF file!
dur = 600;
npt = T*fs;
ffl = '/virgoData/ffl/raw.ffl';
%ch1 = 'V1:Hrec_hoft_20000Hz'; % channel names are read from TF file

datan = getChannel(ffl,char(ch1),gps_current,dur);
data = decimate(datan,dec, 'fir'); 
[d,f] = pwelch(data,hann(npt),[],npt,fs);
hrec = sqrt(d);
datan = getChannel(ffl,char(ch2),gps_current,dur);
data = decimate(datan,dec, 'fir'); 
[d,f] = pwelch(data,hann(npt),[],npt,fs);
MICraw = sqrt(d);

%%-----
MIC=movmean(MICraw,sk_sensor); %  SMOOTHING MIC

TF_measured = TF.*flagUL;
idx = find(TF_measured==0);
TF_measured(idx)=NaN;

projection = TF.*MIC;
projection_measured = TF_measured.*MIC;

%%-----
figure(100)
%loglog(f,hrec,'k'); hold on; 
h(1)=loglog(f,hrec,'k'); hold on; 
h(2)=loglog(f,projection,'s','Color',[0.7,0.7,0.7],'MarkerSize',1,'MarkerFaceColor',[0.7,0.7,0.7]);
h(3)=loglog(f,projection_measured,'sr','MarkerSize',1, 'MarkerFaceColor','r'); 
ylim([1e-26 1e-20]); % axis tight;
xlim([10 2000])
grid
h(4)=loglog(fO4,hO4);
h(5)=loglog(fO5,hO5);
hold off;
ylabel('Hrec ASD[1/sqrt(Hz)]')
xlabel('Frequency [Hz]')
title(['Acoustic projection: ',strrep(char(ch2),'_','\_')]);
legend([h(1), h(3), h(2), h(4), h(5)],[strrep(char(ch1),'_','\_'),' of ',char(datestring_current)],['measured, ',strrep(IDstring,'_',' ')],['upper limit, ',strrep(IDstring,'_',' ')],'O4 design','O5 design');
% h = findobj(gca,'Type','line');
%legend([h(3), h(1), h(2)],'Hrec\_hoft','projection','Upper limit')

%%---- save plots
filename = [OUTdirectory,num2str(gps_current),'_',char(ch2),'_',IDstring,'_projection','.png'];
saveas(gcf,filename,'png');
filename = [OUTdirectory,num2str(gps_current),'_',char(ch2),'_',IDstring,'_projection','.fig'];
saveas(gcf,filename,'fig');
