diff --git a/analysator/calculations/__init__.py b/analysator/calculations/__init__.py index 3257cef4..1aea6eaa 100644 --- a/analysator/calculations/__init__.py +++ b/analysator/calculations/__init__.py @@ -37,7 +37,8 @@ ''' #for usage with "from (package) import *" -__all__=["ids3d","intpol_file","intpol_points","cutthrough","fourier","spectra","variable","timeevolution","pitchangle","gyrophaseangle","themis_observation","cut3d","lineout","fit","fieldtracer","non_maxwellianity","null_lines","interpolator_amr"] +__all__=["ids3d","intpol_file","intpol_points","cutthrough","fourier","spectra","variable","timeevolution","pitchangle","gyrophaseangle","cut3d","lineout","fit","fieldtracer","non_maxwellianity","null_lines","interpolator_amr" + "spacecraft_to_simulation_frame","simulation_to_spacecraft_frame","simulation_to_observation_frame"] # List of functions and classes that should be imported into the interface @@ -53,10 +54,6 @@ from .pitchangle import pitch_angles #from .backstream import extract_velocity_cells_sphere, extract_velocity_cells_non_sphere from .gyrophaseangle import gyrophase_angles_from_file -from .themis_observation import themis_observation_from_file -from .themis_observation import themis_plot_detector -from .themis_observation import themis_plot_phasespace_contour -from .themis_observation import themis_plot_phasespace_helistyle #from .themis_observation import simulation_to_spacecraft_frame from .cut3d import cut3d from .lineout import lineout @@ -66,3 +63,4 @@ from .non_maxwellianity import epsilon_M from .null_lines import LMN_null_lines_FOTE from .interpolator_amr import AMRInterpolator, supported_amr_interpolators +from .virtual_observations import spacecraft_to_simulation_frame,simulation_to_spacecraft_frame,simulation_to_observation_frame diff --git a/analysator/calculations/virtual_observations.py b/analysator/calculations/virtual_observations.py new file mode 100644 index 00000000..9660881a --- /dev/null +++ b/analysator/calculations/virtual_observations.py @@ -0,0 +1,40 @@ + +import numpy as np +def simulation_to_spacecraft_frame(spinvector, detector_axis, phi=0): + ''' Builds a matrix to transform coordinates from simulation frame into spaceraft frame + :param spinvector Spacecraft spin axis, in simulation coordinates + :param detector_axis Detector plane normal axis + :param phi Rotate spacecraft around spin axis after setting up coordinate system + ''' + # TODO: Normalize vectors? + y = np.cross(detector_axis,spinvector) + z = np.cross(spinvector, y) + yr = np.cos(phi)*y - np.sin(phi)*z + zr = np.sin(phi)*y + np.cos(phi)*z + m = np.array([spinvector, yr, zr]) + + return m + +def spacecraft_to_simulation_frame(spinvector, detector_axis, phi=0): + ''' Builds a matrix to transform coordinates from spaceraft frame back to simulation frame + :param spinvector Spacecraft spin axis, in simulation coordinates + :param detector_axis Detector plane normal axis + :param phi Rotate spacecraft around spin axis after setting up coordinate system + ''' + return simulation_to_spacecraft_frame(spinvector,detector_axis,phi).T + +def simulation_to_observation_frame(x_axis,y_axis): + ''' Builds a 3x3 matrix to transform velocities into an observation plane + :param x_axis: x-axis of the observation plane (in simulation coordinates) + :param y_axis: y-axis of the observation plane (gets orthonormalized) + ''' + xn = np.linalg.norm(x_axis) + x_axis /= xn + p = x_axis.dot(y_axis) + y_axis -= p*x_axis + yn = np.linalg.norm(y_axis) + y_axis /= yn + z_axis = np.cross(x_axis,y_axis) + return np.array([x_axis,y_axis,z_axis]) + + diff --git a/analysator/plot/__init__.py b/analysator/plot/__init__.py index 0e71900d..e591a6e8 100644 --- a/analysator/plot/__init__.py +++ b/analysator/plot/__init__.py @@ -36,7 +36,10 @@ #for usage with "from (package) import *" __all__ = ['plot_colormap', 'plot_vdf', 'plot_vdfdiff', 'plot_vdf_profiles', 'plot_colormap3dslice', 'plot_threeslice', 'plot_ionosphere', 'plot_isosurface', 'plot_neutral_sheet', - 'plot_variables','colormaps', 'plot_helpers'] + 'plot_variables','colormaps', 'plot_helpers' + 'themis_plot_phasespace_contour','themis_plot_detector','themis_observation_from_file' + 'themis_plot_phasespace_helistyle',"themis_observation" + ] from .plot_variables import plot_variables, plot_multiple_variables @@ -46,7 +49,7 @@ import matplotlib.pyplot as plt import matplotlib from . import colormaps - +from .plot_themis_observation import themis_plot_detector,themis_plot_phasespace_contour,themis_plot_phasespace_helistyle,themis_observation,themis_observation_from_file from . import plot_helpers from .plot_colormap import plot_colormap from .plot_vdf import plot_vdf diff --git a/analysator/calculations/themis_observation.py b/analysator/plot/plot_themis_observation.py similarity index 84% rename from analysator/calculations/themis_observation.py rename to analysator/plot/plot_themis_observation.py index f93bb4e5..10875ef7 100644 --- a/analysator/calculations/themis_observation.py +++ b/analysator/plot/plot_themis_observation.py @@ -25,11 +25,10 @@ import analysator as pt import matplotlib.pyplot as plt import matplotlib -from .rotation import rotateVectorToVector from scipy.interpolate import griddata from scipy.signal import sepfir2d -from packaging.version import Version import logging +from analysator.calculations import virtual_observations as vsc # Detector data obtained from the Themis ESA instrument paper # http://dx.doi.org/10.1007/s11214-008-9440-2 @@ -48,64 +47,6 @@ themis_colors=[(0,0,0),(.5,0,.5),(0,0,1),(0,1,1),(0,1,0),(1,1,0),(1,0,0)] themis_colormap = matplotlib.colors.LinearSegmentedColormap.from_list("themis",themis_colors) -def get_dv(vlsvReader,pop="proton"): - # Get velocity grid sizes: - vel_mesh_size = vlsvReader.get_velocity_mesh_size(pop=pop) - vel_block_size = vlsvReader.get_velocity_block_size(pop=pop) - vxcells = vel_mesh_size[0]*vel_block_size[0] - vycells = vel_mesh_size[1]*vel_block_size[1] - vzcells = vel_mesh_size[2]*vel_block_size[2] - - vel_mesh_limits = vlsvReader.get_velocity_mesh_extent(pop=pop) - vxmin = vel_mesh_limits[0] - vymin = vel_mesh_limits[1] - vzmin = vel_mesh_limits[2] - vxmax = vel_mesh_limits[3] - vymax = vel_mesh_limits[4] - vzmax = vel_mesh_limits[5] - - dvx = (vxmax - vxmin) / (float)(vxcells) - dvy = (vymax - vymin) / (float)(vycells) - dvz = (vzmax - vzmin) / (float)(vzcells) - return [dvx,dvy,dvz] - -def simulation_to_spacecraft_frame(spinvector, detector_axis, phi=0): - ''' Builds a matrix to transform coordinates from simulation frame into spaceraft frame - :param spinvector Spacecraft spin axis, in simulation coordinates - :param detector_axis Detector plane normal axis - :param phi Rotate spacecraft around spin axis after setting up coordinate system - ''' - # TODO: Normalize vectors? - y = np.cross(detector_axis,spinvector) - z = np.cross(spinvector, y) - yr = np.cos(phi)*y - np.sin(phi)*z - zr = np.sin(phi)*y + np.cos(phi)*z - m = np.array([spinvector, yr, zr]) - - return m - -def spacecraft_to_simulation_frame(spinvector, detector_axis, phi=0): - ''' Builds a matrix to transform coordinates from spaceraft frame back to simulation frame - :param spinvector Spacecraft spin axis, in simulation coordinates - :param detector_axis Detector plane normal axis - :param phi Rotate spacecraft around spin axis after setting up coordinate system - ''' - return simulation_to_spacecraft_frame(spinvector,detector_axis,phi).T - -def simulation_to_observation_frame(x_axis,y_axis): - ''' Builds a 3x3 matrix to transform velocities into an observation plane - :param x_axis: x-axis of the observation plane (in simulation coordinates) - :param y_axis: y-axis of the observation plane (gets orthonormalized) - ''' - xn = np.linalg.norm(x_axis) - x_axis /= xn - p = x_axis.dot(y_axis) - y_axis -= p*x_axis - yn = np.linalg.norm(y_axis) - y_axis /= yn - z_axis = np.cross(x_axis,y_axis) - return np.array([x_axis,y_axis,z_axis]) - def themis_plot_detector(vlsvReader, cellID, detector_axis=np.array([0,1,0]), pop="proton"): ''' Plots a view of the detector countrates using matplotlib :param vlsvReader: Some VlsvReader class with a file open @@ -114,7 +55,7 @@ def themis_plot_detector(vlsvReader, cellID, detector_axis=np.array([0,1,0]), po :param detector_axis: detector axis direction (note: this is not spacecraft spin axis!) ''' - matrix = spacecraft_to_simulation_frame(np.cross(np.array([1.,0,0]),detector_axis),detector_axis) + matrix = vsc.spacecraft_to_simulation_frame(np.cross(np.array([1.,0,0]),detector_axis),detector_axis) logging.info("Getting phasespace data...") angles, energies, vmin, vmax, values = themis_observation_from_file( vlsvReader=vlsvReader, @@ -143,7 +84,7 @@ def themis_plot_phasespace_contour(vlsvReader, cellID, plane_x=np.array([1.,0,0] :param plane_x and plane_y: x and y direction of the resulting plot plane ''' - matrix = simulation_to_observation_frame(plane_x,plane_y) + matrix = vsc.simulation_to_observation_frame(plane_x,plane_y) angles, energies, vmin, vmax, values = themis_observation_from_file( vlsvReader=vlsvReader, cellid=cellID, matrix=matrix,pop=pop) @@ -187,7 +128,7 @@ def themis_plot_phasespace_helistyle(vlsvReader, cellID, plane_x=np.array([1.,0, :param plane_x and plane_y: x and y direction of the resulting plot plane ''' - matrix = simulation_to_observation_frame(plane_x,plane_y) + matrix = vsc.simulation_to_observation_frame(plane_x,plane_y) angles, energies, vmin, vmax, values = themis_observation_from_file( vlsvReader=vlsvReader, cellid=cellID, matrix=matrix, countrates=False) if vmin == 0: @@ -248,7 +189,7 @@ def themis_observation_from_file( vlsvReader, cellid, matrix=np.array([[1,0,0],[ pl.show() ''' # Get velocity space resolution - dvx,dvy,dvz = get_dv(vlsvReader,pop=pop) + dvx,dvy,dvz = vlsvReader.get_velocity_mesh_dv(pop=pop) # Read the velocity cells: velocity_cell_data = vlsvReader.read_velocity_cells(cellid,pop=pop)