# Python file for plotting the results from the cyclic water retention test
# on Hamburg Sand with parallel CT-imaging
# --------------------------------------------------------------------------
# CT scans using the Laboratoire 3SR X-ray tomograph at Univ. Grenoble Alpes
# 21.07.2019 - 22.07.2019, Grenoble
# 
# --------------------------------------------------------------------------
# This Python script visualises measured contact angles on drainage
# and imbibition paths.

# Data file: Contact_angle_and_radii_of_curvature_data.txt

# Import libraries:
# ==========================================================================

import matplotlib.pyplot as plt
from matplotlib import style
style.use('ggplot')

import matplotlib.colors as colors

from matplotlib import rcParams

from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,
                               AutoMinorLocator)

from numpy.polynomial.polynomial import polyfit

import numpy as np
# -----------------------------------------
savefigs = True    # Save figures to file?

# Dimensions of figures:
cm_x = 16
cm_y = 14
# =======================================================================
SMALL_SIZE = 12
MEDIUM_SIZE = 12
BIGGER_SIZE = 14

plt.rc('font', size=MEDIUM_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=MEDIUM_SIZE)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
# ---------------------------
COLOR = 'black'
plt.rcParams['axes.linewidth'] = 0.5
plt.rcParams['axes.edgecolor'] = COLOR
plt.rcParams['text.color'] = COLOR
plt.rcParams['axes.labelcolor'] = COLOR
plt.rcParams['xtick.color'] = COLOR
plt.rcParams['ytick.color'] = COLOR
# ---------------------------
#plt.rcParams["font.family"] = "Times New Roman"
rcParams['font.family'] = 'serif'
# =======================================================================

my_filename_001 = "Contact_angle_and_radii_of_curvature_data.txt"

mm_per_pix = 0.01           # mm/px
surface_tension = 0.072750  # N/m


image_number,step,direction,meniscus_number,angle_left,angle_right,r_curv = np.loadtxt(my_filename_001,                                                  
                                   unpack = True,
                                   delimiter = '\t',
                                   skiprows = 2,
                                   usecols = (0,1,2,3,4,5,8),
                                   dtype = {'names': ('Image number','Step','Direction','Meniscus_number','contact_angle_1','contact_angle_2','radius'),
                                            'formats': (np.float, np.float,'|S10', np.float, np.float, np.float, np.float)}
                                                                                       )

direction = direction.astype('U10')

# radius of curvature in mm:
r_curv = r_curv*mm_per_pix                          # mm

# mean curvature (2D-assumption: second radius is infinity):
mean_curvature = 1/r_curv                           # 1/mm


# capillary pressure (Young-Laplace equation):
pc = 1/(r_curv/1000.0)*surface_tension/1000.0       # kPa

# Find number of drainage data:
count = 0
for i in range(0,len(angle_left)):
    if direction[i] == 'Drainage':
        count = count+1

n_drainage = count
n_imbibition = len(angle_left) - count
#print(n_drainage)

#print(n_imbibition)

# ==============
# Macroscopic suction and saturation values before and after the CT scans:
suction_before_scans_vect =   [0,0.7991,0.8486,1.0902,1.2327,1.3566,0.9354,0.6133,0.4460,0.3531,0.7743,1.0716,
    1.3442,0.7124,0.4894,1.0655,1.2141,0.8301,0.6566,0.9911];
saturation_before_scans_vect = [1.0000,0.8570,0.7642,0.5785,0.3928,0.3000,0.3928,0.5785,0.7642,0.8570,0.7642,
    0.5785,0.3928,0.5785,0.7642,0.5785,0.4857,0.5785,0.6713,0.5785];
# ----------------
suction_after_scans_vect = [0.3035,0.7681,0.8610,1.0283,1.1893,1.3256,0.9725,0.6566,0.4832,0.4274,0.7557,1.0221,
    1.3070,0.7433,0.5265,1.0283,1.2079,0.8734,0.6814,1.0531];
saturation_after_scans_vect = [1.0000,0.8570,0.7642,0.5785,0.3928,0.3000,0.3928,0.5785,0.7642,0.8570,0.7642,
    0.5785,0.3928,0.5785,0.7642,0.5785,0.4857,0.5785,0.6713,0.5785];
step_suction_vect = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] 
# ==============

rows_drainage = np.where(direction=='Drainage')
rows_imbibition = np.where(direction=='Imbibition')

rows_step_1 = np.where(step==1)
rows_step_2 = np.where(step==2)
rows_step_3 = np.where(step==3)
rows_step_4 = np.where(step==4)
rows_step_5 = np.where(step==5)
rows_step_6 = np.where(step==6)
rows_step_7 = np.where(step==7)
rows_step_8 = np.where(step==8)
rows_step_9 = np.where(step==9)
rows_step_10 = np.where(step==10)
rows_step_11 = np.where(step==11)
rows_step_12 = np.where(step==12)
rows_step_13 = np.where(step==13)
rows_step_14 = np.where(step==14)
rows_step_15 = np.where(step==15)
rows_step_16 = np.where(step==16)
rows_step_17 = np.where(step==17)
rows_step_18 = np.where(step==18)
rows_step_19 = np.where(step==19)

angle_left_step_1 = angle_left[rows_step_1]
angle_right_step_1 = angle_right[rows_step_1]
angle_step_1 = np.concatenate((angle_left_step_1,angle_right_step_1))
mean_curvature_step_1 = mean_curvature[rows_step_1]
pc_step_1 = pc[rows_step_1]
r_curv_step_1 = r_curv[rows_step_1]
#
angle_left_step_2 = angle_left[rows_step_2]
angle_right_step_2 = angle_right[rows_step_2]
angle_step_2 = np.concatenate((angle_left_step_2,angle_right_step_2))
mean_curvature_step_2 = mean_curvature[rows_step_2]
pc_step_2 = pc[rows_step_2]
r_curv_step_2 = r_curv[rows_step_2]
#
angle_left_step_3 = angle_left[rows_step_3]
angle_right_step_3 = angle_right[rows_step_3]
angle_step_3 = np.concatenate((angle_left_step_3,angle_right_step_3))
mean_curvature_step_3 = mean_curvature[rows_step_3]
pc_step_3 = pc[rows_step_3]
r_curv_step_3 = r_curv[rows_step_3]
#
angle_left_step_4 = angle_left[rows_step_4]
angle_right_step_4 = angle_right[rows_step_4]
angle_step_4 = np.concatenate((angle_left_step_4,angle_right_step_4))
mean_curvature_step_4 = mean_curvature[rows_step_4]
pc_step_4 = pc[rows_step_4]
r_curv_step_4 = r_curv[rows_step_4]
#
angle_left_step_5 = angle_left[rows_step_5]
angle_right_step_5 = angle_right[rows_step_5]
angle_step_5 = np.concatenate((angle_left_step_5,angle_right_step_5))
mean_curvature_step_5 = mean_curvature[rows_step_5]
pc_step_5 = pc[rows_step_5]
r_curv_step_5 = r_curv[rows_step_5]
#
angle_left_step_6 = angle_left[rows_step_6]
angle_right_step_6 = angle_right[rows_step_6]
angle_step_6 = np.concatenate((angle_left_step_6,angle_right_step_6))
mean_curvature_step_6 = mean_curvature[rows_step_6]
pc_step_6 = pc[rows_step_6]
r_curv_step_6 = r_curv[rows_step_6]
#
angle_left_step_7 = angle_left[rows_step_7]
angle_right_step_7 = angle_right[rows_step_7]
angle_step_7 = np.concatenate((angle_left_step_7,angle_right_step_7))
mean_curvature_step_7 = mean_curvature[rows_step_7]
pc_step_7 = pc[rows_step_7]
r_curv_step_7 = r_curv[rows_step_7]
#
angle_left_step_8 = angle_left[rows_step_8]
angle_right_step_8 = angle_right[rows_step_8]
angle_step_8 = np.concatenate((angle_left_step_8,angle_right_step_8))
mean_curvature_step_8 = mean_curvature[rows_step_8]
pc_step_8 = pc[rows_step_8]
r_curv_step_8 = r_curv[rows_step_8]
#
angle_left_step_9 = angle_left[rows_step_9]
angle_right_step_9 = angle_right[rows_step_9]
angle_step_9 = np.concatenate((angle_left_step_9,angle_right_step_9))
mean_curvature_step_9 = mean_curvature[rows_step_9]
pc_step_9 = pc[rows_step_9]
r_curv_step_9 = r_curv[rows_step_9]
#
angle_left_step_10 = angle_left[rows_step_10]
angle_right_step_10 = angle_right[rows_step_10]
angle_step_10 = np.concatenate((angle_left_step_10,angle_right_step_10))
mean_curvature_step_10 = mean_curvature[rows_step_10]
pc_step_10 = pc[rows_step_10]
r_curv_step_10 = r_curv[rows_step_10]
#
angle_left_step_11 = angle_left[rows_step_11]
angle_right_step_11 = angle_right[rows_step_11]
angle_step_11 = np.concatenate((angle_left_step_11,angle_right_step_11))
mean_curvature_step_11 = mean_curvature[rows_step_11]
pc_step_11 = pc[rows_step_11]
r_curv_step_11 = r_curv[rows_step_11]
#
angle_left_step_12 = angle_left[rows_step_12]
angle_right_step_12 = angle_right[rows_step_12]
angle_step_12 = np.concatenate((angle_left_step_12,angle_right_step_12))
mean_curvature_step_12 = mean_curvature[rows_step_12]
pc_step_12 = pc[rows_step_12]
r_curv_step_12 = r_curv[rows_step_12]
#
angle_left_step_13 = angle_left[rows_step_13]
angle_right_step_13 = angle_right[rows_step_13]
angle_step_13 = np.concatenate((angle_left_step_13,angle_right_step_13))
mean_curvature_step_13 = mean_curvature[rows_step_13]
pc_step_13 = pc[rows_step_13]
r_curv_step_13 = r_curv[rows_step_13]
#
angle_left_step_14 = angle_left[rows_step_14]
angle_right_step_14 = angle_right[rows_step_14]
angle_step_14 = np.concatenate((angle_left_step_14,angle_right_step_14))
mean_curvature_step_14 = mean_curvature[rows_step_14]
pc_step_14 = pc[rows_step_14]
r_curv_step_14 = r_curv[rows_step_14]
#
angle_left_step_15 = angle_left[rows_step_15]
angle_right_step_15 = angle_right[rows_step_15]
angle_step_15 = np.concatenate((angle_left_step_15,angle_right_step_15))
mean_curvature_step_15 = mean_curvature[rows_step_15]
pc_step_15 = pc[rows_step_15]
r_curv_step_15 = r_curv[rows_step_15]
#
angle_left_step_16 = angle_left[rows_step_16]
angle_right_step_16 = angle_right[rows_step_16]
angle_step_16 = np.concatenate((angle_left_step_16,angle_right_step_16))
mean_curvature_step_16 = mean_curvature[rows_step_16]
pc_step_16 = pc[rows_step_16]
r_curv_step_16 = r_curv[rows_step_16]
#
angle_left_step_17 = angle_left[rows_step_17]
angle_right_step_17 = angle_right[rows_step_17]
angle_step_17 = np.concatenate((angle_left_step_17,angle_right_step_17))
mean_curvature_step_17 = mean_curvature[rows_step_17]
pc_step_17 = pc[rows_step_17]
r_curv_step_17 = r_curv[rows_step_17]
#
angle_left_step_18 = angle_left[rows_step_18]
angle_right_step_18 = angle_right[rows_step_18]
angle_step_18 = np.concatenate((angle_left_step_18,angle_right_step_18))
mean_curvature_step_18 = mean_curvature[rows_step_18]
pc_step_18 = pc[rows_step_18]
r_curv_step_18 = r_curv[rows_step_18]
#
angle_left_step_19 = angle_left[rows_step_19]
angle_right_step_19 = angle_right[rows_step_19]
angle_step_19 = np.concatenate((angle_left_step_19,angle_right_step_19))
mean_curvature_step_19 = mean_curvature[rows_step_19]
pc_step_19 = pc[rows_step_19]
r_curv_step_19 = r_curv[rows_step_19]

angle_mean_vect = [np.nanmean(angle_step_1),np.nanmean(angle_step_2),np.nanmean(angle_step_3),np.nanmean(angle_step_4),np.nanmean(angle_step_5),
                    np.nanmean(angle_step_6),np.nanmean(angle_step_7),np.nanmean(angle_step_8),np.nanmean(angle_step_9),np.nanmean(angle_step_10),
                    np.nanmean(angle_step_11),np.nanmean(angle_step_12),np.nanmean(angle_step_13),np.nanmean(angle_step_14),np.nanmean(angle_step_15),
                    np.nanmean(angle_step_16),np.nanmean(angle_step_17),np.nanmean(angle_step_18),np.nanmean(angle_step_19)]
#
mean_curvature_mean_vect = [np.nanmean(mean_curvature_step_1),np.nanmean(mean_curvature_step_2),np.nanmean(mean_curvature_step_3),
                            np.nanmean(mean_curvature_step_4),np.nanmean(mean_curvature_step_5),np.nanmean(mean_curvature_step_6),
                            np.nanmean(mean_curvature_step_7),np.nanmean(mean_curvature_step_8),np.nanmean(mean_curvature_step_9),
                            np.nanmean(mean_curvature_step_10),np.nanmean(mean_curvature_step_11),np.nanmean(mean_curvature_step_12),
                            np.nanmean(mean_curvature_step_13),np.nanmean(mean_curvature_step_14),np.nanmean(mean_curvature_step_15),
                            np.nanmean(mean_curvature_step_16),np.nanmean(mean_curvature_step_17),np.nanmean(mean_curvature_step_18),
                            np.nanmean(mean_curvature_step_19)]
#
pc_mean_vect = [np.nanmean(pc_step_1),np.nanmean(pc_step_2),np.nanmean(pc_step_3),np.nanmean(pc_step_4),np.nanmean(pc_step_5),
                  np.nanmean(pc_step_6),np.nanmean(pc_step_7),np.nanmean(pc_step_8),np.nanmean(pc_step_9),
                  np.nanmean(pc_step_10),np.nanmean(pc_step_11),np.nanmean(pc_step_12),np.nanmean(pc_step_13),
                  np.nanmean(pc_step_14),np.nanmean(pc_step_15),np.nanmean(pc_step_16),np.nanmean(pc_step_17),
                  np.nanmean(pc_step_18),np.nanmean(pc_step_19)]

r_curv_mean_vect = [np.nanmean(r_curv_step_1),np.nanmean(r_curv_step_2),np.nanmean(r_curv_step_3),np.nanmean(r_curv_step_4),
                      np.nanmean(r_curv_step_5),np.nanmean(r_curv_step_6),np.nanmean(r_curv_step_7),np.nanmean(r_curv_step_8),
                      np.nanmean(r_curv_step_9),np.nanmean(r_curv_step_10),np.nanmean(r_curv_step_11),np.nanmean(r_curv_step_12),
                      np.nanmean(r_curv_step_13),np.nanmean(r_curv_step_14),np.nanmean(r_curv_step_15),np.nanmean(r_curv_step_16),
                      np.nanmean(r_curv_step_17),np.nanmean(r_curv_step_18),np.nanmean(r_curv_step_19)]
#
step_vect = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]

# Contact angle statistics:
contact_angles_drainage = np.concatenate((angle_left[rows_drainage],angle_right[rows_drainage]))
contact_angles_imbibition = np.concatenate((angle_left[rows_imbibition],angle_right[rows_imbibition]))
#


cond_drainage_left = np.where(np.logical_and(~np.isnan(angle_left[rows_drainage]),~np.isnan(r_curv[rows_drainage])))
cond_drainage_right = np.where(np.logical_and(~np.isnan(angle_right[rows_drainage]),~np.isnan(r_curv[rows_drainage])))


# For Curve fitting:
cond_imbibition_left = np.where(np.logical_and(~np.isnan(angle_left[rows_imbibition]),~np.isnan(r_curv[rows_imbibition])))
cond_imbibition_right = np.where(np.logical_and(~np.isnan(angle_right[rows_imbibition]),~np.isnan(r_curv[rows_imbibition])))


contact_angles_drainage_not_nan = np.concatenate((angle_left[rows_drainage][cond_drainage_left],angle_right[rows_drainage][cond_drainage_right]))
contact_angles_imbibition_not_nan = np.concatenate((angle_left[rows_imbibition][cond_imbibition_left],angle_right[rows_imbibition][cond_imbibition_right]))

pc_drainage_not_nan = np.concatenate((pc[rows_drainage][cond_drainage_left],pc[rows_drainage][cond_drainage_right]))
pc_imbibition_not_nan = np.concatenate((pc[rows_imbibition][cond_imbibition_left],pc[rows_imbibition][cond_imbibition_right]))

r_curv_drainage_not_nan = np.concatenate((r_curv[rows_drainage][cond_drainage_left],r_curv[rows_drainage][cond_drainage_right]))
r_curv_imbibition_not_nan = np.concatenate((r_curv[rows_imbibition][cond_imbibition_left],r_curv[rows_imbibition][cond_imbibition_right]))


# -------------------------------


mean_contact_angle_drainage = np.nanmean(contact_angles_drainage)
mean_contact_angle_imbibition = np.nanmean(contact_angles_imbibition)

print('==============================')
print('Mean drainage contact angle: ', mean_contact_angle_drainage)
print('Number of measured drainage contact angles (not nan): ', len(contact_angles_drainage[~np.isnan(contact_angles_drainage)]))
print('Mean imbibition contact angle: ', mean_contact_angle_imbibition)
print('Number of measured imbibition contact angles (not nan): ', len(contact_angles_imbibition[~np.isnan(contact_angles_imbibition)]))
print('==============================')

          
# ===============================================
# Histograms:

# Remove nan-entries for histograms:
contact_angles_drainage_hist = contact_angles_drainage[~np.isnan(contact_angles_drainage)]
contact_angles_imbibition_hist = contact_angles_imbibition[~np.isnan(contact_angles_imbibition)]


weights_drainage = np.ones_like(contact_angles_drainage_hist) / (len(contact_angles_drainage_hist))
weights_imbibition = np.ones_like(contact_angles_imbibition_hist) / (len(contact_angles_imbibition_hist))


binwidth = 5

my_bins=np.arange(15, 115 + binwidth, binwidth)

# --------------------------------------------------------
          
plot1 = plt.figure(1)
#
ax = plt.gca()
ax.patch.set_facecolor('white')


plt.plot(step[rows_drainage],angle_left[rows_drainage],
                label = ('Drainage'),
                linestyle='',
                color = 'gray',
                marker = 'o',
                markerfacecolor = 'gray',
                markeredgecolor = 'gray',
                markersize = 5.0,
                linewidth = 1.0)

plt.plot(step[rows_drainage],angle_right[rows_drainage],
                #
                linestyle='',
                color = 'gray',
                marker = 'o',
                markerfacecolor = 'gray',
                markeredgecolor = 'gray',
                markersize = 5.0,
                linewidth = 1.0)
# -----


plt.plot(step[rows_imbibition],angle_left[rows_imbibition],
                label = ('Imbibition'),
                linestyle='',
                color = 'b',
                marker = 's',
                markerfacecolor = 'b',
                markeredgecolor = 'b',
                markersize = 4.0,
                linewidth = 1.0)

plt.plot(step[rows_imbibition],angle_right[rows_imbibition],
                #
                linestyle='',
                color = 'b',
                marker = 's',
                markerfacecolor = 'b',
                markeredgecolor = 'b',
                markersize = 4.0,
                linewidth = 1.0)
# ======
plt.plot(step_vect,angle_mean_vect,
                label = ('Mean'),
                linestyle = '-',
                color = 'k',
                marker = 's',
                markerfacecolor = 'k',
                markeredgecolor = 'k',
                markersize = 4.0,
                linewidth = 1.0)         
          
# -----------------------
plt.xlabel('Hydraulic step [-]')
plt.ylabel(r'Contact angle $\theta$ [$^\circ$]')
plt.legend(loc='lower right',facecolor='white', framealpha=1, ncol=2)
plt.grid(True,color='k',linestyle='--')
plt.xlim((0, 20))
plt.ylim((0, 120))
locs, labels = plt.xticks()            
plt.xticks([0.,5.,10.,15.,20.], [0,5,10,15,20],)  

ax.xaxis.set_minor_locator(MultipleLocator(1))
ax.yaxis.set_minor_locator(MultipleLocator(5))

plot1.set_size_inches(cm_x/2.54,cm_y/2.54)
# ===============================================
plot1.show()
# ===============================================
plot2 = plt.figure(2)
ax = plt.subplot(111)
ax.patch.set_facecolor('white')

plt.hist(contact_angles_drainage_hist,bins=my_bins,histtype='step',weights=weights_drainage,label='Drainage',
         linewidth=2.0)
# ---
plt.hist(contact_angles_imbibition_hist,bins=my_bins,histtype='step',weights=weights_imbibition,label='Imbibition',
         linewidth=2.0)

box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1,
                 box.width, box.height * 0.9])


plt.xlabel(r'Contact angle $\theta$ [$^\circ$]')
plt.ylabel('Probability [-]')
plt.legend(loc = 'upper right', facecolor='white', framealpha=1)
plt.grid(True,color='k',linestyle='--')

plot2.set_size_inches(cm_x/2.54,cm_x/2.54)
plot2.show()
# ===============================================
plot3 = plt.figure(3)
#
ax = plt.gca()
ax.patch.set_facecolor('white')


plt.plot(step[rows_drainage],pc[rows_drainage],
                label = ('Drainage'),
                linestyle='',
                color = 'gray',
                marker = 'o',
                markerfacecolor = 'gray',
                markeredgecolor = 'gray',
                markersize = 5.0,
                linewidth = 1.0)
# -----


plt.plot(step[rows_imbibition],pc[rows_imbibition],
                label = ('Imbibition'),
                linestyle='',
                color = 'b',
                marker = 's',
                markerfacecolor = 'b',
                markeredgecolor = 'b',
                markersize = 4.0,
                linewidth = 1.0)
# =====
plt.plot(step_vect,pc_mean_vect,
                label = ('Mean'),
                linestyle='-',
                color = 'k',
                marker = 's',
                markerfacecolor = 'k',
                markeredgecolor = 'k',
                markersize = 4.0,
                linewidth = 1.0)         

# =====
plt.plot(step_suction_vect,suction_before_scans_vect,
                label = ('Start of CT scan'),
                linestyle='--',
                color = 'k',
                marker = 'o',
                markerfacecolor = 'None',
                markeredgecolor = 'k',
                markersize = 5.0,
                linewidth = 1.0)

plt.plot(step_suction_vect,suction_after_scans_vect,
                label = ('End of CT scan'),
                linestyle='--',
                color = 'k',
                marker = 's',
                markerfacecolor = 'None',
                markeredgecolor = 'k',
                markersize = 5.0,
                linewidth = 1.0)
          
# -----------------------
plt.xlabel('Hydraulic step [-]')
plt.ylabel('Capillary pressure $p_{\mathrm{c}}$ [kPa]')
plt.legend(facecolor='white', framealpha=1)
plt.grid(True,color='k',linestyle='--')

plt.xlim((0, 20))
plt.ylim((-1.0, 1.75))
locs, labels = plt.xticks()            
plt.xticks([0.,5.,10.,15.,20.], [0,5,10,15,20],)  
ax.xaxis.set_minor_locator(MultipleLocator(1))
ax.yaxis.set_minor_locator(MultipleLocator(0.1))

plot3.set_size_inches(cm_x/2.54,cm_y/2.54)
plot3.show()
# ===============================================
# -------------------------------
# Curve fitting:

# Sample data
x_drainage = contact_angles_drainage_not_nan 
y_drainage = pc_drainage_not_nan
#
x_imbibition = contact_angles_imbibition_not_nan 
y_imbibition = pc_imbibition_not_nan

plot4 = plt.figure(4)
#
ax = plt.gca()
ax.patch.set_facecolor('white')

plt.plot(x_drainage, y_drainage, 'o',markersize = 5.0, markerfacecolor = 'gray', markeredgecolor = 'gray', label = 'Drainage')
#
plt.plot(x_imbibition, y_imbibition, 's',markersize = 4.0, markerfacecolor = 'blue', markeredgecolor = 'blue', label = 'Imbibition')
#
# For unsorted data:
plt.plot(np.unique(x_drainage), np.poly1d(np.polyfit(x_drainage, y_drainage, 1))(np.unique(x_drainage)),'k-',
         linewidth=1.0,
         label ='Drainage, linear fit')
#
plt.plot(np.unique(x_imbibition), np.poly1d(np.polyfit(x_imbibition, y_imbibition, 1))(np.unique(x_imbibition)),'k--',
         linewidth=1.0,
         label ='Imbibition, linear fit')


# -----------------------
plt.xlabel(r'Contact angle $\theta$ [$^\circ$]')
plt.ylabel('Capillary pressure $p_{\mathrm{c}}$ [kPa]')
plt.legend(facecolor='white', framealpha=1)
plt.grid(True,color='k',linestyle='--')
plot4.set_size_inches(cm_x/2.54,cm_y/2.54)

plot4.show()

# --------------------------------------------------------
          
plot5 = plt.figure(5)
#
ax = plt.gca()
ax.patch.set_facecolor('white')


plt.plot(step[rows_drainage],r_curv[rows_drainage],
                label = ('Drainage'),
                linestyle='',
                color = 'gray',
                marker = 'o',
                markerfacecolor = 'gray',
                markeredgecolor = 'gray',
                markersize = 5.0,
                linewidth = 1.0)
# -----


plt.plot(step[rows_imbibition],r_curv[rows_imbibition],
                label = ('Imbibition'),
                linestyle='',
                color = 'b',
                marker = 's',
                markerfacecolor = 'b',
                markeredgecolor = 'b',
                markersize = 4.0,
                linewidth = 1.0)
# -----
plt.plot(step_vect,r_curv_mean_vect,
                label = ('Mean'),
                linestyle='-',
                color = 'k',
                marker = 's',
                markerfacecolor = 'k',
                markeredgecolor = 'k',
                markersize = 4.0,
                linewidth = 1.0)
          
# -----------------------
plt.xlabel('Hydraulic step [-]')
plt.ylabel('Radius of curvature $R_{1}$ [mm]')
plt.legend(loc='lower right',facecolor='white', framealpha=1, ncol=2)
plt.grid(True,color='k',linestyle='--')
# ===============================================
plt.xlim((0, 20))

locs, labels = plt.xticks()            
plt.xticks([0.,5.,10.,15.,20.], [0,5,10,15,20],)  
#
ax.xaxis.set_minor_locator(MultipleLocator(1))
ax.yaxis.set_minor_locator(MultipleLocator(0.02))

plot5.set_size_inches(cm_x/2.54,cm_y/2.54)
plot5.show()
# ===============================================


# -------------------------------
# Curve fitting:

# Sample data
x1_drainage = contact_angles_drainage_not_nan 
y1_drainage = r_curv_drainage_not_nan
#
x1_imbibition = contact_angles_imbibition_not_nan 
y1_imbibition = r_curv_imbibition_not_nan

plot6 = plt.figure(6)
#
ax = plt.gca()
ax.patch.set_facecolor('white')

plt.plot(x1_drainage, y1_drainage, 'o',markersize = 5.0, markerfacecolor = 'gray', markeredgecolor = 'gray', label = 'Drainage')
#
plt.plot(x1_imbibition, y1_imbibition, 's',markersize = 4.0, markerfacecolor = 'blue', markeredgecolor = 'blue', label = 'Imbibition')
#
# For unsorted data:
plt.plot(np.unique(x1_drainage), np.poly1d(np.polyfit(x1_drainage, y1_drainage, 1))(np.unique(x1_drainage)),'k-',
         linewidth=1.0,
         label ='Drainage, linear fit')
#
plt.plot(np.unique(x1_imbibition), np.poly1d(np.polyfit(x1_imbibition, y1_imbibition, 1))(np.unique(x1_imbibition)),'k--',
         linewidth=1.0,
         label ='Imbibition, linear fit')


# -----------------------
plt.xlabel(r'Contact angle $\theta$ [$^\circ$]')
plt.ylabel('Radius of curvature $R_{1}$ [mm]')
plt.legend(facecolor='white', framealpha=1)
plt.grid(True,color='k',linestyle='--')
plot6.set_size_inches(cm_x/2.54,cm_y/2.54)

plot6.show()

# ==========================================================================
# Save figures to file if savefigs == True:

if savefigs == True:
    print("Saving figures to file...")
    plot1.savefig('contact_angles_vs_steps_CT.pdf', bbox_inches='tight')
    #
    plot2.savefig('histogram_of_contact_angles_CT.pdf', bbox_inches='tight')
    #
    plot3.savefig('capillary_pressure_vs_steps_CT.pdf', bbox_inches='tight')
    #
    plot4.savefig('capillary_pressure_vs_contact_angle_regression_CT.pdf', bbox_inches='tight')
    #
    plot5.savefig('radius_of_curvature_vs_steps_CT.pdf', bbox_inches='tight') 
    #  
    plot6.savefig('radius_of_curvature_vs_contact_angle_CT.pdf', bbox_inches='tight')

    print('Done.')

else:
    print('Figures not saved to file.')
raise SystemExit