function pitotDataAnalysis

%% Loading data

% CSC vent (14 August 2016)
cd('~/Dropbox/eLab_Chantier01/PitotData');
% [timeCSC_20160814 vData] = localFunction2('CSC_20160814.csv'); % Peut-être pour plus tard...
load('CSC_20160814.mat')
timeCSC_20160814 = Vdiff2010.XData;
dataCSC_20160814 = Vdiff2010.YData;
dataCSC_20160814 = 1000*(dataCSC_20160814-min(dataCSC_20160814))/0.25;
directory = '~/Dropbox/DIAPHANE_Soufriere_Summit_Vents_SlowData/CS_Fumerolles/MatFiles';
filePrefix = 'CS_Fumerolles*';
channelList = {'CSC_Temp'};
[tCSC,tempCSC] = localFunction1(directory,filePrefix,channelList,[2016 08 17 18 20 00],[2016 08 17 18 30 00],90,105);
tempMeanCSC_20160814 = mean(tempCSC);

% CSN vent (14 August 2016)
cd('~/Dropbox/eLab_Chantier01/PitotData');
load('CSN_20160814.mat');
timeCSN_20160814 = Vdiff2010.XData(200:end);
dataCSN_20160814 = Vdiff2010.YData(200:end);
dataCSN_20160814 = 1000*(dataCSN_20160814-min(dataCSN_20160814))/0.25;
directory = '~/Dropbox/DIAPHANE_Soufriere_Summit_Vents_SlowData/CS_Fumerolles/MatFiles';
filePrefix = 'CS_Fumerolles*';
channelList = {'CSN_Temp_2'};
[tCSN2,tempCSN2] = localFunction1(directory,filePrefix,channelList,[2016 08 17 18 20 00],[2016 08 17 18 30 00],90,105);
tempMeanCSN_20160814 = mean(tempCSN2);


% CSN vent (17 August 2016)
load('CSN_20160817.mat');
timeCSN_20160817 = Vpdiff2010.XData;
dataCSN_20160817 = Vpdiff2010.YData;
dataCSN_20160817 = 1000*(dataCSN_20160817-min(dataCSN_20160817))/0.25;
directory = '~/Dropbox/DIAPHANE_Soufriere_Summit_Vents_SlowData/CS_Fumerolles/MatFiles';
filePrefix = 'CS_Fumerolles*';
channelList = {'CSN_Temp_2'};
[tCSN2,tempCSN2] = localFunction1(directory,filePrefix,channelList,[2016 08 17 19 38 00],[2016 08 17 19 42 00],90,105);
tempMeanCSN_20160817 = mean(tempCSN2);

% CSS vent
load('CSS_20160817.mat');
timeCSS_20160817 = Vpdiff2010.XData;
dataCSS_20160817 = Vpdiff2010.YData;
dataCSS_20160817 = 1000*(dataCSS_20160817-min(dataCSS_20160817))/0.25;
directory = '~/Dropbox/DIAPHANE_Soufriere_Summit_Vents_SlowData/CS_Fumerolles/MatFiles';
filePrefix = 'CS_Fumerolles*';
channelList = {'CSS_Temp_3'};
[tCSS3,tempCSS3] = localFunction1(directory,filePrefix,channelList,[2016 08 17 20 23 00],[2016 08 17 20 33 00],90,105);
tempMeanCSS_20160817 = mean(tempCSS3);

%% Dessin des données brutes (converties en hPa)

fontSize = 18;
f1 = figure('color','white');
screenSize = get(0,'ScreenSize');
set(f1,'Position',[0 0 screenSize(3) screenSize(4)]);
subplot1(4,1,'Gap',[0.02 0.06],'XTickL','All');

subplot1(1);
set(gca,'FontSize',fontSize); set(gca,'FontName','Helvetica');
hold on
plot(timeCSC_20160814,dataCSC_20160814,'lineWidth',2,'color','black');
ylabel('CSC 2016-08-14 (hPa)')
axis tight
grid on

subplot1(2);
set(gca,'FontSize',fontSize); set(gca,'FontName','Helvetica');
hold on
plot(timeCSN_20160814,dataCSN_20160814,'lineWidth',2,'color','black');
ylabel('CSN 2016-08-14 (hPa)')
axis tight
grid on

subplot1(3);
set(gca,'FontSize',fontSize); set(gca,'FontName','Helvetica');
hold on
plot(timeCSN_20160817,dataCSN_20160817,'lineWidth',2,'color','black');
ylabel('CSN 2016-08-17 (hPa)')
axis tight
grid on

subplot1(4);
set(gca,'FontSize',fontSize); set(gca,'FontName','Helvetica');
hold on
plot(timeCSS_20160817,dataCSS_20160817,'lineWidth',2,'color','black');
xlabel('Time (s)')
ylabel('CSS 2016-08-17 (hPa)')
axis tight
grid on

export_fig ~/Dropbox/eLab_Chantier01/PitotData/pitotFigure1.png -r300

%% Calculs des histogrammes

fontSize = 18;
f1 = figure('color','white');
screenSize = get(0,'ScreenSize');
set(f1,'Position',[0 0 screenSize(3) screenSize(4)]);
subplot1(2,2,'Gap',[0.06 0.06],'XTickL','All','YTickL','All');

subplot1(1);
set(gca,'FontSize',fontSize); set(gca,'FontName','Helvetica');
hold on
[f,xi] = ksdensity(dataCSC_20160814); 
plot(xi,f,'lineWidth',2,'color','black');
ylabel('CSC 2016-08-14')
axis tight
grid on

subplot1(2);
set(gca,'FontSize',fontSize); set(gca,'FontName','Helvetica');
hold on
[f,xi] = ksdensity(dataCSN_20160814); 
plot(xi,f,'lineWidth',2,'color','black');
ylabel('CSN 2016-08-14')
axis tight
grid on

subplot1(3);
set(gca,'FontSize',fontSize); set(gca,'FontName','Helvetica');
hold on
[f,xi] = ksdensity(dataCSN_20160817); 
plot(xi,f,'lineWidth',2,'color','black');
xlabel('p (hPa)')
ylabel('CSN 2016-08-17')
axis tight
grid on

subplot1(4);
set(gca,'FontSize',fontSize); set(gca,'FontName','Helvetica');
hold on
[f,xi] = ksdensity(dataCSS_20160817); 
plot(xi,f,'lineWidth',2,'color','black');
xlabel('p (hPa)')
ylabel('CSS 2016-08-17')
axis tight
grid on

export_fig ~/Dropbox/eLab_Chantier01/PitotData/pitotFigure2.png -r300

%% Petits calculs

atmPressure = 856*100; % Ambiant pressure (Pa)

% CSS vent
ventDiameter = 0.25; % Vent diameter (m)
rhoSteam = steamProperties(tempMeanCSS_20160817,2,4); % Density (kg/m3)
deltaP = 6500; % Differential pressure (Pa)
vPitot = sqrt(2*deltaP/rhoSteam); % Output velocity (m/s)
innerPressure = atmPressure+deltaP;
qS = steamFlowOrifice(innerPressure/100000,atmPressure/100000,ventDiameter);
vThermo = qS./steamProperties(atmPressure/100000,1,4)/(pi/4*ventDiameter^2)/3600; % Output velocity (m/s)
totalHeat = steamProperties(tempMeanCSS_20160817,2,7);
powerVent = qS.*totalHeat/3600/1000; % Power (MW)

disp(['Thermal power CSS vent   = ' num2str(powerVent) ' MW'])
disp(['velocity Pitot CSS vent  = ' num2str(vPitot) ' m/s'])
disp(['velocity thermo CSS vent = ' num2str(vThermo) ' m/s'])

% CSN vent
ventDiameter = 0.12; % Vent diameter (m)
rhoSteam = steamProperties(tempMeanCSN_20160814,2,4); % Density (kg/m3)
deltaP = 910; % Differential pressure (Pa)
vPitot = sqrt(2*deltaP/rhoSteam); % Output velocity (m/s)
innerPressure = atmPressure+deltaP;
qS = steamFlowOrifice(innerPressure/100000,atmPressure/100000,ventDiameter);
vThermo = qS./steamProperties(atmPressure/100000,1,4)/(pi/4*ventDiameter^2)/3600; % Output velocity (m/s)
totalHeat = steamProperties(tempMeanCSN_20160814,2,7);
powerVent = qS.*totalHeat/3600/1000; % Power (MW)

disp(['Thermal power CSN vent   = ' num2str(powerVent) ' MW'])
disp(['velocity Pitot CSN vent  = ' num2str(vPitot) ' m/s'])
disp(['velocity thermo CSN vent = ' num2str(vThermo) ' m/s'])

%% Thermal power of fumarole vent

ventDiameter = 0.25; % Vent diameter (m)
tempSteam = 106;
atmPressure = 0.860;
pSteam = steamProperties(tempSteam,2,1); % Pressure (bar)
qS = steamFlowOrifice(pSteam,atmPressure,ventDiameter);
vThermo = round(qS./steamProperties(atmPressure,1,4)/(pi/4*ventDiameter^2)/3600); % Output velocity (m/s)
totalHeat = steamProperties(tempSteam,2,7);
powerVent = qS.*totalHeat/3600/1000; % Power (MW)

disp(['Thermal power CSS vent   = ' num2str(powerVent) ' MW'])
disp(['velocity thermo CSS vent = ' num2str(vThermo) ' m/s = ' num2str(3.6*vThermo) ' km/h = Mach ' num2str(vThermo/340)])
disp(['Steam production = ' num2str(round(100*qS/1000)/100) ' t/h'])

end

function [tData,sData] = localFunction1(directory,filePrefix,channelList,dateStart,dateEnd,sMinCorrect,sMaxCorrect)

disp(['Going to directory: ' directory])
cd(directory);
disp(filePrefix)
[tData sData] = selectData(dateStart,dateEnd,filePrefix,channelList);
if ~isempty(tData)
    disp(['Size of data array = ' num2str(length(tData))])
    disp(['Date start = ' datestr(tData(1))])
    disp(['Date end = ' datestr(tData(end))])
    disp(['Channel list: ' channelList])
    
    % median filter
    nMedian = 5;
    sData = filterMedian(sData,nMedian);
    tData = tData((nMedian+1)/2:end-(nMedian-1)/2);
    
    % Remove eventual false data (e.g. maintenance operations on acquisition system)
    [tData ix] = unique(tData);
    sData = sData(ix);
    [tData ix] = sort(tData);
    sData = sData(ix);
    ixFalse = find(sData > sMaxCorrect | sData < sMinCorrect);
    sData(ixFalse) = NaN;
    ixCorrect = find(~(sData > sMaxCorrect | sData < sMinCorrect));
%     sData(ixFalse) = interp1(tData(ixCorrect),sData(ixCorrect),tData(ixFalse));
    disp([channelList ' Min - Moy - Max: ' num2str(0.01*round(100*min(sData(ixCorrect)))) ...
        num2str(0.01*round(100*mean(sData(ixCorrect)))) num2str(0.01*round(100*max(sData(ixCorrect))))])
end

end

function [tData vData] = localFunction2(filename)


txt = textread(filename, '%s', 'delimiter', '\n');
for k = 3:length(txt)
    tData(k) = datenum(txt{k}(1:19),'dd/mm/yyyy,HH:MM:SS');
end
tData = tData(3:end);
vData = csvread(filename,2,2);

end
