diff --git a/examples/script_design_capacity.py b/examples/script_design_capacity.py
new file mode 100644
index 0000000000000000000000000000000000000000..435c9473bd59fa90415c11f2240e62ae3b0a717a
--- /dev/null
+++ b/examples/script_design_capacity.py
@@ -0,0 +1,576 @@
+#******************************************************************************
+#******************************************************************************
+
+# TIP: this script requires the data found in the hyhetra-pipedata and hyhetra-
+# -fluiddata repositories
+
+pipedata_files = ['pipes/isoplus_singlepipes_s1.csv']
+
+fluiddata_file = 'fluids/incropera2006_saturated_water.csv'
+
+# other files
+
+file_best2018 = 'papers/Best2018_design_flow_speeds.csv'
+
+file_sven2013 = 'papers/Sven2013_design_flow_speeds.csv'
+
+file_nussbaumer2016 = 'papers/Nussbaumer2016_design_flow_speeds.csv'
+
+file_euroheat2021 = 'papers/EuroHeat2021_design_flow_speeds.csv'
+
+# import libraries
+import numpy as np
+import matplotlib.pyplot as plt
+# topupheat
+import topupheat.pipes.fic as fic
+from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
+import topupheat.pipes.utils as utils
+from topupheat.common import formulas as core
+
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+
+# specifications
+
+# maximum specific pressure drop
+
+list_max_specific_pressure_loss = [100, 200, 200, 200, 450] # Pa/m
+
+list_pipe_roughness = [0.1/1000, 0.1/1000, 0.05/1000, 0.01/1000, 0.1/1000]
+
+# district heating details
+
+dh_flow_temperature = 273.15+80 # K
+
+dh_return_temperature = 273.15+50 # K
+
+# temperature fluid bulk
+
+temperature_fluid_bulk = 0.5*(dh_return_temperature+dh_flow_temperature)
+
+# ground temperature
+
+ground_temperature = 273.15+10 # K
+
+# pipe absolute/effective roughness
+
+pipe_e_eff = 0.01/1000 # m
+
+#******************************************************************************
+
+# specify the pipes that should be considered using (DN,S) tuples
+
+list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
+ (25, 1, 'isoplus', 'DRE-25-STD'),
+ (32, 1, 'isoplus', 'DRE-32-STD'),
+ (40, 1, 'isoplus', 'DRE-40-STD'),
+ (50, 1, 'isoplus', 'DRE-50-STD'),
+ (65, 1, 'isoplus', 'DRE-65-STD'),
+ (80, 1, 'isoplus', 'DRE-80-STD'),
+ (100, 1, 'isoplus', 'DRE-100-ST'),
+ (125, 1, 'isoplus', 'DRE-125-ST'),
+ (150, 1, 'isoplus', 'DRE-150-ST'),
+ (200, 1, 'isoplus', 'DRE-200-ST'),
+ (250, 1, 'isoplus', 'DRE-250-ST'),
+ (300, 1, 'isoplus', 'DRE-300-ST'),
+ (350, 1, 'isoplus', 'DRE-350-ST'),
+ (400, 1, 'isoplus', 'DRE-400-ST'),
+ (450, 1, 'isoplus', 'DRE-450-ST'),
+ (500, 1, 'isoplus', 'DRE-500-ST'),
+ (600, 1, 'isoplus', 'DRE-600-ST'),
+ (700, 1, 'isoplus', 'DRE-700-ST'),
+ (800, 1, 'isoplus', 'DRE-800-ST'),
+ (900, 1, 'isoplus', 'DRE-900-ST'),
+ (1000, 1, 'isoplus', 'DRE-1000-S')]
+
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+
+# pipe data
+
+pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+# water data
+
+waterdb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+
+#******************************************************************************
+#******************************************************************************
+
+# create the objects for the fluids in the flow and return pipes
+
+dh_flow = fic.Fluid(
+ phase='l',
+ temperature=dh_flow_temperature,
+ pressure=1e5,
+ db=waterdb)
+
+dh_return = fic.Fluid(
+ phase='l',
+ temperature=dh_return_temperature,
+ pressure=1e5,
+ db=waterdb)
+
+fluid_bulk = fic.Fluid(
+ phase='l',
+ temperature=temperature_fluid_bulk,
+ pressure=1e5,
+ db=waterdb)
+
+#******************************************************************************
+#******************************************************************************
+
+# define the pipes under consideration for task 2
+
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+
+# 4) compute the maximum heat flow rate for a given maximum specific pressure drop
+# as a function of pipe diameters
+
+#******************************************************************************
+
+# dictionaries
+
+dict_heat_transfer_rate = {}
+
+dict_friction_factor = {}
+
+dict_fluid_speed = {}
+
+dict_reynolds_number = {}
+
+dict_specific_pressure_loss = {}
+
+# for each specific pressure loss level
+
+for j, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
+
+ # declare the final and intermediate variables
+
+ list_heat_transfer_rate = []
+
+ list_friction_factor = []
+
+ list_fluid_speed = []
+
+ list_reynolds_number = []
+
+ list_specific_pressure_loss = []
+
+ temp_fluid_speed, temp_reynolds_number, temp_friction_factor = 0, 0, 0
+
+ # pipe object list
+
+ list_dh_pipes = [
+ StandardisedPipe(
+ pipe_tuple=pipe_tuple,
+ e_eff=list_pipe_roughness[j],
+ db=pipedb)
+ for pipe_tuple in list_pipe_tuples]
+
+ # for each pipe under consideration
+
+ for i, pipe in enumerate(list_dh_pipes):
+
+ #**********************************************************************
+
+ # get an initial solution
+ # >> should be valid and as close to the actual solution as possible
+ # >> one option is to use explicit methods
+
+ if i == 0:
+
+ x0_fluid_speed = None
+ x0_reynolds_number = None
+ x0_friction_factor = None
+
+ else:
+
+ # use previous solution
+
+ (x0_fluid_speed,
+ x0_reynolds_number,
+ x0_friction_factor) = (list_fluid_speed[i-1],
+ list_reynolds_number[i-1],
+ list_friction_factor[i-1])
+
+ # flow regime
+
+ # compute the flow speed that leads to the specified pressure loss
+
+ temp_fluid_speed, temp_reynolds_number, temp_friction_factor = (
+ utils.find_specific_pressure_loss_matching_flow(
+ pipe,
+ fluid_bulk,
+ max_specific_pressure_loss,
+ dff_model=fic.DFF_COLEBROOK_WHITE,
+ #dff_model=fic.DFF_HAALAND,
+ #dff_model=fic.DFF_NIKURADSE,
+ #dff_model=fic.DFF_PETHUKOV,
+ #dff_model=fic.DFF_BLASIUS,
+ #dff_model=fic.DFF_PRANDTL_NIKURADSE_KARMAN,
+ x0_friction_factor=x0_friction_factor,
+ x0_reynolds_number=x0_reynolds_number,
+ x0_fluid_speed=x0_fluid_speed)
+ )
+
+ temp_specific_pressure_loss = fic.SpecificPressureLossInPipe(
+ fluid_bulk.mass_density,
+ pipe.d_int,
+ temp_fluid_speed,
+ temp_friction_factor)
+
+ # store the results
+
+ list_fluid_speed.append(temp_fluid_speed)
+
+ list_reynolds_number.append(temp_reynolds_number)
+
+ list_friction_factor.append(temp_friction_factor)
+
+ list_specific_pressure_loss.append(temp_specific_pressure_loss)
+
+ # compute the heat flow rate
+
+ temp_heat_flow_rate = fic.HeatTransferRateInPipe(
+ temp_fluid_speed,
+ dh_flow,
+ dh_return,
+ pipe,
+ fluid_bulk=fluid_bulk
+ )
+
+ # store the result
+
+ list_heat_transfer_rate.append(temp_heat_flow_rate)
+
+ #**************************************************************************
+
+ # add lists to the dicts
+
+ dict_specific_pressure_loss.update(
+ {(max_specific_pressure_loss,list_pipe_roughness[j]):
+ list_specific_pressure_loss})
+
+ dict_fluid_speed.update(
+ {(max_specific_pressure_loss,list_pipe_roughness[j]):
+ list_fluid_speed})
+
+ dict_reynolds_number.update(
+ {(max_specific_pressure_loss,list_pipe_roughness[j]):
+ list_reynolds_number})
+
+ dict_friction_factor.update(
+ {(max_specific_pressure_loss,list_pipe_roughness[j]):
+ list_friction_factor})
+
+ dict_heat_transfer_rate.update(
+ {(max_specific_pressure_loss,list_pipe_roughness[j]):
+ list_heat_transfer_rate})
+
+ #**************************************************************************
+
+ # print('******************************************************************')
+ # print('******************************************************************')
+ # print('******************************************************************')
+ # print('******************************************************************')
+
+#******************************************************************************
+
+isoplus_diameters = [20,
+ 25,
+ 32,
+ 40,
+ 50,
+ 65,
+ 80,
+ 100,
+ 125,
+ 150,
+ 200,
+ 250,
+ 300,
+ 350,
+ 400,
+ 450,
+ 500,
+ 600,
+ 700,
+ 800,
+ 900,
+ 1000]
+
+isoplus_mass_flow_rate_min_t_hour = [0.4,
+ 0.8,
+ 1.7,
+ 2.5,
+ 4.7,
+ 9.3,
+ 14.5,
+ 28.5,
+ 50.0,
+ 82.0,
+ 167.0,
+ 300,
+ 472,
+ 610,
+ 862,
+ 1180,
+ 1570,
+ 2520,
+ 3770,
+ 5390,
+ 7400,
+ 9200]
+
+isoplus_mass_flow_rate_max_t_hour = [0.5,
+ 1.0,
+ 2.0,
+ 3.0,
+ 5.5,
+ 11.0,
+ 16.5,
+ 33.0,
+ 58.0,
+ 95.0,
+ 193.0,
+ 348,
+ 547,
+ 705,
+ 1000,
+ 1370,
+ 1820,
+ 2920,
+ 4370,
+ 6240,
+ 9500,
+ np.nan]
+
+isoplus_fluid_speed_min = [
+ core.MeanFluidSpeedViaVolumetricFlowRate(
+ core.VolumetricFlowRateFromMassFlowRate(
+ m*0.27777778, #t/h to kg/s
+ fluid_bulk.mass_density
+ ),
+ fic.Pipe.CrossSectionalArea(fic.Pipe, isoplus_diameters[i]/1e3)
+ )
+ for i, m in enumerate(isoplus_mass_flow_rate_min_t_hour)]
+
+isoplus_fluid_speed_max = [
+ core.MeanFluidSpeedViaVolumetricFlowRate(
+ core.VolumetricFlowRateFromMassFlowRate(
+ m*0.27777778, #t/h to kg/s
+ fluid_bulk.mass_density
+ ),
+ fic.Pipe.CrossSectionalArea(fic.Pipe, isoplus_diameters[i]/1e3)
+ )
+ for i, m in enumerate(isoplus_mass_flow_rate_max_t_hour)]
+
+# get data from Best et al. (2018)
+
+data_best2018 = np.genfromtxt(file_best2018,
+ #dtype=(float,float,float,float),
+ names=True,
+ # skip_header=1,
+ # skip_footer=0,
+ delimiter=',')
+
+data_sven2013 = np.genfromtxt(file_sven2013,
+ #dtype=(float,float,float,float),
+ names=True,
+ # skip_header=1,
+ # skip_footer=0,
+ delimiter=',')
+
+data_nussbaumer2016 = np.genfromtxt(
+ file_nussbaumer2016,
+ #dtype=(float,float,float,float),
+ names=True,
+ # skip_header=1,
+ # skip_footer=0,
+ delimiter=',')
+
+data_euroheat2021 = np.genfromtxt(
+ file_euroheat2021,
+ #dtype=(float,float,float,float),
+ names=True,
+ # skip_header=1,
+ # skip_footer=0,
+ delimiter=',')
+
+# list_pipe_diameters = [t[0] for t in list_pipe_tuples]
+
+list_pipe_diameters = [pipe.d_int*1000 for pipe in list_dh_pipes]
+
+# plot the max heat flow rate as function of pipe diameter (1st pos. in tuple)
+
+fig, ax = plt.subplots()
+
+for i, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
+
+ ax.plot(list_pipe_diameters,
+ [-y/1e6 for y in dict_heat_transfer_rate[
+ (max_specific_pressure_loss,list_pipe_roughness[i])]],
+ label='MSP='+str(max_specific_pressure_loss)+' Pa/m')
+
+ax.set(xlabel='Pipe internal diameter [mm]',
+ ylabel='Heat Transfer Rate [MW]',
+ title='Maximum heat capacity')
+
+ax.legend()
+
+ax.grid()
+
+# fig.savefig("test.png")
+plt.show()
+
+#******************************************************************************
+
+# plot the max fluid speed as a function of pipe diameter (1st pos. in tuple)
+
+fig, ax = plt.subplots()
+
+fig.set_size_inches(12,8) # (15,10) (12,8) # (3,2) (3*2,2*2) (3*3,2*3)
+
+list_linestyle_our = ['rx-','bx-','gx-','cx-','yx-']
+
+list_linestyle_euroheat2021 = ['rs','bs','gs','cs','ys']
+
+list_linestyle_best2018 = ['r>--','b>--','g>--']
+
+list_linestyle_nussbaumer2016 = ['r<-.','b<-.','g<-.']
+
+published_results_labels_euroheat2021 = ['100_Pa_m_01mm',
+ '200_Pa_m_01mm',
+ '200_Pa_m_005mm',
+ '200_Pa_m_001mm',
+ '450_Pa_m_01mm']
+
+published_results_labels_best2018 = ['100_Pa_m','200_Pa_m','300_Pa_m']
+
+published_results_labels_nussbaumer2016 = ['100_Pa_m','200_Pa_m','300_Pa_m']
+
+legend_euroheat2021 = ['100 Pa/m, e=0.1mm',
+ '200 Pa/m, e=0.1mm',
+ '200 Pa/m, e=0.05mm',
+ '200 Pa/m, e=0.01mm',
+ '450 Pa/m, e=0.1mm']
+
+legend_best2018 = ['100 Pa/m, e=0.01mm',
+ '200 Pa/m, e=0.01mm',
+ '300 Pa/m, e=0.01mm']
+
+legend_nussbaumer2016 = ['100 Pa/m, e=0.01mm',
+ '200 Pa/m, e=0.01mm',
+ '300 Pa/m, e=0.01mm']
+
+# x,100_Pa_m_01mm,200_Pa_m_01mm,200_Pa_m_005mm,200_Pa_m_001mm,450_Pa_m_01mm
+
+for i, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
+
+ ax.plot(list_pipe_diameters,
+ [fluid_speed for fluid_speed in dict_fluid_speed[
+ (max_specific_pressure_loss,list_pipe_roughness[i])]],
+ #marker='s', linestyle='-',
+ list_linestyle_our[i],
+ label='topupheat, '+str(max_specific_pressure_loss)+' Pa/m, e='+str(list_pipe_roughness[i]*1e3)+' mm')
+
+for i, label in enumerate(published_results_labels_euroheat2021):
+
+ ax.plot(data_euroheat2021['x']*1e3,
+ data_euroheat2021[label],
+ #marker='x', linestyle='-.',
+ list_linestyle_euroheat2021[i],
+ label='EuroHeat&Power (2021), '+str(legend_euroheat2021[i]))
+
+# for i, label in enumerate(published_results_labels_best2018):
+
+# ax.plot(data_best2018['x']*1e3,
+# data_best2018[label],
+# #marker='x', linestyle='-.',
+# list_linestyle_best2018[i],
+# label='Best et al. (2018), '+str(legend_best2018[i]))
+
+# for i, label in enumerate(published_results_labels_nussbaumer2016):
+
+# ax.plot(data_nussbaumer2016['x']*1e3,
+# data_nussbaumer2016[label],
+# #marker='x', linestyle='-.',
+# list_linestyle_nussbaumer2016[i],
+# label='Nussbaumer et al. (2016), '+str(legend_nussbaumer2016[i]))
+
+# ax.plot(data_sven2013['x']*1e3,
+# data_sven2013['100_Pa_m'],
+# #marker='x', linestyle='-.',
+# 'kd-.',
+# label='Svend and Sven (2013), 100 Pa/m')
+
+#******************************************************************************
+ax.set(xlim=(0, 1000),ylim=(0, 9))
+ax.set(xlabel='Pipe internal diameter [mm]',
+ ylabel='Fluid speed [m/s]',
+ title='Maximum mean fluid speed for a given specific pressure loss')
+
+ax.legend()
+
+ax.grid()
+
+# fig.savefig("test.png")
+plt.show()
+
+#******************************************************************************
+
+# plot the specific pressure loss as a function of pipe diameter
+
+fig, ax = plt.subplots()
+
+for i, max_specific_pressure_loss in enumerate(list_max_specific_pressure_loss):
+
+ ax.plot(list_pipe_diameters,
+ [spl for spl in dict_specific_pressure_loss[
+ (max_specific_pressure_loss,list_pipe_roughness[i])]],
+ label='MSP='+str(max_specific_pressure_loss)+' Pa/m')
+
+ # ax.plot(data_best2018['x']*1e3,
+ # data_best2018[published_results_labels[i]],
+ # #marker='o', linestyle='--',
+ # list_linestyle_best2018[i],
+ # label='Best et al. (2018), '+str(max_specific_pressure_loss)+' Pa/m')
+
+ # ax.plot(data_nussbaumer2016['x']*1e3,
+ # data_nussbaumer2016[published_results_labels[i]],
+ # #marker='x', linestyle='-.',
+ # list_linestyle_nussbaumer2016[i],
+ # label='Nussbaumer et al. (2016), '+str(max_specific_pressure_loss)+' Pa/m')
+
+
+ax.set(xlabel='Pipe internal diameter [mm]',
+ ylabel='Specific pressure loss [Pa/m]',
+ title='...')
+
+ax.legend()
+
+ax.grid()
+
+# fig.savefig("test.png")
+plt.show()
+
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
+#******************************************************************************
\ No newline at end of file
diff --git a/examples/script_moody_diagram.py b/examples/script_moody_diagram.py
new file mode 100644
index 0000000000000000000000000000000000000000..76486e25563d9b07beb419f427abea5983435e8d
--- /dev/null
+++ b/examples/script_moody_diagram.py
@@ -0,0 +1,208 @@
+#******************************************************************************
+#******************************************************************************
+
+# TIP: this script requires the data found in the hyhetra-pipedata and hyhetra-
+# -fluiddata repositories
+
+pipedata_files = ['pipes/isoplus_singlepipes_s1.csv']
+
+fluiddata_file = 'fluids/incropera2006_saturated_water.csv'
+
+# import libraries
+import numpy as np
+import matplotlib.pyplot as plt
+import topupheat.pipes.fic as fic
+from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
+
+#******************************************************************************
+#******************************************************************************
+
+# specify the pipes that should be considered
+
+list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
+ (25, 1, 'isoplus', 'DRE-25-STD'),
+ (32, 1, 'isoplus', 'DRE-32-STD'),
+ (40, 1, 'isoplus', 'DRE-40-STD'),
+ (50, 1, 'isoplus', 'DRE-50-STD'),
+ (65, 1, 'isoplus', 'DRE-65-STD'),
+ (80, 1, 'isoplus', 'DRE-80-STD'),
+ (100, 1, 'isoplus', 'DRE-100-ST'),
+ (125, 1, 'isoplus', 'DRE-125-ST'),
+ (150, 1, 'isoplus', 'DRE-150-ST'),
+ (200, 1, 'isoplus', 'DRE-200-ST')]
+
+# relative pipe roughness
+
+list_relative_pipe_roughness = [
+ 5e-2,
+ 2.5e-2,
+ 1e-2,
+ 5e-3,
+ 2e-3,
+ 1e-3,
+ 5e-4,
+ 2e-4,
+ 1e-4,
+ 5e-5,
+ 2e-5,
+ 1e-5,
+ 0
+ ]
+
+labels_relative_pipe_roughness = [
+ '5e-2',
+ '2.5e-2',
+ '1e-2',
+ '5e-3',
+ '2e-3',
+ '1e-3',
+ '5e-4',
+ '2e-4',
+ '1e-4',
+ '5e-5',
+ '2e-5',
+ '1e-5',
+ '0']
+
+# specify the pipe that is to be considered individually
+
+selected_pipe_index = 5
+
+# moody diagram
+
+number_points_laminar = 5
+
+number_points_turbulent = 100
+
+#******************************************************************************
+#******************************************************************************
+
+# pipe data
+
+pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+# water data
+
+waterdb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+
+#******************************************************************************
+#******************************************************************************
+
+# 1) plot the moody diagram
+
+#******************************************************************************
+
+# define the reynolds numbers' vector
+
+list_reynolds_number = (
+ np.logspace(0,
+ 8,
+ number_points_turbulent+number_points_laminar)
+ )
+
+# define the models to be used
+
+list_dff_models = [fic.DFF_COLEBROOK_WHITE]
+
+#******************************************************************************
+#******************************************************************************
+
+# variables
+
+dict_friction_factor = {}
+
+dict_pipe = {}
+
+# for each roughness value
+
+for relative_pipe_roughness in list_relative_pipe_roughness:
+
+ # create the pipe instance
+
+ pipe = StandardisedPipe(
+ pipe_tuple=list_pipe_tuples[selected_pipe_index],
+ e_rel=relative_pipe_roughness,
+ db=pipedb)
+
+ dict_pipe[relative_pipe_roughness] = pipe
+
+ # for each dff model
+
+ for model in list_dff_models:
+
+ # get the friction factors
+
+ list_friction_factor = [fic.DarcyFrictionFactor(
+ re,
+ pipe,
+ model=model)
+ for re in list_reynolds_number]
+
+ # store it
+
+ dict_friction_factor[
+ (model, relative_pipe_roughness)] = list_friction_factor
+
+ #**********************************************************************
+ #**********************************************************************
+
+ #**************************************************************************
+ #**************************************************************************
+
+#******************************************************************************
+#******************************************************************************
+
+list_linestyle_our = ['ro-','go-','bx-',
+ 'rv-','g^-','bd-',
+ 'r<-','g>-','b.-',
+ 'ro--','go--','bx--',
+ 'rv--','g^--','bd--',
+ 'r<--','g>--','b.--']
+
+# plot
+
+fig, ax = plt.subplots()
+
+fig.set_size_inches(12,8)
+
+for j, model in enumerate(list_dff_models):
+
+ for i, relative_pipe_roughness in enumerate(list_relative_pipe_roughness):
+
+ # plot them using a log scale
+
+ linestyle_index = i+(j)*len(list_relative_pipe_roughness)
+
+ ax.loglog(list_reynolds_number,
+ dict_friction_factor[(model, relative_pipe_roughness)],
+ list_linestyle_our[linestyle_index],
+ #label=model+', e/D='+labels_relative_pipe_roughness[i],
+ label='CWE, e/D='+labels_relative_pipe_roughness[i])
+
+ #**********************************************************************
+
+ #**************************************************************************
+
+#******************************************************************************
+#******************************************************************************
+
+ax.set(xlabel='Reynolds number [-]',
+ ylabel='Friction factor [-]',
+ title='Moody diagram')
+
+# plt.autoscale(enable=True, axis='x', tight=True)
+
+ax.grid(which='both')
+
+ax.legend()
+
+# left, right = xlim() # return the current xlim
+
+plt.xlim((0.6e3, 1e8)) # set the xlim to left, right
+plt.ylim((0.006, 0.1)) # set the xlim to left, right
+
+# fig.savefig("test.png")
+plt.show()
+
+#******************************************************************************
+#******************************************************************************
\ No newline at end of file
diff --git a/examples/script_nikuradse_experiments.py b/examples/script_nikuradse_experiments.py
new file mode 100644
index 0000000000000000000000000000000000000000..cae7390478e1deb07ef00369c4a2fc610a0aa84d
--- /dev/null
+++ b/examples/script_nikuradse_experiments.py
@@ -0,0 +1,240 @@
+#******************************************************************************
+#******************************************************************************
+
+# TIP: this script requires the data found in the hyhetra-pipedata and hyhetra-
+# -fluiddata repositories
+
+pipedata_files = ['pipes/isoplus_singlepipes_s1.csv']
+
+fluiddata_file = 'fluids/incropera2006_saturated_water.csv'
+
+# Note: this script uses data published by Li and Huai (2016) about Nikuradse's
+# experiments. DOI: https://doi.org/10.1371/journal.pone.0154408
+
+# The xls file included in the paper as an appendix should then be converted to
+# csv format and its two first lines (headers) replaced by the following one:
+# A_run,A_f,A_1/e,A_Re,B_run,B_f,B_1/e,B_Re,C_run,C_f,C_1/e,C_Re,D_run,D_f,D_1/e,D_Re,E_run,E_f,E_1/e,E_Re,F_run,F_f,F_1/e,F_Re
+# How? copy the line above (minus the leading '# ') and paste it using a text editor
+
+nikuradse_file = 'papers/Li2016_edited.csv'
+
+# import libraries
+import numpy as np
+import matplotlib.pyplot as plt
+import topupheat.pipes.fic as fic
+from topupheat.pipes.single import StandardisedPipe, StandardisedPipeDatabase
+
+#******************************************************************************
+#******************************************************************************
+
+# specify the pipes that should be considered using (DN,S) tuples
+
+# 20,25,32,40,50,65,80,100,125,150,200
+
+list_pipe_tuples = [(20, 1, 'isoplus', 'DRE-20-STD'),
+ (25, 1, 'isoplus', 'DRE-25-STD'),
+ (32, 1, 'isoplus', 'DRE-32-STD'),
+ (40, 1, 'isoplus', 'DRE-40-STD'),
+ (50, 1, 'isoplus', 'DRE-50-STD'),
+ (65, 1, 'isoplus', 'DRE-65-STD'),
+ (80, 1, 'isoplus', 'DRE-80-STD'),
+ (100, 1, 'isoplus', 'DRE-100-ST'),
+ (125, 1, 'isoplus', 'DRE-125-ST'),
+ (150, 1, 'isoplus', 'DRE-150-ST'),
+ (200, 1, 'isoplus', 'DRE-200-ST')]
+
+# relative pipe roughness
+
+list_relative_pipe_roughness = [
+ 0.5/507,
+ 0.5/252,
+ 0.5/126,
+ 0.5/60,
+ 0.5/30.6,
+ 0.5/15]
+
+labels_relative_pipe_roughness = [
+ '1/1014',
+ '1/504',
+ '1/252',
+ '1/120',
+ '1/61.2',
+ '1/30']
+
+# specify the pipe that is to be considered individually
+
+selected_pipe_index = 5
+
+# moody diagram
+
+number_points_laminar = 5
+
+number_points_turbulent = 100
+
+#******************************************************************************
+#******************************************************************************
+
+# pipe data
+
+pipedb = StandardisedPipeDatabase(source=pipedata_files)
+
+# water data
+
+waterdb = fic.FluidDatabase(fluid='water', phase='l', source=fluiddata_file)
+
+#******************************************************************************
+#******************************************************************************
+
+nikuradse_data = np.genfromtxt(nikuradse_file,
+ #dtype=(float,float,float,float),
+ names=True,
+ #skip_header=0,
+ # skip_header=1,
+ # skip_footer=0,
+ delimiter=',')
+
+nikuradse_experiments = ['A','B','C','D','E','F']
+
+nikuradse_data_length = [48, 38, 71, 72, 63, 70]
+
+nikuradse_data_f = {nikuradse_experiments[i]:
+ nikuradse_data[value+'_f'][0:nikuradse_data_length[i]]
+ for i, value in enumerate(nikuradse_experiments)}
+
+nikuradse_data_re = {nikuradse_experiments[i]:
+ nikuradse_data[value+'_Re'][0:nikuradse_data_length[i]]
+ for i, value in enumerate(nikuradse_experiments)}
+
+nikuradse_data_1e = {nikuradse_experiments[i]:
+ nikuradse_data[value+'_1e'][1]
+ for i, value in enumerate(nikuradse_experiments)}
+
+#******************************************************************************
+#******************************************************************************
+
+# 1) plot the moody diagram
+
+#******************************************************************************
+
+# define the reynolds numbers' vector
+
+list_reynolds_number = (
+ np.logspace(0,
+ 8,
+ number_points_turbulent+number_points_laminar)
+ )
+
+# define the models to be used
+
+list_dff_models = [fic.DFF_COLEBROOK_WHITE]
+
+#******************************************************************************
+#******************************************************************************
+
+# variables
+
+dict_friction_factor = {}
+
+dict_pipe = {}
+
+# for each roughness value
+
+for relative_pipe_roughness in list_relative_pipe_roughness:
+
+ # create the pipe instance
+
+ pipe = StandardisedPipe(
+ pipe_tuple=list_pipe_tuples[selected_pipe_index],
+ e_rel=relative_pipe_roughness,
+ db=pipedb)
+
+ dict_pipe[relative_pipe_roughness] = pipe
+
+ # for each dff model
+
+ for model in list_dff_models:
+
+ # get the friction factors
+
+ list_friction_factor = [fic.DarcyFrictionFactor(
+ re,
+ pipe,
+ model=model)
+ for re in list_reynolds_number]
+
+ # store it
+
+ dict_friction_factor[
+ (model, relative_pipe_roughness)] = list_friction_factor
+
+ #**********************************************************************
+ #**********************************************************************
+
+ #**************************************************************************
+ #**************************************************************************
+
+#******************************************************************************
+#******************************************************************************
+
+list_linestyle_our = ['rs-','rd-','rx-',
+ 'r^-','r>-','r<-']
+
+list_linestyle_nikuradse = ['bs','bd','bx',
+ 'b^','b>','b<']
+
+# plot
+
+fig, ax = plt.subplots()
+
+fig.set_size_inches(12,8)
+
+for j, model in enumerate(list_dff_models):
+
+ for i, relative_pipe_roughness in enumerate(list_relative_pipe_roughness):
+
+ # plot them using a log scale
+
+ linestyle_index = i+(j)*len(list_relative_pipe_roughness)
+
+ ax.loglog(list_reynolds_number,
+ dict_friction_factor[(model, relative_pipe_roughness)],
+ list_linestyle_our[linestyle_index],
+ #label=model+', e/D='+labels_relative_pipe_roughness[i],
+ label='CWE, e/D='+labels_relative_pipe_roughness[i])
+
+ #**********************************************************************
+
+ #**************************************************************************
+
+#******************************************************************************
+
+for e, exp in enumerate(nikuradse_experiments):
+
+ ax.loglog(nikuradse_data_re[exp],
+ nikuradse_data_f[exp],
+ list_linestyle_nikuradse[e],
+ label='N33, e/D=1/'+str(round(2*nikuradse_data_1e[exp])))
+
+#******************************************************************************
+
+ax.set(xlabel='Reynolds number [-]',
+ ylabel='Friction factor [-]',
+ title='Colebrook-White equation (CWE) vs results from Nikuradse (1933)')
+
+# plt.autoscale(enable=True, axis='x', tight=True)
+
+ax.grid(which='both')
+
+ax.legend()
+
+# left, right = xlim() # return the current xlim
+
+# plt.xlim((1e3, 1e8)) # set the xlim to left, right
+plt.xlim((1e3, 1e6)) # set the xlim to left, right
+plt.ylim((1.5e-2, 8e-2)) # set the xlim to left, right
+
+# fig.savefig("test.png")
+plt.show()
+
+#******************************************************************************
+#******************************************************************************
\ No newline at end of file