#!/usr/bin/python import matplotlib matplotlib.use('TkAgg') import calendar import datetime import HeliosE1Plots import glob import numpy as np import matplotlib.pyplot as plt import matplotlib.patches from matplotlib.ticker import AutoMinorLocator from mpl_toolkits.mplot3d.axes3d import Axes3D import operators import pylib from scipy.interpolate import Rbf from scipy.interpolate import griddata import time import TransformCoordSys import scipy.constants import sys class HeliosE1(object): """ ------------------------------------------------------------------------------------------------------------------ description data type HeliosE1.data: data dictionary (shorthand - HeliosE1.d) dict key: property value: array with a length of HeliosE1.dps (np.ndarray) HeliosE1.data[key]: 'fn': filename string array(type='|S24') 'sc': spacecraft (1 or 2 for Helios 1/2) array(type=int8) 'time': seconds since the Unix epoch (1970-01-01 00:00:00), array(type=int32) assuming that 'time' and thus the dataset timestamp represent a time in UTC calculated with calendar.timegm() 'hdm': high data mode (True/False for on/off) array(type=bool) hdm off means normal data mode on, obviously 'shift': shift (True/False for on/off) array(type=bool) 'perihelion_shift': ? (True/False for on/off) array(type=bool) 'minus': indicates a hdm file which contained bad data array(type=bool) (frames), but could be handled as a ndm file 'i3_on': i3 on instead of i1a (True/False for on/off) array(type=bool) 'hel_dist_sun': distance Helios S/C to sun / AU array(type=np.float_) 'hel_car_long': Carrington longitude at Helios S/C / deg array(type=float_) 'hel_car_lat': Carrington latitude at Helios S/C / deg array(type=float_) 'hel_car_rot': Carrington rotation number at Helios S/C array(type=int16) 'ear_car_long': Carrington longitude at earth / deg array(type=float_) 'ear_car_lat': Carrington latitude at earth / deg array(type=float_) 'ear_car_rot': Carrington rotation number at earth array(type=int16) 'angle': Helios S/C - sun - earth angle / deg array(type=float_) 'hel_vel_rad': Helios S/C radial velocity / (AU/day) array(type=float_) 'hel_vel_norm': Helios S/C (in ecliptic) tangential velocity / (AU/day) array(type=float_) 'i1a_pro_den': proton density / (1/cm^3) array(type=float_) 'i1a_pro_vel': proton velocity / (km/s) array(type=float_) 'i1a_pro_temp': proton temperature / K array(type=float_) 'i1a_alp_den': alpha density / (1/cm^3) array(type=float_) 'i1a_alp_vel': alpha velocity / (km/s) array(type=float_) 'i1a_alp_temp': alpha temperature / K array(type=float_) 'i1b_pro_den': proton density / (1/cm^3) array(type=float_) 'i1b_pro_vel': proton velocity / (km/s) array(type=float_) 'i1b_pro_temp': proton temperature / K array(type=float_) 'Bx': x-component of magnetic field / nT array(type=float_) 'By': y-component of magnetic field / nT array(type=float_) 'Bz': z-component of magnetic field / nT array(type=float_) 'sig_Bx': standard deviation of the x-component / nT array(type=float_) 'sig_By': standard deviation of the y-component / nT array(type=float_) 'sig_Bz': standard deviation of the z-component / nT array(type=float_) 'pizzo_corr': ??? pizzo correction / deg array(type=float_) 'eff_sol_diameter': ??? correction for effective solar diameter |v| / ??? array(type=float_) 'eff_sol_azimuth': ??? correction for effective solar azimuth / ??? array(type=float_) '1D_i1a_int': 32 onboard integrated energy channels of proton array(type=np.float_) instrument i1a '1D_i1a_vel': corresponding velocities to 1D_i1a_int array(type=np.float_) '1D_i1b': 32 onboard integrated energy channels of proton array(type=np.float_) instrument i1b '1D_i1b_vel': corresponding velocities to 1D_i1b array(type=np.float_) '3D': '2D': 'ear_dist_sun': distance earth to sun (always zero, historical) array(type=float_) generally disabled, initialize with unknown=True to enable ------------------------------------------------------------------------------------------------------------------ """ def __init__(self, *args, **kwargs): self.load_data(*args, **kwargs) return def load_data(self, path_fn_list, load=['1D', '3D', ], depth=['default', ], init_multiplier=200): """ """ if isinstance(load, str): load = [load, ] self.load = [dim for dim in load if dim in ['1D', '3D', '2D']] datatype = {'fn': '|S24', 'sc': np.int8, 'time': np.int32, 'hdm': bool, 'shift': bool, 'perihelion_shift': bool, 'minus': bool, 'i3_on': bool, 'hel_dist_sun': np.float_, 'hel_car_long': np.float_, 'hel_car_lat': np.float_, 'hel_car_rot': np.int16, 'ear_car_long': np.float_, 'ear_car_lat': np.float_, 'angle': np.float_, 'ear_car_rot': np.int16, 'hel_vel_rad': np.float_, 'hel_vel_norm': np.float_, 'i1a_pro_den': np.float_, 'i1a_pro_vel': np.float_, 'i1a_pro_temp': np.float_, 'i1a_alp_den': np.float_, 'i1a_alp_vel': np.float_, 'i1a_alp_temp': np.float_, 'i1b_pro_den': np.float_, 'i1b_pro_vel': np.float_, 'i1b_pro_temp': np.float_, 'Bx': np.float_, 'By': np.float_, 'Bz': np.float_, 'sig_Bx': np.float_, 'sig_By': np.float_, 'sig_Bz': np.float_, 'ear_dist_sun': np.float_, 'pizzo_corr': np.float_, 'eff_sol_diameter': np.float_, 'eff_sol_azimuth': np.float_, 'i1a_pro_azi': np.float_, 'i1a_pro_ele': np.float_, '1D_i1a_int': np.float_, '1D_i1a_vel': np.float_, '1D_i1b': np.float_, '1D_i1b_vel': np.float_, '2D_azimuth_ch': np.int8, '2D_energy_ch': np.int8, '2D_phasespace_density': np.float_, '2D_counts': np.int32, '2D_vel_x': np.float_, '2D_vel_y': np.float_, '3D_azimuth_ch': np.int8, '3D_elevation_ch': np.int8, '3D_energy_ch': np.int8, '3D_phasespace_density': np.float_, '3D_counts': np.int32, '3D_vel_x': np.float_, '3D_vel_y': np.float_, '3D_vel_z': np.float_, '3D_velocity': np.float_, '3D_azimuth_ang': np.float_, '3D_elevation_ang': np.float_} essential = ['sc', 'time', 'hdm', 'shift', 'minus', 'i3_on', 'hel_dist_sun', 'hel_car_long', 'hel_car_lat', 'hel_car_rot', 'hel_vel_rad', 'hel_vel_norm', 'i1a_pro_den', 'i1a_pro_vel', 'i1a_pro_temp', 'i1b_pro_den', 'i1b_pro_vel', 'i1b_pro_temp', 'Bx', 'By', 'Bz', 'sig_Bx', 'sig_By', 'sig_Bz'] default = ['sc', 'time', 'hdm', 'shift', 'minus', 'i3_on', 'hel_dist_sun', 'hel_car_long', 'hel_car_lat', 'hel_car_rot', 'ear_car_long', 'ear_car_lat', 'angle', 'ear_car_rot', 'hel_vel_rad', 'hel_vel_norm', 'i1a_pro_den', 'i1a_pro_vel', 'i1a_pro_temp', 'i1a_pro_azi', 'i1a_pro_ele', 'i1a_alp_den', 'i1a_alp_vel', 'i1a_alp_temp', 'i1b_pro_den', 'i1b_pro_vel', 'i1b_pro_temp', 'Bx', 'By', 'Bz', 'sig_Bx', 'sig_By', 'sig_Bz'] debug = ['fn', 'sc', 'time', 'hdm', 'shift', 'perihelion_shift', 'minus', 'i3_on', 'hel_dist_sun', 'hel_car_long', 'hel_car_lat', 'hel_car_rot', 'ear_car_long', 'ear_car_lat', 'angle', 'ear_car_rot', 'hel_vel_rad', 'hel_vel_norm', 'i1a_pro_den', 'i1a_pro_vel', 'i1a_pro_temp', 'i1a_alp_den', 'i1a_alp_vel', 'i1a_alp_temp', 'i1b_pro_den', 'i1b_pro_vel', 'i1b_pro_temp', 'Bx', 'By', 'Bz', 'sig_Bx', 'sig_By', 'sig_Bz', 'ear_dist_sun', 'pizzo_corr', 'eff_sol_diameter', 'eff_sol_azimuth', 'i1a_pro_azi', 'i1a_pro_ele'] keys_1D = ['1D_i1a_int', '1D_i1a_vel', '1D_i1b', '1D_i1b_vel'] keys_2D = ['2D_azimuth_ch', '2D_energy_ch', '2D_phasespace_density', '2D_counts', '2D_vel_x', '2D_vel_y'] keys_3D = ['3D_azimuth_ch', '3D_elevation_ch', '3D_energy_ch', '3D_phasespace_density', '3D_counts', '3D_vel_x', '3D_vel_y', '3D_vel_z', '3D_velocity', '3D_azimuth_ang', '3D_elevation_ang'] self.data = pylib.DataDict() self.d = self.data self.TwoFilesSameTime = [] self.TwoFilesInOne = [] self.hdm_fn_file_neq = [] self.shift_fn_file_neq = [] fn_path_list = [] #list of format [[fn_0, path_0], [fn_1, path_1], ... ] if isinstance(path_fn_list, str): path_fn_list = [path_fn_list] for path_fn in path_fn_list: if path_fn.find('/') != -1: fn_path_list.append(path_fn.rsplit('/', 1)[::-1]) fn_path_list[-1][1] = fn_path_list[-1][1] + '/' else: fn_path_list.append([path_fn, './']) ext = fn_path_list[-1][0].rsplit('.', 1)[-1] if len(fn_path_list[-1][0]) != 24 or (ext != '0' and ext != '1'): print "\'{}\' does not have the expected format, its data has not been loaded.".format(fn_path_list[-1][0]) del fn_path_list[-1] fn_path_list.sort() metadata_length = np.int_(len(fn_path_list)) assumed_data_length = np.int_(metadata_length * init_multiplier) load_keys = [] if 'essential' in depth: load_keys = essential if 'default' in depth: load_keys = default if 'debug' in depth: load_keys = debug if '1D' in self.load: for key in keys_1D: self.d[key] = np.zeros(assumed_data_length, dtype=datatype[key]) if '3D' in self.load: for key in keys_3D: self.d[key] = np.zeros(assumed_data_length, dtype=datatype[key]) if '2D' in self.load: for key in keys_2D: self.d[key] = np.zeros(assumed_data_length, dtype=datatype[key]) for key in depth: if key in datatype and key not in load_keys: load_keys.append(key) for key in load_keys: self.d[key] = np.empty(metadata_length, dtype=datatype[key]) expansion = np.empty(assumed_data_length, dtype=np.int_) n_metadata = 0 n_data = 0 skip_next = False for i in np.arange(metadata_length): fn, path = fn_path_list[i] if i < metadata_length - 1: next_fn = fn_path_list[i+1][0] if fn[:18] == next_fn[:18]: self.TwoFilesSameTime.append([fn, next_fn]) skip_next = True continue elif skip_next: skip_next = False continue elif skip_next: skip_next = False continue sc = fn[1:2] timestamp = calendar.timegm(time.strptime(fn[2:18], 'y%yd%jh%Hm%Ms%S')) if fn[19:22] == 'hdm': hdm_fn = '1' else: hdm_fn = '0' shift_fn = fn.rsplit('.', 1)[-1] of = open(path + fn, 'r') pizzo_corr = of.readline().split()[0] eff_sol_diameter, eff_sol_azimuth = of.readline().split() of.readline() hdm_file, shift_file, perihelion_shift, minus, i3_on = of.readline().split() hel_dist_sun, hel_car_long, hel_car_lat = of.readline().split() hel_car_rot = of.readline() ear_dist_sun, ear_car_long, ear_car_lat = of.readline().split() angle, ear_car_rot = of.readline().split() hel_vel_rad, hel_vel_norm = of.readline().split() i1a_pro_den, i1a_pro_vel, i1a_pro_temp = of.readline().split() i1a_pro_azi, i1a_pro_ele = of.readline().split() #values remain unidentified, what do they mean? i1a_alp_den, i1a_alp_vel, i1a_alp_temp = of.readline().split() i1b_pro_den, i1b_pro_vel, i1b_pro_temp = of.readline().split() Bx, By, Bz = of.readline().split() sig_Bx, sig_By, sig_Bz = of.readline().split() i1a_int_1D = [] i1a_vel_1D = [] i1b_1D = [] i1b_vel_1D = [] azimuth_ch_3D = [] elevation_ch_3D = [] energy_ch_3D = [] phasespace_density_3D = [] counts_3D = [] vel_x_3D = [] vel_y_3D = [] vel_z_3D = [] velocity_3D = [] elevation_ang_3D = [] azimuth_ang_3D = [] azimuth_ch_2D = [] energy_ch_2D = [] phasespace_density_2D = [] counts_2D = [] vel_x_2D = [] vel_y_2D = [] if '1D' in self.load: of.readline() i1a_int_1D = of.readline().split() of.readline() i1a_vel_1D = of.readline().split() of.readline() i1b_1D = of.readline().split() of.readline() i1b_vel_1D = of.readline().split() data_found_3D = 0 data_found_2D = 0 two_files_in_one = False for line in of: if line.find('Pizzo') != -1: two_files_in_one = True self.TwoFilesInOne.append(fn) break if '3D' in self.load and data_found_3D != -1: if line.find('2-D') != -1: data_found_3D = -1 data_found_2D = 1 continue if data_found_3D == 3: t = line.split() azimuth_ch_3D.append(t[0]) elevation_ch_3D.append(t[1]) energy_ch_3D.append(t[2]) phasespace_density_3D.append(t[3]) counts_3D.append(t[4]) #don't know if '*****', i.e. values above 99999 occur, therefore no #exception handling here vel_x_3D.append(t[5]) vel_y_3D.append(t[6]) vel_z_3D.append(t[7]) velocity_3D.append(t[8]) azimuth_ang_3D.append(t[9]) elevation_ang_3D.append(t[10]) continue if data_found_3D == 2 and line.find('Max') != -1: data_found_3D = 3 continue if data_found_3D == 1 and line.find('pon') != -1: data_found_3D = 2 continue if line.find('3-D') != -1: data_found_3D = 1 continue if '2D' in self.load and data_found_2D != -1: if line.find('no el') != -1: data_found_2D = -1 continue if data_found_2D == 3: t = line.split() azimuth_ch_2D.append(t[0]) energy_ch_2D.append(t[1]) phasespace_density_2D.append(t[2]) if t[3] == '*****': counts_2D.append('100001') #if '*****', i.e. values above 99999 occur, set to dummy 100001 else: counts_2D.append(t[3]) vel_x_2D.append(t[4]) vel_y_2D.append(t[5]) continue if data_found_2D == 1 and line.find('Max') != -1: data_found_2D = 3 continue if line.find('2-D') != -1: data_found_2D = 1 continue values = {'fn': fn, 'sc': sc, 'time': timestamp, 'hdm': np.int8(hdm_fn), 'shift': np.int8(shift_fn), 'perihelion_shift': np.int8(perihelion_shift), 'minus': np.int8(minus), 'i3_on': np.int8(i3_on) - 1, 'hel_dist_sun': hel_dist_sun, 'hel_car_long': hel_car_long, 'hel_car_lat': hel_car_lat, 'hel_car_rot': hel_car_rot, 'ear_car_long': ear_car_long, 'ear_car_lat': ear_car_lat, 'angle': angle, 'ear_car_rot': ear_car_rot, 'hel_vel_rad': hel_vel_rad, 'hel_vel_norm': hel_vel_norm, 'i1a_pro_den': i1a_pro_den, 'i1a_pro_vel': i1a_pro_vel, 'i1a_pro_temp': i1a_pro_temp, 'i1a_alp_den': i1a_alp_den, 'i1a_alp_vel': i1a_alp_vel, 'i1a_alp_temp': i1a_alp_temp, 'i1b_pro_den': i1b_pro_den, 'i1b_pro_vel': i1b_pro_vel, 'i1b_pro_temp': i1b_pro_temp, 'Bx': Bx, 'By': By, 'Bz': Bz, 'sig_Bx': sig_Bx, 'sig_By': sig_By, 'sig_Bz': sig_Bz, 'ear_dist_sun': ear_dist_sun, 'pizzo_corr': pizzo_corr, 'eff_sol_diameter': eff_sol_diameter, 'eff_sol_azimuth': eff_sol_azimuth, 'i1a_pro_azi': i1a_pro_azi, 'i1a_pro_ele': i1a_pro_ele, '1D_i1a_int': i1a_int_1D, '1D_i1a_vel': i1a_vel_1D, '1D_i1b': i1b_1D, '1D_i1b_vel': i1b_vel_1D, '2D_azimuth_ch': azimuth_ch_2D, '2D_energy_ch': energy_ch_2D, '2D_phasespace_density': phasespace_density_2D, '2D_counts': counts_2D, '2D_vel_x': vel_x_2D, '2D_vel_y': vel_y_2D, '3D_azimuth_ch': azimuth_ch_3D, '3D_elevation_ch': elevation_ch_3D, '3D_energy_ch': energy_ch_3D, '3D_phasespace_density': phasespace_density_3D, '3D_counts': counts_3D, '3D_vel_x': vel_x_3D, '3D_vel_y': vel_y_3D, '3D_vel_z': vel_z_3D, '3D_velocity': velocity_3D, '3D_azimuth_ang': azimuth_ang_3D, '3D_elevation_ang': elevation_ang_3D} if not two_files_in_one: if self.load: add_length = max([len(values[key]) for key in ['1D_i1a_vel', '3D_counts', '2D_counts'] if key[:2] in self.load]) else: add_length = 0 if '1D' in self.load: for key in keys_1D: self.d[key][n_data:n_data+len(values[key])] = values[key] if '3D' in self.load: for key in keys_3D: self.d[key][n_data:n_data+len(values[key])] = values[key] if '2D' in self.load: for key in keys_2D: self.d[key][n_data:n_data+len(values[key])] = values[key] for key in load_keys: self.d[key][n_metadata] = values[key] expansion[n_data:n_data+add_length] = n_metadata if add_length: n_metadata += 1 n_data += add_length if hdm_fn != hdm_file: self.hdm_fn_file_neq.append(fn) if shift_fn != shift_file: self.shift_fn_file_neq.append(fn) print fn of.close() if '1D' in self.load: for key in keys_1D: self.d[key] = self.d[key][:n_data] if '3D' in self.load: for key in keys_3D: self.d[key] = self.d[key][:n_data] if '2D' in self.load: for key in keys_2D: self.d[key] = self.d[key][:n_data] expansion = expansion[:n_data] first_key = load_keys.pop() self.d.add_reduced_data(first_key, self.d[first_key][:n_metadata], 'metadata expansion', expansion) for key in load_keys: self.d.add_reduced_data(key, self.d[key][:n_metadata], 'metadata expansion') if self.hdm_fn_file_neq: print "For the files:" for fn in self.hdm_fn_file_neq: print "\t{}".format(fn) print "the hdm flag of the filename is different than the hdm flag in the file itself, please check." if self.shift_fn_file_neq: print "For the files:" for fn in self.shift_fn_file_neq: print "\t{}".format(fn) print "the shift flag of the filename is different than the shift flag in the file itself, please check." return