diff --git a/auto/Examples/varia/RoughSurface.py b/auto/Examples/varia/RoughSurface.py
index e941e69426deada0f76d05c97cb96f80992fd204..c14f8adc7d44356d689966c504ebc48b40ad4671 100755
--- a/auto/Examples/varia/RoughSurface.py
+++ b/auto/Examples/varia/RoughSurface.py
@@ -6,10 +6,12 @@ heights statistics and lateral roughness spectrum.
 
 import numpy as np
 import bornagain as ba
-from bornagain import nm
+from bornagain import nm, nm2, nm3
 from matplotlib import pyplot as plt, gridspec
 
 def plot(h):
+    rms = np.std(surface)
+    
     fig = plt.figure(figsize=(10, 6))
     gs = gridspec.GridSpec(3, 2)
     gs.update(left=0.1, right=0.48, wspace=0.15, hspace=0.5)
@@ -34,7 +36,7 @@ def plot(h):
 
     # hist
     fig.add_subplot(gs[-1, 1])
-    plt.hist(h.flatten(), range=[-4*sigma, 4*sigma],
+    plt.hist(h.flatten(), range=[-4*rms, 4*rms],
              bins=300, color='black')
     plt.title('Height histogram')
     plt.xlabel('h, nm')
@@ -42,27 +44,68 @@ def plot(h):
 
     plt.show()
 
-# roughness parameters
-sigma = 1*nm
-alpha = 0.5
-xi = 35*nm
+def get_sample():
+    """
+    A sample with one layer on a substrate, with partially 
+    correlated roughnesses.
+    """
+    # defining materials
+    vacuum = ba.RefractiveMaterial("ambience", 0, 0)
+    material_layer = ba.RefractiveMaterial("layer", 5e-6, 0)
+    material_substrate = ba.RefractiveMaterial("substrate", 15e-6, 0)
 
-# sample size
-Lx = 1000*nm
-Ly = 1000*nm
-X_points = 512
-Y_points = 512
+    # defining layers
+    l_ambience = ba.Layer(vacuum)
+    l_layer = ba.Layer(material_layer, 25*nm)
+    l_substrate = ba.Layer(material_substrate)
 
-# create roughness model
-autocorr = ba.K_CorrelationModel(sigma, alpha, xi)
-height_distribution = ba.ErfInterlayer()
-roughness = ba.LayerRoughness(autocorr, height_distribution)
+    # defining roughness
+    height_distribution = ba.ErfInterlayer()
+    max_spat_freq = 0.5*nm
 
-# generate roughness map
-roughness_map = ba.RoughnessMap(X_points, Y_points, Lx, Ly, roughness)
-surface = roughness_map.generate()
+    # for layer
+    omega = 1*nm3
+    a1 = 1
+    a2 = 10*nm
+    a3 = 100*nm2
+    a4 = 1000*nm3
+    autocorr_layer = ba.LinearGrowthModel(omega, a1, a2, a3, a4, max_spat_freq)
+    roughness_layer = ba.LayerRoughness(autocorr_layer, height_distribution)
+
+    # for substrate   
+    sigma = 0.9*nm
+    alpha = 0.5
+    xi = 35*nm
+    autocorr_sub = ba.K_CorrelationModel(sigma, alpha, xi, max_spat_freq)    
+    roughness_sub = ba.LayerRoughness(autocorr_sub, height_distribution)
+
+    # adding layers
+    my_sample = ba.MultiLayer()
+    my_sample.addLayer(l_ambience)
+    my_sample.addLayerWithTopRoughness(l_layer, roughness_layer)
+    my_sample.addLayerWithTopRoughness(l_substrate, roughness_sub)
+
+    return my_sample
+ 
+if __name__ == '__main__':    
+
+    # sample size
+    Lx = 1000*nm
+    Ly = 1000*nm
+    X_points = 512
+    Y_points = 512
+    
+    sample = get_sample()
+    interface_index = 0 # top interface
+    
+    # generate roughness map
+    roughness_map = ba.RoughnessMap(X_points, Y_points, Lx, Ly, 
+                                    sample, interface_index)
+    surface = roughness_map.generate()
+    
+    print("rms = {:.3}".format(np.std(surface)),"nm")
+    
+    plot(surface)
+    plt.show()
 
-print("rms = {:.3}".format(np.std(surface)),"nm")
 
-plot(surface)
-plt.show()
diff --git a/auto/FigExamples/varia/RoughSurface.py b/auto/FigExamples/varia/RoughSurface.py
index 64ac3864635b3b020b0a97daff807940c1b506d4..753819d7feda052bce516219ff4b9f5ed63f579b 100755
--- a/auto/FigExamples/varia/RoughSurface.py
+++ b/auto/FigExamples/varia/RoughSurface.py
@@ -6,10 +6,12 @@ heights statistics and lateral roughness spectrum.
 
 import numpy as np
 import bornagain as ba
-from bornagain import nm
+from bornagain import nm, nm2, nm3
 from matplotlib import pyplot as plt, gridspec
 
 def plot(h):
+    rms = np.std(surface)
+    
     fig = plt.figure(figsize=(10, 6))
     gs = gridspec.GridSpec(3, 2)
     gs.update(left=0.1, right=0.48, wspace=0.15, hspace=0.5)
@@ -34,7 +36,7 @@ def plot(h):
 
     # hist
     fig.add_subplot(gs[-1, 1])
-    plt.hist(h.flatten(), range=[-4*sigma, 4*sigma],
+    plt.hist(h.flatten(), range=[-4*rms, 4*rms],
              bins=300, color='black')
     plt.title('Height histogram')
     plt.xlabel('h, nm')
@@ -42,28 +44,69 @@ def plot(h):
 
     plt.show()
 
-# roughness parameters
-sigma = 1*nm
-alpha = 0.5
-xi = 35*nm
+def get_sample():
+    """
+    A sample with one layer on a substrate, with partially 
+    correlated roughnesses.
+    """
+    # defining materials
+    vacuum = ba.RefractiveMaterial("ambience", 0, 0)
+    material_layer = ba.RefractiveMaterial("layer", 5e-6, 0)
+    material_substrate = ba.RefractiveMaterial("substrate", 15e-6, 0)
 
-# sample size
-Lx = 1000*nm
-Ly = 1000*nm
-X_points = 512
-Y_points = 512
+    # defining layers
+    l_ambience = ba.Layer(vacuum)
+    l_layer = ba.Layer(material_layer, 25*nm)
+    l_substrate = ba.Layer(material_substrate)
 
-# create roughness model
-autocorr = ba.K_CorrelationModel(sigma, alpha, xi)
-height_distribution = ba.ErfInterlayer()
-roughness = ba.LayerRoughness(autocorr, height_distribution)
+    # defining roughness
+    height_distribution = ba.ErfInterlayer()
+    max_spat_freq = 0.5*nm
 
-# generate roughness map
-roughness_map = ba.RoughnessMap(X_points, Y_points, Lx, Ly, roughness)
-surface = roughness_map.generate()
+    # for layer
+    omega = 1*nm3
+    a1 = 1
+    a2 = 10*nm
+    a3 = 100*nm2
+    a4 = 1000*nm3
+    autocorr_layer = ba.LinearGrowthModel(omega, a1, a2, a3, a4, max_spat_freq)
+    roughness_layer = ba.LayerRoughness(autocorr_layer, height_distribution)
+
+    # for substrate   
+    sigma = 0.9*nm
+    alpha = 0.5
+    xi = 35*nm
+    autocorr_sub = ba.K_CorrelationModel(sigma, alpha, xi, max_spat_freq)    
+    roughness_sub = ba.LayerRoughness(autocorr_sub, height_distribution)
+
+    # adding layers
+    my_sample = ba.MultiLayer()
+    my_sample.addLayer(l_ambience)
+    my_sample.addLayerWithTopRoughness(l_layer, roughness_layer)
+    my_sample.addLayerWithTopRoughness(l_substrate, roughness_sub)
+
+    return my_sample
+ 
+if __name__ == '__main__':    
+
+    # sample size
+    Lx = 1000*nm
+    Ly = 1000*nm
+    X_points = 512
+    Y_points = 512
+    
+    sample = get_sample()
+    interface_index = 0 # top interface
+    
+    # generate roughness map
+    roughness_map = ba.RoughnessMap(X_points, Y_points, Lx, Ly, 
+                                    sample, interface_index)
+    surface = roughness_map.generate()
+    
+    print("rms = {:.3}".format(np.std(surface)),"nm")
+    
+    from bornagain import ba_plot as bp
+    plotargs = bp.parse_commandline()
+    bp.export(**plotargs)
 
-print("rms = {:.3}".format(np.std(surface)),"nm")
 
-from bornagain import ba_plot as bp
-plotargs = bp.parse_commandline()
-bp.export(**plotargs)
diff --git a/auto/MiniExamples/varia/RoughSurface.py b/auto/MiniExamples/varia/RoughSurface.py
index f04346dbc59b630aa99382054195579609362e0a..8dcf9bed08f3945659ff55f2ae8425ee439a9172 100755
--- a/auto/MiniExamples/varia/RoughSurface.py
+++ b/auto/MiniExamples/varia/RoughSurface.py
@@ -6,10 +6,12 @@ heights statistics and lateral roughness spectrum.
 
 import numpy as np
 import bornagain as ba
-from bornagain import nm
+from bornagain import nm, nm2, nm3
 from matplotlib import pyplot as plt, gridspec
 
 def plot(h):
+    rms = np.std(surface)
+    
     fig = plt.figure(figsize=(10, 6))
     gs = gridspec.GridSpec(3, 2)
     gs.update(left=0.1, right=0.48, wspace=0.15, hspace=0.5)
@@ -34,7 +36,7 @@ def plot(h):
 
     # hist
     fig.add_subplot(gs[-1, 1])
-    plt.hist(h.flatten(), range=[-4*sigma, 4*sigma],
+    plt.hist(h.flatten(), range=[-4*rms, 4*rms],
              bins=300, color='black')
     plt.title('Height histogram')
     plt.xlabel('h, nm')
@@ -42,28 +44,69 @@ def plot(h):
 
     plt.show()
 
-# roughness parameters
-sigma = 1*nm
-alpha = 0.5
-xi = 35*nm
+def get_sample():
+    """
+    A sample with one layer on a substrate, with partially 
+    correlated roughnesses.
+    """
+    # defining materials
+    vacuum = ba.RefractiveMaterial("ambience", 0, 0)
+    material_layer = ba.RefractiveMaterial("layer", 5e-6, 0)
+    material_substrate = ba.RefractiveMaterial("substrate", 15e-6, 0)
 
-# sample size
-Lx = 1000*nm
-Ly = 1000*nm
-X_points = 10
-Y_points = 10
+    # defining layers
+    l_ambience = ba.Layer(vacuum)
+    l_layer = ba.Layer(material_layer, 25*nm)
+    l_substrate = ba.Layer(material_substrate)
 
-# create roughness model
-autocorr = ba.K_CorrelationModel(sigma, alpha, xi)
-height_distribution = ba.ErfInterlayer()
-roughness = ba.LayerRoughness(autocorr, height_distribution)
+    # defining roughness
+    height_distribution = ba.ErfInterlayer()
+    max_spat_freq = 0.5*nm
 
-# generate roughness map
-roughness_map = ba.RoughnessMap(X_points, Y_points, Lx, Ly, roughness, 0)
-surface = roughness_map.generate()
+    # for layer
+    omega = 1*nm3
+    a1 = 1
+    a2 = 10*nm
+    a3 = 100*nm2
+    a4 = 1000*nm3
+    autocorr_layer = ba.LinearGrowthModel(omega, a1, a2, a3, a4, max_spat_freq)
+    roughness_layer = ba.LayerRoughness(autocorr_layer, height_distribution)
+
+    # for substrate   
+    sigma = 0.9*nm
+    alpha = 0.5
+    xi = 35*nm
+    autocorr_sub = ba.K_CorrelationModel(sigma, alpha, xi, max_spat_freq)    
+    roughness_sub = ba.LayerRoughness(autocorr_sub, height_distribution)
+
+    # adding layers
+    my_sample = ba.MultiLayer()
+    my_sample.addLayer(l_ambience)
+    my_sample.addLayerWithTopRoughness(l_layer, roughness_layer)
+    my_sample.addLayerWithTopRoughness(l_substrate, roughness_sub)
+
+    return my_sample
+ 
+if __name__ == '__main__':    
+
+    # sample size
+    Lx = 1000*nm
+    Ly = 1000*nm
+    X_points = 10
+    Y_points = 10
+    
+    sample = get_sample()
+    interface_index = 0 # top interface
+    
+    # generate roughness map
+    roughness_map = ba.RoughnessMap(X_points, Y_points, Lx, Ly, 
+                                    sample, interface_index)
+    surface = roughness_map.generate()
+    
+    print("rms = {:.3}".format(np.std(surface)),"nm")
+    
+    from bornagain import ba_plot as bp
+    plotargs = bp.parse_commandline()
+    bp.export(**plotargs)
 
-print("rms = {:.3}".format(np.std(surface)),"nm")
 
-from bornagain import ba_plot as bp
-plotargs = bp.parse_commandline()
-bp.export(**plotargs)
diff --git a/rawEx/varia/RoughSurface.py b/rawEx/varia/RoughSurface.py
index 1e22299eaa0fd73d51dd2f8e52777f27c85972e4..c12911350febe3a931e1283231d2ad4292b8c45b 100755
--- a/rawEx/varia/RoughSurface.py
+++ b/rawEx/varia/RoughSurface.py
@@ -6,10 +6,12 @@ heights statistics and lateral roughness spectrum.
 
 import numpy as np
 import bornagain as ba
-from bornagain import nm
+from bornagain import nm, nm2, nm3
 from matplotlib import pyplot as plt, gridspec
 
 def plot(h):
+    rms = np.std(surface)
+    
     fig = plt.figure(figsize=(10, 6))
     gs = gridspec.GridSpec(3, 2)
     gs.update(left=0.1, right=0.48, wspace=0.15, hspace=0.5)
@@ -34,7 +36,7 @@ def plot(h):
 
     # hist
     fig.add_subplot(gs[-1, 1])
-    plt.hist(h.flatten(), range=[-4*sigma, 4*sigma],
+    plt.hist(h.flatten(), range=[-4*rms, 4*rms],
              bins=300, color='black')
     plt.title('Height histogram')
     plt.xlabel('h, nm')
@@ -42,33 +44,74 @@ def plot(h):
 
     plt.show()
 
-# roughness parameters
-sigma = 1*nm
-alpha = 0.5
-xi = 35*nm
-
-# sample size
-Lx = 1000*nm
-Ly = 1000*nm
-X_points = <%= test_mode ? 10 : 512 %>
-Y_points = <%= test_mode ? 10 : 512 %>
-
-# create roughness model
-autocorr = ba.K_CorrelationModel(sigma, alpha, xi)
-height_distribution = ba.ErfInterlayer()
-roughness = ba.LayerRoughness(autocorr, height_distribution)
-
-# generate roughness map
-roughness_map = ba.RoughnessMap(X_points, Y_points, Lx, Ly, roughness<%= test_mode ? ", 0" : "" %>)
-surface = roughness_map.generate()
-
-print("rms = {:.3}".format(np.std(surface)),"nm")
-
-<%- if test_mode or figure_mode -%>
-from bornagain import ba_plot as bp
-plotargs = bp.parse_commandline()
-bp.export(**plotargs)
-<%- else -%>
-plot(surface)
-plt.show()
-<%- end -%>
+def get_sample():
+    """
+    A sample with one layer on a substrate, with partially 
+    correlated roughnesses.
+    """
+    # defining materials
+    vacuum = ba.RefractiveMaterial("ambience", 0, 0)
+    material_layer = ba.RefractiveMaterial("layer", 5e-6, 0)
+    material_substrate = ba.RefractiveMaterial("substrate", 15e-6, 0)
+
+    # defining layers
+    l_ambience = ba.Layer(vacuum)
+    l_layer = ba.Layer(material_layer, 25*nm)
+    l_substrate = ba.Layer(material_substrate)
+
+    # defining roughness
+    height_distribution = ba.ErfInterlayer()
+    max_spat_freq = 0.5*nm
+
+    # for layer
+    omega = 1*nm3
+    a1 = 1
+    a2 = 10*nm
+    a3 = 100*nm2
+    a4 = 1000*nm3
+    autocorr_layer = ba.LinearGrowthModel(omega, a1, a2, a3, a4, max_spat_freq)
+    roughness_layer = ba.LayerRoughness(autocorr_layer, height_distribution)
+
+    # for substrate   
+    sigma = 0.9*nm
+    alpha = 0.5
+    xi = 35*nm
+    autocorr_sub = ba.K_CorrelationModel(sigma, alpha, xi, max_spat_freq)    
+    roughness_sub = ba.LayerRoughness(autocorr_sub, height_distribution)
+
+    # adding layers
+    my_sample = ba.MultiLayer()
+    my_sample.addLayer(l_ambience)
+    my_sample.addLayerWithTopRoughness(l_layer, roughness_layer)
+    my_sample.addLayerWithTopRoughness(l_substrate, roughness_sub)
+
+    return my_sample
+ 
+if __name__ == '__main__':    
+
+    # sample size
+    Lx = 1000*nm
+    Ly = 1000*nm
+    X_points = <%= test_mode ? 10 : 512 %>
+    Y_points = <%= test_mode ? 10 : 512 %>
+    
+    sample = get_sample()
+    interface_index = 0 # top interface
+    
+    # generate roughness map
+    roughness_map = ba.RoughnessMap(X_points, Y_points, Lx, Ly, 
+                                    sample, interface_index)
+    surface = roughness_map.generate()
+    
+    print("rms = {:.3}".format(np.std(surface)),"nm")
+    
+    <%- if test_mode or figure_mode -%>
+    from bornagain import ba_plot as bp
+    plotargs = bp.parse_commandline()
+    bp.export(**plotargs)
+    <%- else -%>
+    plot(surface)
+    plt.show()
+    <%- end -%>
+
+