← Back to branch summary

~emmanuel-lambert/python-meep/intec

« back to all changes in this revision

Viewing changes to tests-intec/.svn/text-base/known_results.py.svn-base

  • Committer: emmanuel.lambert at ugent
  • Date: 2009-11-23 10:23:17 UTC
  • Revision ID: emmanuel.lambert@intec.ugent.be-20091123102317-m00vrwyw0m6kpcmm
further merging with NS / replacing tests dir with default tests and creating new 'tests-intec' subdir

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests-intec/.svn/text-base/known_results.py.svn-base
 
1
#!/usr/bin/env python 
 
2
# passed
 
3
from meep import *
 
4
from math import fabs,sqrt
 
5
 
 
6
dpml = 1.0
 
7
 
 
8
class epsilon(Callback):
 
9
    def __init__(self):
 
10
        Callback.__init__(self)
 
11
    def double_vec(self,r):
 
12
        return self.eps(r)
 
13
 
 
14
    def eps(self,r):
 
15
        p = r
 
16
        while (p.x() < -0.5):
 
17
            p.set_direction(X, p.x() + 1.0)
 
18
        while (p.x() >  0.5):
 
19
            p.set_direction(X, p.x() - 1.0)
 
20
        while (p.y() < -0.5):
 
21
            p.set_direction(Y, p.y() + 1.0)
 
22
        while (p.y() >  0.5):
 
23
            p.set_direction(Y, p.y() - 1.0)
 
24
        if (p.x()*p.x() + p.y()*p.y() < 0.3):
 
25
            return 12.0
 
26
        return 1.0
 
27
 
 
28
 
 
29
class sigma(Callback):
 
30
    def __init__(self, sig):
 
31
        Callback.__init__(self)
 
32
        self.sig=sig
 
33
        
 
34
    def double_vec(self,r):
 
35
        return self.sig
 
36
   
 
37
 
 
38
rods=epsilon()
 
39
set_EPS_Callback(rods.__disown__())
 
40
 
 
41
 
 
42
def compare(a, b, n):
 
43
    if (fabs(a-b) > fabs(b)*1e-5): 
 
44
        master_printf("%s differs by\t%g out of\t%g\n", n, a-b, b)
 
45
        master_printf("This gives a fractional error of %g\n", fabs(a-b)/fabs(b))
 
46
        return False
 
47
    else: 
 
48
        master_printf("Completed %s \n", n)
 
49
        return True 
 
50
 
 
51
def using_pml_ez(v, eps):
 
52
    ttot = 30.0
 
53
    s= structure(v, eps, pml(dpml))
 
54
    f= fields (s)
 
55
    f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(), complex(0,-2*pi*0.2))
 
56
    while (f.time() < ttot):
 
57
        f.step()
 
58
    p=monitor_point()
 
59
    f.get_point(p, v.center())
 
60
    return p.get_component(Ez).real
 
61
    
 
62
 
 
63
def x_periodic_y_pml(v, eps):
 
64
    ttot = 30.0
 
65
    s= structure(v, eps, pml(dpml, Y))
 
66
    f= fields_complex(s)
 
67
    f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(), complex(0,-2*pi*0.2))
 
68
    f.use_bloch(X, 0.1)
 
69
    while (f.time() < ttot):
 
70
        f.step()
 
71
    p=monitor_point()
 
72
    f.get_point(p, v.center())
 
73
    return p.get_component(Ez).real
 
74
 
 
75
def x_periodic(v,eps):
 
76
    ttot = 30.0
 
77
    s= structure(v, eps)
 
78
    f= fields_complex(s)
 
79
    f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(), complex(0,-2*pi*0.2))
 
80
    f.use_bloch(X, 0.1)
 
81
    while (f.time() < ttot):
 
82
        f.step()
 
83
    p= monitor_point()
 
84
    f.get_point(p, v.center())
 
85
    return p.get_component(Ez).real
 
86
 
 
87
def periodic_ez(v, eps):
 
88
    ttot = 30.0
 
89
    s= structure(v, eps)
 
90
    f= fields_complex(s)
 
91
    f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(), complex(0,-2*pi*0.2))
 
92
  
 
93
    if v.dim == D1:
 
94
        k = vec(0.3)
 
95
    elif v.dim ==  D2:
 
96
        k = vec(0.3,0.4)
 
97
    elif v.dim == D3:
 
98
        k = vec(0.3,0.5,0.8)
 
99
    else:
 
100
        k = veccyl(0.3,0.2)
 
101
  
 
102
    f.use_bloch(k)
 
103
    while (f.time() < ttot):
 
104
        f.step()
 
105
    p=monitor_point()
 
106
    f.get_point(p, v.center())
 
107
    return p.get_component(Ez).real
 
108
 
 
109
def metallic_ez(v,eps):
 
110
    ttot = 10.0
 
111
    s = structure(v, eps)
 
112
    f = fields_complex(s)
 
113
    f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(), complex(0,-2*pi*0.2))
 
114
    while (f.time() < ttot):
 
115
        f.step()
 
116
    p=monitor_point()
 
117
    f.get_point(p, v.center())
 
118
    return p.get_component(Ez).real
 
119
 
 
120
 
 
121
 
 
122
def polariton_ex(v, eps):
 
123
    ttot = 10.0
 
124
    s= structure(v, eps)
 
125
 
 
126
    pol1 = sigma(7.63)
 
127
    set_DBL1_Callback(pol1.__disown__())
 
128
    
 
129
    s.add_polarizability(DBL1, 0.3, 0.1)
 
130
    
 
131
    f=fields_complex(s)
 
132
    f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center(), complex(0,-2*pi*0.2))
 
133
    while (f.time() < ttot):
 
134
        f.step()
 
135
    p=monitor_point()
 
136
    f.get_point(p, v.center())
 
137
    del_DBL1_Callback()
 
138
    return p.get_component(Ex).real
 
139
 
 
140
def polariton_energy(v, eps):
 
141
    ttot = 10.0
 
142
    s= structure(v, eps)
 
143
    
 
144
    pol1 = sigma(7.63)
 
145
    set_DBL1_Callback(pol1.__disown__())
 
146
    s.add_polarizability(DBL1, 0.3, 0.1)
 
147
 
 
148
    #s.add_polarizability(one, 0.3, 0.1, 7.63)
 
149
    f=fields_complex(s, 0, 1)
 
150
    f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center(), complex(0,-2*pi*0.2))
 
151
    while (f.time() < ttot):
 
152
        f.step()
 
153
 
 
154
    del_DBL1_Callback()
 
155
    return f.total_energy()
 
156
 
 
157
def saturated_polariton_ex(v, eps):
 
158
    ttot = 10.0
 
159
    s=structure(v, eps)
 
160
 
 
161
    pol1 = sigma(7.63)
 
162
    set_DBL1_Callback(pol1.__disown__())
 
163
    s.add_polarizability(DBL1, 0.3, 0.1)
 
164
    
 
165
#    thep = s.add_polarizability(one, 0.3, 0.1, -0.063, 0.1)
 
166
    f= fields (s)
 
167
    f.use_real_fields()
 
168
    f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center(),  1/sqrt(4*pi) * complex(0,-2*pi*0.2))
 
169
    while (f.time() < ttot):
 
170
        f.step()
 
171
    p=monitor_point()
 
172
    f.get_point(p, v.center())
 
173
    del_DBL1_Callback()
 
174
 
 
175
    return p.get_component(Ex).real * sqrt(4*pi)
 
176
 
 
177
master_printf("Testing with some known results...\n")
 
178
a = 10.0
 
179
compare(-0.0894851, polariton_ex(volone(1.0, a), ONE), "1D polariton")
 
180
compare(0.32617294, polariton_energy(volone(1.0, a), ONE), "1D polariton energy")
 
181
compare(5.20605, metallic_ez(voltwo(1.0, 1.0, a), ONE), "1x1 metallic 2D TM")
 
182
compare(0.883776, using_pml_ez(voltwo(1.0+2*dpml, 1.0+2*dpml, a), ONE), "1x1 PML 2D TM")
 
183
compare(0.110425, x_periodic(voltwo(1.0, 1.0, a), ONE), "1x1 X periodic 2D TM")
 
184
compare(-4.78767, periodic_ez(voltwo(1.0, 3.0, a), EPS),"1x1 fully periodic 2D TM rods")
 
185
compare(1.12502, periodic_ez(voltwo(1.0, 3.0, a), ONE), "1x1 fully periodic 2D TM")
 
186
compare(0.608815, x_periodic_y_pml(voltwo(1.0, 1.0+2*dpml, a), ONE), "1x1 X periodic Y PML 2D TM")
 
187
compare(-41.8057, metallic_ez(vol3d(1.0, 1.0, 1.0, a), ONE), "1x1x1 metallic 3D")
 
188
compare(-100.758, x_periodic(vol3d(1.0, 1.0, 1.0, a), ONE), "1x1x1 X periodic 3D")
 
189
compare(-101.398, x_periodic_y_pml(vol3d(1.0, 1.0+2*dpml, 1.0, a), ONE), "1x1x1 X periodic Y PML 3D")
 
190
compare(-103.844, periodic_ez(vol3d(1.0, 1.0, 1.0, a), EPS), "1x1x1 fully periodic 3D rods")
 
191
compare(-99.1618, periodic_ez(vol3d(1.0, 1.0, 1.0, a), ONE), "1x1x1 fully periodic 3D")
 
192