Skip to content
Snippets Groups Projects
Commit 78d70d7c authored by Jan Caron's avatar Jan Caron
Browse files

scripts: Splitted the main in several smaller scripts 

test: changed the dataloader TestCase
pyramex: first experiments with Cython and SCons
dataloader: changed the order in which the dimensions are saved,
        z_len is no longer part of the magnetization (this was the reason why
        we had to divide by res, which is equal to z_len for calculating coeff)
magcreator: shapes are now separate functions, mag_shape is now an argument for
        one unified create_hom_mag-function
phasemap: unification of the phasemap-functions in real space (just the pixel-
        field is calculated different now)
parent 0afb984a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 514 additions and 770 deletions
SConstruct 0 → 100644
+ 25
− 0
View file @ 78d70d7c
import os
import distutils.sysconfig
import sys
PYTHON_LIBPATH = distutils.sysconfig.get_python_lib()
PYTHON_INCPATH = distutils.sysconfig.get_python_inc()
PYTHON_LIBRARY = "python" + sys.version[:3]
print 'PYTHON_LIBPATH: ' + PYTHON_LIBPATH
print 'PYTHON_INCPATH: ' + PYTHON_INCPATH
print 'PYTHON_LIBRARY: ' + PYTHON_LIBRARY
env = Environment(ENV = os.environ)
if ARGUMENTS.get('VERBOSE') != '1':
env['CCCOMSTR'] = "Compiling $TARGET"
env['LINKCOMSTR'] = "Linking $TARGET"
Progress('$TARGET\r', overwrite=True)
env.AppendUnique(LIBPATH=[PYTHON_LIBPATH], CPPPATH=[PYTHON_INCPATH])
env.Program('hello', 'hellotest.c', CPPPATH='.')
env.Decider('MD5-timestamp')
\ No newline at end of file
scripts/compare_methods.py
+ 37
− 90
View file @ 78d70d7c
......@@ -25,75 +25,42 @@ def phase_from_mag():
None
'''
'''INPUT'''
# TODO: Input via GUI
filename = '../output/output.txt'
b_0 = 10.0 # in T
b_0 = 1.0 # in T
v_0 = 0 # TODO: units?
v_acc = 30000 # in V
padding = 20
density = 100
density = 10
dim = (50, 50) # in px (y,x)
res = 1.0 # in nm
res = 10.0 # in nm
beta = pi/4
plot_mag_distr = True
# Slab:
shape_fun = mc.slab
center = (24, 24) # in px (y,x) index starts with 0!
width = (25, 25) # in px (y,x)
params = (center, width)
# # Disc:
# shape_fun = mc.disc
# center = (4, 4) # in px (y,x)
# radius = 2.5
# params = (center, radius)
# # Filament:
# shape_fun = mc.filament
# pos = 5
# x_or_y = 'y'
# params = (pos, x_or_y)
# # Alternating Filaments:
# shape_fun = mc.alternating_filaments
# spacing = 5
# x_or_y = 'y'
# params = (spacing, x_or_y)
# # Single Pixel:
# shape_fun = mc.single_pixel
# pixel = (5, 5)
# params = pixel
# '''CREATE LOGO'''
# mc.create_logo(128, res, beta, filename, plot_mag_distr)
# mag_data = dl.MagDataLLG(filename)
# phase= pm.real_space_slab(mag_data, b_0)
# holo = hi.holo_image(phase, mag_data.res, density)
# hi.display_holo(holo, '')
radius = 12.5
filename = '../output/output.txt'
geometry = 'slab'
if geometry == 'slab':
mag_shape = mc.slab(dim, center, width)
phase_ana = an.phasemap_slab(dim, res, beta, center, width, b_0)
elif geometry == 'slab':
mag_shape = mc.disc(dim, center, radius)
phase_ana = an.phasemap_disc(dim, res, beta, center, radius, b_0)
holo_ana = hi.holo_image(phase_ana, res, density)
display_combined(phase_ana, res, holo_ana, 'Analytic Solution')
'''CREATE MAGNETIC DISTRIBUTION'''
mc.create_hom_mag(dim, res, beta, shape_fun, params,
filename, plot_mag_distr)
mc.create_hom_mag(dim, res, beta, mag_shape, filename)
'''LOAD MAGNETIC DISTRIBUTION'''
mag_data = dl.MagDataLLG(filename)
# TODO: get it to work:
phase_el = pm.phase_elec(mag_data, v_0=1, v_acc=200000)
phi_cos_real_slab = np.zeros((2*dim[0]-1, 2*dim[1]-1))
pm.real_space_slab(mag_data, b_0, pixel_map=phi_cos_real_slab)
pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
phi_cos_real_disc = np.zeros((2*dim[0]-1, 2*dim[1]-1))
pm.real_space_disc(mag_data, b_0, pixel_map=phi_cos_real_disc)
pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
phi_cos_diff = phi_cos_real_slab - phi_cos_real_disc
pm.display_phase(phi_cos_diff, mag_data.res, 'Difference: One Pixel Slab - Disc')
phase_el = pm.phase_elec(mag_data, v_0, v_acc)
'''NUMERICAL SOLUTION'''
# numerical solution Fourier Space:
......@@ -102,54 +69,34 @@ def phase_from_mag():
toc = time.clock()
print 'Time for Fourier Space Approach: ' + str(toc - tic)
holo_fft = hi.holo_image(phase_fft, mag_data.res, density)
display_combined(phase_fft, mag_data.res, holo_fft,
'Fourier Space Approach')
display_combined(phase_fft, mag_data.res, holo_fft, 'Fourier Space Approach')
# numerical solution Real Space (Slab):
jacobi = np.zeros((2*dim[0]*dim[1], dim[0]*dim[1]))
tic = time.clock()
phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
phase_slab = pm.real_space(mag_data, 'slab', b_0)
toc = time.clock()
np.savetxt('../output/jacobi.npy', jacobi)
print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
holo_real_slab = hi.holo_image(phase_real_slab, mag_data.res, density)
display_combined(phase_real_slab, mag_data.res, holo_real_slab,
'Real Space Approach (Slab)')
holo_slab = hi.holo_image(phase_slab, mag_data.res, density)
display_combined(phase_slab, mag_data.res, holo_slab, 'Real Space Approach (Slab)')
# numerical solution Real Space (Disc):
tic = time.clock()
phase_real_disc = pm.real_space_disc(mag_data, b_0)
phase_disc = pm.real_space(mag_data, 'disc', b_0)
toc = time.clock()
print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
holo_real_disc = hi.holo_image(phase_real_disc, mag_data.res, density)
display_combined(phase_real_disc, mag_data.res, holo_real_disc,
'Real Space Approach (Disc)')
'''ANALYTIC SOLUTION'''
# analytic solution slab:
phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
holo = hi.holo_image(phase, res, density)
display_combined(phase, res, holo, 'Analytic Solution - Slab')
# # analytic solution disc:
# phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
# holo = hi.holo_image(phase, res, density)
# display_combined(phase, res, holo, 'Analytic Solution - Disc')
# # analytic solution sphere:
# phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
# holo = hi.holo_image(phase, res, density)
# display_combined(phase, res, holo, 'Analytic Solution - Sphere')
holo_disc = hi.holo_image(phase_disc, mag_data.res, density)
display_combined(phase_disc, mag_data.res, holo_disc, 'Real Space Approach (Disc)')
'''DIFFERENCES'''
diff_fft_to_ana = phase_fft - phase
diff_real_slab_to_ana = phase_real_slab - phase
diff_real_disc_to_ana = phase_real_disc - phase
diff_slab_to_disc = phase_real_disc - phase_real_slab
pm.display_phase(diff_fft_to_ana, res, 'Difference: FFT - Analytic')
pm.display_phase(diff_real_slab_to_ana, res, 'Difference: Slab - Analytic')
pm.display_phase(diff_real_disc_to_ana, res, 'Difference: Disc - Analytic')
pm.display_phase(diff_slab_to_disc, res, 'Difference: Disc - Slab')
# TODO: Delete
# import pdb; pdb.set_trace()
diff_fft = phase_fft - phase_ana
diff_slab = phase_slab - phase_ana
diff_disc = phase_disc - phase_ana
rms_fft = np.sqrt((diff_fft**2).mean())
rms_slab = np.sqrt((diff_slab**2).mean())
rms_disc = np.sqrt((diff_disc**2).mean())
pm.display_phase(diff_fft, res, 'FFT - Analytic (RMS = ' + '{:3.2e}'.format(rms_fft) + ')')
pm.display_phase(diff_slab, res, 'Slab - Analytic (RMS = ' +'{:3.2e}'.format(rms_slab) + ')')
pm.display_phase(diff_disc, res, 'Disc - Analytic (RMS = ' + '{:3.2e}'.format(rms_disc) + ')')
def display_combined(phase, res, holo, title):
......
scripts/compare_pixel_fields.py 0 → 100644
+ 44
− 0
View file @ 78d70d7c
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 03 11:15:38 2013
@author: Jan
"""
import numpy as np
import pyramid.phasemap as pm
import pdb, traceback, sys
def compare_pixel_fields():
'''Calculate and display the phase map from a given magnetization.
Arguments:
None
Returns:
None
'''
# TODO: Input via GUI
b_0 = 1.0 # in T
res = 10.0
dim = (10, 10)
x_big = np.linspace(-(dim[1]-1), dim[1]-1, num=2*dim[1]-1)
y_big = np.linspace(-(dim[0]-1), dim[0]-1, num=2*dim[0]-1)
xx_big, yy_big = np.meshgrid(x_big, y_big)
phi_cos_real_slab = pm.phi_pixel('slab', xx_big, yy_big, res, b_0)
pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
phi_cos_real_disc = pm.phi_pixel('disc', xx_big, yy_big, res, b_0)
pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
phi_cos_diff = phi_cos_real_disc - phi_cos_real_slab
pm.display_phase(phi_cos_diff, res, 'Phase of one Pixel-Disc (Cos - Part)')
if __name__ == "__main__":
try:
compare_pixel_fields()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
\ No newline at end of file
scripts/create_logo.py
+ 20
− 6
View file @ 78d70d7c
......@@ -10,6 +10,7 @@ import pyramid.dataloader as dl
import pyramid.phasemap as pm
import pyramid.holoimage as hi
import pdb, traceback, sys
import numpy as np
from numpy import pi
......@@ -21,21 +22,34 @@ def create_logo():
None
'''
filename = '../output/logo_magnetization.txt'
filename = '../output/mag_distr_logo.txt'
b_0 = 1.0 # in T
res = 10.0 # in nm
beta = pi/2
beta = pi/2 # in rad
density = 10
dim = (128, 128)
x = range(dim[1])
y = range(dim[0])
xx, yy = np.meshgrid(x, y)
bottom = (yy >= 0.25*dim[0])
left = (yy <= 0.75/0.5 * dim[0]/dim[1] * xx)
right = np.fliplr(left)
mag_shape = np.logical_and(np.logical_and(left, right), bottom)
mc.create_logo(128, res, beta, filename)
mc.create_hom_mag(dim, res, beta, mag_shape, filename)
mag_data = dl.MagDataLLG(filename)
phase= pm.real_space_slab(mag_data, b_0)
phase= pm.real_space(mag_data, 'slab', b_0)
holo = hi.holo_image(phase, mag_data.res, density)
hi.display_holo(holo, '')
hi.display_holo(holo, 'PYRAMID - LOGO')
if __name__ == "__main__":
# parser = argparse.ArgumentParser(description='Create the PYRAMID logo')
# parser.add_argument('-d','--dimensions', help='Logo dimensions.', required=False)
# args = parser.parse_args()
# if args.dimensions is None:
# args.dimensions = (128,128)
try:
create_logo()
except:
......
scripts/create_sample.py 0 → 100644
+ 59
− 0
View file @ 78d70d7c
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 03 11:15:38 2013
@author: Jan
"""
import pyramid.magcreator as mc
import pdb, traceback, sys
from numpy import pi
def create_sample():
'''Calculate and display the phase map from a given magnetization.
Arguments:
None
Returns:
None
'''
# TODO: Input via GUI
key = 'slab'
filename = '../output/mag_distr_' + key + '.txt'
dim = (50, 50) # in px (y,x)
res = 1.0 # in nm
beta = pi/4
plot_mag_distr = True
center = (24, 24) # in px (y,x) index starts with 0!
width = (25, 25) # in px (y,x)
radius = 12.5 # in px
pos = 24 # in px
spacing = 5 # in px
x_or_y = 'y'
pixel = (24, 24) # in px
if key == 'slab':
mag_shape = mc.slab(dim, center, width)
elif key == 'disc':
mag_shape = mc.disc(dim, center, radius)
elif key == 'filament':
mag_shape = mc.filament(dim, pos, x_or_y)
elif key == 'alternating_filaments':
mag_shape = mc.alternating_filaments(dim, spacing, x_or_y)
elif key == 'pixel':
mag_shape = mc.single_pixel(dim, pixel)
mc.create_hom_mag(dim, res, beta, mag_shape, filename, plot_mag_distr)
if __name__ == "__main__":
try:
create_sample()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
\ No newline at end of file
scripts/create_samples.py deleted 100644 → 0
+ 0
− 172
View file @ 0afb984a
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 03 11:15:38 2013
@author: Jan
"""
import numpy as np
import matplotlib.pyplot as plt
import pyramid.magcreator as mc
import pyramid.dataloader as dl
import pyramid.phasemap as pm
import pyramid.holoimage as hi
import pyramid.analytic as an
import time
import pdb, traceback, sys
from numpy import pi
def phase_from_mag():
'''Calculate and display the phase map from a given magnetization.
Arguments:
None
Returns:
None
'''
'''INPUT'''
# TODO: Input via GUI
filename = '../output/output.txt'
b_0 = 10.0 # in T
v_0 = 0 # TODO: units?
v_acc = 30000 # in V
padding = 20
density = 100
dim = (50, 50) # in px (y,x)
res = 1.0 # in nm
beta = pi/4
plot_mag_distr = True
# Slab:
shape_fun = mc.slab
center = (24, 24) # in px (y,x) index starts with 0!
width = (25, 25) # in px (y,x)
params = (center, width)
# # Disc:
# shape_fun = mc.disc
# center = (4, 4) # in px (y,x)
# radius = 2.5
# params = (center, radius)
# # Filament:
# shape_fun = mc.filament
# pos = 5
# x_or_y = 'y'
# params = (pos, x_or_y)
# # Alternating Filaments:
# shape_fun = mc.alternating_filaments
# spacing = 5
# x_or_y = 'y'
# params = (spacing, x_or_y)
# # Single Pixel:
# shape_fun = mc.single_pixel
# pixel = (5, 5)
# params = pixel
# '''CREATE LOGO'''
# mc.create_logo(128, res, beta, filename, plot_mag_distr)
# mag_data = dl.MagDataLLG(filename)
# phase= pm.real_space_slab(mag_data, b_0)
# holo = hi.holo_image(phase, mag_data.res, density)
# hi.display_holo(holo, '')
'''CREATE MAGNETIC DISTRIBUTION'''
mc.create_hom_mag(dim, res, beta, shape_fun, params,
filename, plot_mag_distr)
'''LOAD MAGNETIC DISTRIBUTION'''
mag_data = dl.MagDataLLG(filename)
# TODO: get it to work:
phase_el = pm.phase_elec(mag_data, v_0=1, v_acc=200000)
phi_cos_real_slab = np.zeros((2*dim[0]-1, 2*dim[1]-1))
pm.real_space_slab(mag_data, b_0, pixel_map=phi_cos_real_slab)
pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
phi_cos_real_disc = np.zeros((2*dim[0]-1, 2*dim[1]-1))
pm.real_space_disc(mag_data, b_0, pixel_map=phi_cos_real_disc)
pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
phi_cos_diff = phi_cos_real_slab - phi_cos_real_disc
pm.display_phase(phi_cos_diff, mag_data.res, 'Difference: One Pixel Slab - Disc')
'''NUMERICAL SOLUTION'''
# numerical solution Fourier Space:
tic = time.clock()
phase_fft = pm.fourier_space(mag_data, b_0, padding)
toc = time.clock()
print 'Time for Fourier Space Approach: ' + str(toc - tic)
holo_fft = hi.holo_image(phase_fft, mag_data.res, density)
display_combined(phase_fft, mag_data.res, holo_fft,
'Fourier Space Approach')
# numerical solution Real Space (Slab):
jacobi = np.zeros((2*dim[0]*dim[1], dim[0]*dim[1]))
tic = time.clock()
phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
toc = time.clock()
np.savetxt('../output/jacobi.npy', jacobi)
print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
holo_real_slab = hi.holo_image(phase_real_slab, mag_data.res, density)
display_combined(phase_real_slab, mag_data.res, holo_real_slab,
'Real Space Approach (Slab)')
# numerical solution Real Space (Disc):
tic = time.clock()
phase_real_disc = pm.real_space_disc(mag_data, b_0)
toc = time.clock()
print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
holo_real_disc = hi.holo_image(phase_real_disc, mag_data.res, density)
display_combined(phase_real_disc, mag_data.res, holo_real_disc,
'Real Space Approach (Disc)')
'''ANALYTIC SOLUTION'''
# analytic solution slab:
phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
holo = hi.holo_image(phase, res, density)
display_combined(phase, res, holo, 'Analytic Solution - Slab')
# # analytic solution disc:
# phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
# holo = hi.holo_image(phase, res, density)
# display_combined(phase, res, holo, 'Analytic Solution - Disc')
# # analytic solution sphere:
# phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
# holo = hi.holo_image(phase, res, density)
# display_combined(phase, res, holo, 'Analytic Solution - Sphere')
'''DIFFERENCES'''
diff_fft_to_ana = phase_fft - phase
diff_real_slab_to_ana = phase_real_slab - phase
diff_real_disc_to_ana = phase_real_disc - phase
diff_slab_to_disc = phase_real_disc - phase_real_slab
pm.display_phase(diff_fft_to_ana, res, 'Difference: FFT - Analytic')
pm.display_phase(diff_real_slab_to_ana, res, 'Difference: Slab - Analytic')
pm.display_phase(diff_real_disc_to_ana, res, 'Difference: Disc - Analytic')
pm.display_phase(diff_slab_to_disc, res, 'Difference: Disc - Slab')
# TODO: Delete
# import pdb; pdb.set_trace()
def display_combined(phase, res, holo, title):
fig = plt.figure(figsize=(14, 7))
fig.suptitle(title, fontsize=20)
holo_axis = fig.add_subplot(1,2,1, aspect='equal')
hi.display_holo(holo, 'Holography Image', holo_axis)
phase_axis = fig.add_subplot(1,2,2, aspect='equal')
pm.display_phase(phase, res, 'Phasemap', phase_axis)
if __name__ == "__main__":
try:
phase_from_mag()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
\ No newline at end of file
scripts/cython_example.py 0 → 100644
+ 55
− 0
View file @ 78d70d7c
# -*- coding: utf-8 -*-
"""
Created on Fri May 03 09:20:51 2013
@author: Jan
"""
import cython
import numpy as np
n= 100
A = np.array(range(n**2)).reshape((n,n))
def naive_convolve(f, g):
# f is an image and is indexed by (v, w)
# g is a filter kernel and is indexed by (s, t),
# it needs odd dimensions
# h is the output image and is indexed by (x, y),
# it is not cropped
if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
raise ValueError("Only odd dimensions on filter supported")
# smid and tmid are number of pixels between the center pixel
# and the edge, ie for a 5x5 filter they will be 2.
#
# The output size is calculated by adding smid, tmid to each
# side of the dimensions of the input image.
vmax = f.shape[0]
wmax = f.shape[1]
smax = g.shape[0]
tmax = g.shape[1]
smid = smax // 2
tmid = tmax // 2
xmax = vmax + 2*smid
ymax = wmax + 2*tmid
# Allocate result image.
h = np.zeros([xmax, ymax], dtype=f.dtype)
# Do convolution
for x in range(xmax):
for y in range(ymax):
# Calculate pixel value for h at (x,y). Sum one component
# for each pixel (s, t) of the filter g.
s_from = max(smid - x, -smid)
s_to = min((xmax - x) - smid, smid + 1)
t_from = max(tmid - y, -tmid)
t_to = min((ymax - y) - tmid, tmid + 1)
value = 0
for s in range(s_from, s_to):
for t in range(t_from, t_to):
v = x - smid + s
w = y - tmid + t
value += g[smid - s, tmid - t] * f[v, w]
h[x, y] = value
return h
\ No newline at end of file
scripts/get_jacobi.py
+ 13
− 125
View file @ 78d70d7c
......@@ -6,12 +6,9 @@ Created on Wed Apr 03 11:15:38 2013
"""
import numpy as np
import matplotlib.pyplot as plt
import pyramid.magcreator as mc
import pyramid.dataloader as dl
import pyramid.phasemap as pm
import pyramid.holoimage as hi
import pyramid.analytic as an
import time
import pdb, traceback, sys
from numpy import pi
......@@ -25,147 +22,38 @@ def phase_from_mag():
None
'''
'''INPUT'''
# TODO: Input via GUI
filename = '../output/output.txt'
b_0 = 1.0 # in T
v_0 = 0 # TODO: units?
v_acc = 30000 # in V
padding = 20
density = 100
dim = (50, 50) # in px (y,x)
res = 1.0 # in nm
beta = pi/4
plot_mag_distr = True
b_0 = 1.0 # in T
dim = (3, 3) # in px (y,x)
res = 10.0 # in nm
beta = pi/4
# Slab:
shape_fun = mc.slab
center = (24, 24) # in px (y,x) index starts with 0!
width = (25, 25) # in px (y,x)
params = (center, width)
# # Disc:
# shape_fun = mc.disc
# center = (4, 4) # in px (y,x)
# radius = 2.5
# params = (center, radius)
# # Filament:
# shape_fun = mc.filament
# pos = 5
# x_or_y = 'y'
# params = (pos, x_or_y)
# # Alternating Filaments:
# shape_fun = mc.alternating_filaments
# spacing = 5
# x_or_y = 'y'
# params = (spacing, x_or_y)
# # Single Pixel:
# shape_fun = mc.single_pixel
# pixel = (5, 5)
# params = pixel
# '''CREATE LOGO'''
# mc.create_logo(128, res, beta, filename, plot_mag_distr)
# mag_data = dl.MagDataLLG(filename)
# phase= pm.real_space_slab(mag_data, b_0)
# holo = hi.holo_image(phase, mag_data.res, density)
# hi.display_holo(holo, '')
center = (1, 1) # in px (y,x) index starts with 0!
width = (1, 1) # in px (y,x)
mag_shape = mc.slab(dim, center, width)
'''CREATE MAGNETIC DISTRIBUTION'''
mc.create_hom_mag(dim, res, beta, shape_fun, params,
filename, plot_mag_distr)
mc.create_hom_mag(dim, res, beta, mag_shape, filename)
'''LOAD MAGNETIC DISTRIBUTION'''
mag_data = dl.MagDataLLG(filename)
# TODO: get it to work:
phase_el = pm.phase_elec(mag_data, v_0=1, v_acc=200000)
phi_cos_real_slab = np.zeros((2*dim[0]-1, 2*dim[1]-1))
pm.real_space_slab(mag_data, b_0, pixel_map=phi_cos_real_slab)
pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
phi_cos_real_disc = np.zeros((2*dim[0]-1, 2*dim[1]-1))
pm.real_space_disc(mag_data, b_0, pixel_map=phi_cos_real_disc)
pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
phi_cos_diff = phi_cos_real_slab - phi_cos_real_disc
pm.display_phase(phi_cos_diff, mag_data.res, 'Difference: One Pixel Slab - Disc')
'''NUMERICAL SOLUTION'''
# numerical solution Fourier Space:
tic = time.clock()
phase_fft = pm.fourier_space(mag_data, b_0, padding)
toc = time.clock()
print 'Time for Fourier Space Approach: ' + str(toc - tic)
holo_fft = hi.holo_image(phase_fft, mag_data.res, density)
display_combined(phase_fft, mag_data.res, holo_fft,
'Fourier Space Approach')
# numerical solution Real Space (Slab):
jacobi = np.zeros((2*dim[0]*dim[1], dim[0]*dim[1]))
jacobi = np.zeros((dim[0]*dim[1], 2*dim[0]*dim[1]))
tic = time.clock()
phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
pm.real_space(mag_data, 'slab', b_0, jacobi=jacobi)
toc = time.clock()
np.savetxt('../output/jacobi.npy', jacobi)
print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
holo_real_slab = hi.holo_image(phase_real_slab, mag_data.res, density)
display_combined(phase_real_slab, mag_data.res, holo_real_slab,
'Real Space Approach (Slab)')
# numerical solution Real Space (Disc):
tic = time.clock()
phase_real_disc = pm.real_space_disc(mag_data, b_0)
toc = time.clock()
print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
holo_real_disc = hi.holo_image(phase_real_disc, mag_data.res, density)
display_combined(phase_real_disc, mag_data.res, holo_real_disc,
'Real Space Approach (Disc)')
'''ANALYTIC SOLUTION'''
# analytic solution slab:
phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
holo = hi.holo_image(phase, res, density)
display_combined(phase, res, holo, 'Analytic Solution - Slab')
# # analytic solution disc:
# phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
# holo = hi.holo_image(phase, res, density)
# display_combined(phase, res, holo, 'Analytic Solution - Disc')
# # analytic solution sphere:
# phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
# holo = hi.holo_image(phase, res, density)
# display_combined(phase, res, holo, 'Analytic Solution - Sphere')
'''DIFFERENCES'''
diff_fft_to_ana = phase_fft - phase
diff_real_slab_to_ana = phase_real_slab - phase
diff_real_disc_to_ana = phase_real_disc - phase
diff_slab_to_disc = phase_real_disc - phase_real_slab
pm.display_phase(diff_fft_to_ana, res, 'Difference: FFT - Analytic')
pm.display_phase(diff_real_slab_to_ana, res, 'Difference: Slab - Analytic')
pm.display_phase(diff_real_disc_to_ana, res, 'Difference: Disc - Analytic')
pm.display_phase(diff_slab_to_disc, res, 'Difference: Disc - Slab')
# TODO: Delete
# import pdb; pdb.set_trace()
def display_combined(phase, res, holo, title):
fig = plt.figure(figsize=(14, 7))
fig.suptitle(title, fontsize=20)
holo_axis = fig.add_subplot(1,2,1, aspect='equal')
hi.display_holo(holo, 'Holography Image', holo_axis)
print 'Time for Real Space Approach with Jacobi-Matrix (Slab): ' + str(toc - tic)
phase_axis = fig.add_subplot(1,2,2, aspect='equal')
pm.display_phase(phase, res, 'Phasemap', phase_axis)
return jacobi
if __name__ == "__main__":
try:
phase_from_mag()
jacobi = phase_from_mag()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
......
scripts/main.py deleted 100644 → 0
+ 0
− 172
View file @ 0afb984a
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 03 11:15:38 2013
@author: Jan
"""
import numpy as np
import matplotlib.pyplot as plt
import pyramid.magcreator as mc
import pyramid.dataloader as dl
import pyramid.phasemap as pm
import pyramid.holoimage as hi
import pyramid.analytic as an
import time
import pdb, traceback, sys
from numpy import pi
def phase_from_mag():
'''Calculate and display the phase map from a given magnetization.
Arguments:
None
Returns:
None
'''
'''INPUT'''
# TODO: Input via GUI
filename = '../output/output.txt'
b_0 = 1.0 # in T
v_0 = 0 # TODO: units?
v_acc = 30000 # in V
padding = 20
density = 10
dim = (50, 50) # in px (y,x)
res = 10.0 # in nm
beta = pi/4
plot_mag_distr = True
# Slab:
shape_fun = mc.slab
center = (24, 24) # in px (y,x) index starts with 0!
width = (25, 25) # in px (y,x)
params = (center, width)
# # Disc:
# shape_fun = mc.disc
# center = (4, 4) # in px (y,x)
# radius = 2.5
# params = (center, radius)
# # Filament:
# shape_fun = mc.filament
# pos = 5
# x_or_y = 'y'
# params = (pos, x_or_y)
# # Alternating Filaments:
# shape_fun = mc.alternating_filaments
# spacing = 5
# x_or_y = 'y'
# params = (spacing, x_or_y)
# # Single Pixel:
# shape_fun = mc.single_pixel
# pixel = (5, 5)
# params = pixel
# '''CREATE LOGO'''
# mc.create_logo(128, res, beta, filename, plot_mag_distr)
# mag_data = dl.MagDataLLG(filename)
# phase= pm.real_space_slab(mag_data, b_0)
# holo = hi.holo_image(phase, mag_data.res, density)
# hi.display_holo(holo, '')
'''CREATE MAGNETIC DISTRIBUTION'''
mc.create_hom_mag(dim, res, beta, shape_fun, params,
filename, plot_mag_distr)
'''LOAD MAGNETIC DISTRIBUTION'''
mag_data = dl.MagDataLLG(filename)
# TODO: get it to work:
phase_el = pm.phase_elec(mag_data, v_0=1, v_acc=200000)
phi_cos_real_slab = np.zeros((2*dim[0]-1, 2*dim[1]-1))
pm.real_space_slab(mag_data, b_0, pixel_map=phi_cos_real_slab)
pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
phi_cos_real_disc = np.zeros((2*dim[0]-1, 2*dim[1]-1))
pm.real_space_disc(mag_data, b_0, pixel_map=phi_cos_real_disc)
pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
phi_cos_diff = phi_cos_real_slab - phi_cos_real_disc
pm.display_phase(phi_cos_diff, mag_data.res, 'Difference: One Pixel Slab - Disc')
'''NUMERICAL SOLUTION'''
# numerical solution Fourier Space:
tic = time.clock()
phase_fft = pm.fourier_space(mag_data, b_0, padding)
toc = time.clock()
print 'Time for Fourier Space Approach: ' + str(toc - tic)
holo_fft = hi.holo_image(phase_fft, mag_data.res, density)
display_combined(phase_fft, mag_data.res, holo_fft,
'Fourier Space Approach')
# numerical solution Real Space (Slab):
jacobi = np.zeros((2*dim[0]*dim[1], dim[0]*dim[1]))
tic = time.clock()
phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
toc = time.clock()
np.savetxt('../output/jacobi.npy', jacobi)
print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
holo_real_slab = hi.holo_image(phase_real_slab, mag_data.res, density)
display_combined(phase_real_slab, mag_data.res, holo_real_slab,
'Real Space Approach (Slab)')
# numerical solution Real Space (Disc):
tic = time.clock()
phase_real_disc = pm.real_space_disc(mag_data, b_0)
toc = time.clock()
print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
holo_real_disc = hi.holo_image(phase_real_disc, mag_data.res, density)
display_combined(phase_real_disc, mag_data.res, holo_real_disc,
'Real Space Approach (Disc)')
'''ANALYTIC SOLUTION'''
# analytic solution slab:
phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
holo = hi.holo_image(phase, res, density)
display_combined(phase, res, holo, 'Analytic Solution - Slab')
# # analytic solution disc:
# phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
# holo = hi.holo_image(phase, res, density)
# display_combined(phase, res, holo, 'Analytic Solution - Disc')
# # analytic solution sphere:
# phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
# holo = hi.holo_image(phase, res, density)
# display_combined(phase, res, holo, 'Analytic Solution - Sphere')
'''DIFFERENCES'''
diff_fft_to_ana = phase_fft - phase
diff_real_slab_to_ana = phase_real_slab - phase
diff_real_disc_to_ana = phase_real_disc - phase
diff_slab_to_disc = phase_real_disc - phase_real_slab
pm.display_phase(diff_fft_to_ana, res, 'Difference: FFT - Analytic')
pm.display_phase(diff_real_slab_to_ana, res, 'Difference: Slab - Analytic')
pm.display_phase(diff_real_disc_to_ana, res, 'Difference: Disc - Analytic')
pm.display_phase(diff_slab_to_disc, res, 'Difference: Disc - Slab')
# TODO: Delete
# import pdb; pdb.set_trace()
def display_combined(phase, res, holo, title):
fig = plt.figure(figsize=(14, 7))
fig.suptitle(title, fontsize=20)
holo_axis = fig.add_subplot(1,2,1, aspect='equal')
hi.display_holo(holo, 'Holography Image', holo_axis)
phase_axis = fig.add_subplot(1,2,2, aspect='equal')
pm.display_phase(phase, res, 'Phasemap', phase_axis)
if __name__ == "__main__":
try:
phase_from_mag()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
\ No newline at end of file
src/pyramex/.sconsign.dblite 0 → 100644
+ 0
− 0
View file @ 78d70d7c
File added
src/pyramex/SConstruct 0 → 100644
+ 16
− 0
View file @ 78d70d7c
import distutils.sysconfig
import os
import sys
PYTHON_LIBPATH = distutils.sysconfig.get_python_lib()
PYTHON_INCPATH = distutils.sysconfig.get_python_inc()
PYTHON_LIBRARY = "python" + sys.version[:3]
cythonBuilder = Builder(action ='cython $SOURCE',
suffix = '.c',
src_suffix = '.pyx')
env = Environment(ENV = os.environ,
BUILDERS = {'Py2C':cythonBuilder})
env.Py2C('c1.c', 'c1.pyx')
\ No newline at end of file
src/pyramex/SConstructBACKUP 0 → 100644
+ 70
− 0
View file @ 78d70d7c
import distutils.sysconfig
import os
import sys
PYTHON_LIBPATH = distutils.sysconfig.get_python_lib()
PYTHON_INCPATH = distutils.sysconfig.get_python_inc()
PYTHON_LIBRARY = "python" + sys.version[:3]
print 'PYTHON_LIBPATH: ' + PYTHON_LIBPATH
print 'PYTHON_INCPATH: ' + PYTHON_INCPATH
print 'PYTHON_LIBRARY: ' + PYTHON_LIBRARY
env = Environment(ENV = os.environ)
env.AppendUnique(LIBPATH=[PYTHON_LIBPATH], CPPPATH=[PYTHON_INCPATH])
if ARGUMENTS.get('VERBOSE') != '1':
env['CCCOMSTR'] = "Compiling $TARGET"
env['LINKCOMSTR'] = "Linking $TARGET"
Progress('$TARGET\r', overwrite=True)
env.Decider('MD5-timestamp')
# This builder takes the cython source code and generates c code.
cythonBuilder = Builder(action ='cython $SOURCE',
suffix = '.c',
src_suffix = '.pyx')
# Setup the compilation environment.
env = Environment(PYEXT_USE_DISTUTILS=True,
ENV = os.environ,
tools=['mingw'],
BUILDERS={'Py2C':cythonBuilder},
CCFLAGS='-LC:\Python27\libs -IC:\Python27\include',
LINKFLAGS='-lpython2.7 -shared -o $SOURCE')
# Method that calls routines neccessary to create shared library.
def cythonPseudoBuilder(env,source):
cCode = env.Py2C(source)
#env.SharedLibrary(source,cCode,LIBPREFIX='')
env.AddMethod(cythonPseudoBuilder,'Cython')
env.Cython('c1')
#env.SharedLibrary('c1.pyd', ['c1.c'])
c1_module = env.LoadableModule('c1.pyd', ['c1.c'])
#import setup
#for ext in setup.get_extensions():
# print e.sources[:]
#env.Command('hello.pyd', '', 'python setup.py build_ext --compiler=mingw32 --inplace')
\ No newline at end of file
src/pyramex/c1.pyx 0 → 100644
+ 13
− 0
View file @ 78d70d7c
import math
def great_circle(float lon1,float lat1,float lon2,float lat2):
cdef float radius = 3956.0
cdef float pi = 3.14159265
cdef float x = pi/180.0
cdef float a,b,theta,c
a = (90.0-lat1)*(x)
b = (90.0-lat2)*(x)
theta = (lon2-lon1)*(x)
c = math.acos((math.cos(a)*math.cos(b)) + (math.sin(a)*math.sin(b)*math.cos(theta)))
return radius*c
\ No newline at end of file
src/pyramex/c2.pyx 0 → 100644
+ 16
− 0
View file @ 78d70d7c
cdef extern from "math.h":
float cosf(float theta)
float sinf(float theta)
float acosf(float theta)
def great_circle(float lon1,float lat1,float lon2,float lat2):
cdef float radius = 3956.0
cdef float pi = 3.14159265
cdef float x = pi/180.0
cdef float a,b,theta,c
a = (90.0-lat1)*(x)
b = (90.0-lat2)*(x)
theta = (lon2-lon1)*(x)
c = acosf((cosf(a)*cosf(b)) + (sinf(a)*sinf(b)*cosf(theta)))
return radius*c
\ No newline at end of file
src/pyramex/c3.pyx 0 → 100644
+ 23
− 0
View file @ 78d70d7c
cdef extern from "math.h":
float cosf(float theta)
float sinf(float theta)
float acosf(float theta)
cdef float _great_circle(float lon1,float lat1,float lon2,float lat2):
cdef float radius = 3956.0
cdef float pi = 3.14159265
cdef float x = pi/180.0
cdef float a,b,theta,c
a = (90.0-lat1)*(x)
b = (90.0-lat2)*(x)
theta = (lon2-lon1)*(x)
c = acosf((cosf(a)*cosf(b)) + (sinf(a)*sinf(b)*cosf(theta)))
return radius*c
def great_circle(float lon1,float lat1,float lon2,float lat2,int num):
cdef int i
cdef float x
for i from 0 <= i < num:
x = _great_circle(lon1,lat1,lon2,lat2)
return x
\ No newline at end of file
src/pyramex/hello.pyx 0 → 100644
+ 2
− 0
View file @ 78d70d7c
def say_hello_to(name):
print 'Hello %s!' % name
\ No newline at end of file
src/pyramex/setup.py 0 → 100644
+ 36
− 0
View file @ 78d70d7c
# -*- coding: utf-8 -*-
"""
Created on Fri May 03 10:27:04 2013
@author: Jan
"""
#from distutils.core import setup, Extension
#
#module1 = Extension('c1',
# sources = ['c1.pyx'])
#
#setup (name = 'pyramex',
# version = '1.0',
# description = 'This is a demo package',
# ext_modules = [module1])
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
def get_extensions():
return [
Extension('hello', ['hello.pyx']),
Extension('c1', ['c1.pyx']),
Extension('c2', ['c2.pyx']),
Extension('c3', ['c3.pyx'])
]
setup(
name = 'pyramex',
cmdclass = {'build_ext': build_ext},
ext_modules = get_extensions()
)
\ No newline at end of file
src/pyramid/dataloader.py
+ 6
− 6
View file @ 78d70d7c
......@@ -18,21 +18,21 @@ class MagDataLLG:
'''
scale = 1.0E-9 / 1.0E-2 #from cm to nm
data = np.genfromtxt(filename, skip_header=2)
data = np.genfromtxt(filename, skip_header=2)
x_dim, y_dim, z_dim = np.genfromtxt(filename, dtype=int,
skip_header=1,
skip_footer=len(data[:, 0]))
res = (data[1, 0] - data[0, 0]) / scale
x_len, y_len, z_len = [data[-1, i]/scale+res/2 for i in range(3)]
x_mag, y_mag, z_mag = [data[:, i].reshape(z_dim, y_dim, x_dim).mean(0)
*z_len for i in range(3,6)]
x_mag, y_mag, z_mag = [data[:, i].reshape(z_dim, y_dim, x_dim)
for i in range(3,6)]
#Reshape in Python and Igor is different,
#Python fills rows first, Igor columns!
self.filename = filename
self.res = res
self.dim = (x_dim, y_dim, z_dim)
self.length = (x_len, y_len, z_len)
self.magnitude = (x_mag, y_mag, z_mag)
self.dim = (z_dim, y_dim, x_dim)
self.length = (z_len, y_len, x_len)
self.magnitude = (z_mag, y_mag, x_mag)
def __str__(self):
'''Return the filename of the loaded file.
......
src/pyramid/magcreator.py
+ 27
− 80
View file @ 78d70d7c
# -*- coding: utf-8 -*-
"""Create simple LLG Files which describe magnetization in 2D (z-Dim=1)."""
import numpy as np
import matplotlib.pyplot as plt
from numpy import pi
def slab(dim, params):
def slab(dim, center, width):
'''Get the magnetic shape of a slab.
Arguments:
dim - the dimensions of the grid, shape(y, x)
......@@ -17,14 +15,13 @@ def slab(dim, params):
the magnetic shape as a 2D-boolean-array.
'''
center, width = params
shape_mag = np.array([[abs(x - center[1]) <= width[1] / 2
mag_shape = np.array([[abs(x - center[1]) <= width[1] / 2
and abs(y - center[0]) <= width[0] / 2
for x in range(dim[1])] for y in range(dim[0])])
return shape_mag
return mag_shape
def disc(dim, params):
def disc(dim, center, radius):
'''Get the magnetic shape of a disc.
Arguments:
dim - the dimensions of the grid, shape(y, x)
......@@ -34,13 +31,12 @@ def disc(dim, params):
the magnetic shape as a 2D-boolean-array.
'''
center, radius = params
shape_mag = np.array([[(np.hypot(x-center[1], y-center[0]) <= radius)
mag_shape = np.array([[(np.hypot(x-center[1], y-center[0]) <= radius)
for x in range(dim[1])] for y in range(dim[0])])
return shape_mag
return mag_shape
def filament(dim, params):
def filament(dim, pos, x_or_y):
'''Get the magnetic shape of a single filament in x or y direction.
Arguments:
dim - the dimensions of the grid, shape(y, x)
......@@ -51,16 +47,15 @@ def filament(dim, params):
the magnetic shape as a 2D-boolean-array.
'''
pos, x_or_y = params
shape_mag = np.zeros(dim)
mag_shape = np.zeros(dim)
if x_or_y == 'y':
shape_mag[:, pos] = 1
mag_shape[:, pos] = 1
else:
shape_mag[pos, :] = 1
return shape_mag
mag_shape[pos, :] = 1
return mag_shape
def alternating_filaments(dim, params):
def alternating_filaments(dim, spacing, x_or_y):
'''Get the magnetic shape of alternating filaments in x or y direction.
Arguments:
dim - the dimensions of the grid, shape(y, x)
......@@ -71,19 +66,18 @@ def alternating_filaments(dim, params):
the magnetic shape as a 2D-boolean-array.
'''
spacing, x_or_y = params
shape_mag = np.zeros(dim)
mag_shape = np.zeros(dim)
if x_or_y == 'y':
# TODO: List comprehension:
for i in range(0, dim[1], spacing):
shape_mag[:, i] = 1 - 2 * (int(i / spacing) % 2) # 1 or -1
mag_shape[:, i] = 1 - 2 * (int(i / spacing) % 2) # 1 or -1
else:
for i in range(0, dim[0], spacing):
shape_mag[i, :] = 1 - 2 * (int(i / spacing) % 2) # 1 or -1
return shape_mag
mag_shape[i, :] = 1 - 2 * (int(i / spacing) % 2) # 1 or -1
return mag_shape
def single_pixel(dim, params):
def single_pixel(dim, pixel):
'''Get the magnetic shape of a single magnetized pixel.
Arguments:
dim - the dimensions of the grid, shape(y, x)
......@@ -92,14 +86,12 @@ def single_pixel(dim, params):
the magnetic shape as a 2D-boolean-array.
'''
pixel = params
shape_mag = np.zeros(dim)
shape_mag[pixel[0], pixel[1]] = 1
return shape_mag
mag_shape = np.zeros(dim)
mag_shape[pixel[0], pixel[1]] = 1
return mag_shape
def create_hom_mag(dim, res, beta, shape_fun, params,
filename='output.txt', plot_mag_distr=False):
def create_hom_mag(dim, res, beta, mag_shape, filename='output.txt', plot_mag_distr=False):
'''Create homog. magnetization data, saved in a file with LLG-format.
Arguments:
dim - the dimensions of the grid, shape(y, x)
......@@ -107,23 +99,20 @@ def create_hom_mag(dim, res, beta, shape_fun, params,
(real space distance between two points)
beta - the angle of the magnetization
filename - the name of the file in which to save the data
shape_fun - the function to determine the shape of the magnetization,
several are defined in this module
params - a tuple of parameters used in shape_fun
mag_shape - an array of shape dim, representing the shape of the magnetic object,
a few are supplied in this module
Returns:
the the magnetic shape as a 2D-boolean-array.
the the magnetic distribution as a 2D-boolean-array.
'''
res *= 1.0E-9 / 1.0E-2 # from nm to cm
x = np.linspace(res / 2, dim[1] * res - res / 2, num=dim[1])
y = np.linspace(res / 2, dim[0] * res - res / 2, num=dim[0])
xx, yy = np.meshgrid(x, y)
shape_mag = shape_fun(dim, params)
xx, yy = np.meshgrid(x, y)
x_mag = np.array(np.ones(dim)) * np.cos(beta) * shape_mag
y_mag = np.array(np.ones(dim)) * np.sin(beta) * shape_mag
x_mag = np.array(np.ones(dim)) * np.cos(beta) * mag_shape
y_mag = np.array(np.ones(dim)) * np.sin(beta) * mag_shape
z_mag = np.array(np.zeros(dim))
if (plot_mag_distr):
......@@ -138,49 +127,7 @@ def create_hom_mag(dim, res, beta, shape_fun, params,
x_mag = np.reshape(x_mag,(-1))
y_mag = np.reshape(y_mag,(-1))
z_mag = np.array(np.zeros(dim[0] * dim[1]))
data = np.array([xx, yy, zz, x_mag, y_mag, z_mag]).T
with open(filename,'w') as mag_file:
mag_file.write('LLGFileCreator2D: %s\n' % filename.replace('.txt', ''))
mag_file.write(' %d %d %d\n' % (dim[1], dim[0], 1))
mag_file.writelines('\n'.join(' '.join('{:7.6e}'.format(cell)
for cell in row) for row in data) )
def create_logo(edge, res, beta = pi/2, filename='logo.txt', plot_mag_distr=False):
# TODO: rewrite so that every possible shape_mag can be given as argument
# TODO: write a script for creating, displaying and saving the logo
dim = (edge, edge)
res *= 1.0E-9 / 1.0E-2 # from nm to cm
x = np.linspace(res / 2, dim[1] * res - res / 2, num=dim[1])
y = np.linspace(res / 2, dim[0] * res - res / 2, num=dim[0])
xx, yy = np.meshgrid(x, y)
bottom = (yy >= 0.25*edge*res)
left = (yy <= 0.75/0.5 * xx)
right = np.fliplr(left)#(yy <= (edge-1)*res - 0.75/0.5 * (xx - 0.5*(edge-1)*res))
shape_mag = np.logical_and(np.logical_and(left, right), bottom)
x_mag = np.array(np.ones(dim)) * np.cos(beta) * shape_mag
y_mag = np.array(np.ones(dim)) * np.sin(beta) * shape_mag
z_mag = np.array(np.zeros(dim))
if (plot_mag_distr):
fig = plt.figure()
fig.add_subplot(111, aspect='equal')
plt.quiver(x_mag, y_mag, pivot='middle', angles='xy', scale_units='xy',
scale=1, headwidth=6, headlength=7)
xx = np.reshape(xx,(-1))
yy = np.reshape(yy,(-1))
zz = np.array(np.ones(dim[0] * dim[1]) * res / 2)
x_mag = np.reshape(x_mag,(-1))
y_mag = np.reshape(y_mag,(-1))
z_mag = np.array(np.zeros(dim[0] * dim[1]))
data = np.array([xx, yy, zz, x_mag, y_mag, z_mag]).T
with open(filename,'w') as mag_file:
mag_file.write('LLGFileCreator2D: %s\n' % filename.replace('.txt', ''))
......
src/pyramid/phasemap.py
+ 52
− 119
View file @ 78d70d7c
......@@ -14,7 +14,7 @@ Q_E = 1.602E-19 # TODO: unit
C = 2.998E8 # TODO: unit
def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
def fourier_space(mag_data, b_0=1, padding=0):
'''Calculate phasemap from magnetization data (fourier approach).
Arguments:
mag_data - MagDataLLG object (from magneticimaging.dataloader) storing
......@@ -31,8 +31,8 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
'''
res = mag_data.res
x_dim, y_dim, z_dim = mag_data.dim
x_mag, y_mag, z_mag = mag_data.magnitude
z_dim, y_dim, x_dim = mag_data.dim
z_mag, y_mag, x_mag = mag_data.magnitude
# TODO: interpolation (redefine xDim,yDim,zDim) thickness ramp
......@@ -40,6 +40,7 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
y_pad = y_dim/2 * padding
x_mag_big = np.zeros(((1 + padding) * y_dim, (1 + padding) * x_dim))
y_mag_big = np.zeros(((1 + padding) * y_dim, (1 + padding) * x_dim))
# TODO: padding so that x_dim and y_dim = 2^n
x_mag_big[y_pad : y_pad + y_dim, x_pad : x_pad + x_dim] = x_mag
y_mag_big[y_pad : y_pad + y_dim, x_pad : x_pad + x_dim] = y_mag
......@@ -50,7 +51,7 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
y = np.linspace( -nyq/2, nyq/2, x_mag_fft.shape[0]+1)[:-1]
xx, yy = np.meshgrid(x, y)
phase_fft = (1j * b_0) / (2 * PHI_0) * ((x_mag_fft * yy - y_mag_fft * xx)
phase_fft = (1j * res * b_0) / (2 * PHI_0) * ((x_mag_fft * yy - y_mag_fft * xx)
/ (xx ** 2 + yy ** 2 + 1e-18))
phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
......@@ -59,56 +60,45 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
# TODO: Electrostatic Component
return phase
def real_space_slab(mag_data, b_0=1, v_0=0, v_acc=30000,
pixel_map=None, jacobi=None):
def phi_pixel(method, n, m, res, b_0): # TODO: rename
if method == 'slab':
def F_h(n, m):
a = np.log(res**2 * (n**2 + m**2))
b = np.arctan(n / m)
return n*a - 2*n + 2*m*b
coeff = -b_0 * res**2 / ( 2 * PHI_0 )
return coeff * 0.5 * ( F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
-F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5) )
elif method == 'disc':
coeff = - b_0 * res**2 / ( 2 * PHI_0 )
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
def real_space(mag_data, method, b_0=1, jacobi=None):
'''Calculate phasemap from magnetization data (real space approach).
Arguments:
mag_data - MagDataLLG object (from magneticimaging.dataloader) storing
the filename, coordinates and magnetization in x, y and z
mag_data - MagDataLLG object (from magneticimaging.dataloader) storing the filename,
coordinates and magnetization in x, y and z
Returns:
the phasemap as a 2 dimensional array
'''
# TODO: combine with disc and sphere approach, pixel field as function
# TODO: Expand docstring!
'''
# TODO: Expand docstring!
# TODO: Delete
# import pdb; pdb.set_trace()
res = mag_data.res
x_dim, y_dim, z_dim = mag_data.dim
x_mag, y_mag, z_mag = mag_data.magnitude
z_dim, y_dim, x_dim = mag_data.dim
z_mag, y_mag, x_mag = mag_data.magnitude
# TODO: proper projection algorithm
x_mag, y_mag, z_mag = x_mag.mean(0), y_mag.mean(0), z_mag.mean(0)
beta = np.arctan2(y_mag, x_mag)
mag = np.hypot(x_mag, y_mag)
coeff = - b_0 * res**2 / ( 2 * PHI_0 ) / res # TODO: why / res
def F_h(n, m):
a = np.log(res**2 * (n**2 + m**2))
b = np.arctan(n / m)
return n*a - 2*n + 2*m*b
def phi_pixel(n, m): # TODO: rename
return coeff * 0.5 * ( F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
-F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5) )
mag = np.hypot(x_mag, y_mag)
def phi_mag(i, j): # TODO: rename
return (np.cos(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i]
-np.sin(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i])
def phi_mag_deriv(i, j): # TODO: rename
return -(np.sin(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i]
+ np.cos(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i])
'''CREATE COORDINATE GRIDS'''
x = np.linspace(0,(x_dim-1),num=x_dim)
y = np.linspace(0,(y_dim-1),num=y_dim)
......@@ -118,95 +108,38 @@ def real_space_slab(mag_data, b_0=1, v_0=0, v_acc=30000,
y_big = np.linspace(-(y_dim-1), y_dim-1, num=2*y_dim-1)
xx_big, yy_big = np.meshgrid(x_big, y_big)
phi_cos = phi_pixel(xx_big, yy_big)
phi_sin = phi_pixel(yy_big, xx_big)
if pixel_map is not None:
pixel_map[:] = phi_cos
phi_cos = phi_pixel(method, xx_big, yy_big, res, b_0)
phi_sin = phi_pixel(method, yy_big, xx_big, res, b_0)
def phi_mag(i, j): # TODO: rename
return (np.cos(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i]
-np.sin(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i])
def phi_mag_deriv(i, j): # TODO: rename
return -(np.sin(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i]
+np.cos(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i])
'''CALCULATE THE PHASE'''
phase = np.zeros((y_dim, x_dim))
# TODO: only iterate over pixels that have a magn. > threshold (first >0)
# import pdb; pdb.set_trace()
for j in range(y_dim):
for i in range(x_dim):
if (mag[j, i] != 0): # TODO: same result with or without?
if (mag[j, i] != 0 ):#or jacobi is not None): # TODO: same result with or without?
phi_mag_cache = phi_mag(i, j)
phase += mag[j,i] * phi_mag_cache
if jacobi is not None:
jacobi[i+x_dim*j] = phi_mag_cache.reshape(-1)
jacobi[x_dim*y_dim+i+x_dim*j] = (mag[j,i]*phi_mag_deriv(i, j)).reshape(-1)
jacobi[:,i+x_dim*j] = phi_mag_cache.reshape(-1)
jacobi[:,x_dim*y_dim+i+x_dim*j] = (mag[j,i]*phi_mag_deriv(i, j)).reshape(-1)
return phase
def real_space_disc(mag_data, b_0=1, v_0=0, v_acc=30000,
pixel_map=None, jacobi=None):
'''Calculate phasemap from magnetization data (real space approach).
Arguments:
mag_data - MagDataLLG object (from magneticimaging.dataloader) storing
the filename, coordinates and magnetization in x, y and z
Returns:
the phasemap as a 2 dimensional array
'''
# TODO: Expand docstring!
# TODO: Delete
# import pdb; pdb.set_trace()
res = mag_data.res
x_dim, y_dim, z_dim = mag_data.dim
x_mag, y_mag, z_mag = mag_data.magnitude
beta = np.arctan2(y_mag, x_mag)
mag = np.hypot(x_mag, y_mag)
coeff = - b_0 * res**2 / ( 2 * PHI_0 ) / res # TODO: why / res
def phi_pixel(n, m):
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
result = coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
return result
# r = np.hypot(n, m)
# r[y_dim-1,x_dim-1] = 1E-18 # Prevent div zero error at disc center
# return coeff / r**2 * m * (r != 0) # (r > R) = 0 for disc center
def phi_mag(i, j):
return mag[j,i]*(np.cos(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i]
-np.sin(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j,
x_dim-1-i:(2*x_dim-1)-i])
'''CREATE COORDINATE GRIDS'''
x = np.linspace(0,(x_dim-1),num=x_dim)
y = np.linspace(0,(y_dim-1),num=y_dim)
xx, yy = np.meshgrid(x,y)
x_big = np.linspace(-(x_dim-1), x_dim-1, num=2*x_dim-1)
y_big = np.linspace(-(y_dim-1), y_dim-1, num=2*y_dim-1)
xx_big, yy_big = np.meshgrid(x_big, y_big)
phi_cos = phi_pixel(xx_big, yy_big)
phi_sin = phi_pixel(yy_big, xx_big)
if pixel_map is not None:
pixel_map[:] = phi_cos
'''CALCULATE THE PHASE'''
phase = np.zeros((y_dim, x_dim))
# TODO: only iterate over pixels that have a magn. > threshold (first >0)
for j in range(y_dim):
for i in range(x_dim):
# if (mag[j, i] != 0): # TODO: same result with or without?
phase += phi_mag(i, j)
return phase
def phase_elec(mag_data, b_0=1, v_0=0, v_acc=30000):
......@@ -216,8 +149,8 @@ def phase_elec(mag_data, b_0=1, v_0=0, v_acc=30000):
# import pdb; pdb.set_trace()
res = mag_data.res
x_dim, y_dim, z_dim = mag_data.dim
x_mag, y_mag, z_mag = mag_data.magnitude
z_dim, y_dim, x_dim = mag_data.dim
z_mag, y_mag, x_mag = mag_data.magnitude
phase = np.logical_or(x_mag, y_mag, z_mag)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment