summaryrefslogtreecommitdiff
path: root/openEMS/python/openEMS/nf2ff.py
blob: b6d38d03114198277b4973994ae1e17c5a20d626 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program 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 3 of the License, or
# (at your option) any later version.
#
# This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
#

import os
import numpy as np
import h5py
from openEMS import _nf2ff
from openEMS import utilities

class nf2ff:
    """
    Create an nf2ff recording box. The nf2ff can either record in time-domain
    or frequency-domain. Further more certain directions and boundary condition
    mirroring can be enabled/disabled.

    :param name: str -- Name for this recording box.
    :param start/stop: (3,) array -- Box start/stop coordinates.
    :param directions: (6,) bool array -- Enable/Disables directions.
    :param mirror: (6,) int array -- 0 (Off), 1 (PEC) or 2 (PMC) boundary mirroring
    :param frequency: array like -- List of frequencies (FD-domain recording)
    """
    def __init__(self, CSX, name, start, stop, **kw):
        self.CSX   = CSX
        self.name  = name
        self.start = start
        self.stop  = stop

        self.freq  = None
        self.theta = None
        self.phi   = None
        self.center = None

        self.directions = [True]*6 # all directions by default
        if 'directions' in kw:
            self.directions = kw['directions']
            del kw['directions']
            assert len(self.directions)==6

        self.mirror = [0]*6
        if 'mirror' in kw:
            self.mirror = kw['mirror']
            del kw['mirror']
            assert len(self.mirror)==6

        self.dump_type = 0  # default Et/Ht
        self.dump_mode = 1  # default cell interpolated

        self.freq = None   # broadband recording by defualt
        if 'frequency' in kw:
            self.freq = kw['frequency']
            del kw['frequency']
            self.dump_type = 10  # Ef/Hf

            if np.isscalar(self.freq):
                self.freq = [self.freq]

        self.e_file = '{}_E'.format(self.name)
        self.h_file = '{}_H'.format(self.name)

        self.e_dump = CSX.AddDump(self.e_file, dump_type=self.dump_type  , dump_mode=self.dump_mode, file_type=1, **kw)
        self.h_dump = CSX.AddDump(self.h_file, dump_type=self.dump_type+1, dump_mode=self.dump_mode, file_type=1, **kw)
        if self.freq is not None:
            self.e_dump.SetFrequency(self.freq)
            self.h_dump.SetFrequency(self.freq)

#        print(self.directions)
        for ny in range(3):
            pos = 2*ny
            if self.directions[pos]:
                l_start = np.array(start)
                l_stop  = np.array(stop)
                l_stop[ny] = l_start[ny]
                self.e_dump.AddBox(l_start, l_stop)
                self.h_dump.AddBox(l_start, l_stop)
            if self.directions[pos+1]:
                l_start = np.array(start)
                l_stop  = np.array(stop)
                l_start[ny] = l_stop[ny]
                self.e_dump.AddBox(l_start, l_stop)
                self.h_dump.AddBox(l_start, l_stop)

    def CalcNF2FF(self, sim_path, freq, theta, phi, radius=1, center=[0,0,0], outfile=None, read_cached=False, verbose=0):
        """ CalcNF2FF(sim_path, freq, theta, phi, center=[0,0,0], outfile=None, read_cached=True, verbose=0):

        Calculate the far-field after the simulation is done.

        :param sim_path: str -- Simulation path
        :param freq: array like -- list of frequency for transformation
        :param theta/phi: array like -- Theta/Phi angles to calculate the far-field
        :param radius: float -- Radius to calculate the far-field (default is 1m)
        :param center: (3,) array -- phase center, must be inside the recording box
        :param outfile: str -- File to save results in. (defaults to recording name)
        :param read_cached: bool -- enable/disable read already existing results (default off)
        :param verbose: int -- set verbose level (default 0)

        :returns: nf2ff_results class instance
        """
        if np.isscalar(freq):
            freq = [freq]
        self.freq  = freq
        if np.isscalar(theta):
            theta = [theta]
        self.theta = theta
        if np.isscalar(phi):
            phi = [phi]
        self.phi   = phi
        self.center = center

        if outfile is None:
            fn = os.path.join(sim_path, self.name + '.h5')
        else:
            fn = os.path.join(sim_path,  outfile)
        if  not read_cached or not os.path.exists(fn):
            nfc = _nf2ff._nf2ff(self.freq, np.deg2rad(theta), np.deg2rad(phi), center, verbose=verbose)

            for ny in range(3):
                nfc.SetMirror(self.mirror[2*ny]  , ny, self.start[ny])
                nfc.SetMirror(self.mirror[2*ny+1], ny, self.stop[ny])

            nfc.SetRadius(radius)

            for n in range(6):
                fn_e = os.path.join(sim_path, self.e_file + '_{}.h5'.format(n))
                fn_h = os.path.join(sim_path, self.h_file + '_{}.h5'.format(n))
                if os.path.exists(fn_e) and os.path.exists(fn_h):
                    assert nfc.AnalyseFile(fn_e, fn_h)

            nfc.Write2HDF5(fn)

        result = nf2ff_results(fn)
        if result.phi is not None:
            assert np.abs((result.r-radius)/radius)<1e-6, 'Radius does not match. Did you read an invalid chached result? Try "read_cached=False"'
            assert utilities.Check_Array_Equal(np.rad2deg(result.theta), self.theta, 1e-4), 'Theta array does not match. Did you read an invalid chached result? Try "read_cached=False"'
            assert utilities.Check_Array_Equal(np.rad2deg(result.phi), self.phi, 1e-4), 'Phi array does not match. Did you read an invalid chached result? Try "read_cached=False"'
            assert utilities.Check_Array_Equal(result.freq, self.freq, 1e-6, relative=True), 'Frequency array does not match. Did you read an invalid chached result? Try "read_cached=False"'
        return result

class nf2ff_results:
    """
    nf2ff result class containing all results obtained by the nf2ff calculation.
    Usueally returned from nf2ff.CalcNF2FF

    Available attributes:

    * `fn`   : file name
    * `theta`: theta angles
    * `phi`  : phi angles
    * `r`    : radius
    * `freq` : frequencies
    * `Dmax` : directivity over frequency
    * `Prad` : total radiated power over frequency

    * `E_theta` : theta component of electric field over frequency/theta/phi
    * `E_phi`   : phi   component of electric field over frequency/theta/phi
    * `E_norm`  : abs   component of electric field over frequency/theta/phi
    * `E_cprh`  : theta component of electric field over frequency/theta/phi
    * `E_cplh`  : theta component of electric field over frequency/theta/phi
    * `P_rad`   : radiated power (S) over frequency/theta/phi
    """
    def __init__(self, fn):
        self.fn  = fn
        h5_file  = h5py.File(fn, 'r')
        mesh_grp = h5_file['Mesh']
        self.phi   = np.array(mesh_grp['phi'])
        self.theta = np.array(mesh_grp['theta'])
        self.r     = np.array(mesh_grp['r'])

        data  = h5_file['nf2ff']
        self.freq = np.array(data.attrs['Frequency'])

        self.Dmax = np.array(data.attrs['Dmax'])
        self.Prad = np.array(data.attrs['Prad'])

        THETA, PHI = np.meshgrid(self.theta, self.phi, indexing='ij')
        cos_phi = np.cos(PHI)
        sin_phi = np.sin(PHI)

        self.E_theta = []
        self.E_phi   = []
        self.P_rad   = []
        self.E_norm  = []
        self.E_cprh  = []
        self.E_cplh  = []
        for n in range(len(self.freq)):
            E_theta = np.array(h5_file['/nf2ff/E_theta/FD/f{}_real'.format(n)]) + 1j*np.array(h5_file['/nf2ff/E_theta/FD/f{}_imag'.format(n)])
            E_theta = np.swapaxes(E_theta, 0, 1)
            E_phi   = np.array(h5_file['/nf2ff/E_phi/FD/f{}_real'.format(n)])   + 1j*np.array(h5_file['/nf2ff/E_phi/FD/f{}_imag'.format(n)])
            E_phi   = np.swapaxes(E_phi, 0, 1)
            self.P_rad  .append(np.swapaxes(np.array(h5_file['/nf2ff/P_rad/FD/f{}'.format(n)]), 0, 1))

            self.E_theta.append(E_theta)
            self.E_phi  .append(E_phi)
            self.E_norm .append(np.sqrt(np.abs(E_theta)**2 + np.abs(E_phi)**2))
            self.E_cprh .append((cos_phi+1j*sin_phi) * (E_theta+1j*E_phi)/np.sqrt(2.0))
            self.E_cplh .append((cos_phi-1j*sin_phi) * (E_theta-1j*E_phi)/np.sqrt(2.0))