2525import analysator as pt
2626import matplotlib .pyplot as plt
2727import matplotlib
28- from .rotation import rotateVectorToVector
2928from scipy .interpolate import griddata
3029from scipy .signal import sepfir2d
31- from packaging .version import Version
3230import logging
31+ from analysator .calculations import virtual_observations as vsc
3332
3433# Detector data obtained from the Themis ESA instrument paper
3534# http://dx.doi.org/10.1007/s11214-008-9440-2
4847themis_colors = [(0 ,0 ,0 ),(.5 ,0 ,.5 ),(0 ,0 ,1 ),(0 ,1 ,1 ),(0 ,1 ,0 ),(1 ,1 ,0 ),(1 ,0 ,0 )]
4948themis_colormap = matplotlib .colors .LinearSegmentedColormap .from_list ("themis" ,themis_colors )
5049
51- def get_dv (vlsvReader ,pop = "proton" ):
52- # Get velocity grid sizes:
53- vel_mesh_size = vlsvReader .get_velocity_mesh_size (pop = pop )
54- vel_block_size = vlsvReader .get_velocity_block_size (pop = pop )
55- vxcells = vel_mesh_size [0 ]* vel_block_size [0 ]
56- vycells = vel_mesh_size [1 ]* vel_block_size [1 ]
57- vzcells = vel_mesh_size [2 ]* vel_block_size [2 ]
58-
59- vel_mesh_limits = vlsvReader .get_velocity_mesh_extent (pop = pop )
60- vxmin = vel_mesh_limits [0 ]
61- vymin = vel_mesh_limits [1 ]
62- vzmin = vel_mesh_limits [2 ]
63- vxmax = vel_mesh_limits [3 ]
64- vymax = vel_mesh_limits [4 ]
65- vzmax = vel_mesh_limits [5 ]
66-
67- dvx = (vxmax - vxmin ) / (float )(vxcells )
68- dvy = (vymax - vymin ) / (float )(vycells )
69- dvz = (vzmax - vzmin ) / (float )(vzcells )
70- return [dvx ,dvy ,dvz ]
71-
72- def simulation_to_spacecraft_frame (spinvector , detector_axis , phi = 0 ):
73- ''' Builds a matrix to transform coordinates from simulation frame into spaceraft frame
74- :param spinvector Spacecraft spin axis, in simulation coordinates
75- :param detector_axis Detector plane normal axis
76- :param phi Rotate spacecraft around spin axis after setting up coordinate system
77- '''
78- # TODO: Normalize vectors?
79- y = np .cross (detector_axis ,spinvector )
80- z = np .cross (spinvector , y )
81- yr = np .cos (phi )* y - np .sin (phi )* z
82- zr = np .sin (phi )* y + np .cos (phi )* z
83- m = np .array ([spinvector , yr , zr ])
84-
85- return m
86-
87- def spacecraft_to_simulation_frame (spinvector , detector_axis , phi = 0 ):
88- ''' Builds a matrix to transform coordinates from spaceraft frame back to simulation frame
89- :param spinvector Spacecraft spin axis, in simulation coordinates
90- :param detector_axis Detector plane normal axis
91- :param phi Rotate spacecraft around spin axis after setting up coordinate system
92- '''
93- return simulation_to_spacecraft_frame (spinvector ,detector_axis ,phi ).T
94-
95- def simulation_to_observation_frame (x_axis ,y_axis ):
96- ''' Builds a 3x3 matrix to transform velocities into an observation plane
97- :param x_axis: x-axis of the observation plane (in simulation coordinates)
98- :param y_axis: y-axis of the observation plane (gets orthonormalized)
99- '''
100- xn = np .linalg .norm (x_axis )
101- x_axis /= xn
102- p = x_axis .dot (y_axis )
103- y_axis -= p * x_axis
104- yn = np .linalg .norm (y_axis )
105- y_axis /= yn
106- z_axis = np .cross (x_axis ,y_axis )
107- return np .array ([x_axis ,y_axis ,z_axis ])
108-
109- def themis_plot_detector (vlsvReader , cellID , detector_axis = np .array ([0 ,1 ,0 ]), pop = "proton" ):
50+ def themis_plot_detector (vlsvReader , cellID ,outputfile = "./themis_plot_detector.png" ,nooverwrite = False ,draw = True , detector_axis = np .array ([0 ,1 ,0 ]), pop = "proton" ):
11051 ''' Plots a view of the detector countrates using matplotlib
11152 :param vlsvReader: Some VlsvReader class with a file open
11253 :type vlsvReader: :class:`vlsvfile.VlsvReader`
11354 :param cellid: The cell id where the distribution is supposet to be sampled NOTE: The cell id must have a velocity distribution!
11455 :param detector_axis: detector axis direction (note: this is not spacecraft spin axis!)
56+ :param draw: Set to false to save to file instead of drawing on screen
57+ :kward outputfile: File to output the image to if Draw=False
58+ :kward nooverwrite: Whether to allow overwriting of files when saving, Default False
59+
11560 '''
11661
117- matrix = spacecraft_to_simulation_frame (np .cross (np .array ([1. ,0 ,0 ]),detector_axis ),detector_axis )
62+ matrix = vsc . spacecraft_to_simulation_frame (np .cross (np .array ([1. ,0 ,0 ]),detector_axis ),detector_axis )
11863
11964 logging .info ("Getting phasespace data..." )
12065 angles , energies , vmin , vmax , values = themis_observation_from_file ( vlsvReader = vlsvReader ,
@@ -133,17 +78,24 @@ def themis_plot_detector(vlsvReader, cellID, detector_axis=np.array([0,1,0]), po
13378 cax = ax .pcolormesh (grid_theta ,grid_r ,values , norm = matplotlib .colors .LogNorm (vmin = vmin ,vmax = vmax ), cmap = themis_colormap )
13479 ax .grid (True )
13580 fig .colorbar (cax )
136- plt .show ()
81+ if not draw :
82+ outputpath = pt .plot .output_path (outputfile ,None ,None ,nooverwrite )
83+ plt .savefig (outputpath )
84+ else :
85+ plt .show ()
13786
138- def themis_plot_phasespace_contour (vlsvReader , cellID , plane_x = np .array ([1. ,0 ,0 ]), plane_y = np .array ([0 ,0 ,1. ]), smooth = False , xlabel = "Vx" , ylabel = "Vy" , pop = "proton" ):
87+ def themis_plot_phasespace_contour (vlsvReader , cellID ,outputfile = './themis_plot_phasespace_contour.png' , nooverwrite = False , draw = True , plane_x = np .array ([1. ,0 ,0 ]), plane_y = np .array ([0 ,0 ,1. ]), smooth = False , xlabel = "Vx" , ylabel = "Vy" , pop = "proton" ):
13988 ''' Plots a contour view of phasespace, as seen by a themis detector, at the given cellID
14089 :param vlsvReader: Some VlsvReader class with a file open
14190 :type vlsvReader: :class:`vlsvfile.VlsvReader`
14291 :param cellid: The cell id where the distribution is supposet to be sampled NOTE: The cell id must have a velocity distribution!
92+ :param draw: Set to false to save to file instead of drawing on screen
93+ :kward outputfile: File to output the image to if Draw=False
94+ :kward nooverwrite: Whether to allow overwriting of files when saving, Default False
14395 :param plane_x and plane_y: x and y direction of the resulting plot plane
14496 '''
14597
146- matrix = simulation_to_observation_frame (plane_x ,plane_y )
98+ matrix = vsc . simulation_to_observation_frame (plane_x ,plane_y )
14799
148100 angles , energies , vmin , vmax , values = themis_observation_from_file ( vlsvReader = vlsvReader , cellid = cellID , matrix = matrix ,pop = pop )
149101
@@ -176,18 +128,25 @@ def themis_plot_phasespace_contour(vlsvReader, cellID, plane_x=np.array([1.,0,0]
176128 cax = ax .contour (xi ,yi ,vi .T , levels = np .logspace (np .log10 (vmin ),np .log10 (vmax ),20 ), norm = matplotlib .colors .LogNorm (vmin = vmin ,vmax = vmax ))
177129 ax .grid (True )
178130 fig .colorbar (cax )
179- plt .show ()
131+ if not draw :
132+ outputpath = pt .plot .output_path (outputfile ,None ,None ,nooverwrite )
133+ plt .savefig (outputpath )
134+ else :
135+ plt .show ()
180136
181- def themis_plot_phasespace_helistyle (vlsvReader , cellID , plane_x = np .array ([1. ,0 ,0 ]), plane_y = np .array ([0 ,0 ,1. ]), smooth = True , xlabel = "Vx" , ylabel = "Vy" ):
137+ def themis_plot_phasespace_helistyle (vlsvReader , cellID ,outputfile = './themis_plot_phasespace_helistyle' , plane_x = np .array ([1. ,0 ,0 ]), plane_y = np .array ([0 ,0 ,1. ]), smooth = True , xlabel = "Vx" , ylabel = "Vy" , draw = True , nooverwrite = False ):
182138 ''' Plots a view of phasespace, as seen by a themis detector, at the given cellID, in the style that heli likes.
183139 :param vlsvReader: Some VlsvReader class with a file open
184140 :type vlsvReader: :class:`vlsvfile.VlsvReader`
185141 :param cellid: The cell id where the distribution is supposet to be sampled NOTE: The cell id must have a velocity distribution!
186142 :param smooth: Smooth re-gridded phasespace before plotting
143+ :param draw: Set to false to save to file instead of drawing on screen
144+ :kward outputfile: File to output the image to if Draw=False
145+ :kward nooverwrite: Whether to allow overwriting of files when saving, Default False
187146 :param plane_x and plane_y: x and y direction of the resulting plot plane
188147 '''
189148
190- matrix = simulation_to_observation_frame (plane_x ,plane_y )
149+ matrix = vsc . simulation_to_observation_frame (plane_x ,plane_y )
191150
192151 angles , energies , vmin , vmax , values = themis_observation_from_file ( vlsvReader = vlsvReader , cellid = cellID , matrix = matrix , countrates = False )
193152 if vmin == 0 :
@@ -218,12 +177,16 @@ def themis_plot_phasespace_helistyle(vlsvReader, cellID, plane_x=np.array([1.,0,
218177 ax .set_ylabel (ylabel + " (km/s)" )
219178 cmapuse = pt .plot .get_cmap ("Blues" )
220179
221- cax = ax .pcolormesh (xi ,yi ,vi .T , norm = matplotlib .colors .LogNorm (vmin = vmin ,vmax = vmax ), vmin = vmin , vmax = vmax , cmap = cmapuse , shading = 'flat' )
180+ cax = ax .pcolormesh (xi ,yi ,vi .T , norm = matplotlib .colors .LogNorm (vmin = vmin ,vmax = vmax ), cmap = cmapuse )
222181 cax2 = ax .contourf (xi ,yi ,vi .T , levels = np .logspace (np .log10 (vmin ),np .log10 (vmax ),20 ), norm = matplotlib .colors .LogNorm (vmin = vmin ,vmax = vmax ), vmin = vmin , vmax = vmax , cmap = cmapuse )
223182 #cax3 = ax.contour(xi,yi,vi.T, levels=np.logspace(np.log10(vmin),np.log10(vmax),20), norm=matplotlib.colors.LogNorm(vmin=vmin,vmax=vmax), cmap=pl.get_cmap("binary"))
224183 ax .grid (True )
225184 fig .colorbar (cax )
226- plt .show ()
185+ if not draw :
186+ outputpath = pt .plot .output_path (outputfile ,None ,None ,nooverwrite )
187+ fig .savefig (outputpath )
188+ else :
189+ plt .show ()
227190def themis_observation_from_file ( vlsvReader , cellid , matrix = np .array ([[1 ,0 ,0 ],[0 ,1 ,0 ],[0 ,0 ,1 ]]), countrates = True , interpolate = True ,binOffset = [0. ,0. ],pop = 'proton' ):
228191 ''' Calculates artificial THEMIS EMS observation from the given cell
229192 :param vlsvReader: Some VlsvReader class with a file open
@@ -248,7 +211,7 @@ def themis_observation_from_file( vlsvReader, cellid, matrix=np.array([[1,0,0],[
248211 pl.show()
249212 '''
250213 # Get velocity space resolution
251- dvx ,dvy ,dvz = get_dv ( vlsvReader , pop = pop )
214+ dvx ,dvy ,dvz = vlsvReader . get_velocity_mesh_dv ( pop = pop )
252215
253216 # Read the velocity cells:
254217 velocity_cell_data = vlsvReader .read_velocity_cells (cellid ,pop = pop )
0 commit comments