%% Alexander Penn, 10.2019, ETH Zurich
% Open data to the manuscript: Penn et al., CEJ 2019, https://doi.org/10.1016/j.cej.2019.12318
% This code reads in MRI images time series and performs a frequency
% analysis of the break-off of bubbles presented in Figure 7
% For further information or assistance please contact Alexander Penn at his email (current email address: apenn@ethz.ch)
%% Frequency analysis 
clear all
close all

load('MRImages.mat')

%% Relevant scans of the bubble detaching regimes
scanIDs = [3:15,16:30,32:45,51:64,66:79, 84:89,99:112];
lineCenter(scanIDs) = 25; % Corresponds to a vertical distance of 150 mm between the the ROI and the distributor.
columnCenter(scanIDs) = 49;

CollapsePoint = 65;  % For the pulsating jet with bubble collapse regime, the vertical position of the ROI in the bed is 40 mm lower compared to the pulsating jet regime and the puslating jet and freely bubbling regime as described in the manuscript. 
lineCenter(109) = CollapsePoint;
lineCenter(110) = CollapsePoint;
lineCenter(14) = CollapsePoint;
lineCenter(15) = CollapsePoint;
lineCenter(30) = CollapsePoint;

lineNumbers = [lineCenter-4 : lineCenter+4]; % Defines the width of the ROI in which data is analyzed (mm)
rowNumbers = [columnCenter-3 : columnCenter+3]; % Defines the height of the ROI in which data is analyzed (mm)


%% Frequency analysis
fs=3000;
fBubble = zeros(2,106);
for scanID = scanIDs
    clear kLine dataLine peakLocs
    lineNumbers = [lineCenter(scanID)-4 : lineCenter(scanID)+4]; 
    rowNumbers = [columnCenter(scanID)-3 : columnCenter(scanID)+3]; 
    dataROI = squeeze(sum(sum(data(lineNumbers,rowNumbers,:,scanID),1),2)); 
    DCOffset = mean(dataROI);
    dataNorm = dataROI - DCOffset;
    kLine = abs(fftshift(fft(dataNorm,fs)));
    kNorm = smooth(kLine,3);
    [peakVals,peakLocs,peakWidth, peakProminence] = findpeaks(kNorm,'sortstr','descend','NPeaks',2,'MinPeakProminence',3e3);
    findpeaks(kNorm,'sortstr','descend','NPeaks',2,'MinPeakDistance',fs/100,'MinPeakProminence',3e3);
    fBubble(1,scanID)=0.5*(fs/abs(peakLocs(2)-peakLocs(1))*Tdyn)^-1; %The factor 0.5 comes from the fact that we are looking at the separation between the two strongest peaks in the symmetric spectrum
    fBubble(2,scanID)=(fs/abs(peakWidth(2))*Tdyn)^-1;
   end
%% FroudeNumberConversion 
dV_mf = 442; %[lpm] Gas flow rate of minimum fluidization
dV_mf_SI = dV_mf*10^(-3)/60; % [m^3/s] Gas flow rate of minimum fluidization in SI units
dV_or = 15*linspace(1,14,14);% [lpm] Gas flow rate through orifice
dV_or_SI = dV_or*10^(-3)/60; % [m^3/s] Gas flow rate through orifice in SI units
d_or = 0.006; % [m] Diameter of orifice
d_bed = 0.190; % [m] Diameter of bed
A_or = (d_or/2)^2*pi; %[m^2] Surface area of orifice 
A_bed = (d_bed/2)^2*pi;% [m^2] Surface area of fluidized bed
U_or = dV_or_SI/A_or; % [m/2] Gas velocity through the orifice
U_mf = dV_mf_SI/A_bed; % [m/s]
Fr_or = (U_or).^2/(9.81*d_or);% [-] Froude number of orifice flow
U_or_norm = (U_or/U_mf)*(A_or/A_bed);% [-] Gas flow thhrough orifice nrmalized by surface area ratio between the fluidized bed and the orifice
Fr_flipped = fliplr(Fr_or)/10^5;

%% Plot Figure 7

markerSize = 14;
lineWidth = 1.5;
fig1=figure(1);
fontSize = 17.5;
set(fig1,'Units', 'pixels', 'Position', [0, 300, 550, 550]);

plot(nan, nan,'k s','markersize',markerSize,'linewidth',lineWidth), hold on % These blank plots are placedholders to get the legend right.
plot(nan, nan,'k o','markersize',markerSize,'linewidth',lineWidth), hold on
plot(nan, nan,'k x','markersize',markerSize,'linewidth',lineWidth), hold on
plot(nan, nan,'k d','markersize',markerSize,'linewidth',lineWidth), hold on
plot(nan, nan,'k +','markersize',markerSize,'linewidth',lineWidth), hold on
plot(nan, nan,'k v','markersize',markerSize,'linewidth',lineWidth), hold on
plot(nan, nan,'k >','markersize',markerSize,'linewidth',lineWidth), hold on

plot(Fr_flipped(1:5),fBubble(1,[84,85,87:89]),'s','markersize',markerSize,'MarkerEdgeColor',[0,1,0],'linewidth',lineWidth)
plot(Fr_flipped(1:8),fBubble(1,99:106),'o','markersize',markerSize,'MarkerEdgeColor',	[0,1,0],'linewidth',lineWidth)
plot(Fr_flipped(1:11),fBubble(1,3:13),'x','markersize',markerSize,'MarkerEdgeColor',	[0,1,0],'linewidth',lineWidth)
plot(Fr_flipped(1:13), fBubble(1,17:29),'d','markersize',markerSize,'MarkerEdgeColor',	[0,1,0],'linewidth',lineWidth)
plot(Fr_flipped(1:14), fBubble(1,32:45),'+','markersize',markerSize,'MarkerEdgeColor',[0,1,0],'linewidth',lineWidth)
plot(Fr_flipped(1:14), fBubble(1,66:79),'v','markersize',markerSize,'MarkerEdgeColor',	[1,0.5,0],'linewidth',lineWidth)
plot(Fr_flipped(1:14), fBubble(1,51:64),'>','markersize',markerSize,'MarkerEdgeColor',[1,0.5,0],'linewidth',lineWidth)
plot(Fr_flipped(11:12), fBubble(1,109:110),'o','markersize',markerSize,'MarkerEdgeColor',[1, 0, 1],'linewidth',lineWidth)
plot(Fr_flipped(12:13), fBubble(1,14:15),'x','markersize',markerSize,'MarkerEdgeColor',[1, 0, 1],'linewidth',lineWidth)
plot(Fr_flipped(14), fBubble(1,30),'d','markersize',markerSize,'MarkerEdgeColor',[1, 0, 1],'linewidth',lineWidth)
text(1.65,13,'Pulsating jet','Color',[0,1,0], 'FontSize',fontSize)
text(-0.05,13.7,{'Pulsating jet with', 'bubble collapse'},'Color',[1, 0, 1], 'FontSize',fontSize)
text(0.85,5,{'Pulsating jet and', 'freely bubbling'},'Color',[1,0.5,0], 'FontSize',fontSize)

h=legend('0.57','0.68','0.79','0.90','1.02','1.13','1.24','location',[0.74 0.24 0.1374 0.2464]);
[hleg,att] = legend('show');
title(hleg,'{\itU}_{bg}/{\itU}_{mf}')
ylabel('{\itf}_{breakoff} (Hz)');
xlabel('{\itFr} × 10^5 (-)');
set(gca,'XDir','Reverse');
set(gca,'XTickMode','manual');
set(gca,'fontsize',fontSize)
set(gca,'Xdir','normal')
ylim([0 15.1])
xlim([-0.1 2.9])
yticks([0 2 4 6 8 10 12 14])
set(gcf,'renderer','painters')
set(gcf,'color','w');

print('-dtiff','Figure7FrequencyAnalysis','-r450')