Added Scripts for Paper 1

scripts/paper 1/ch5-0-evaluation_and_comparison.py

+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jul 26 14:37:20 2013
+
+@author: Jan
+"""
+import sys
+import traceback
+import pdb
+import os
+
+from numpy import pi
+
+import shelve
+
+import pyramid.magcreator as mc
+from pyramid.magdata import MagData
+
+import matplotlib.pyplot as plt
+
+
+def run():
+    
+    print '\nACCESS SHELVE'
+    # Create / Open databank:
+    directory = '../../output/paper 1'
+    if not os.path.exists(directory):
+        os.makedirs(directory)
+    data_shelve = shelve.open(directory + '/paper_1_shelve')
+
+    ###############################################################################################
+    print 'CH5-0 MAGNETIC DISTRIBUTIONS'
+
+    key = 'ch5-0-magnetic_distributions'
+    if key in data_shelve:
+        print '--LOAD MAGNETIC DISTRIBUTIONS'
+        (mag_data_disc, mag_data_vort) = data_shelve[key]
+    else:
+        print '--CREATE MAGNETIC DISTRIBUTIONS'
+        # Input parameters:
+        res = 0.5  # in nm
+        phi = pi/2
+        dim = (32, 256, 256)  # in px (z, y, x)
+        # Create magnetic shape:
+        center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts with 0!
+        radius = dim[1]/4  # in px
+        height = dim[0]/2  # in px
+        mag_shape = mc.Shapes.disc(dim, center, radius, height)
+        print '--CREATE MAGN. DISTR. OF HOMOG. MAG. DISC'
+        mag_data_disc = MagData(res, mc.create_mag_dist(mag_shape, phi))
+        mag_data_disc.scale_down(3)
+        print '--CREATE MAGN. DISTR. OF VORTEX STATE DISC'
+        mag_data_vort = MagData(res, mc.create_mag_dist_vortex(mag_shape, center))
+        mag_data_vort.scale_down(3)
+        # Mayavi-Plots    
+        mag_data_disc.quiver_plot3d()    
+        mag_data_vort.quiver_plot3d()
+        print '--SHELVE MAGNETIC DISTRIBUTIONS'
+        data_shelve[key] = (mag_data_disc, mag_data_vort)
+
+    print '--PLOT/SAVE HOMOG. MAGN. DISC'
+    # Plot and save MagData (Disc):
+    mag_data_disc.quiver_plot('Homog. magn. disc')
+    plt.savefig(directory + '/ch5-0-mag_data_disc.png', bbox_inches='tight')
+
+    
+    print '--PLOT/SAVE VORTEX STATE DISC'
+    # Plot and save MagData (Vortex):
+    mag_data_vort.quiver_plot('Vortex state disc')
+    plt.savefig(directory + '/ch5-0-mag_data_vort.png', bbox_inches='tight')
+    
+    ###############################################################################################
+    print 'CLOSING SHELVE\n'
+    # Close shelve:
+    data_shelve.close()
+    
+    ###############################################################################################
+
+if __name__ == "__main__":
+    try:
+        run()
+    except:
+        type, value, tb = sys.exc_info()
+        traceback.print_exc()
+        pdb.post_mortem(tb)

scripts/paper 1/ch5-1-evaluation_real_space.py

+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jul 26 14:37:30 2013
+
+@author: Jan
+"""
+import pdb
+import traceback
+import sys
+import os
+
+import numpy as np
+from numpy import pi
+
+import shelve
+
+import pyramid.magcreator as mc
+import pyramid.projector as pj
+import pyramid.phasemapper as pm
+import pyramid.holoimage as hi
+import pyramid.analytic as an
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+
+import matplotlib.pyplot as plt
+from matplotlib.colors import BoundaryNorm
+from matplotlib.ticker import MaxNLocator
+from matplotlib.cm import RdBu
+from matplotlib.patches import Rectangle
+
+
+PHI_0 = -2067.83  # magnetic flux in T*nm²
+
+
+def run():
+
+    print '\nACCESS SHELVE'
+    # Create / Open databank:
+    directory = '../../output/paper 1'
+    if not os.path.exists(directory):
+        os.makedirs(directory)
+    data_shelve = shelve.open(directory + '/paper_1_shelve')
+
+    ###############################################################################################
+    print 'CH5-1 ANALYTIC SOLUTIONS'
+
+    # Input parameters:
+    res = 0.125  # in nm
+    phi = pi/2
+    dim = (128, 1024, 1024)  # in px (z, y, x)
+    # Create magnetic shape:
+    center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts with 0!
+    radius = dim[1]/4  # in px
+    height = dim[0]/2  # in px
+    print '--CALCULATE ANALYTIC SOLUTIONS'
+    # Get analytic solution:
+    phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
+    phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+    phase_map_ana_disc = PhaseMap(res, phase_ana_disc)
+    phase_map_ana_vort = PhaseMap(res, phase_ana_vort)
+    print '--PLOT/SAVE ANALYTIC SOLUTIONS'
+    hi.display_combined(phase_map_ana_disc, 100, 'Analytic solution: hom. magn. disc', 'bilinear')
+    axis = plt.gcf().add_subplot(1, 2, 2, aspect='equal')
+    axis.axhline(y=512, linewidth=3, linestyle='--', color='r')
+    plt.savefig(directory + '/ch5-1-analytic_solution_disc.png', bbox_inches='tight')
+    hi.display_combined(phase_map_ana_vort, 100, 'Analytic solution: Vortex state', 'bilinear')
+    axis = plt.gcf().add_subplot(1, 2, 2, aspect='equal')
+    axis.axhline(y=512, linewidth=3, linestyle='--', color='r')
+    plt.savefig(directory + '/ch5-1-analytic_solution_vort.png', bbox_inches='tight')
+    # Get colorwheel:
+    hi.make_color_wheel()
+    plt.savefig(directory + '/ch5-1-colorwheel.png', bbox_inches='tight')
+
+    ###############################################################################################
+    print 'CH5-1 PHASE SLICES REAL SPACE'  # TODO: Verschieben
+    
+    # Input parameters:
+    res = 0.25  # in nm
+    phi = pi/2
+    density = 20
+    dim = (64, 512, 512)  # in px (z, y, x)
+    # Create magnetic shape:
+    center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts with 0!
+    radius = dim[1]/4  # in px
+    height = dim[0]/2  # in px
+    
+    key = 'ch5-1-phase_slice_mag_dist'
+    if key in data_shelve:
+        print '--LOAD MAGNETIC DISTRIBUTION'
+        mag_shape = data_shelve[key]
+    else:
+        print '--CREATE MAGNETIC DISTRIBUTION'
+        mag_shape = mc.Shapes.disc(dim, center, radius, height)
+        print '--SAVE MAGNETIC DISTRIBUTION'
+        data_shelve[key] = mag_shape
+
+    key = 'ch5-1-phase_slice_disc'
+    if key in data_shelve:
+        print '--LOAD PHASE SLICES HOMOG. MAGN. DISC'
+        (x, y) = data_shelve[key]
+    else:
+        print '--CREATE PHASE SLICES HOMOG. MAGN. DISC'
+        # Arrays for plotting:
+        x = []
+        y = []
+        # Analytic solution:
+        L = dim[1] * res  # in px/nm
+        Lz = 0.5 * dim[0] * res  # in px/nm
+        R = 0.25 * L  # in px/nm
+        x0 = L / 2  # in px/nm
+    
+        def F_disc(x):
+            coeff = - pi * Lz / (2*PHI_0)
+            result = coeff * (- (x - x0) * np.sin(phi))
+            result *= np.where(np.abs(x - x0) <= R, 1, (R / (x - x0)) ** 2)
+            return result
+        
+        x.append(np.linspace(0, L, 5000))
+        y.append(F_disc(x[0]))
+        # Create and plot MagData (Disc):
+        mag_data_disc = MagData(res, mc.create_mag_dist(mag_shape, phi))
+        for i in range(5):
+            mag_data_disc.scale_down()
+            print '----res =', mag_data_disc.res, 'nm', 'dim =', mag_data_disc.dim
+            projection = pj.simple_axis_projection(mag_data_disc)
+            phase_map = PhaseMap(mag_data_disc.res, 
+                                 pm.phase_mag_real(mag_data_disc.res, projection, 'slab'))
+            hi.display_combined(phase_map, density, 'Disc, res = {}'.format(res))
+            x.append(np.linspace(0, mag_data_disc.dim[1]*mag_data_disc.res, mag_data_disc.dim[1]))
+            y.append(phase_map.phase[int(mag_data_disc.dim[1]/2), :])
+        # Shelve x and y:
+        print '--SAVE PHASE SLICES HOMOG. MAGN. DISC'
+        data_shelve[key] = (x, y)
+    
+    print '--PLOT/SAVE PHASE SLICES HOMOG. MAGN. DISC'
+    # Plot phase slices:
+    fig = plt.figure()
+    axis = fig.add_subplot(1, 1, 1)
+    axis.plot(x[0], y[0], 'k', label='analytic')
+    axis.plot(x[1], y[1], 'r', label='0.5 nm')
+    axis.plot(x[2], y[2], 'm', label='1 nm')
+    axis.plot(x[3], y[3], 'y', label='2 nm')
+    axis.plot(x[4], y[4], 'g', label='4 nm')
+    axis.plot(x[5], y[5], 'c', label='8 nm')
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    axis.set_title('DISC', fontsize=18)
+    axis.set_xlabel('x [nm]', fontsize=15)
+    axis.set_ylabel('phase [rad]', fontsize=15)
+    axis.set_xlim(0, 128)
+    axis.set_ylim(-0.22, 0.22)
+    axis.legend()
+    # Plot Zoombox and Arrow:
+    rect1 = Rectangle((23.5, 0.16), 15, 0.04, fc='w', ec='k')
+    axis.add_patch(rect1)
+    plt.arrow(32.5, 0.16, 0.0, -0.16, length_includes_head=True, 
+              head_width=2, head_length=0.02, fc='k', ec='k')
+    # Plot zoom inset:
+    ins_axis = plt.axes([0.2, 0.2, 0.3, 0.3])
+    ins_axis.plot(x[0], y[0], 'k', label='analytic')
+    ins_axis.plot(x[1], y[1], 'r', label='0.5 nm')
+    ins_axis.plot(x[2], y[2], 'm', label='1 nm')
+    ins_axis.plot(x[3], y[3], 'y', label='2 nm')
+    ins_axis.plot(x[4], y[4], 'g', label='4 nm')
+    ins_axis.plot(x[5], y[5], 'c', label='8 nm')
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    ins_axis.set_xlim(23.5, 38.5)
+    ins_axis.set_ylim(0.16, 0.2)
+    ins_axis.xaxis.set_major_locator(MaxNLocator(nbins=4, integer= True))
+    ins_axis.yaxis.set_major_locator(MaxNLocator(nbins=3))
+    plt.show()
+    plt.savefig(directory + '/ch5-1-disc_slice_comparison.png', bbox_inches='tight')
+    
+    key = 'ch5-1-phase_slice_vort'
+    if key in data_shelve:
+        print '--LOAD PHASE SLICES VORTEX STATE DISC'
+        (x, y) = data_shelve[key]
+    else:
+        print '--CREATE PHASE SLICES VORTEX STATE DISC'
+        x = []
+        y = []
+        # Analytic solution:
+        L = dim[1] * res  # in px/nm
+        Lz = 0.5 * dim[0] * res  # in px/nm
+        R = 0.25 * L  # in px/nm
+        x0 = L / 2  # in px/nm
+    
+        def F_vort(x):
+            coeff = pi*Lz/PHI_0
+            result = coeff * np.where(np.abs(x - x0) <= R, (np.abs(x-x0)-R), 0)
+            return result
+        
+        x.append(np.linspace(0, L, 5001))
+        y.append(F_vort(x[0]))
+        # Create and plot MagData (Vortex):
+        mag_data_vort = MagData(res, mc.create_mag_dist_vortex(mag_shape))
+        for i in range(5):
+            mag_data_vort.scale_down()
+            print '----i =', i, 'dim =', mag_data_vort.dim, 'res =', mag_data_vort.res, 'nm'
+            projection = pj.simple_axis_projection(mag_data_vort)
+            phase_map = PhaseMap(mag_data_vort.res, 
+                                 pm.phase_mag_real(mag_data_vort.res, projection, 'slab'))
+            hi.display_combined(phase_map, density, 'Disc, res = {}'.format(res))
+            x.append(np.linspace(0, mag_data_vort.dim[1]*mag_data_vort.res, mag_data_vort.dim[1]))
+            y.append(phase_map.phase[int(mag_data_vort.dim[1]/2), :])
+        # Shelve x and y:
+        print '--SAVE PHASE SLICES VORTEX STATE DISC'
+        data_shelve[key] = (x, y)
+
+    print '--PLOT/SAVE PHASE SLICES VORTEX STATE DISC'
+    # Plot:
+    fig = plt.figure()
+    axis = fig.add_subplot(1, 1, 1)
+    axis.plot(x[0], y[0], 'k', label='analytic')
+    axis.plot(x[1], y[1], 'r', label='0.5 nm')
+    axis.plot(x[2], y[2], 'm', label='1 nm')
+    axis.plot(x[3], y[3], 'y', label='2 nm')
+    axis.plot(x[4], y[4], 'g', label='4 nm')
+    axis.plot(x[5], y[5], 'c', label='8 nm')
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    axis.set_title('VORTEX', fontsize=18)
+    axis.set_xlabel('x [nm]', fontsize=15)
+    axis.set_ylabel('phase [rad]', fontsize=15)
+    axis.set_xlim(0, 128)
+    axis.legend()
+    # Plot Zoombox and Arrow:
+    zoom = (59, 0.34, 10, 0.055)
+    rect1 = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
+    axis.add_patch(rect1)
+    plt.arrow(zoom[0]+zoom[2]/2, zoom[1], 0, -0.19, length_includes_head=True, 
+              head_width=2, head_length=0.02, fc='k', ec='k')
+    # Plot zoom inset:
+    ins_axis = plt.axes([0.47, 0.15, 0.15, 0.3])
+    ins_axis.plot(x[0], y[0], 'k', label='analytic')
+    ins_axis.plot(x[1], y[1], 'r', label='0.5 nm')
+    ins_axis.plot(x[2], y[2], 'm', label='1 nm')
+    ins_axis.plot(x[3], y[3], 'y', label='2 nm')
+    ins_axis.plot(x[4], y[4], 'g', label='4 nm')
+    ins_axis.plot(x[5], y[5], 'c', label='8 nm')
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    ins_axis.set_xlim(zoom[0], zoom[0]+zoom[2])
+    ins_axis.set_ylim(zoom[1], zoom[1]+zoom[3])
+    ins_axis.xaxis.set_major_locator(MaxNLocator(nbins=4, integer= True))
+    ins_axis.yaxis.set_major_locator(MaxNLocator(nbins=4))
+    plt.savefig(directory + '/ch5-1-vortex_slice_comparison.png', bbox_inches='tight')
+    
+    ###############################################################################################
+    print 'CH5-1 PHASE DIFFERENCES REAL SPACE'
+
+    # Input parameters:
+    res = 1.0  # in nm
+    phi = pi/2
+    dim = (16, 128, 128)  # in px (z, y, x)
+    center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts with 0!
+    radius = dim[1]/4  # in px
+    height = dim[0]/2  # in px
+
+    key = 'ch5-1-phase_diff_mag_dist'
+    if key in data_shelve:
+        print '--LOAD MAGNETIC DISTRIBUTIONS'
+        (mag_data_disc, mag_data_vort) = data_shelve[key]
+    else:
+        print '--CREATE MAGNETIC DISTRIBUTIONS'
+        # Create magnetic shape (4 times the size):
+        res_big = res / 2
+        dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
+        center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
+        radius_big = dim_big[1]/4  # in px
+        height_big = dim_big[0]/2  # in px
+        mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
+        # Create MagData (4 times the size):
+        mag_data_disc = MagData(res_big, mc.create_mag_dist(mag_shape, phi))
+        mag_data_vort = MagData(res_big, mc.create_mag_dist_vortex(mag_shape, center_big))
+        # Scale mag_data, resolution and dimensions:
+        mag_data_disc.scale_down()
+        mag_data_vort.scale_down()
+        print '--SAVE MAGNETIC DISTRIBUTIONS'
+        # Shelve magnetic distributions:
+        data_shelve[key] = (mag_data_disc, mag_data_vort)
+
+    print '--PLOT/SAVE PHASE DIFFERENCES'
+    # Create projections along z-axis:
+    projection_disc = pj.simple_axis_projection(mag_data_disc)
+    projection_vort = pj.simple_axis_projection(mag_data_vort)
+    # Get analytic solutions:
+    phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
+    phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+    # Create norm for the plots:
+    bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
+    norm = BoundaryNorm(bounds, RdBu.N)
+    # Disc:
+    phase_num_disc = pm.phase_mag_real(res, projection_disc, 'slab')
+    phase_diff_disc = PhaseMap(res, (phase_num_disc-phase_ana_disc)*1E3)  # in mrad -> *1000)
+    RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
+    phase_diff_disc.display('Deviation (homog. magn. disc), RMS = {:3.2f} mrad'.format(RMS_disc),
+                            labels=('x-axis [nm]', 'y-axis [nm]', '$\Delta$phase [mrad]'),
+                            limit=np.max(bounds), norm=norm)
+    plt.savefig(directory + '/ch5-1-disc_phase_diff.png', bbox_inches='tight')
+    # Vortex:
+    phase_num_vort = pm.phase_mag_real(res, projection_vort, 'slab')
+    phase_diff_vort = PhaseMap(res, (phase_num_vort-phase_ana_vort)*1E3)  # in mrad -> *1000
+    RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
+    phase_diff_vort.display('Deviation (vortex state disc), RMS = {:3.2f} mrad'.format(RMS_vort),
+                            labels=('x-axis [nm]', 'y-axis [nm]', '$\Delta$phase [mrad]'),
+                            limit=np.max(bounds), norm=norm)
+    plt.savefig(directory + '/ch5-1-vortex_phase_diff.png', bbox_inches='tight')
+
+    ###############################################################################################
+    print 'CLOSING SHELVE\n'
+    # Close shelve:
+    data_shelve.close()
+    
+    ###############################################################################################
+
+
+if __name__ == "__main__":
+    try:
+        run()
+    except:
+        type, value, tb = sys.exc_info()
+        traceback.print_exc()
+        pdb.post_mortem(tb)
+

scripts/paper 1/ch5-2-evaluation_fourier_space.py

+"""Compare the different methods to create phase maps."""
+
+
+import time
+import pdb
+import traceback
+import sys
+import os
+
+import numpy as np
+from numpy import pi
+
+import shelve
+
+import pyramid.magcreator as mc
+import pyramid.projector as pj
+import pyramid.phasemapper as pm
+import pyramid.analytic as an
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+
+import matplotlib.pyplot as plt
+from matplotlib.colors import BoundaryNorm
+from matplotlib.ticker import MaxNLocator
+from matplotlib.cm import RdBu
+
+
+def run():
+
+    print '\nACCESS SHELVE'
+    # Create / Open databank:
+    directory = '../../output/paper 1'
+    if not os.path.exists(directory):
+        os.makedirs(directory)
+    data_shelve = shelve.open(directory + '/paper_1_shelve')
+
+    ###############################################################################################
+    print 'CH5-2 PHASE DIFFERENCES FOURIER SPACE'
+
+    # Input parameters:
+    res = 1.0  # in nm
+    phi = pi/2
+    dim = (16, 128, 128)  # in px (z, y, x)
+    center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts with 0!
+    radius = dim[1]/4  # in px
+    height = dim[0]/2  # in px
+
+    key = 'ch5-1-phase_diff_mag_dist'
+    if key in data_shelve:
+        print '--LOAD MAGNETIC DISTRIBUTIONS'
+        (mag_data_disc, mag_data_vort) = data_shelve[key]
+    else:
+        print '--CREATE MAGNETIC DISTRIBUTIONS'
+        # Create magnetic shape (4 times the size):
+        res_big = res / 2
+        dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
+        center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
+        radius_big = dim_big[1]/4  # in px
+        height_big = dim_big[0]/2  # in px
+        mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
+        # Create MagData (4 times the size):
+        mag_data_disc = MagData(res_big, mc.create_mag_dist(mag_shape, phi))
+        mag_data_vort = MagData(res_big, mc.create_mag_dist_vortex(mag_shape, center_big))
+        # Scale mag_data, resolution and dimensions:
+        mag_data_disc.scale_down()
+        mag_data_vort.scale_down()
+        print '--SAVE MAGNETIC DISTRIBUTIONS'
+        # Shelve magnetic distributions:
+        data_shelve[key] = (mag_data_disc, mag_data_vort)
+
+    print '--PLOT/SAVE PHASE DIFFERENCES'
+    # Create projections along z-axis:
+    projection_disc = pj.simple_axis_projection(mag_data_disc)
+    projection_vort = pj.simple_axis_projection(mag_data_vort)
+    # Get analytic solutions:
+    phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
+    phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+    # Create norm for the plots:
+    bounds = np.array([-100, -50, -25, -5, 0, 5, 25, 50, 100])
+    norm = BoundaryNorm(bounds, RdBu.N)
+    # Disc:
+    phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=0)
+    phase_diff_disc = PhaseMap(res, (phase_num_disc-phase_ana_disc)*1E3)  # in mrad -> *1000)
+    RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
+    phase_diff_disc.display('Deviation (homog. magn. disc), RMS = {:3.2f} mrad'.format(RMS_disc),
+                            labels=('x-axis [nm]', 'y-axis [nm]', 
+                                    '$\Delta$phase [mrad] (padding = 0)'),
+                            limit=np.max(bounds), norm=norm)
+    plt.savefig(directory + '/ch5-2-disc_phase_diff_no_padding.png', bbox_inches='tight')
+    # Vortex:
+    phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=0)
+    phase_diff_vort = PhaseMap(res, (phase_num_vort-phase_ana_vort)*1E3)  # in mrad -> *1000
+    RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
+    phase_diff_vort.display('Deviation (vortex state disc), RMS = {:3.2f} mrad'.format(RMS_vort),
+                            labels=('x-axis [nm]', 'y-axis [nm]', 
+                                    '$\Delta$phase [mrad] (padding = 0)'),
+                            limit=np.max(bounds), norm=norm)
+    plt.savefig(directory + '/ch5-2-vortex_phase_diff_no_padding.png', bbox_inches='tight')
+
+    # Create norm for the plots:
+    bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
+    norm = BoundaryNorm(bounds, RdBu.N)
+    # Disc:
+    phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=10)
+    phase_diff_disc = PhaseMap(res, (phase_num_disc-phase_ana_disc)*1E3)  # in mrad -> *1000)
+    RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
+    phase_diff_disc.display('Deviation (homog. magn. disc), RMS = {:3.2f} mrad'.format(RMS_disc),
+                            labels=('x-axis [nm]', 'y-axis [nm]', 
+                                    '$\Delta$phase [mrad] (padding = 10)'),
+                            limit=np.max(bounds), norm=norm)
+    plt.savefig(directory + '/ch5-2-disc_phase_diff_padding_10.png', bbox_inches='tight')
+    # Vortex:
+    phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=10)
+    phase_diff_vort = PhaseMap(res, (phase_num_vort-phase_ana_vort)*1E3)  # in mrad -> *1000
+    RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
+    phase_diff_vort.display('Deviation (vortex state disc), RMS = {:3.2f} mrad'.format(RMS_vort),
+                            labels=('x-axis [nm]', 'y-axis [nm]', 
+                                    '$\Delta$phase [mrad] (padding = 10)'),
+                            limit=np.max(bounds), norm=norm)
+    plt.savefig(directory + '/ch5-2-vortex_phase_diff_padding_10.png', bbox_inches='tight')
+
+    ###############################################################################################
+    print 'CH5-2 FOURIER PADDING VARIATION'
+
+    # Input parameters:
+    res = 1.0  # in nm
+    phi = pi/2
+    dim = (16, 128, 128)  # in px (z, y, x)
+    center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts with 0!
+    radius = dim[1]/4  # in px
+    height = dim[0]/2  # in px
+    
+    key = 'ch5-2-fourier_padding_mag_dist'
+    if key in data_shelve:
+        print '--LOAD MAGNETIC DISTRIBUTIONS'
+        (mag_data_disc, mag_data_vort) = data_shelve[key]
+    else:
+        print '--CREATE MAGNETIC DISTRIBUTIONS'
+        # Create magnetic shape (4 times the size):
+        res_big = res / 2
+        dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
+        center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
+        radius_big = dim_big[1]/4  # in px
+        height_big = dim_big[0]/2  # in px
+        mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
+        # Create MagData (4 times the size):
+        mag_data_disc = MagData(res_big, mc.create_mag_dist(mag_shape, phi))
+        mag_data_vort = MagData(res_big, mc.create_mag_dist_vortex(mag_shape, center_big))
+        # Scale mag_data, resolution and dimensions:
+        mag_data_disc.scale_down()
+        mag_data_vort.scale_down()
+        print '--SAVE MAGNETIC DISTRIBUTIONS'
+        # Shelve magnetic distributions:
+        data_shelve[key] = (mag_data_disc, mag_data_vort)
+
+    # Create projections along z-axis:
+    projection_disc = pj.simple_axis_projection(mag_data_disc)
+    projection_vort = pj.simple_axis_projection(mag_data_vort)
+    # Get analytic solutions:
+    phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
+    phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+
+    # List of applied paddings:
+    padding_list = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
+
+    print '--LOAD/CREATE PADDING SERIES OF HOMOG. MAGN. DISC'
+    data_disc = np.zeros((3, len(padding_list)))
+    data_disc[0, :] = padding_list
+    for (i, padding) in enumerate(padding_list):
+        key = ', '.join(['Padding->RMS|duration', 'Fourier', 'padding={}'.format(padding_list[i]),
+                        'res={}'.format(res), 'dim={}'.format(dim), 'phi={}'.format(phi), 'disc'])
+        if key in data_shelve:
+            data_disc[:, i] = data_shelve[key]
+        else:
+            print '----calculate and save padding =', padding_list[i]
+            start_time = time.time()
+            phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding_list[i])
+            data_disc[2, i] = time.time() - start_time
+            phase_diff = (phase_num_disc-phase_ana_disc) * 1E3  # in mrad -> *1000)
+            phase_map_diff = PhaseMap(res, phase_diff)
+            phase_map_diff.display(labels=('x-axis [nm]', 'y-axis [nm]', 'phase [mrad]'))
+            data_disc[1, i] = np.sqrt(np.mean(phase_diff**2))
+            data_shelve[key] = data_disc[:, i]
+
+    print '--PLOT/SAVE PADDING SERIES OF HOMOG. MAGN. DISC'
+    # Plot RMS against padding:
+    fig = plt.figure()
+    axis = fig.add_subplot(1, 1, 1)
+    axis.axhline(y=0.18, linestyle='--', color='g', label='RMS [mrad] (real space)')
+    axis.plot(data_disc[0], data_disc[1], 'go-', label='RMS [mrad] (Fourier space)')
+    axis.set_title('Variation of the Padding (homog. magn. disc)', fontsize=18)
+    axis.set_xlabel('padding', fontsize=15)
+    axis.set_ylabel('RMS [mrad]', fontsize=15)
+    axis.set_xlim(-0.5, 16.5)
+    axis.set_ylim(-5, 45)
+    axis.xaxis.set_major_locator(MaxNLocator(nbins=10, integer= True))
+    axis.legend()
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    # Plot zoom inset:
+    ins_axis = plt.axes([0.3, 0.3, 0.55, 0.4])
+    ins_axis.axhline(y=0.18, linestyle='--', color='g')
+    ins_axis.plot(data_disc[0], data_disc[1], 'go-')
+    ins_axis.set_yscale('log')
+    ins_axis.set_xlim(5.5, 16.5)
+    ins_axis.set_ylim(0.1, 1.1)
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    plt.savefig(directory + '/ch5-2-disc_padding_RMS.png', bbox_inches='tight')
+    # Plot duration against padding:
+    fig = plt.figure()
+    axis = fig.add_subplot(1, 1, 1)
+    axis.plot(data_disc[0], data_disc[2], 'bo-')
+    axis.set_title('Variation of the Padding (homog. magn. disc)', fontsize=18)
+    axis.set_xlabel('padding', fontsize=15)
+    axis.set_ylabel('duration [s]', fontsize=15)
+    axis.set_xlim(-0.5, 16.5)
+    axis.set_ylim(-0.05, 1.5)
+    axis.xaxis.set_major_locator(MaxNLocator(nbins=10, integer= True))
+    axis.yaxis.set_major_locator(MaxNLocator(nbins=10))
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    plt.savefig(directory + '/ch5-2-disc_padding_duration.png', bbox_inches='tight')
+
+
+    print '--LOAD/CREATE PADDING SERIES OF VORTEX STATE DISC'
+    data_vort = np.zeros((3, len(padding_list)))
+    data_vort[0, :] = padding_list
+    for (i, padding) in enumerate(padding_list):
+        key = ', '.join(['Padding->RMS|duration', 'Fourier', 'padding={}'.format(padding_list[i]),
+                        'res={}'.format(res), 'dim={}'.format(dim), 'phi={}'.format(phi), 'vort'])
+        if key in data_shelve:
+            data_vort[:, i] = data_shelve[key]
+        else:
+            print '----calculate and save padding =', padding_list[i]
+            start_time = time.time()
+            phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding_list[i])
+            data_vort[2, i] = time.time() - start_time
+            phase_diff = (phase_num_vort-phase_ana_vort) * 1E3  # in mrad -> *1000)
+            phase_map_diff = PhaseMap(res, phase_diff)
+            phase_map_diff.display(labels=('x-axis [nm]', 'y-axis [nm]', 'phase [mrad]'))
+            data_vort[1, i] = np.sqrt(np.mean(phase_diff**2))
+            data_shelve[key] = data_vort[:, i]
+
+    print '--PLOT/SAVE PADDING SERIES OF VORTEX STATE DISC'
+    # Plot RMS against padding:
+    fig = plt.figure()
+    axis = fig.add_subplot(1, 1, 1)
+    axis.axhline(y=0.22, linestyle='--', color='g', label='RMS [mrad] (real space)')
+    axis.plot(data_vort[0], data_vort[1], 'go-', label='RMS [mrad] (Fourier space)')
+    axis.set_title('Variation of the Padding (vortex state disc)', fontsize=18)
+    axis.set_xlabel('padding', fontsize=15)
+    axis.set_ylabel('RMS [mrad]', fontsize=15)
+    axis.set_xlim(-0.5, 16.5)
+    axis.set_ylim(-5, 45)
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    axis.xaxis.set_major_locator(MaxNLocator(nbins=10, integer= True))
+    axis.legend()
+    # Plot zoom inset:
+    ins_axis = plt.axes([0.3, 0.3, 0.55, 0.4])
+    ins_axis.axhline(y=0.22, linestyle='--', color='g')
+    ins_axis.plot(data_vort[0], data_vort[1], 'go-')
+    ins_axis.set_yscale('log')
+    ins_axis.set_xlim(5.5, 16.5)
+    ins_axis.set_ylim(0.1, 1.1)
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    plt.savefig(directory + '/ch5-2-vortex_padding_RMS.png', bbox_inches='tight')
+    # Plot duration against padding:
+    fig = plt.figure()
+    axis = fig.add_subplot(1, 1, 1)
+    axis.plot(data_vort[0], data_vort[2], 'bo-')
+    axis.set_title('Variation of the Padding (vortex state disc)', fontsize=18)
+    axis.set_xlabel('padding', fontsize=15)
+    axis.set_ylabel('duration [s]', fontsize=15)
+    axis.set_xlim(-0.5, 16.5)
+    axis.set_ylim(-0.05, 1.5)
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    plt.savefig(directory + '/ch5-2-vortex_padding_duration.png', bbox_inches='tight')
+
+    ###############################################################################################
+    print 'CLOSING SHELVE\n'
+    # Close shelve:
+    data_shelve.close()
+    
+    ###############################################################################################
+
+
+if __name__ == "__main__":
+    try:
+        run()
+    except:
+        type, value, tb = sys.exc_info()
+        traceback.print_exc()
+        pdb.post_mortem(tb)

scripts/paper 1/ch5-3-comparison_of_methods.py

+#! python
+# -*- coding: utf-8 -*-
+"""Compare the different methods to create phase maps."""
+
+
+import time
+import pdb
+import traceback
+import sys
+import os
+import numpy as np
+from numpy import pi
+import matplotlib.pyplot as plt
+
+import pyramid.magcreator as mc
+import pyramid.projector as pj
+import pyramid.phasemapper as pm
+import pyramid.analytic as an
+from pyramid.magdata import MagData
+import shelve
+
+
+def run():
+
+    print '\nACCESS SHELVE'
+    # Create / Open databank:
+    directory = '../../output/paper 1'
+    if not os.path.exists(directory):
+        os.makedirs(directory)
+    data_shelve = shelve.open(directory + '/paper_1_shelve')
+
+    ###############################################################################################
+    print 'CH5-3 METHOD COMPARISON'
+    
+    key = 'ch5-3-compare_method_data'
+    if key in data_shelve:
+        print '--LOAD METHOD DATA'
+        (data_disc_fourier0, data_vort_fourier0,
+         data_disc_fourier1, data_vort_fourier1,
+         data_disc_fourier10, data_vort_fourier10,
+         data_disc_real_s, data_vort_real_s,
+         data_disc_real_d, data_vort_real_d) = data_shelve[key]
+    else:
+        # Input parameters:
+        steps = 6 
+        res = 0.25  # in nm
+        phi = pi/2
+        dim = (64, 512, 512)  # in px (z, y, x)
+        center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts with 0!
+        radius = dim[1]/4  # in px
+        height = dim[0]/2  # in px
+        
+        print '--CREATE MAGNETIC SHAPE'
+        mag_shape = mc.Shapes.disc(dim, center, radius, height)
+        # Create MagData (4 times the size):
+        print '--CREATE MAG. DIST. HOMOG. MAGN. DISC'
+        mag_data_disc = MagData(res, mc.create_mag_dist(mag_shape, phi))
+        print '--CREATE MAG. DIST. VORTEX STATE DISC'
+        mag_data_vort = MagData(res, mc.create_mag_dist_vortex(mag_shape, center))
+    
+        # Create Data Arrays
+        res_list = [res*2**i for i in np.linspace(1, steps, steps)]
+        data_disc_fourier0 = np.vstack((res_list, np.zeros((2, steps))))
+        data_vort_fourier0 = np.vstack((res_list, np.zeros((2, steps))))
+        data_disc_fourier1 = np.vstack((res_list, np.zeros((2, steps))))
+        data_vort_fourier1 = np.vstack((res_list, np.zeros((2, steps))))
+        data_disc_fourier10 = np.vstack((res_list, np.zeros((2, steps))))
+        data_vort_fourier10 = np.vstack((res_list, np.zeros((2, steps))))
+        data_disc_real_s = np.vstack((res_list, np.zeros((2, steps))))
+        data_vort_real_s = np.vstack((res_list, np.zeros((2, steps))))
+        data_disc_real_d = np.vstack((res_list, np.zeros((2, steps))))
+        data_vort_real_d = np.vstack((res_list, np.zeros((2, steps))))
+
+        for i in range(steps):
+            # Scale mag_data, resolution and dimensions:
+            mag_data_disc.scale_down()
+            mag_data_vort.scale_down()
+            dim = mag_data_disc.dim
+            res = mag_data_disc.res
+            center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts with 0!
+            radius = dim[1]/4  # in px
+            height = dim[0]/2  # in px
+            
+            print '--res =', res, 'nm', 'dim =', dim
+            
+            print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
+            # Create projections along z-axis:
+            projection_disc = pj.simple_axis_projection(mag_data_disc)
+            # Analytic solution:
+            phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
+            # Fourier unpadded:
+            padding = 0
+            start_time = time.clock()
+            phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+            data_disc_fourier0[2, i] = time.clock() - start_time
+            print '------time (disc, fourier0) =', data_disc_fourier0[2, i]
+            phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3  # in mrad -> *1000
+            data_disc_fourier0[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+            # Fourier padding 1:
+            padding = 1
+            start_time = time.clock()
+            phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+            data_disc_fourier1[2, i] = time.clock() - start_time
+            print '------time (disc, fourier1) =', data_disc_fourier1[2, i]
+            phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3  # in mrad -> *1000
+            data_disc_fourier1[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+            # Fourier padding 10:
+            padding = 10
+            start_time = time.clock()
+            phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+            data_disc_fourier10[2, i] = time.clock() - start_time
+            print '------time (disc, fourier10) =', data_disc_fourier10[2, i]
+            phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3  # in mrad -> *1000
+            data_disc_fourier10[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+            # Real space slab:
+            start_time = time.clock()
+            phase_num_disc = pm.phase_mag_real(res, projection_disc, 'slab')
+            data_disc_real_s[2, i] = time.clock() - start_time
+            print '------time (disc, real slab) =', data_disc_real_s[2, i]
+            phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3  # in mrad -> *1000
+            data_disc_real_s[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+            # Real space disc:
+            start_time = time.clock()
+            phase_num_disc = pm.phase_mag_real(res, projection_disc, 'disc')
+            data_disc_real_d[2, i] = time.clock() - start_time
+            print '------time (disc, real disc) =', data_disc_real_d[2, i]
+            phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3  # in mrad -> *1000
+            data_disc_real_d[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+
+            print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
+            # Create projections along z-axis:
+            projection_vort = pj.simple_axis_projection(mag_data_vort)
+            # Analytic solution:
+            phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+            # Fourier unpadded:
+            padding = 0
+            start_time = time.clock()
+            phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+            data_vort_fourier0[2, i] = time.clock() - start_time
+            print '------time (vortex, fourier0) =', data_vort_fourier0[2, i]
+            phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3  # in mrad -> *1000
+            data_vort_fourier0[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+            # Fourier padding 1:
+            padding = 1
+            start_time = time.clock()
+            phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+            data_vort_fourier1[2, i] = time.clock() - start_time
+            print '------time (vortex, fourier1) =', data_vort_fourier1[2, i]
+            phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3  # in mrad -> *1000
+            data_vort_fourier1[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+            # Fourier padding 10:
+            padding = 10
+            start_time = time.clock()
+            phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+            data_vort_fourier10[2, i] = time.clock() - start_time
+            print '------time (vortex, fourier10) =', data_vort_fourier10[2, i]
+            phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3  # in mrad -> *1000
+            data_vort_fourier10[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+            # Real space slab:
+            start_time = time.clock()
+            phase_num_vort = pm.phase_mag_real(res, projection_vort, 'slab')
+            data_vort_real_s[2, i] = time.clock() - start_time
+            print '------time (vortex, real slab) =', data_vort_real_s[2, i]
+            phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3  # in mrad -> *1000
+            data_vort_real_s[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+            # Real space disc:
+            start_time = time.clock()
+            phase_num_vort = pm.phase_mag_real(res, projection_vort, 'disc')
+            data_vort_real_d[2, i] = time.clock() - start_time
+            print '------time (vortex, real disc) =', data_vort_real_d[2, i]
+            phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3  # in mrad -> *1000
+            data_vort_real_d[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+            
+        print '--SHELVE METHOD DATA'
+        data_shelve[key] = (data_disc_fourier0, data_vort_fourier0,
+                            data_disc_fourier1, data_vort_fourier1,
+                            data_disc_fourier10, data_vort_fourier10,
+                            data_disc_real_s, data_vort_real_s,
+                            data_disc_real_d, data_vort_real_d)
+    
+    print '--PLOT/SAVE METHOD DATA'
+    
+    # row and column sharing
+    fig, axes = plt.subplots(2, 2, sharex='col', sharey='row', figsize=(16, 7))
+#    ax1.plot(x, y)
+#    ax1.set_title('Sharing x per column, y per row')
+#    ax2.scatter(x, y)
+#    ax3.scatter(x, 2 * y ** 2 - 1, color='r')
+#    ax4.plot(x, 2 * y ** 2 - 1, color='r')
+    
+    
+    
+    # Plot duration against res (disc):
+#    fig = plt.figure()
+#    axis = fig.add_subplot(1, 1, 1)
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    axes[1, 0].set_yscale('log')
+    axes[1, 0].plot(data_disc_fourier0[0], data_disc_fourier0[1], '--b+', label='Fourier padding=0')
+    axes[1, 0].plot(data_disc_fourier1[0], data_disc_fourier1[1], '--bx', label='Fourier padding=1')
+    axes[1, 0].plot(data_disc_fourier10[0], data_disc_fourier10[1], '--b*', label='Fourier padding=10')
+    axes[1, 0].plot(data_disc_real_s[0], data_disc_real_s[1], '--rs', label='Real space (slab)')
+    axes[1, 0].plot(data_disc_real_d[0], data_disc_real_d[1], '--ro', label='Real space (disc)')
+    axes[1, 0].set_title('Variation of the resolution (disc)', fontsize=18)
+#    axes[1, 0].set_xlabel('res [nm]', fontsize=15)
+#    axes[1, 0].set_ylabel('RMS [mrad]', fontsize=15)
+#    axes[1, 0].set_xlim(-0.5, 16.5)
+#    axes[1, 0].legend(loc=4)
+#    plt.tick_params(axis='both', which='major', labelsize=14)
+#    plt.savefig(directory + '/ch5-3-disc_RMS_against_res.png', bbox_inches='tight')
+    # Plot RMS against res (disc):
+#    fig = plt.figure()
+#    axis = fig.add_subplot(1, 1, 1)
+    axes[0, 0].set_yscale('log')
+    axes[0, 0].plot(data_disc_fourier0[0], data_disc_fourier0[2], '--b+', label='Fourier padding=0')
+    axes[0, 0].plot(data_disc_fourier1[0], data_disc_fourier1[2], '--bx', label='Fourier padding=1')
+    axes[0, 0].plot(data_disc_fourier10[0], data_disc_fourier10[2], '--b*', label='Fourier padding=10')
+    axes[0, 0].plot(data_disc_real_s[0], data_disc_real_s[2], '--rs', label='Real space (slab)')
+    axes[0, 0].plot(data_disc_real_d[0], data_disc_real_d[2], '--ro', label='Real space (disc)')
+    axes[0, 0].set_title('Variation of the resolution (disc)', fontsize=18)
+#    axes[0, 0].set_xlabel('res [nm]', fontsize=15)
+    axes[0, 0].set_ylabel('duration [s]', fontsize=15)
+#    axes[0, 0].set_xlim(-0.5, 16.5)
+#    axes[0, 0].legend(loc=1)
+#    plt.tick_params(axis='both', which='major', labelsize=14)
+#    plt.savefig(directory + '/ch5-3-disc_duration_against_res.png', bbox_inches='tight')
+
+    # Plot duration against res (vortex):
+#    fig = plt.figure()
+#    axis = fig.add_subplot(1, 1, 1)
+    axes[1, 1].set_yscale('log')
+    axes[1, 1].plot(data_vort_fourier0[0], data_vort_fourier0[1], '--b+', label='Fourier padding=0')
+    axes[1, 1].plot(data_vort_fourier1[0], data_vort_fourier1[1], '--bx', label='Fourier padding=1')
+    axes[1, 1].plot(data_vort_fourier10[0], data_vort_fourier10[1], '--b*', label='Fourier padding=10')
+    axes[1, 1].plot(data_vort_real_s[0], data_vort_real_s[1], '--rs', label='Real space (slab)')
+    axes[1, 1].plot(data_vort_real_d[0], data_vort_real_d[1], '--ro', label='Real space (disc)')
+#    axes[1, 1].set_title('Variation of the resolution (vortex)', fontsize=18)
+    axes[1, 1].set_xlabel('res [nm]', fontsize=15)
+#    axes[1, 1].set_ylabel('RMS [mrad]', fontsize=15)
+    axes[1, 1].set_xlim(-0.5, 16.5)
+#    axes[1, 1].legend(loc=4)
+#    plt.tick_params(axis='both', which='major', labelsize=14)
+#    plt.savefig(directory + '/ch5-3-vortex_RMS_against_res.png', bbox_inches='tight')
+    # Plot RMS against res (vort):
+#    fig = plt.figure()
+#    axis = fig.add_subplot(1, 1, 1)
+    axes[0, 1].set_yscale('log')
+    axes[0, 1].plot(data_vort_fourier0[0], data_vort_fourier0[2], '--b+', label='Fourier padding=0')
+    axes[0, 1].plot(data_vort_fourier1[0], data_vort_fourier1[2], '--bx', label='Fourier padding=1')
+    axes[0, 1].plot(data_vort_fourier10[0], data_vort_fourier10[2], '--b*', label='Fourier padding=10')
+    axes[0, 1].plot(data_vort_real_s[0], data_vort_real_s[2], '--rs', label='Real space (slab)')
+    axes[0, 1].plot(data_vort_real_d[0], data_vort_real_d[2], '--ro', label='Real space (disc)')
+    axes[0, 1].set_title('Variation of the resolution (vortex)', fontsize=18)
+#    axes[0, 1].set_xlabel('res [nm]', fontsize=15)
+#    axes[0, 1].set_ylabel('duration [s]', fontsize=15)
+    axes[0, 1].set_xlim(-0.5, 16.5)
+    axes[0, 1].legend(loc=1)
+#    plt.tick_params(axis='both', which='major', labelsize=14)
+    plt.savefig(directory + '/ch5-3-vortex_duration_against_res.png', bbox_inches='tight')
+
+
+
+
+
+
+
+
+
+
+    ###############################################################################################
+    print 'CLOSING SHELVE\n'
+    # Close shelve:
+    data_shelve.close()
+    
+    ###############################################################################################
+
+
+if __name__ == "__main__":
+    try:
+        run()
+    except:
+        type, value, tb = sys.exc_info()
+        traceback.print_exc()
+        pdb.post_mortem(tb)