diff --git a/demo/comparison-oskar/generate_basefunction_plots.py b/demo/comparison-oskar/generate_basefunction_plots.py index 8b2b42abd3b6adfc6ec1b4541cc5b71c587efcf8..8da6a0167fd52a296a921f9e23733b7e8c3fe22b 100755 --- a/demo/comparison-oskar/generate_basefunction_plots.py +++ b/demo/comparison-oskar/generate_basefunction_plots.py @@ -33,47 +33,25 @@ for em_idx in range(2): run_oskar.main() 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.imshow(v[:,:,0].T, clim=(0,.25), origin='lower') - #plt.imshow(v[:,:,0].T, origin='lower') + plt.imshow(np.abs(A[:,:,0,0]).T, clim=(0,.25), origin='lower') plt.colorbar() plt.title('abs(Etheta)') plt.ylabel('OSKAR') 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.title('abs(Ephi)') 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.title('angle(Etheta)') 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.title('angle(Ephi)') @@ -96,53 +74,26 @@ for em_idx in range(2): if apply_transpose: 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.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.title('abs(Etheta)') plt.ylabel('EveryBeam') 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.title('abs(Ephi)') 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.title('angle(Ex)') + plt.title('angle(Etheta)') 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.title('angle(Ey)') - - #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.title('angle(Ephi)') plt.gcf().suptitle('l = {}, m = {}, s = {}'.format(l,m,s)) diff --git a/demo/comparison-oskar/read_oskar_beams.py b/demo/comparison-oskar/read_oskar_beams.py index d57b3c52b87b46771319d6ac216cf119ae0d7518..c6b53de3c05ab8c992dfcfa7a67f4a1a0abc9b57 100644 --- a/demo/comparison-oskar/read_oskar_beams.py +++ b/demo/comparison-oskar/read_oskar_beams.py @@ -15,21 +15,10 @@ def read_oskar_beams(): N = d.shape[-1] A = np.zeros((N,N,2,2), dtype=np.complex128) A[:,:,i,j] = d - - return A - - - -if __name__ == "__main__": - - - A = read_oskar_beams() - N = A.shape[0] + 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): @@ -37,40 +26,47 @@ if __name__ == "__main__": 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)]) + # Need to swap the sign of phi to get the transformation right + # 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 - 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.mod(np.angle(A[i,j,0,0]),2*np.pi) - v[i,j,3] = np.angle(A[i,j,0,1]) + # This matrix applies the transformation defined in + # https://github.com/OxfordSKA/OSKAR/blob/master/oskar/convert/define_convert_ludwig3_to_theta_phi_components.h + T = np.concatenate((e_theta, e_phi), axis=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.imshow(v[:,:,0], clim=(0,.25)) + plt.imshow(np.abs(A[:,:,0,0]).T, clim=(0,.25), origin='lower') plt.colorbar() plt.title('abs(Etheta)') 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.title('abs(Ephi)') 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.title('angle(Ex)') + plt.title('angle(Etheta)') 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.title('angle(Ey)') + plt.title('angle(Ephi)') plt.show() - - - - - - -