← Back to branch summary

~emmanuel-lambert/python-meep/intec

« back to all changes in this revision

Viewing changes to tests-intec/.svn/text-base/harmonics.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/harmonics.py.svn-base
 
1
#!/usr/bin/env python 
 
2
# passed
 
3
 
 
4
from meep import *
 
5
from math import fabs,sqrt
 
6
 
 
7
class epsilon(Callback):
 
8
    def __init__(self, the_value = 1.0):
 
9
        Callback.__init__(self)
 
10
        self.eps = the_value 
 
11
    def double_vec(self,r):
 
12
        return self.eps 
 
13
 
 
14
rods=epsilon()
 
15
set_EPS_Callback(rods.__disown__())
 
16
 
 
17
def harmonics(freq, chi2, chi3, J):
 
18
    dpml = 5.0
 
19
    res = 20
 
20
    sz = 100+2*dpml
 
21
    v = vol1d(sz, res)
 
22
    v.center_origin()
 
23
 
 
24
    rods.eps= 1.0
 
25
    s=structure(v, EPS, pml(dpml))
 
26
    rods.eps = chi2
 
27
    s.set_chi2(EPS)
 
28
    rods.eps = chi3
 
29
    s.set_chi3(EPS)
 
30
 
 
31
    f=fields(s)
 
32
    f.use_real_fields()
 
33
  
 
34
    src=gaussian_src_time(freq, freq / 20)
 
35
    f.add_point_source(Ex, src, vec(-0.5 * sz + dpml), J)
 
36
 
 
37
    fpt=vec(0.5 * sz - dpml - 0.5)
 
38
    d1 = f.add_dft_flux(Z, volume(fpt), freq, freq, 1)
 
39
    d2 = f.add_dft_flux(Z, volume(fpt), 2*freq, 2*freq, 1)
 
40
    d3 = f.add_dft_flux(Z, volume(fpt), 3*freq, 3*freq, 1)
 
41
 
 
42
    emax = 0
 
43
 
 
44
    while (f.time() < f.last_source_time()):
 
45
        emax = max(emax, abs(f.get_field(Ex, fpt)))
 
46
        f.step()
 
47
 
 
48
    while True  :
 
49
        emaxcur = 0;
 
50
        T = f.time() + 50
 
51
        while (f.time() < T):
 
52
            e = abs(f.get_field(Ex, fpt))
 
53
            emax = max(emax, e)
 
54
            emaxcur = max(emaxcur, e)
 
55
            f.step()
 
56
    
 
57
        if (emaxcur < 1e-6 * emax):
 
58
            break
 
59
  
 
60
    #dft_flux.flux() returns a pointer  !
 
61
    # to obtain an array of fluxes use get_flux()
 
62
    d1f = pdouble_value(d1.flux())
 
63
    d2f = pdouble_value(d2.flux())
 
64
    d3f = pdouble_value(d3.flux())
 
65
 
 
66
    A2 = d2f / d1f
 
67
    A3 = d3f / d1f
 
68
    master_printf("harmonics(%g,%g,%g) = %g, %g\n", chi2, chi3, J, A2, A3)
 
69
    return ( A2, A3 )
 
70
 
 
71
def different(a, a0, thresh, msg):
 
72
    if (fabs(a - a0) > thresh * fabs(a0)):
 
73
        master_printf("error: %s\n --- %g vs. %g (%g error > %g)\n",msg, a, a0, fabs(a - a0)/fabs(a0), thresh)
 
74
        return 1
 
75
    return 0
 
76
 
 
77
freq = 1.0 / 3.0
 
78
a2=0
 
79
a3=0
 
80
a2_2=0
 
81
a3_2=0
 
82
thresh = 1e-4
 
83
 
 
84
a2,a3 = harmonics(freq, 0.27e-4, thresh, 1.0)
 
85
if (different(a2, 9.80298e-07, thresh, "2nd harmonic mismatches known value")):
 
86
    abort('error 1')    
 
87
if (different(a3, 9.97759e-07, thresh, "3rd harmonic mismatches known value")):
 
88
    abort('error 2')    
 
89
 
 
90
a2_2,a3_2 = harmonics(freq, 0.54e-4, 2e-4, 1.0)
 
91
master_printf("doubling chi2, chi3 = %g x 2nd harmonic, %g x 3rd\n", a2_2 / a2, a3_2 / a3)
 
92
if (different(a2_2 / a2, 4.0, 0.01, "incorrect chi2 scaling")):
 
93
   abort('error 3')    
 
94
if (different(a3_2 / a3, 4.0, 0.01, "incorrect chi3 scaling")):
 
95
   abort('error 4')    
 
96
 
 
97
a2_2, a3_2 = harmonics(freq, 0.27e-4, 1e-4, 2.0)
 
98
master_printf("doubling J = %g x 2nd harmonic, %g x 3rd\n", a2_2 / a2, a3_2 / a3)
 
99
if (different(a2_2 / a2, 4.0, 0.01, "incorrect J scaling for 2nd harm.")):
 
100
   abort('error 5')    
 
101
if (different(a3_2 / a3, 16.0, 0.01, "incorrect J scaling for 3rd harm.")):
 
102
   abort('error 6')    
 
103
 
 
104
a2_2, a3_2 = harmonics(freq, 0.27e-4, 0.0, 1.0)
 
105
if (different(a2, a2_2, 1e-2, "chi3 has too big effect on 2nd harmonic")):
 
106
   abort('error 7')    
 
107
if (a3_2 / a3 > 1e-4):
 
108
   master_printf("error: too much 3rd harmonic without chi3\n");
 
109
   abort('error 8')    
 
110
 
 
111
a2_2, a3_2 = harmonics(freq, 0.0, 1e-4, 1.0)
 
112
if (different(a3, a3_2, 1e-3, "chi2 has too big effect on 3rd harmonic")):
 
113
   abort('error 9')    
 
114
if (a2_2 / a2 > 1e-5):
 
115
   master_printf("error: too much 2nd harmonic without chi3\n");
 
116
   abort('error 10')
 
117