-
Notifications
You must be signed in to change notification settings - Fork 1
/
main.py
208 lines (173 loc) · 8.71 KB
/
main.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
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
## -*- coding: utf-8 -*- heloooooooooooooooooo
import rayTracing as rt_line
import radPat as rp
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import input as I
import pandas as pd
import rayTubes as rtube
import warnings
warnings.filterwarnings("ignore")
import time
start = time.time()
############### Input import ##################
h2 = I.h2
p = I.p
n2 = I.n_diec #dielectric refractive index
n1 = I.n1 #air refractive indeix
nML = I.nML
nML2 = I.nML2 #matching layer dielectric
N = I.N #number of rays
L = I.L
Array = I.Array #x-values indicating where the array is
D = I.D #x-space
s1 = I.s1
s2 = I.s2
s0 = I.s0
matchingLayer1 = I.matchingLayer1
matchingLayer2 = I.matchingLayer2
surface1_arr = I.surface1
MLayer1_arr = I.MLayer1
surface2_arr = I.surface2
MLayer2_arr = I.MLayer2
aperture_plane = I.aperture_plane
################ end input import###########
# df = pd.DataFrame(surface1_arr*1e3, p*1e3)
# df.to_excel('surface1.xlsx', sheet_name='Sheet1')
# df = pd.DataFrame(surface2_arr*1e3, p*1e3)
# df.to_excel('surface2.xlsx', sheet_name='Sheet1')
# df = pd.DataFrame(MLayer1_arr*1e3, p*1e3)
# df.to_excel('MLayer1_arr.xlsx', sheet_name='Sheet1')
# df = pd.DataFrame(MLayer2_arr*1e3, p*1e3)
# df.to_excel('MLayer2_arr.xlsx', sheet_name='Sheet1')
############ Create Segments ####################
segments =[]
#function calling: discretize_function(f, n_segments, n1, n2, isLast, isFirst):
segments = np.append(segments, rt_line.discretize_function(s0, 1, 1,1, False, True))
segments = np.append(segments, rt_line.discretize_function(aperture_plane, 1, 1, 1, True, False))
if I.matchingLayers: #if there are matching layers
segments = np.append(segments, rt_line.discretize_function(matchingLayer1, 70, n1, nML, False, False))
segments = np.append(segments, rt_line.discretize_function(s1, 70, nML, n2, False, False))
segments = np.append(segments, rt_line.discretize_function(s2, 70, n2, nML2, False, False))
segments = np.append(segments, rt_line.discretize_function(matchingLayer2, 70, nML2, n1, False, False))
else:
segments = np.append(segments, rt_line.discretize_function(s1, 7, n1, n2, False, False))
segments = np.append(segments, rt_line.discretize_function(s2, 7, n2, n1, False, False))
################# end create segments #######################
# if 0:
# #if we want to import an aritrary shape from a file
# surface1 = np.loadtxt('surface1.csv', delimiter=',')
# surface2 = np.loadtxt('surface2.csv', delimiter=',')
################################# REVERSE ##########################################3
angle_in = []
angle_position = []
angle_in, angle_position = rt_line.reverseRayTracing_segments(I.output_angle, segments)
f = interp1d(angle_position, angle_in, kind='cubic') #interpolate to find function f that fits angle position and angle in
angles_for_direct = f(Array) #Find angles corresponding to defined points on x-axis
#angles_for_direct = np.deg2rad(np.ones(N)*I.output_angle)
plot = 1
################################# DIRECT ##########################################3
#Create plot for direct raytracing
if plot:
fig = plt.figure()
fig.set_dpi(700)
ax = fig.add_subplot(111)
csfont = {'fontname':'Times New Roman', 'fontsize':'10'}
font = {'family' : 'Times New Roman',
'weight' : 'normal',
'size' : 12}
plt.rc('font', **font)
plt.title('Direct RT')
plt.ylim([-0,1])
plt.xlim([-1.5, 1.5])
plt.ylabel('z (m)')
plt.xlabel('x (m)')
plt.plot(p, surface1_arr, color='gray', linewidth = 1)
plt.plot(p, surface2_arr, color='gray', linewidth = 1)
ax.fill_between(p, surface2_arr, surface1_arr, color = '#c3c5e2')
# if I.nSurfaces == 4:
# plt.plot(p, MLayer1_arr, color = 'chocolate', linewidth = 0.1)
# plt.plot(p, MLayer2_arr, color = 'chocolate', linewidth = 0.1)
# ax.fill_between(p, MLayer1_arr, surface1_arr, color = 'orange')
# ax.fill_between(p, MLayer2_arr, surface2_arr, color = 'orange')
ax.set_aspect(1, adjustable='box')
for i in range(0, len(segments)):
plt.plot([segments[i].A[0], segments[i].B[0]], [segments[i].A[1], segments[i].B[1]], color = 'red', linewidth = 0.5)
Ak_ap = [] #E field amplitude on the aperture
phi_a = np.zeros(N) #phase distribution on the array
dck = [] #infinitesimal arc length on the aperture
rays = [] #the ray objects
x_ap = np.zeros(N) #x value on aperture
y_ap = np.zeros(N) #y value on aperture
rays = rt_line.directRayTracing_segments(angles_for_direct, segments)
for i in range(0, len(rays)):
Pk_np = rays[i].Pk
for j in range(0, len(Pk_np)-3):
if j % 2 == 0: #to get each point, Pk_np = [x1 y1 x2 y2...]
if plot: plt.plot([Pk_np[j], Pk_np[j+2]], [Pk_np[j+1], Pk_np[j+3]], color='black', linewidth = 0.5)
x_ap[i] = Pk_np[j]
y_ap[i] = Pk_np[j+1]
################################# END DIRECT, END GO ##########################################3
############################## RAY TUBE THEORY ##########################################3
path_length = np.zeros(N, dtype=np.complex_) #length that ray has travelled
nk = np.zeros([N,2]) #normal of the aperture
sk = np.zeros([N,2]) #pointying vector
ts_cascade = np.ones(N, dtype=np.complex_) #reflection coefficient normal to plane of incident
ts_aggr = np.ones(N, dtype=np.complex_)
#all these functions can be optimized
path_length = rtube.getPathLength(rays, segments) #length that ray have travelled, including material parameters
nk, sk = rtube.getLastNormal(rays) #get normal and poynting vector at aperture
ts_cascade, ts_aggr = rtube.getTransmissionCoef(rays, segments) #transmission coefficients
Ak_ap, dck = rtube.getAmplitude(rays, segments, angles_for_direct) #electric field amplitude over aperture and infinitesimal arc length over aperture
############################## RADIATION PATTERN - KIRCHOFF ##########################################3
# Etotal, theta, Ap_field_casc, dck_casc = rp.getRadiationPattern(Ak_ap, path_length[1:N-1], nk[1:N-1], sk[1:N-1], dck, x_ap[1:N-1], y_ap[1:N-1], ts_cascade[1:N-1]) #does not include the outer rays
# Etotal_dB = 20*np.log10(abs(Etotal))
# print('Cascade: '+ str(max(Etotal_dB)))
# #plot the radiation pattern
# fig2 = plt.figure()
# fig2.set_dpi(400)
# ax2 = fig2.add_subplot(111)
# ax2.set_aspect(1.5, adjustable='box')
# plt.plot(-theta*180/np.pi+90, 20*np.log10(abs(Etotal)/max(abs(Etotal))), linewidth=1, color = 'blue')
# plt.ylabel('Normalized Pattern, dB')
# plt.title('Cascade')
# plt.xlim([-70, 70])
# plt.ylim([-35, 0])
# plt.xlabel('$\u03B8 $, degrees')
# plt.xticks(range(-90, 91, 10))
# plt.yticks(range(-35, 10, 5))
# plt.rcParams["font.family"] = "Times New Roman"
# ax2.xaxis.label.set_fontsize(10)
# ax2.yaxis.label.set_fontsize(10)
# plt.grid()
# plt.show()
# # ##saving the radiation pattern results in an excel
# df = pd.DataFrame(Etotal_dB, theta)
# df.to_excel('RT_radpat_' + str(I.output_angle) + 'deg.xlsx', sheet_name='Sheet1')
Etotal, theta, Ap_field_aggr, dck_agg = rp.getRadiationPattern(Ak_ap, path_length[1:N-1], nk[1:N-1], sk[1:N-1], dck, x_ap[1:N-1], y_ap[1:N-1], ts_aggr[1:N-1]) #does not include the outer rays
Etotal_dB = 20*np.log10(abs(Etotal))
print('Aggregate: '+ str(max(Etotal_dB)))
#plot the radiation pattern
fig2 = plt.figure()
fig2.set_dpi(400)
ax2 = fig2.add_subplot(111)
ax2.set_aspect(1.5, adjustable='box')
plt.plot(-theta*180/np.pi+90, 20*np.log10(abs(Etotal)/max(abs(Etotal))), linewidth=1, color = 'red')
plt.ylabel('Normalized Pattern, dB')
plt.title('Aggregate')
plt.xlim([-70, 70])
plt.ylim([-35, 0])
plt.xlabel('$\u03B8 $, degrees')
plt.xticks(range(-90, 91, 10))
plt.yticks(range(-35, 10, 5))
plt.rcParams["font.family"] = "Times New Roman"
ax2.xaxis.label.set_fontsize(10)
ax2.yaxis.label.set_fontsize(10)
plt.grid()
plt.show()
# #saving the radiation pattern results in an excel
df = pd.DataFrame(Etotal_dB, theta)
df.to_excel('RT_radpat_aggr' + str(I.output_angle) + 'deg.xlsx', sheet_name='Sheet1')
end = time.time()
print("Execution time : ", (end-start))