← Back to branch summary

~emmanuel-lambert/python-meep/intec

« back to all changes in this revision

Viewing changes to tests-intec/bragg_transmission.py

  • 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/bragg_transmission.py
 
1
#!/usr/bin/env python 
 
2
#passed
 
3
from meep import *
 
4
from math import fabs,sqrt, atan2, cos, sin, modf
 
5
import cmath 
 
6
 
 
7
nhi = 3.0
 
8
nlo = 1.0
 
9
wlo = nhi / (nlo + nhi)
 
10
Nperiods = 4
 
11
zsize = 10
 
12
 
 
13
class eps_nlo(Callback):
 
14
    def __init__(self):
 
15
        Callback.__init__(self)
 
16
        self.nlo=1
 
17
    def double_vec(self,v):
 
18
        return self.nlo*self.nlo
 
19
 
 
20
class eps_bragg(Callback):
 
21
    def __init__(self):
 
22
        Callback.__init__(self)
 
23
        self.nlo=1
 
24
 
 
25
    def double_vec(self,v):
 
26
        return self.eps(v)
 
27
 
 
28
    def eps(self,v):
 
29
        z = v.z() - zsize * 0.5
 
30
        if (fabs(z)*2 > Nperiods):
 
31
            return nlo*nlo
 
32
        else:
 
33
            #zi=0
 
34
            (zf,zi) = modf(z)
 
35
            if (zf < 0):
 
36
                zf += 1
 
37
            if (zf < wlo):
 
38
                return (nlo*nlo)
 
39
            else:
 
40
                return (nhi*nhi)
 
41
 
 
42
def byT12(m, n1, n2):
 
43
    n12 = n1 / n2
 
44
 
 
45
    td = 0.5 * (1 + n12)
 
46
    tod = 0.5 * (1 - n12)
 
47
 
 
48
    m00 = m[0][0]
 
49
    m01 = m[0][1]
 
50
    m10 = m[1][0]
 
51
    m11 = m[1][1]
 
52
 
 
53
    m[0][0] = m00 * td + m01 * tod
 
54
    m[0][1] = m00 * tod + m01 * td
 
55
    m[1][0] = m10 * td + m11 * tod
 
56
    m[1][1] = m10 * tod + m11 * td
 
57
 
 
58
 
 
59
def byP(m, n, w, dz):
 
60
    p = cmath.exp(complex(0, n * w * dz))
 
61
    pc = p.conjugate()
 
62
   
 
63
    m[0][0] *= p
 
64
    m[0][1] *= pc
 
65
    m[1][0] *= p
 
66
    m[1][1] *= pc
 
67
 
 
68
 
 
69
def abs2(x):
 
70
    ax = abs(x)
 
71
    return ax*ax
 
72
 
 
73
def bragg_transmission_analytic(freq_min, freq_max, nfreq, T, R):
 
74
    for i in range (0,nfreq):
 
75
        omega = 2*pi * (freq_min + i * (freq_max - freq_min) / (nfreq - 1))
 
76
        Tm = [ [ 1, 0 ], [ 0, 1 ] ]
 
77
        for j in range (0, Nperiods):
 
78
 
 
79
            byT12(Tm, nlo, nhi)
 
80
            byP(Tm, nhi, omega, 1 - wlo)
 
81
            byT12(Tm, nhi, nlo)
 
82
            byP(Tm, nlo, omega, wlo)
 
83
 
 
84
        refl = - Tm[1][0] / Tm[1][1];
 
85
        T[i] = abs2(Tm[0][0] + refl * Tm[0][1]);
 
86
        R[i] = abs2(refl)
 
87
 
 
88
  
 
89
 
 
90
def bragg_transmission(a, freq_min, freq_max, nfreq, T, R, use_hdf5):
 
91
    v = volone(zsize, a)
 
92
    # eps_bragg
 
93
    s=structure(v, MU, pml(0.5))
 
94
    f=fields(s)
 
95
    f.use_real_fields()
 
96
 
 
97
    #eps_nlo
 
98
    s0=structure(v, EPS, pml(0.5))
 
99
    f0=fields(s0)
 
100
    f0.use_real_fields()
 
101
 
 
102
    srcpt=vec(0.1)
 
103
    Tfluxpt=volume(vec(zsize - 0.1))
 
104
    Rfluxpt=volume(vec(0.1))
 
105
 
 
106
    src=gaussian_src_time((freq_min + freq_max) * 0.5, 0.5 / fabs(freq_max - freq_min), 0, 5 / fabs(freq_max - freq_min))
 
107
    f.add_point_source(Ex, src, srcpt)
 
108
    f0.add_point_source(Ex, src, srcpt)
 
109
 
 
110
    ft = f.add_dft_flux_plane(Tfluxpt,freq_min, freq_max, nfreq)
 
111
    fr = f.add_dft_flux_plane(Rfluxpt,freq_min, freq_max, nfreq)
 
112
    ft0 = f0.add_dft_flux_plane(Tfluxpt, freq_min, freq_max, nfreq)
 
113
    fr0 = f0.add_dft_flux_plane(Rfluxpt, freq_min, freq_max, nfreq)
 
114
 
 
115
    while (f0.time() < nfreq / fabs(freq_max - freq_min) / 2):
 
116
        f0.step()
 
117
 
 
118
    if (use_hdf5):
 
119
        fr0.save_hdf5(f, "flux", "reflection")
 
120
        fr.load_hdf5(f, "flux", "reflection")
 
121
        fr.scale_dfts(-1.0)
 
122
        
 
123
        ff = f.open_h5file("flux", h5file.READONLY)
 
124
        ff.remove();
 
125
        del ff
 
126
    else:
 
127
        #fr -=fr0
 
128
        # operator -= cannot be wrapped !
 
129
        # use subtract() function as shown below  
 
130
        fr.subtract(fr0)
 
131
    while (f.time() < nfreq / fabs(freq_max - freq_min) / 2):
 
132
        f.step()
 
133
 
 
134
    flux = get_flux(ft)
 
135
    flux0 = get_flux(ft0)
 
136
 
 
137
    for i in range(0,nfreq):
 
138
        T[i] = flux[i] / flux0[i]
 
139
        
 
140
    flux = get_flux(fr)
 
141
    for i in range (0,nfreq):
 
142
        R[i] = -flux[i] / flux0[i]
 
143
 
 
144
    
 
145
def max2(a, b):
 
146
    return a if (a > b) else b
 
147
 
 
148
def min2(a, b):
 
149
    return a if (a < b) else b
 
150
 
 
151
def max2a(a, b):
 
152
    return max2(abs(a), abs(b))
 
153
 
 
154
def sqr( x):
 
155
    return x*x
 
156
 
 
157
def distance_from_curve(n, dx, ys, x, y):
 
158
    d = infinity
 
159
    for i in range (1, n):
 
160
        theta = atan2(ys[i] - ys[i-1], dx)
 
161
        L = sqrt(sqr(dx) + sqr(ys[i]-ys[i-1]))
 
162
        x0 = x - (i-1) * dx
 
163
        y0 = y - ys[i-1]
 
164
        x0p = x0 * cos(theta) + y0 * sin(theta)
 
165
        y0p = y0 * cos(theta) - x0 * sin(theta)
 
166
 
 
167
        if (x0p < 0):
 
168
            d = min2(sqrt(sqr(x0) + sqr(y0)), d)
 
169
        elif (x0p > L):
 
170
            d = min2(sqrt(sqr(x-i*dx) + sqr(y-ys[i])), d)
 
171
        else:
 
172
            d = min2(abs(y0p), d)
 
173
 
 
174
    return d
 
175
 
 
176
 
 
177
def doit(use_hdf5):
 
178
    nfreq = 100
 
179
    freq_min = 0.1
 
180
    freq_max = 0.5
 
181
 
 
182
 
 
183
    T = d_array(nfreq)
 
184
    R = d_array(nfreq)
 
185
    bragg_transmission(40.0, freq_min, freq_max, nfreq, T, R, use_hdf5)
 
186
    T0 = d_array(nfreq)
 
187
    R0 = d_array(nfreq)
 
188
    bragg_transmission_analytic(freq_min, freq_max, nfreq, T0, R0)
 
189
 
 
190
    dfreq = (freq_max - freq_min) / (nfreq - 1)
 
191
 
 
192
    maxerrT = 0
 
193
    maxerrR = 0
 
194
 
 
195
    for i in range(0,nfreq):
 
196
        errT = distance_from_curve(nfreq, dfreq, T0, i * dfreq, T[i])
 
197
        errR = distance_from_curve(nfreq, dfreq, R0, i * dfreq, R[i])
 
198
        if (errT > maxerrT):
 
199
            maxerrT = errT
 
200
        if (errR > maxerrR):
 
201
            maxerrR = errR
 
202
        if (errT * sqr(freq_min / (freq_min + i*dfreq)) > 0.01):
 
203
            abort("large error %g at freq = %g: T = %g instead of %g\n", errT, freq_min + i*dfreq, T[i], T0[i])
 
204
        if (errR * sqr(freq_min / (freq_min + i*dfreq)) > 0.01):
 
205
            abort("large error %g at freq = %g: R = %g instead of %g\n", errR, freq_min + i*dfreq, R[i], R0[i])
 
206
 
 
207
  
 
208
    master_printf("Done (max. err in T = %e, in R = %e)\n", maxerrT, maxerrR)
 
209
 
 
210
 
 
211
ep=eps_nlo()
 
212
set_EPS_Callback(ep.__disown__())
 
213
 
 
214
mu=eps_bragg()
 
215
set_MU_Callback(mu.__disown__())
 
216
 
 
217
doit(True)
 
218
doit(False)
 
219