Skip to content
Snippets Groups Projects
Commit 2115d08a authored by Bas van der Tol's avatar Bas van der Tol
Browse files

Cleanup comparison-oskar demo scripts

parent 4baab9f0
No related branches found
No related tags found
No related merge requests found
Pipeline #2375 passed
Stage: build-base
Stage: format-and-test
Stage: build-everybeam
Stage: build-doc
Show whitespace changes
Inline Side-by-side
demo/comparison-oskar/generate_basefunction_plots.py
+ 10
− 59
View file @ 2115d08a
...@@ -33,47 +33,25 @@ for em_idx in range(2): ...@@ -33,47 +33,25 @@ for em_idx in range(2):
run_oskar.main() run_oskar.main()
A = read_oskar_beams() A = read_oskar_beams()
N = A.shape[0]
print(N)
v = np.zeros((N,N,4))
for i in range(N):
x = 2.0*i/(N-1) - 1.0
for j in range(N):
y = 2.0*j/(N-1) - 1.0
theta = np.arcsin(np.sqrt(x*x + y*y));
phi = np.arctan2(y,x);
e_theta = np.array([np.cos(phi), np.sin(phi)])
e_phi = np.array([-np.sin(phi), np.cos(phi)])
T = np.stack((e_theta, e_phi))
v[i,j,0] = np.abs(np.dot(A[i,j,:,:], T))[0,0]
v[i,j,1] = np.abs(np.dot(A[i,j,:,:], T))[0,1]
v[i,j,2] = np.angle(np.dot(A[i,j,:,:],T))[0,0]
v[i,j,3] = np.angle(np.dot(A[i,j,:,:],T))[0,1]
plt.subplot(2,4,1) plt.subplot(2,4,1)
plt.imshow(v[:,:,0].T, clim=(0,.25), origin='lower') plt.imshow(np.abs(A[:,:,0,0]).T, clim=(0,.25), origin='lower')
#plt.imshow(v[:,:,0].T, origin='lower')
plt.colorbar() plt.colorbar()
plt.title('abs(Etheta)') plt.title('abs(Etheta)')
plt.ylabel('OSKAR') plt.ylabel('OSKAR')
plt.subplot(2,4,2) plt.subplot(2,4,2)
plt.imshow(v[:,:,1].T, clim=(0,.25), origin='lower') plt.imshow(np.abs(A[:,:,0,1]).T, clim=(0,.25), origin='lower')
plt.colorbar() plt.colorbar()
plt.title('abs(Ephi)') plt.title('abs(Ephi)')
plt.subplot(2,4,3) plt.subplot(2,4,3)
plt.imshow(v[:,:,2].T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower') plt.imshow(np.angle(A[:,:,0,0]).T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
plt.colorbar() plt.colorbar()
plt.title('angle(Etheta)') plt.title('angle(Etheta)')
plt.subplot(2,4,4) plt.subplot(2,4,4)
plt.imshow(v[:,:,3].T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower') plt.imshow(np.angle(A[:,:,0,1]).T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
plt.colorbar() plt.colorbar()
plt.title('angle(Ephi)') plt.title('angle(Ephi)')
...@@ -96,53 +74,26 @@ for em_idx in range(2): ...@@ -96,53 +74,26 @@ for em_idx in range(2):
if apply_transpose: if apply_transpose:
A *= (-1)**(m+1) A *= (-1)**(m+1)
N = A.shape[0]
print(N)
v = np.zeros((N,N,4))
for i in range(N):
x = 2.0*i/(N-1) - 1.0
for j in range(N):
y = 2.0*j/(N-1) - 1.0
theta = np.arcsin(np.sqrt(x*x + y*y));
phi = np.arctan2(y,x);
e_theta = np.array([np.cos(phi), np.sin(phi)])
e_phi = np.array([-np.sin(phi), np.cos(phi)])
T = np.stack((e_theta, e_phi))
T = np.eye(2)
v[i,j,0] = np.abs(np.dot(A[i,j,:,:], T))[0,0]
v[i,j,1] = np.abs(np.dot(A[i,j,:,:], T))[0,1]
v[i,j,2] = np.angle(np.dot(A[i,j,:,:],T))[0,0]
v[i,j,3] = np.angle(np.dot(A[i,j,:,:],T))[0,1]
plt.subplot(2,4,5) plt.subplot(2,4,5)
plt.imshow(v[:,:,0].T, clim=(0,.25), origin='lower') plt.imshow(np.abs(A[:,:,0,0]).T, clim=(0,.25), origin='lower')
plt.colorbar() plt.colorbar()
plt.title('abs(Etheta)') plt.title('abs(Etheta)')
plt.ylabel('EveryBeam') plt.ylabel('EveryBeam')
plt.subplot(2,4,6) plt.subplot(2,4,6)
plt.imshow(v[:,:,1].T, clim=(0,.25), origin='lower') plt.imshow(np.abs(A[:,:,0,1]).T, clim=(0,.25), origin='lower')
plt.colorbar() plt.colorbar()
plt.title('abs(Ephi)') plt.title('abs(Ephi)')
plt.subplot(2,4,7) plt.subplot(2,4,7)
plt.imshow(v[:,:,2].T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower') plt.imshow(np.angle(A[:,:,0,0]).T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
plt.colorbar() plt.colorbar()
plt.title('angle(Ex)') plt.title('angle(Etheta)')
plt.subplot(2,4,8) plt.subplot(2,4,8)
plt.imshow(v[:,:,3].T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower') plt.imshow(np.angle(A[:,:,0,1]).T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
plt.colorbar() plt.colorbar()
plt.title('angle(Ey)') plt.title('angle(Ephi)')
#with open('telescope.tm/element_pattern_spherical_wave_x_te_re_0_50.txt') as f:
#coeff_line = f.readline()
#plt.gcf().suptitle('telescope.tm/element_pattern_spherical_wave_x_te_re_0_50.txt\n{}'.format(coeff_line))
plt.gcf().suptitle('l = {}, m = {}, s = {}'.format(l,m,s)) plt.gcf().suptitle('l = {}, m = {}, s = {}'.format(l,m,s))
......
demo/comparison-oskar/read_oskar_beams.py
+ 28
− 32
View file @ 2115d08a
...@@ -15,21 +15,10 @@ def read_oskar_beams(): ...@@ -15,21 +15,10 @@ def read_oskar_beams():
N = d.shape[-1] N = d.shape[-1]
A = np.zeros((N,N,2,2), dtype=np.complex128) A = np.zeros((N,N,2,2), dtype=np.complex128)
A[:,:,i,j] = d A[:,:,i,j] = d
return A
if __name__ == "__main__":
A = read_oskar_beams()
N = A.shape[0] N = A.shape[0]
print(N) print(N)
v = np.zeros((N,N,4))
for i in range(N): for i in range(N):
x = 2.0*i/(N-1) - 1.0 x = 2.0*i/(N-1) - 1.0
for j in range(N): for j in range(N):
...@@ -37,40 +26,47 @@ if __name__ == "__main__": ...@@ -37,40 +26,47 @@ if __name__ == "__main__":
theta = np.arcsin(np.sqrt(x*x + y*y)); theta = np.arcsin(np.sqrt(x*x + y*y));
phi = np.arctan2(y,x); phi = np.arctan2(y,x);
e_theta = np.array([np.cos(phi), np.sin(phi)]) # Need to swap the sign of phi to get the transformation right
e_phi = np.array([-np.sin(phi), np.cos(phi)]) # Still need to figure out why
phi = -phi
e_theta = np.array([[np.cos(phi)], [np.sin(phi)]])
e_phi = np.array([[-np.sin(phi)], [np.cos(phi)]])
T = np.stack((e_theta, e_phi)).T # This matrix applies the transformation defined in
v[i,j,0] = np.abs(np.dot(A[i,j,:,:], T))[0,0] # https://github.com/OxfordSKA/OSKAR/blob/master/oskar/convert/define_convert_ludwig3_to_theta_phi_components.h
v[i,j,1] = np.abs(np.dot(A[i,j,:,:], T))[0,1] T = np.concatenate((e_theta, e_phi), axis=1)
v[i,j,2] = np.mod(np.angle(A[i,j,0,0]),2*np.pi)
v[i,j,3] = np.angle(A[i,j,0,1]) # Transformation is applied to the right side of A
# so the rows a of A are tranformed
A[i,j,:,:] = np.dot(A[i,j,:,:], T)
return A
if __name__ == "__main__":
A = read_oskar_beams()
plt.subplot(2,2,1) plt.subplot(2,2,1)
plt.imshow(v[:,:,0], clim=(0,.25)) plt.imshow(np.abs(A[:,:,0,0]).T, clim=(0,.25), origin='lower')
plt.colorbar() plt.colorbar()
plt.title('abs(Etheta)') plt.title('abs(Etheta)')
plt.subplot(2,2,2) plt.subplot(2,2,2)
plt.imshow(v[:,:,1], clim=(0,.25)) plt.imshow(np.abs(A[:,:,0,1]).T, clim=(0,.25), origin='lower')
plt.colorbar() plt.colorbar()
plt.title('abs(Ephi)') plt.title('abs(Ephi)')
plt.subplot(2,2,3) plt.subplot(2,2,3)
plt.imshow(v[:,:,2]) plt.imshow(np.angle(A[:,:,0,0]).T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
plt.colorbar() plt.colorbar()
plt.title('angle(Ex)') plt.title('angle(Etheta)')
plt.subplot(2,2,4) plt.subplot(2,2,4)
plt.imshow(v[:,:,3]) plt.imshow(np.angle(A[:,:,0,1]).T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
plt.colorbar() plt.colorbar()
plt.title('angle(Ey)') plt.title('angle(Ephi)')
plt.show() plt.show()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment