magpar_scripts.py

Go to the documentation of this file.

#################################
#
# magpar standard Python library
#
#################################

# This file is part of magpar.

# Copyright (C) 2009 Werner Scholz

# www:   http://www.magpar.net/
# email: magpar(at)magpar.net

# magpar is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.

# magpar is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with magpar; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

import math, random, unittest, cProfile

# Vector class inspired by
# http://www.gamedev.net/community/forums/topic.asp?topic_id=486122

# see also:
#http://code.activestate.com/recipes/52272/

class Vector(object):
    """
    """
    def __init__(self, p = [0, 0, 0]):
        self.x = float(p[0])
        self.y = float(p[1])
        self.z = float(p[2])

    def __getitem__(self, key):
        if(key == 0):
            return self.x
        elif(key == 1):
            return self.y
        elif(key == 2):
            return self.z
        else:
            raise Exception("Invalid key to Vector")

    def __setitem__(self, key, value):
        if(key == 0):
            self.x = value
        elif(key == 1):
            self.y = value
        elif(key == 2):
            self.z = value
        else:
            raise Exception("Invalid key to Vector")

    def copy(self):
        return Vector([self.x,self.y,self.z])

    def __neg__(self):
        return Vector([-self.x,-self.y,-self.z])

    def __abs__(self):
        return math.sqrt(self.x*self.x+self.y*self.y+self.z*self.z)

    def __pow__(self, val):
        return abs(self)**val

    def __add__(self, val):
        return Vector([self.x + val.x, self.y + val.y, self.z + val.z])

    def __iadd__(self, val):
        self.x += val.x
        self.y += val.y
        self.z += val.z
        return self

    def __sub__(self,val):
        return Vector([self.x - val.x, self.y - val.y, self.z - val.z])

    def __isub__(self, val):
        self.x -= val.x
        self.y -= val.y
        self.z -= val.z
        return self

    def __mul__(self, val):
        if type(val) == Vector:
          return self.x * val.x + self.y * val.y + self.z * val.z
        else:
          return Vector([self.x * val, self.y * val, self.z * val])

    def __rmul__(self, val):
        if type(val) == Vector:
          return self.x * val.x + self.y * val.y + self.z * val.z
        else:
          return Vector([self.x * val, self.y * val, self.z * val])

    def __imul__(self, val):
        self.x *= val
        self.y *= val
        self.z *= val
        return self

    def __div__(self, val):
        return Vector([self.x / val, self.y / val, self.z / val])

    def __idiv__(self, val):
        self.x /= val
        self.y /= val
        self.z /= val
        return self

    def outer(self,b):
        return Vector([self.y*b.z-self.z*b.y,
                       self.z*b.x-self.x*b.z,
                       self.x*b.y-self.y*b.x])

    def normalize(self):
        if abs(self)<1e-9:
            self.x = 0.0
            self.y = 0.0
            self.z = 0.0
        else:
            l = abs(self)
            self.x /= l
            self.y /= l
            self.z /= l
        return self

    def __eq__(self,other):
        if (self.x==other.x) & (self.y==other.y) & (self.z==other.z):
            return True
        else:
            return False

    def __ne__(self,other):
        return not(self==other)

    def __repr__(self):
        return "[" + str(self.x) + "," + str(self.y) + "," + str(self.z) + "]"

    def __str__(self):
        return "(" + str(self.x) + "," + str(self.y) + "," + str(self.z) + ")"

class TestVectorClass(unittest.TestCase):
    def setUp(self):
        self.x=Vector([1,0,0])
        self.y=Vector([0,1,0])
        self.z=Vector([0,0,1])
        self.a=Vector([3,0,0])
        self.b=Vector([0,4,0])
        self.c=Vector([3,4,0])

    def testNeg(self):
        t = self.c.copy()
        t = -t
        self.assertEqual(t,Vector([-3,-4,0]))

    def testAbs(self):
        t=abs(self.a)
        self.assertEqual(t,3.0)
        t=abs(self.c)
        self.assertEqual(t,5.0)

    def testPow(self):
        t=self.a**2.0
        self.assertEqual(t,abs(self.a)**2.0)
        t=abs(self.c)
        self.assertEqual(t,5.0)

    def testSimpleArithmetics(self):
        t = self.x+self.y
        self.assertEqual(t,Vector([1,1,0]))
        t -= self.y
        self.assertEqual(t,self.x)
        t = self.x - self.y
        self.assertEqual(t,Vector([1,-1,0]))
        t += self.y
        self.assertEqual(t,self.x)
        t = self.x * 3.0
        self.assertEqual(t,self.a)
        t /= 3.0
        self.assertEqual(t,self.x)
        t = self.a / 3.0
        self.assertEqual(t,self.x)
        t *= 3.0
        self.assertEqual(t,self.a)
        t = self.x * self.a
        self.assertEqual(t,3.0)

    def testOuter(self):
        t=self.x.outer(self.y)
        self.assertEqual(t,self.z)

    def testNormalize(self):
        t=self.a.normalize()
        self.assertEqual(t,self.x)

def initialize_magnetization_mx(r):
    """
    """
    return Vector([1.0, 0.0, 0.0])

def initialize_magnetization_my(r):
    """
    """
    return Vector([0.0, 1.0, 0.0])

def initialize_magnetization_mz(r):
    """
    """
    return Vector([0.0, 0.0, 1.0])

def initialize_magnetization_111(r):
    """
    """
    a=math.sqrt(1.0/3.0)
    return Vector([a, a, a])

def initialize_magnetization_theta_rad(r, theta = 0):
    """Set magnetization in x-z plane to desired angle theta (from z-axis)
    """

    return Vector([math.sin(theta), 0.0, math.cos(theta)])

def initialize_magnetization_theta_deg(r, theta_deg = 0):
    m = initialize_magnetization_theta_rad(r,math.radians(theta_deg))

    return m

def initialize_magnetization_phi_rad(r, phi = 0):
    """Set magnetization in x-y plane to desired angle
    """

    return Vector([math.cos(theta), math.sin(theta), 0.0])

def initialize_magnetization_phi_deg(r,phi_deg = 0):
    m = initialize_magnetization_phi_rad(r,math.radians(phi_deg))

    return m

def initialize_magnetization_random(r):
    """Random magnetization - inaccurate: not uniform on unit sphere!
    """

    mx = random.uniform(-1, 1)
    my = random.uniform(-1, 1)
    mz = random.uniform(-1, 1)

    return Vector([mx, my, mz])

def initialize_magnetization_flower(r, c = Vector([0.0,0.0,0.0]), f = 1.0):
    """Approximate flower state (around z-axis)
    """

    m = Vector([r[0]-c[0],r[1]-c[1],1.0])
    if (r[2]<0):
        m *= -1.0

    return m

def vortex(p, R):
    """Approximate vortex state

    vortex centered at (cx, cy) with radius R:
    if R>0, then vortex is counter-clockwise in xy plane.
    if R<0, then vortex is clockwise in xy plane.
    center of vortex points in pol*z direction.
    """

    r=p.copy()
    r[2]=0
    absr=abs(r)

    # inside vortex core
    if absr <= math.fabs(R):
        A = 2.0*R*absr/(R*R+absr*absr)
        # avoid division by zero
        if (absr > 1e-8):
            m = Vector([A*r[1]/absr, -A*r[0]/absr, pol*math.sqrt(1.0-A*A)])
        else:
            m = Vector([0.0, 0.0, pol])
    # outside vortex core
    else:
        m = Vector([r[1]/absr, -r[0]/absr, 0.0])

    if R<0:
        m[0]=-m[0]
        m[1]=-m[1]

    return m

def initialize_magnetization_vortex(r, c = Vector([0.0,0.0,0.0])):
    """Approximate vortex state
    """

    # vortex radius
    R = 0.1

    # polarity (+1.0 or -1.0)
    pol = 1.0

    m = vortex(r-c, R)
    m[2] *= pol

    return Vector(m)

def initialize_magnetization_blochwall(r, xpos = 0.0, width = 0.1):
    """Bloch wall (parallel to y-z plane)
    """

    if (r[0] < xpos-width/2.0):
        m = Vector([0.0, 0.0, -1.0])
    elif (r[0] > xpos+width/2.0):
        m = Vector([0.0, 0.0, 1.0])
    else:
        m = Vector([0.0,
                    math.cos((r[0]-xpos)/(width/2.0)*(math.pi/2.0)),
                    math.sin((r[0]-xpos)/(width/2.0)*(math.pi/2.0))])

    return m

def initialize_magnetization_head_to_head_transverse(x, y, z, xpos = 0.0, width = 0.1):
    """Transverse head-to-head wall
    """

    if (r[0] < xpos-width/2.0):
        m = Vector([1.0, 0.0, 0.0])
    elif (r[0] > xpos+width/2.0):
        m = Vector([-1.0, 0.0, 0.0])
    else:
        m = Vector([-math.sin((x-xpos)/(width/2.0)*(math.pi/2.0)),
                    math.cos((x-xpos)/(width/2.0)*(math.pi/2.0)),
                    0.0])

    return m

def initialize_magnetization_head_to_head_vortex(r, xpos = 0.0, width = 0.1):
    """Approximate vortex head-to-head wall

    Head-to-head vortex wall at (0, 0) with radius R
    if R>0, then vortex is clockwise in xy plane.
    if R<0, then vortex is counter-clockwise in xy plane.
    center of vortex points in -z direction.
    """

    if (x < xpos-width/2.0):
        m = Vector([1.0, 0.0, 0.0])
    elif (x > xpos+width/2.0):
        m = Vector([-1.0, 0.0, 0.0])
    else:
        m = vortex(r-c, R)
        m[2] *= pol

    return m

def external_field_homogeneous(r, theta_deg = 0.0, phi_deg = 0.0, amp_T = 0.0):
    """
    """

    theta = math.radians(theta_deg)
    phi   = math.radians(phi_deg)

    h = Vector([amp_T*math.sin(theta)*math.cos(phi),
                amp_T*math.sin(theta)*math.sin(phi),
                amp_T*math.cos(theta)])

    return h

def external_field_time_scaling(t_ns):
    """
    """

    if (t_ns<0):
        amp_T=0
    else:
        amp_T=t_ns

    return amp_T

class Line:
    def __init__(self,a,b):
      self.start_point = a
      self.end_point = b
      self.vector = b-a
      self.length = abs(b-a)

    def __str__(self):
        return "start_point: " + str(self.start_point) + "\n" + \
               "end_point:   " + str(self.end_point) + "\n"  + \
               "vector:      " + str(self.vector) + "\n"  + \
               "length:      " + str(self.length) + "\n"

class TestLineClass(unittest.TestCase):
    def setUp(self):
        self.a=Vector([3,0,0])
        self.b=Vector([0,4,0])
        self.line=Line(self.a,self.b)

    def testvector(self):
        self.assertEqual(self.line.vector,Vector([-3,4,0]))

    def testlength(self):
        self.assertEqual(self.line.length,5)

class Segment:
    def __init__(self,s,a,b):
      self.start = s
      self.line = Line(a,b)

    def __str__(self):
        return "start: " + str(self.start) + "\n" + str(self.line)

class Polygon:
    def __init__(self):
        self.segments = []
        self.length = 0.0

    def append_segment(self,a,b):
        self.segments.append(Segment(self.length,a,b))
        self.length += abs(b-a)

    def point(self,u):
        for s in self.segments:
            if (u>=s.start/self.length) & (u<(s.start+s.line.length)/self.length):
                p = s.line.start_point + (u-s.start/self.length)*s.line.vector/s.line.length*self.length
                return p

    def vec(self,u,length):
        for s in self.segments:
            if (u>=s.start/self.length) & (u<(s.start+s.line.length)/self.length):
                v = Vector(s.line.vector)
                v.normalize()
                return v*length

    def __str__(self):
        string = "Polygon of length " + str(self.length) + \
                 " with " + str(len(self.segments)) + " segments\n"
        for s in self.segments:
            string += str(s)
        return string

def test_straight_poly():
    poly_points=[]
    poly_points.append(Vector([-100.0,0.0,0.0]))
    poly_points.append(Vector([ -50.0,0.0,0.0]))
    poly_points.append(Vector([ -10.0,0.0,0.0]))
    poly_points.append(Vector([   0.0,0.0,0.0]))
    poly_points.append(Vector([   5.0,0.0,0.0]))
    poly_points.append(Vector([  30.0,0.0,0.0]))
    poly_points.append(Vector([ 100.0,0.0,0.0]))
    polygon = Polygon()
    for s in range(len(poly_points)-1):
        polygon.append_segment(poly_points[s],poly_points[s+1])
    return polygon

class TestPolygonClass(unittest.TestCase):
    def setUp(self):
        self.polygon = test_straight_poly()

    def testpoints(self):
        self.assertEqual(self.polygon.point(0),self.polygon.segments[0].line.start_point)
        self.assertEqual(self.polygon.length,200)
        self.assertEqual(self.polygon.point(0.5),Vector([0,0,0]))
        self.assertEqual(self.polygon.point(0.99),Vector([98,0,0]))
        self.assertEqual(self.polygon.vec(0,1),Vector([1,0,0]))

def biot_savart(current_A, r, ds):
    """Calculate magnetic field according to Biot-Savart law

    dB=mu0*I/(4 pi)* (ds x r)/r^3

    http://en.wikipedia.org/wiki/Biot-Savart_law
    http://maxwell.byu.edu/~spencerr/websumm122/node70.html
    http://hyperphysics.phy-astr.gsu.edu/HBASE/magnetic/magcur.html#c2
    """

#    if abs(r) < abs(ds)*5.0:
#        print "Warning: distance r ", abs(r), r, " too short compared to segment ds ", abs(ds), ds

    if abs(r)<1e-9:
        # catch 1/r divergence
        h = Vector([0,0,0])
    else:
        h = 1e-7*current_A*ds.outer(r)/abs(r)**3.0
    return h

class TestBiotSavart(unittest.TestCase):
    def setUp(self):
        self.polygon = test_straight_poly()

    def testfield(self):
        h=biot_savart(1.0, Vector([0,0,1]), Vector([0.01,0,0]))
        self.assertAlmostEqual(abs(h)*1e9,1,5)
        h.normalize()
        self.assertEqual(h,Vector([0,-1,0]))

def field_of_wire_slow(wire,current_A,r):
    npoints=10000
    h = Vector()
    for u in range(0,npoints):
        #print wire.point(u*1.0/npoints), wire.vec(u*1.0/npoints,wire.length/npoints), r
        h += biot_savart(current_A, \
                    wire.point(u*1.0/npoints)-r, \
                    wire.vec(u*1.0/npoints,wire.length/npoints))
    return h

def field_of_wire(wire,current_A,r):
    npoints=1000
    h = Vector()
    for s in wire.segments:
        for u in range(0,npoints):
            h += biot_savart(current_A, \
                        s.line.start_point+s.line.vector*u/npoints-r, \
                        s.line.vector/npoints)
    return h

class Test_external_field_wire(unittest.TestCase):
    def setUp(self):
        self.polygon = test_straight_poly()

    # http://hyperphysics.phy-astr.gsu.edu/HBASE/magnetic/magcur.html#c2
    def testfield(self):
        h = field_of_wire_slow(self.polygon,1.0,Vector([0,0,0.1]))

        self.assertAlmostEqual(abs(h)*1e7,20,4)
        h.normalize()
        self.assertEqual(h,Vector([0,1,0]))

def initmag(x, y, z):
    """Initialize magnetization with desired function called in here.

    This function must be define after all the initialize_magnetization_* functions!
    """

    r = Vector([x,y,z])
    m = initialize_magnetization_mx(r)

    m /= abs(m)

    return m[0],m[1],m[2]

def hext_xyz2(x, y, z):
    """Set external field with desired function called in here.

    This function must be define after all the external_field_* functions!
    """

    h_T = external_field_homogeneous(x, y, z)

    return h_T[0], h_T[1], h_T[2]

def hext_xyzt(x, y, z, t_ns):
    """Set external field with desired function called in here.

    This function must be define after all the external_field_* functions!
    """

    h_T = external_field_homogeneous(x, y, z)
    f = external_field_time_scaling(t_ns)

    return h_T[0]*f, h_T[1]*f, h_T[2]*f

# Finally, import user defined functions.
# These can shadow/overload the standard functions above!

#from magpar_user_scripts import *

# from: http://www.python.org/doc/essays/list2str.html
import time
    def timing(f, n, a):
        print f.__name__,
        r = range(n)
        t1 = time.clock()
        for i in r:
            f(a); f(a); f(a); f(a); f(a); f(a); f(a); f(a); f(a); f(a)
        t2 = time.clock()
        print round(t2-t1, 3)

def test_function_2col():
    astart = -3.0
    astop = 3.0
    nstep = 100

    for i in range(nstep+1):
        a = astart+i*(astop-astart)/nstep

        # v1, v2, v3 = initialize_magnetization_flower(a, 0, 0)
        v1 = hext_xyzt(0,0,0,a)

        print a, v1

def misc3():
    import matplotlib.pyplot as plt
    plt.plotfile("x",cols=(0,1),delimiter=" ")

def wire_z():
    poly_points=[]
    poly_points.append(Vector([0.0,0.0, -10.0]))
    poly_points.append(Vector([0.0,0.0,  10.0]))
    polygon = Polygon()
    for s in range(len(poly_points)-1):
        polygon.append_segment(poly_points[s],poly_points[s+1])
    return polygon

def hext_f(t_ns):
    """
    """

    tstart = 0.0
    fstart = 0.0
    tend   = 50.0
    fend   = 200.0

    if t_ns < t_start:
        f=fstart
    elif t_ns < fend:
        f=fstart+(fend-fstart)*(t_ns-tstart)/(tend-tstart)

    return f

coils = wire_z()

def hext_xyz(x, y, z):
    """Set external field with desired function called in here.

    This function must be defined after all the external_field_* functions!
    """

    h_T = field_of_wire(coils,0.001,Vector([x,y,z]))

    return h_T[0], h_T[1], h_T[2]

# http://docs.python.org/tutorial/modules.html#executing-modules-as-scripts
if __name__ == "__main__":
    #test_function_2col()

    #print hext_xyz(1.0,0,0)
    cProfile.run('print hext_xyz(1.0,0,0)')

    unittest.main()