-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrasolver.py
162 lines (131 loc) · 5.22 KB
/
rasolver.py
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
from hornpattern import PyramidalHorn, get_default_pyramidal_horn, get_horn_input_power
from arrayinfo import RAInfo
from sources import PlaneWave, Source
from rautils import gsinc, ideal_ref_unit, dB, sph2car_mtx, make_v_mtx, farR, Fresnel
from ratasks import Gain2D, Gain3D, EfieldResult, NearFieldResult
import numpy as np
import matplotlib.pyplot as plt
import time
import concurrent.futures
import threading
"""
1. array info (cell size scale | horns)
2. get phase
3. get tetm
4. add task
5. calc
"""
class RASolver:
def __init__(self, rainfo, type='far'):
self.rainfo = rainfo
self.tasks = []
self.lock = threading.Lock()
self.type = type
def __erxy_fft(self, u, v):
px, py = self.rainfo.get_pxy()
Nx, Ny = self.rainfo.get_Nxy()
dx, dy = self.rainfo.get_dxy()
k0 = self.rainfo.get_k0()
K1 = np.exp(-1j*k0/2.0 * (u*(Nx-1)*dx + v*(Ny-1)*dy))
d_sumx, d_sumy = 0.0, 0.0
for (i, (Exmn, Eymn)) in list(enumerate(self.rainfo)):
n, m = int(i / Nx), int(i % Nx)
ejk = np.exp(1j*k0*(u*m*dx+v*n*dy))
d_sumx += (Exmn * ejk)
d_sumy += (Eymn * ejk)
A = K1 * px *py * gsinc(k0*u*px/2.0) * gsinc(k0*v*py/2.0)
return A*d_sumx, A*d_sumy
def __calc_one_point(self, r, t, p):
u = np.sin(t) * np.cos(p)
v = np.sin(t) * np.sin(p)
Erx, Ery = self.__erxy_fft(u, v)
k0 = self.rainfo.get_k0()
R = r
E_phi = -1j*k0*np.exp(-1j*k0*R)/(2*np.pi*R)*np.cos(t) * (Erx * np.sin(p) - Ery * np.cos(p))
E_theta = 1j*k0*np.exp(-1j*k0*R)/(2*np.pi*R) * (Erx * np.cos(p) + Ery * np.sin(p))
return EfieldResult(E_theta, E_phi, (r, t, p))
def __calc_fresnel_point(self, x, y, z):
Ex, Ey, Ez = 0j, 0j, 0j
px, py = self.rainfo.get_pxy()
Nx, Ny = self.rainfo.get_Nxy()
k0 = self.rainfo.get_k0()
integ_sum = 0j
for (i, (Exmn, Eymn)) in list(enumerate(self.rainfo)):
m, n = int(i / Nx), int(i % Nx)
A = m*px - (Nx-1)*px/2
B = n*py - (Ny-1)*py/2
t2 = (-px/2+x-A) * np.sqrt(k0/(np.pi*z))
t1 = (px/2+x-A) * np.sqrt(k0/(np.pi*z))
t2p = (-py/2+y-B) * np.sqrt(k0/(np.pi*z))
t1p = (py/2+y-B) * np.sqrt(k0/(np.pi*z))
Ix = Fresnel(t2) - Fresnel(t1)
Iy = Fresnel(t2p) - Fresnel(t1p)
integ_sum += Eymn * Ix * Iy
Ey = 1j/2.0 * np.exp(-1j*k0*z) * integ_sum
return NearFieldResult(Ex, Ey, Ez, (x, y, z))
# add far zone task
def append_task(self, tsk):
self.tasks.append(tsk)
# calc a task in a thread
def __calc_one_task__(self, tsk):
#with self.lock:
for (a, b, c) in tsk:
if self.type == 'far':
tsk.set_current_result(self.__calc_one_point(a, b, c))
elif self.type == 'fresnel':
tsk.set_current_result(self.__calc_fresnel_point(a, b, c)) # this is x y z actually
else:
raise ValueError('type error')
return tsk
# iterate on tasks, each task use multi-threading
def run(self):
for task in self.tasks:
field_task = task.assign_task()
start_time = time.clock()
for tsk in field_task:
ret = self.__calc_one_task__(tsk)
task.set_results(ret)
print(time.clock()-start_time, "sec")
# FixMe bug with concurrency
def run_concurrency(self):
for task in self.tasks:
field_task = task.assign_task()
start_time = time.clock()
#t_num = int(len(field_task))
t_num = 4
with concurrent.futures.ThreadPoolExecutor(max_workers=t_num) as e:
futs = [e.submit(self.__calc_one_task__, tsk) for tsk in field_task]
for fut in concurrent.futures.as_completed(futs):
#task.set_results(futs[fut])
task.set_results(fut.result())
print(time.clock()-start_time, "sec")
def test_multi():
freq = 5e9
cell_sz = 30. / 1000.
scale = 30
# pre calculated horn power integration at 5GHz
horn_integ = 22732.769823328235
abg = (np.deg2rad(180), np.deg2rad(180), np.deg2rad(0))
src = Source()
horn = get_default_pyramidal_horn(freq)
src.append(horn, abg, (0., 0., cell_sz*scale))
tp = [(np.deg2rad(20), np.deg2rad(0)), (np.deg2rad(20), np.deg2rad(90)),
(np.deg2rad(20), np.deg2rad(180)), (np.deg2rad(20), np.deg2rad(270))]
#tp = [(np.deg2rad(0), np.deg2rad(0))]
tpm = [(np.deg2rad(0), np.deg2rad(0), 1)]
arr = RAInfo(src, cell_sz, (scale, scale), ('pencil', (tp)), ideal_ref_unit)
#arr = RAInfo(src, cell_sz, (scale, scale), ('oam', (tpm, np.deg2rad(0))), ideal_ref_unit)
solver = RASolver(arr)
tsk1 = Gain2D(np.deg2rad(0), 999)
tsk2 = Gain2D(np.deg2rad(90), 600)
tsk3 = Gain3D(100, 100)
solver.append_task(tsk1)
#solver.append_task(tsk2)
#solver.append_task(tsk3)
#solver.run()
solver.run_concurrency()
tsk1.post_process(horn_integ, True)
#tsk2.post_process(horn_integ, True)
#tsk3.post_process(horn_integ, True)
if __name__ == '__main__':
test_multi()