diff --git a/NuRadioReco/examples/LOFAR/run_lofar_example.py b/NuRadioReco/examples/LOFAR/run_lofar_example.py new file mode 100644 index 00000000..3257cf08 --- /dev/null +++ b/NuRadioReco/examples/LOFAR/run_lofar_example.py @@ -0,0 +1,16 @@ +import matplotlib.pyplot as plt + +import NuRadioReco.modules.io.lofar.readLOFAR +from NuRadioReco.utilities import units + +import logging +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger('readARAexample') + +# Initialize the reader with the file you want to read +LOFARreader = NuRadioReco.modules.io.lofar.readLOFAR.readLOFAR(['/Users/anelles/Experiments/LOFAR/A_NuRadioReco_Test/L78862_D20121205T051644.039Z_CS002_R000_tbb.h5'],DAL=1) + +for iE, event in enumerate(LOFARreader.run()): + print(event.get_id()) + + diff --git a/NuRadioReco/modules/io/lofar/README.md b/NuRadioReco/modules/io/lofar/README.md new file mode 100644 index 00000000..d9fa8af1 --- /dev/null +++ b/NuRadioReco/modules/io/lofar/README.md @@ -0,0 +1,16 @@ +First attempt in writing an LOFAR reader to look into using NuRadioReco for LOFAR + +The files: +metadata.py, raw_tbb_IO.py, readLOFAR.py, utilities.py +have been copied from Brian Hare +https://github.com/Bhare8972/LOFAR-LIM/blob/master/LIM_scripts/ + +His python package is very useful, but does not follow python convention +for modules, which it is copied rather than imported. This should be worked +on in the future. + +The files also need LOFAR metadata to be downloaded: +https://github.com/Bhare8972/LOFAR-LIM/blob/master/LIM_scripts/data.tar.gz +and un-tarred in the folder. + + diff --git a/NuRadioReco/modules/io/lofar/__init__.py b/NuRadioReco/modules/io/lofar/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/NuRadioReco/modules/io/lofar/metadata.py b/NuRadioReco/modules/io/lofar/metadata.py new file mode 100755 index 00000000..938c716a --- /dev/null +++ b/NuRadioReco/modules/io/lofar/metadata.py @@ -0,0 +1,736 @@ +#!/usr/bin/env python3 + +"""This module reads in calibration metadata from file in the early fases of LOFAR. In the future this should be replaced by reading the metadata from the files. + +.. moduleauthor:: Sander ter Veen + +Modified by Brian Hare for use with LOFAR for Lightning Imaging +""" + +## Imports +import numpy as np +import struct + +from NuRadioReco.modules.io.lofar.utilities import SId_to_Sname, latlonCS002, RTD, MetaData_directory + +#### first some simple utility functions ### + +def mapAntennasetKeyword(antennaset): + """Ugly fix to map correct antenna names in input to wrong antenna names + for metadata module. + """ + + # Strip whitespace + antennaset = antennaset.strip() + + allowed = ["LBA_OUTER", "LBA_INNER", "LBA_X", "LBA_Y", "HBA", "HBA_0", "HBA_1"] + incorrect = {'LBA_INNER': 'LBA_INNER', + 'LBA_OUTER': 'LBA_OUTER', + 'HBA_ZERO': 'HBA_0', + 'HBA_ONE': 'HBA_1', + 'HBA_DUAL': 'HBA', + 'HBA_JOINED': 'HBA', + 'HBA_ZERO_INNER': 'HBA_0', # Only true for core stations + 'HBA_ONE_INNER': 'HBA_1', # Only true for core stations + 'HBA_DUAL_INNER': 'HBA', # Only true for core stations + 'HBA_JOINED_INNER': 'HBA'} # Only true for core stations + + if antennaset in incorrect: + antennaset = incorrect[antennaset] + elif antennaset == "HBA_BOTH": + # This keyword is also wrong but present in file headers + print( "Keyword " + antennaset + " does not comply with ICD, mapping...") + antennaset = "HBA" + + assert antennaset in allowed + + return antennaset + +def make_antennaID_filter(antennaIDs): + """For a list of antennaIDs, return a filter to filter data by antenna. + example use: + getStationPhaseCalibration("CS001", "LBA_OUTER")[ make_antennaID_filter(["002000001"]) ] + + note: Only works for one station at a time. Assumes that the array you want to filter includes ALL antennas in the appropriate antenna set""" + + RCU_id = np.array([int(ID[-3:]) for ID in antennaIDs]) + return RCU_id + + + + + + +#### read callibration data ### + +def getStationPhaseCalibration(station, antennaset, file_location=None): + """Read phase calibration data for a station. + + Required arguments: + + ================== ==================================================== + Parameter Description + ================== ==================================================== + *station* station name (as str) or ID. + *mode* observation mode. + ================== ==================================================== + + returns weights for 512 subbands. + + Examples:: + + >>> metadata.getStationPhaseCalibration("CS002","LBA_OUTER") + array([[ 1.14260161 -6.07397622e-18j, 1.14260161 -6.05283530e-18j, + 1.14260161 -6.03169438e-18j, ..., 1.14260161 +4.68675289e-18j, + 1.14260161 +4.70789381e-18j, 1.14260161 +4.72903474e-18j], + [ 0.95669876 +2.41800591e-18j, 0.95669876 +2.41278190e-18j, + 0.95669876 +2.40755789e-18j, ..., 0.95669876 -2.41017232e-19j, + 0.95669876 -2.46241246e-19j, 0.95669876 -2.51465260e-19j], + [ 0.98463207 +6.80081617e-03j, 0.98463138 +6.89975906e-03j, + 0.98463069 +6.99870187e-03j, ..., 0.98299670 +5.71319125e-02j, + 0.98299096 +5.72306908e-02j, 0.98298520 +5.73294686e-02j], + ..., + [ 1.03201290 +7.39535744e-02j, 1.03144532 +8.14880844e-02j, + 1.03082273 +8.90182487e-02j, ..., -0.82551740 -6.23731331e-01j, + -0.82094046 -6.29743206e-01j, -0.81631975 -6.35721497e-01j], + [ 1.12370332 -1.15296909e-01j, 1.12428451 -1.09484545e-01j, + 1.12483564 -1.03669252e-01j, ..., -0.92476286 +6.48703460e-01j, + -0.92810503 +6.43912711e-01j, -0.93142239 +6.39104744e-01j], + [ 1.10043006 -6.18995646e-02j, 1.10075250 -5.58731668e-02j, + 1.10104193 -4.98450938e-02j, ..., -1.01051042 +4.40052904e-01j, + -1.01290481 +4.34513198e-01j, -1.01526883 +4.28960464e-01j]]) + + >>> metadata.getStationPhaseCalibration(122,"LBA_OUTER") + Calibration data not yet available. Returning 1 + array([[ 1.+0.j, 1.+0.j, 1.+0.j, ..., 1.+0.j, 1.+0.j, 1.+0.j], + [ 1.+0.j, 1.+0.j, 1.+0.j, ..., 1.+0.j, 1.+0.j, 1.+0.j], + [ 1.+0.j, 1.+0.j, 1.+0.j, ..., 1.+0.j, 1.+0.j, 1.+0.j], + ..., + [ 1.+0.j, 1.+0.j, 1.+0.j, ..., 1.+0.j, 1.+0.j, 1.+0.j], + [ 1.+0.j, 1.+0.j, 1.+0.j, ..., 1.+0.j, 1.+0.j, 1.+0.j], + [ 1.+0.j, 1.+0.j, 1.+0.j, ..., 1.+0.j, 1.+0.j, 1.+0.j]]) + + """ + + # Return mode nr depending on observation mode + antennasetToMode = {"LBA_OUTER": "1", + "LBA_INNER": "3", + "HBA": "5", + "HBA_0": "5", + "HBA_1": "5", + } + + antennaset = mapAntennasetKeyword( antennaset ) + + if antennaset not in antennasetToMode.keys(): + raise KeyError("Not a valid antennaset " + antennaset) + + modenr = antennasetToMode[antennaset] + if not isinstance(station, str): + # Convert a station id to a station name + station = SId_to_Sname[station] + + stationNr = station[2:] + + # filename + if file_location is None: + file_location = MetaData_directory + '/lofar/StaticMetaData/CalTables' + + filename = file_location + '/CalTable_' + stationNr + '_mode' + modenr + '.dat' + with open(filename, 'rb') as fin: + # Test for header record above raw data - present in newer caltables (starting 2012) + line = fin.readline().decode() + if 'HeaderStart' in line: + while not 'HeaderStop' in line: + line = fin.readline().decode() + else: # no header present, seek to starting position + fin.seek(0) + + data = np.fromfile(fin, dtype=np.double) + + data.resize(512, 96, 2) + + complexdata = np.empty(shape=(512, 96), dtype=complex) + complexdata.real = data[:, :, 0] + complexdata.imag = data[:, :, 1] + + return complexdata.transpose() + +def convertPhase_to_Timing(phase_calibration, sample_time=5.0e-9): + """Given the phase calibration of the 512 LOFAR subbands, such as the output of getStationPhaseCalibration, return the timing callibration of each antenna. + Not sure how well this works with HBA antennas. Sample time should be seconds per sample. Default is 5 ns""" + phases = np.angle(phase_calibration) + delays = (phases[:, 1] - phases[:, 0]) * (1024 / (2*np.pi)) * sample_time ## this just finds the slope from the first two points. Are there better methods? + ### TODO: add a conditional that takes different points if the slope is too large + return delays + + + +#def getStationGainCalibration(station, antennaset, file_location=None): +# """Read phase calibration data for a station. +# +# Required arguments: +# +# ================== ==================================================== +# Parameter Description +# ================== ==================================================== +# *station* station name or ID. +# *mode* observation mode. +# ================== ==================================================== +# +# Optional arguments: +# +# ================== ==================================================== +# Parameter Description +# ================== ==================================================== +# *return_as_hArray* Default False +# ================== ==================================================== +# +# returns one gain per RCU. This gain is calculated using the absolute +# value from the CalTables assuming these are not frequency dependent. +# This seems to be true in current (2013-08) tables. +# """ +# +# cal = getStationPhaseCalibration(station, antennaset, file_location) +# +# gain = np.abs(cal[:,0]) +# return gain + + + +#### information about cable lengths and delays #### + +def getCableDelays(station, antennaset): + """ Get cable delays in seconds. + + Required arguments: + + ================== ==================================================== + Parameter Description + ================== ==================================================== + *station* Station name or ID e.g. "CS302", 142 + *antennaset* Antennaset used for this station. Options: + + * LBA_INNER + * LBA_OUTER + * LBA_X + * LBA_Y + * LBA_SPARSE0 + * LBA_SPARSE1 + * HBA_0 + * HBA_1 + * HBA + + ================== ==================================================== + + returns "array of (rcus * cable delays ) for all dipoles in a station" + + """ + + # Check station id type + if not isinstance(station, str): + # Convert a station id to a station name + station = SId_to_Sname[station] + + antennaset = mapAntennasetKeyword( antennaset ) + + if "LBA_OUTER" == antennaset: + rcu_connection = "LBL" + elif "LBA_INNER" == antennaset: + rcu_connection = "LBH" + elif antennaset in ['HBA', "HBA_1", "HBA_0"]: + rcu_connection = "HBA" + else: + raise KeyError("Not a valid antennaset " + antennaset) + + cabfilename = MetaData_directory + '/lofar/StaticMetaData/CableDelays/' + station + '-CableDelays.conf' + cabfile = open(cabfilename) + + cable_delays = np.zeros(96) + + str_line = '' + while "RCUnr" not in str_line: + str_line = cabfile.readline() + if len(str_line) == 0: + # end of file reached, no data available + assert False + + str_line = cabfile.readline() + for i in range(96): + str_line = cabfile.readline() + sep_line = str_line.split() + if rcu_connection == "LBL": + cable_delays[int(sep_line[0])] = float(sep_line[2]) * 1e-9 + elif rcu_connection == "LBH": + cable_delays[int(sep_line[0])] = float(sep_line[4]) * 1e-9 + elif rcu_connection == "HBA": + cable_delays[int(sep_line[0])] = float(sep_line[6]) * 1e-9 + + return cable_delays + +def getCableLength(station,antennaset): + + # Check station id type + if not isinstance(station, str): + # Convert a station id to a station name + station = SId_to_Sname[station] + + antennaset = mapAntennasetKeyword( antennaset ) + + if "LBA_OUTER" == antennaset: + rcu_connection = "LBL" + elif "LBA_INNER" == antennaset: + rcu_connection = "LBH" + elif antennaset in ['HBA', "HBA_1", "HBA_0"]: + rcu_connection = "HBA" + else: + raise KeyError("Not a valid antennaset " + antennaset) + + cabfilename = MetaData_directory + '/lofar/StaticMetaData/CableDelays/' + station + '-CableDelays.conf' + cabfile = open(cabfilename) + + cable_length = np.zeros(96) + + str_line = '' + while "RCUnr" not in str_line: + str_line = cabfile.readline() + if len(str_line) == 0: + # end of file reached, no data available + assert False + + str_line = cabfile.readline() + for i in range(96): + str_line = cabfile.readline() + sep_line = str_line.split() + if rcu_connection == "LBL": + cable_length[int(sep_line[0])] = float(sep_line[1]) + elif rcu_connection == "LBH": + cable_length[int(sep_line[0])] = float(sep_line[3]) + elif rcu_connection == "HBA": + cable_length[int(sep_line[0])] = float(sep_line[5]) + + return cable_length + +def antennaset2rcumode(antennaset,filter): + antennaset = mapAntennasetKeyword( antennaset ) + rcumode=dict() + rcumode[('LBA_INNER','LBA_10_90')]=1 + rcumode[('LBA_OUTER','LBA_10_90')]=2 + rcumode[('LBA_INNER','LBA_30_90')]=3 + rcumode[('LBA_OUTER','LBA_30_90')]=4 + rcumode[('HBA','HBA_110_190')]=5 + rcumode[('HBA','HBA_170_230')]=6 + rcumode[('HBA','HBA_210_250')]=7 + if 'HBA' in antennaset: + antennaset='HBA' + return rcumode[(antennaset,filter)] + +def getCableAttenuation(station,antennaset,filter=None): + + cable_length = getCableLength(station, antennaset) + + attenuationFactor=dict() + attenuationFactor[1]=-0.0414#{50:-2.05,80:-3.32,85:-3.53,115:-4.74,130:-5.40} + attenuationFactor[2]=attenuationFactor[1] + attenuationFactor[3]=attenuationFactor[1] + attenuationFactor[4]=attenuationFactor[1] + attenuationFactor[5]=-0.0734#{50:-3.64,80:-5.87,85:-6.22,115:-8.53,130:-9.52} + attenuationFactor[6]=-0.0848#{50:-4.24,80:-6.82,85:-7.21,115:-9.70,130:-11.06} + attenuationFactor[7]=-0.0892#{50:-4.46,80:-7.19,85:-7.58,115:-10.18,130:-11.61} + if filter==None: + if 'LBA' in antennaset: + filter='LBA_30_90' + else: + print( "Please specify the filter!") + filter='HBA_110_190' + rcumode=antennaset2rcumode(antennaset,filter) + att=attenuationFactor[rcumode] + return cable_length*att + + + + + + + + + + + +#### functions for antenna and station location ##### + +def getItrfAntennaPosition(station, antennaset): + """Returns the antenna positions of all the antennas in the station + in ITRF coordinates for the specified antennaset. + station can be the name or id of the station. + + Required arguments: + + =================== ============================================== + Parameter Description + =================== ============================================== + *station* Name or id of the station. e.g. "CS302" or 142 + *antennaset* Antennaset used for this station. Options: + + * LBA_INNER + * LBA_OUTER + * LBA_X + * LBA_Y + * LBA_SPARSE0 + * LBA_SPARSE1 + * HBA_0 + * HBA_1 + * HBA + + =================== ============================================== + + """ + # Check station id type + if isinstance(station, int): + # Convert a station id to a station name + station = SId_to_Sname[station] + + antennaset = mapAntennasetKeyword( antennaset ) + + if "LBA" in antennaset: + antennatype = "LBA" + elif "HBA" in antennaset: + antennatype = "HBA" + + + # Obtain filename of antenna positions + filename = MetaData_directory + "/lofar/StaticMetaData/AntennaFields/" + station + "-AntennaField.conf" + + # Open file + f = open(filename, 'r') + + if station[0:2] != "CS": + if "HBA" in antennaset: + antennaset = "HBA" + + # Find position of antennaset in file + str_line = '' + while antennatype != str_line.strip(): + str_line = f.readline() + if len(str_line) == 0: + # end of file reached, no data available + assert False + + # Find the location of the station. Antenna locations are relative to this + str_line = f.readline() + str_split = str_line.split() + stationX = float(str_split[2]) + stationY = float(str_split[3]) + stationZ = float(str_split[4]) + + str_line = f.readline() + + # Get number of antennas and the number of directions + nrantennas = int(str_line.split()[0]) + nrdir = int(str_line.split()[4]) + + antenna_positions = np.empty((2*nrantennas, nrdir), dtype=np.double) + for i in range(nrantennas): + line = f.readline().split() + + antenna_positions[2*i, 0] = float(line[0]) + stationX + antenna_positions[2*i, 1] = float(line[1]) + stationY + antenna_positions[2*i, 2] = float(line[2]) + stationZ + + antenna_positions[2*i+1, 0] = float(line[3]) + stationX + antenna_positions[2*i+1, 1] = float(line[4]) + stationY + antenna_positions[2*i+1, 2] = float(line[5]) + stationZ + + if antennatype == "LBA": + # There are three types of feed + # H for HBA + # h for lbh + # l for lbl + feed = {} + feed["CS"] = {} + feed["RS"] = {} + feed["DE"] = {} + feed["CS"]["LBA_SPARSE_EVEN"] = "24llhh" + feed["CS"]["LBA_SPARSE_ODD"] = "24hhll" + feed["CS"]["LBA_X"] = "48hl" + feed["CS"]["LBA_Y"] = "48lh" + feed["CS"]["LBA_INNER"] = "96h" + feed["CS"]["LBA_OUTER"] = "96l" + feed["RS"]["LBA_SPARSE_EVEN"] = "24llhh" + feed["RS"]["LBA_SPARSE_ODD"] = "24hhll" + feed["RS"]["LBA_X"] = "48hl" + feed["RS"]["LBA_Y"] = "48lh" + feed["RS"]["LBA_INNER"] = "96h" + feed["RS"]["LBA_OUTER"] = "96l" + feed["DE"]["LBA"] = "192h" + if station[0:2] == "CS" or "RS": + feedsel = feed[station[0:2]][antennaset] + nrset = int(feedsel.split('l')[0].split('h')[0].split('H')[0]) + feeds = '' + feedsel = feedsel[len(str(nrset)):] + for i in range(nrset): + feeds += feedsel + + indexselection = [] + for i in range(len(feeds)): + if feeds[i] == 'l': + # The 'l' feeds are the last 96 numbers of the total list + indexselection.append(i + 96) + elif feeds[i] == 'h': + # The 'h' feeds are the first 96 numbers of the total list + indexselection.append(i) + else: + # This selection is not yet supported + assert False + antenna_positions = antenna_positions[indexselection] + + return antenna_positions + +def getStationPositions(station, antennaset, coordinatesystem): + """Returns the antenna positions of all the antennas in the station + relative to the station center for the specified antennaset. + station can be the name or id of the station. + + Required arguments: + + ================== ============================================== + Argument Description + ================== ============================================== + *station* Name or id of the station. e.g. "CS302" or 142 + *antennaset* Antennaset used for this station. Options: + + * LBA_INNER + * LBA_OUTER + * LBA_X + * LBA_Y + * LBA_SPARSE0 + * LBA_SPARSE1 + * HBA_0 + * HBA_1 + * HBA + + *coordinatesystem* WGS84 or ITRF + + output: + if coordinatesystem == "WGS84", then return [lat, lon, alt] as a numpy array + else if coordinatesystem=="ITRF", then return [X, Y, Z] as a numpy array + + """ + + # Check if requested antennaset is known + assert coordinatesystem in ["WGS84", 'ITRF'] + + # Check station id type + if isinstance(station, int): + # Convert a station id to a station name + station = SId_to_Sname[station] + + antennaset = mapAntennasetKeyword( antennaset ) + + # Obtain filename of antenna positions + if 'WGS84' in coordinatesystem: + filename = MetaData_directory + "/lofar/StaticMetaData/AntennaArrays/" + station + "-AntennaArrays.conf" + else: + filename = MetaData_directory + "/lofar/StaticMetaData/AntennaFields/" + station + "-AntennaField.conf" + + # Open file + f = open(filename, 'r') + + if "LBA" in antennaset: + antennaset = "LBA" + + if station[0:2] != "CS": + if "HBA" in antennaset: + antennaset = "HBA" + + # Find position of antennaset in file + str_line = '' + while antennaset != str_line.strip(): + str_line = f.readline() + if len(str_line) == 0: + # end of file reached, no data available + print( "Antenna set not found in calibration file", filename) + return None + + # Skip name and station reference position + str_line = f.readline().split() + + A = float(str_line[2]) ## lon in WGS84, X in ITRF + B = float(str_line[3]) ## lat in WGS84, Y in ITRF + C = float(str_line[4]) ## alt in WGS84, Z in ITRF + + return np.array([A,B,C]) + +ITRFCS002 = getStationPositions('CS002', 'LBA_OUTER', coordinatesystem='ITRF') # ($LOFARSOFT/data/lofar/StaticMetaData/AntennaFields/CS002-AntennaField.conf) + +def convertITRFToLocal(itrfpos, phase_center=ITRFCS002, reflatlon=latlonCS002, out=None): + """ + ================== ============================================== + Argument Description + ================== ============================================== + *itrfpos* an ITRF position as 1D numpy array, or list of positions as a 2D array + *phase_center* the origin of the coordinate system, in ITRF. Default is CS002. + *reflatlon* the rotation of the coordinate system. Is the [lat, lon] (in degrees) on the Earth which defines "UP" + + Function returns a 2D numpy array (even if input is 1D). + Out cannot be same array as itrfpos + """ + if out is itrfpos: + print("out cannot be same as itrfpos in convertITRFToLocal. TODO: make this a real error") + quit() + + lat = reflatlon[0]/RTD + lon = reflatlon[1]/RTD + arg0 = np.array([-np.sin(lon), -np.sin(lat) * np.cos(lon), np.cos(lat) * np.cos(lon)]) + arg1 = np.array([np.cos(lon) , -np.sin(lat) * np.sin(lon), np.cos(lat) * np.sin(lon)]) + arg2 = np.array([ 0.0, np.cos(lat), np.sin(lat)]) + + if out is None: + ret = np.empty(itrfpos.shape, dtype=np.double ) + else: + ret = out + + ret[:] = np.outer(itrfpos[...,0]-phase_center[0], arg0 ) + ret += np.outer(itrfpos[...,1]-phase_center[1], arg1 ) + ret += np.outer(itrfpos[...,2]-phase_center[2], arg2 ) + + return ret + +def geoditic_to_ITRF(latLonAlt): + """for a latLonAlt, (can be list of three numpy arrays), convert to ITRF coordinates. Using information at: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#Geodetic_to/from_ENU_coordinates and + http://itrf.ensg.ign.fr/faq.php?type=answer""" + + a = 6378137.0 #### semi-major axis, m + e2 = 0.00669438002290 ## eccentricity squared + + def N(lat): + return a/np.sqrt( 1 - e2*(np.sin(lat)**2) ) + + lat = latLonAlt[0]/RTD + lon = latLonAlt[1]/RTD + + X = ( N(lat) + latLonAlt[2] ) *np.cos(lat) *np.cos(lon) + Y = ( N(lat) + latLonAlt[2] ) *np.cos(lat) *np.sin(lon) + + b2_a2 = 1-e2 + Z = ( b2_a2*N(lat) + latLonAlt[2] ) *np.sin(lat) + + return np.array( [X,Y,Z] ) + + + +#### previously known clock offsets. Only used for compatibility with past data #### +def getClockCorrectionFromParsetAddition(): + parsetFilename = MetaData_directory + '/lofar/station_clock_offsets/StationCalibration.parset' + + offsetDictX = {} + offsetDictY = {} + + infile = open(parsetFilename, 'r') + for line in infile: + s = line.split('=') + value = s[1] + params = s[0].split('.') + thisStation = params[2][0:5] + thisAntennaSet = params[3] + thisFilter = params[4] + thisValueType = params[5] + thisPolarization = params[6][0] + + if thisAntennaSet == 'LBA_OUTER' and thisFilter == 'LBA_30_90' and thisValueType == 'delay': + if thisPolarization == 'X': + offsetDictX[thisStation] = float(value) + elif thisPolarization == 'Y': + offsetDictY[thisStation] = float(value) + else: + raise ValueError('Wrong!') + infile.close() + + offsetDictCombined = {} + + for key in offsetDictX.keys(): + combined = 0.5 * (offsetDictX[key] + offsetDictY[key]) + offsetDictCombined[key] = combined + + return offsetDictCombined + + +def getClockCorrections( antennaset="LBA", time=1383264000+1000): + """Get clock correction for superterp stations in seconds. Currently static values. + + *station* Station name or number for which to get the correction. + *time* Optional. Linux time of observation. As clocks drift the value from the correct time should be given. Not yet implemented. + """ + + clockcorrection = {} + if "LBA" in antennaset: + if time < (1383264000): + # Values before 1 Nov 2013, eventID-time 120960000, Unix time: add 1262304000. + clockcorrection["CS002"] = 8.32233e-06 # definition, global offset + # Addition is the finetuning using Smilde from 1 or 2 random events, to about +/- 0.2 ns. Need to check constancy over time. + clockcorrection["CS003"] = 6.921444e-06 + 0.35e-9 + clockcorrection["CS004"] = 7.884847e-06 + 1.0e-9 + clockcorrection["CS005"] = 8.537828e-06 + 0.14e-9 + clockcorrection["CS006"] = 7.880705e-06 - 0.24e-9 + clockcorrection["CS007"] = 7.916458e-06 - 0.22e-9 + + clockcorrection["CS001"] = 4.755947e-06 + clockcorrection["CS011"] = 7.55500e-06 - 0.3e-9 + clockcorrection["CS013"] = 9.47910e-06 + clockcorrection["CS017"] = 1.540812e-05 - 0.87e-9 + clockcorrection["CS021"] = 6.044335e-06 + 1.12e-9 + clockcorrection["CS024"] = 4.66335e-06 - 1.24e-9 + clockcorrection["CS026"] = 1.620482e-05 - 1.88e-9 + clockcorrection["CS028"] = 1.6967048e-05 + 1.28e-9 + clockcorrection["CS030"] = 9.7110576e-06 + 3.9e-9 + clockcorrection["CS031"] = 6.375533e-06 + 1.87e-9 + clockcorrection["CS032"] = 8.541675e-06 + 1.1e-9 + clockcorrection["CS101"] = 1.5155471e-05 + clockcorrection["CS103"] = 3.5503206e-05 + clockcorrection["CS201"] = 1.745439e-05 + clockcorrection["CS301"] = 7.685249e-06 + clockcorrection["CS302"] = 1.2317004e-05 + clockcorrection["CS401"] = 8.052200e-06 + clockcorrection["CS501"] = 1.65797e-05 + else: + clockcorrection = getClockCorrectionFromParsetAddition() + clockcorrection["CS003"] = clockcorrection["CS003"] - 1.7e-9 + 2.0e-9 + clockcorrection["CS004"] = clockcorrection["CS004"] - 9.5e-9 + 4.2e-9 + clockcorrection["CS005"] = clockcorrection["CS005"] - 6.9e-9 + 0.4e-9 + clockcorrection["CS006"] = clockcorrection["CS006"] - 8.3e-9 + 3.8e-9 + clockcorrection["CS007"] = clockcorrection["CS007"] - 3.6e-9 + 3.4e-9 + clockcorrection["CS011"] = clockcorrection["CS011"] - 18.7e-9 + 0.6e-9 + +# Old values were + elif "HBA" in antennaset: + # Correct to 2013-03-26 values from parset L111421 + clockcorrection["CS001"] = 4.759754e-06 + clockcorrection["CS002"] = 8.318834e-06 + clockcorrection["CS003"] = 6.917926e-06 + clockcorrection["CS004"] = 7.889961e-06 + clockcorrection["CS005"] = 8.542093e-06 + clockcorrection["CS006"] = 7.882892e-06 + clockcorrection["CS007"] = 7.913020e-06 + clockcorrection["CS011"] = 7.55852e-06 + clockcorrection["CS013"] = 9.47910e-06 + clockcorrection["CS017"] = 1.541095e-05 + clockcorrection["CS021"] = 6.04963e-06 + clockcorrection["CS024"] = 4.65857e-06 + clockcorrection["CS026"] = 1.619948e-05 + clockcorrection["CS028"] = 1.6962571e-05 + clockcorrection["CS030"] = 9.7160576e-06 + clockcorrection["CS031"] = 6.370090e-06 + clockcorrection["CS032"] = 8.546255e-06 + clockcorrection["CS101"] = 1.5157971e-05 + clockcorrection["CS103"] = 3.5500922e-05 + clockcorrection["CS201"] = 1.744924e-05 + clockcorrection["CS301"] = 7.690431e-06 + clockcorrection["CS302"] = 1.2321604e-05 + clockcorrection["CS401"] = 8.057504e-06 + clockcorrection["CS501"] = 1.65842e-05 + + else: + print( "ERROR: no clock offsets available for this antennaset: ", antennaset) + return 0 + + return clockcorrection diff --git a/NuRadioReco/modules/io/lofar/raw_tbb_IO.py b/NuRadioReco/modules/io/lofar/raw_tbb_IO.py new file mode 100755 index 00000000..ee192cdc --- /dev/null +++ b/NuRadioReco/modules/io/lofar/raw_tbb_IO.py @@ -0,0 +1,756 @@ +#!/usr/bin/env python3 + +"""This module implements an interface for reading LOFAR TBB data. + +This module is strongly based on pyCRtools module tbb.py by Pim Schellart, Tobias Winchen, and others. +However, it has been completely re-written for use with LOFAR-LIM + +Author: Brian Hare + +Definitions: +LOFAR is split into a number of different stations. There are three main types: Core Stations (CS), Remote Stations (RS), and international stations +Each station contains 96 low band antennas (LBA) and 48 high band antennas (HBA). Each antenna is dual polarized. + +Each station is refered to by its name (e.g. "CS001"), which is a string, or its ID (e.g. 1), which is an integer. In general, these are different! +The mapping, however, is unique and is given in utilities.py + +There are a few complications with reading the data. +1) The data from each station is often spread over multiple files + There is a class below that can combine multiple files (even from different stations) +2) It is entirely possible that one file could contain multiple stations + This feature is not used, so I assume that it isn't a problem (for now) +3) Each Station has unknown clock offsets. Technically the core stations are all on one clock, but there are some unknown cable delays + This is a difficult problem, not handled here +4) Each Antenna doesn't necisarily start reading data at precisely the same time. + The code below picks the latest start time so this problem can be ignored by the end user +5) The software that inserts metadata (namely antenna positions and callibrations) sometimes "forgets" to do its job + The code below will automatically read the metadata from other files when necisary +6) LOFAR is constantly changing + So..keeping code up-to-date and still backwards compatible will be an interesting challange + +7) LOFAR only has 96 RCUs (reciever control units) per station (at the moment). + Each RCU is essentually one digitizer. Each antenna needs two RCS to record both polarizations. The result is only 1/3 of the antennas can + be read out each time. + + LOFAR keeps track of things with two ways. First, the data is all refered to by its RCUid. 0 is the 0th RCU, ect... However, which antenna + corresponds to which RCU depends on the antennaSet. For LOFAR-LIM the antenna set will generally be "LBA_OUTER". This could change, and sometimes + the antenna set is spelled wrong in the data files. (This should be handeld here though) + + In the code below each RCU is refered to by ANTENNA_NAME or antennaID. These are the same thing (I think). They are, however, a misnomer, as they + actually refer to the RCU, not antenna. The specific antenna depends on the antenna set. For the same antenna set, however, the ANTENNA_NAME will + always refer to the same antenna. + + Each ANTENNA_NAME is a string of 9 digits. First three is the station ID (not name!), next three is the group (no idea, don't ask), final 3 is the RCU id + + For LBA_INNER data set, even RCU ids refer to X-polarized dipoles and odd RCU ids refer to Y-polarized dipoles. This is flipped for LBA_OUTER antenna set. + X-polarization is NE-SW, and Y-polarization is NW-SE. antenna_responce.py, which handles the antenna function, assumes the data is LBA_OUTER. + +""" + +##### TODO: +## add a way to combine event that is spread across close timeID (is this necisary?) +## add proper fitting phase vs frequency to lines, adn func to return the frequency independand phase offset + + + + +import os + +import numpy as np +import h5py + +import NuRadioReco.modules.io.lofar.metadata as md +import NuRadioReco.modules.io.lofar.utilities as util + + +#nyquist_zone = {'LBA_10_90' : 1, 'LBA_30_90' : 1, 'HBA_110_190' : 2, 'HBA_170_230' : 3, 'HBA_210_250' : 3} +conversiondict = {"": 1.0, "kHz": 1000.0, "MHz": 10.0 ** 6, "GHz": 10.0 ** 9, "THz": 10.0 ** 12} + +#### helper functions #### + +def filePaths_by_stationName(timeID, raw_data_loc=None): + """Given a timeID, and a location of raw data (default set in utilities.py), return a dictionary. + The keys of the dictionary are antenna names, the values are lists of file paths to data files that contain that station.""" + + data_file_path = util.raw_data_dir(timeID, raw_data_loc) + h5_files = [f for f in os.listdir(data_file_path) if f[-6:] == 'tbb.h5'] + + ret = {} + for fname in h5_files: + Fpath = data_file_path + '/' + fname + junk, sname, junk, junk = util.Fname_data(Fpath) + + if sname not in ret: + ret[sname] = [] + + ret[sname].append( Fpath ) + + return ret + +def eventData_filePaths(timeID, raw_data_loc=None): + """Given a timeID, and a location of raw data (default set in utilities.py), return a list of file paths of data files""" + data_file_path = util.raw_data_dir(timeID, raw_data_loc) + return [f for f in os.listdir(data_file_path) if f[-6:] == 'tbb.h5'] + +######## The following four functions read what I call "correction files" these are corrections made to improve the data ########## +def read_antenna_pol_flips(fname): + antennas_to_flip = [] + with open(fname) as fin: + for line in fin: + ant_name = line.split()[0] + antennas_to_flip.append( ant_name ) + return antennas_to_flip + +def read_bad_antennas(fname): + bad_antenna_data = [] + with open(fname) as fin: + for line in fin: + ant_name, pol = line.split()[0:2] + bad_antenna_data.append((ant_name,int(pol))) + return bad_antenna_data + +def read_antenna_delays(fname): + additional_ant_delays = {} + with open(fname) as fin: + for line in fin: + ant_name, pol_E_delay, pol_O_delay = line.split()[0:3] + additional_ant_delays[ant_name] = [float(pol_E_delay), float(pol_O_delay)] + return additional_ant_delays + +def read_station_delays(fname): + station_delays = {} + with open(fname) as fin: + for line in fin: + sname, delay = line.split()[0:2] + station_delays[sname] = float(delay) + return station_delays + +#### data reading class #### +# Note: ASTRON will "soon" release a new DAL (data ) + +class TBBData_Dal1: + """a class for reading one station from one file. However, since one station is often spread between different files, + use filePaths_by_stationName combined with MultiFile_Dal1 below""" + + def __init__(self, filename, force_metadata_ant_pos=False, forcemetadata_delays=True): + self.filename = filename + self.force_metadata_ant_pos = force_metadata_ant_pos + self.forcemetadata_delays = forcemetadata_delays + + #### open file and set some basic info#### + self.file = h5py.File(filename, "r") + + stationKeys = [s for s in self.file.keys() if s.startswith('Station')] + ## assume there is only one station in the file + + if len(stationKeys) != 1: + print("WARNING! file", self.filename, "has more then one station") + self.stationKey = stationKeys[0] + + self.antennaSet = self.file.attrs['ANTENNA_SET'][0].decode() + self.dipoleNames = list( self.file[ self.stationKey ].keys() ) + self.StationID = self.file[ self.stationKey ][ self.dipoleNames[0] ].attrs["STATION_ID"][0] + self.StationName = util.SId_to_Sname[ self.StationID ] + ## assume all antennas have the same sample frequency + + self.SampleFrequency = self.file[ self.stationKey ][ self.dipoleNames[0] ].attrs["SAMPLE_FREQUENCY_VALUE" + ][0]*conversiondict[ self.file[ self.stationKey ][ self.dipoleNames[0] ].attrs["SAMPLE_FREQUENCY_UNIT"][0].decode() ] + ## filter selection is typically "LBA_10_90" + self.FilterSelection = self.file.attrs['FILTER_SELECTION'][0].decode() + + + #### check that all antennas start in the same second, and record the same number of samples #### + self.Time = None + self.DataLengths = np.zeros(len(self.dipoleNames), dtype=int) + self.SampleNumbers = np.zeros(len(self.dipoleNames), dtype=int) + for dipole_i, dipole in enumerate(self.dipoleNames): + + if self.Time is None: + self.Time = self.file[ self.stationKey ][ dipole ].attrs["TIME"][0] + else: + if self.Time != self.file[ self.stationKey ][ dipole ].attrs["TIME"][0]: + raise IOError("antennas do not start at same time in "+self.filename) + + self.DataLengths[dipole_i] = self.file[ self.stationKey ][ dipole ].attrs["DATA_LENGTH"][0] + self.SampleNumbers[dipole_i] = self.file[ self.stationKey ][ dipole ].attrs["SAMPLE_NUMBER"][0] + + + #### get position and delay metadata...maybe#### + self.have_metadata = 'DIPOLE_CALIBRATION_DELAY_VALUE' in self.file[ self.stationKey ][ self.dipoleNames[0] ].attrs + self.antenna_filter = md.make_antennaID_filter(self.dipoleNames) + + + # load antenna locations from metadata and from file. IF they are too far apart, then give warning, and use metadata + self.ITRF_dipole_positions = md.getItrfAntennaPosition(self.StationName, self.antennaSet)[ self.antenna_filter ] ## load positions from metadata file + if self.have_metadata and not self.force_metadata_ant_pos: + + use_TBB_positions = True + TBB_ITRF_dipole_positions = np.empty((len(self.dipoleNames), 3), dtype=np.double) + for i,dipole in enumerate(self.dipoleNames): + TBB_ITRF_dipole_positions[i] = self.file[ self.stationKey ][dipole].attrs['ANTENNA_POSITION_VALUE'] + + dif = np.linalg.norm( TBB_ITRF_dipole_positions[i]-self.ITRF_dipole_positions[i] ) + if dif > 1: + print("WARNING: station", self.StationName, "has suspicious antenna locations. Using metadata instead") + use_TBB_positions = False + + if use_TBB_positions: + self.ITRF_dipole_positions = TBB_ITRF_dipole_positions + + + + if self.have_metadata and not self.forcemetadata_delays: + self.calibrationDelays = np.empty( len(self.dipoleNames), dtype=np.double ) + + for i,dipole in enumerate(self.dipoleNames): + self.calibrationDelays[i] = self.file[ self.stationKey ][dipole].attrs['DIPOLE_CALIBRATION_DELAY_VALUE'] + + + + #### get the offset, in number of samples, needed so that each antenna starts at the same time #### + self.nominal_sample_number = np.max( self.SampleNumbers ) + self.sample_offsets = self.nominal_sample_number - self.SampleNumbers + self.nominal_DataLengths = self.DataLengths - self.sample_offsets + + #### GETTERS #### + def needs_metadata(self): + """return true if this file does not have metadata""" + return not self.have_metadata + + def get_station_name(self): + """returns the name of the station, as a string""" + return self.StationName + + def get_station_ID(self): + """returns the ID of the station, as an integer. This is not the same as StationName. Mapping is givin in utilities""" + return self.StationID + + def get_antenna_names(self): + """return name of antenna as a list of strings. This is really the RCU id, and the physical antenna depends on the antennaSet""" + return self.dipoleNames + + def get_antenna_set(self): + """return the antenna set as a string. Typically "LBA_OUTER" """ + return self.antennaSet + + def get_sample_frequency(self): + """gets samples per second. Typically 200 MHz.""" + return self.SampleFrequency + + def get_filter_selection(self): + """return a string that represents the frequency filter used. Typically "LBA_10_90""" + return self.FilterSelection + + def get_timestamp(self): + """return the POSIX timestamp of the first data point""" + return self.Time + + def get_full_data_lengths(self): + """get the number of samples stored for each antenna. Note that due to the fact that the antennas do not start recording + at the exact same instant (in general), this full data length is not all usable + returns array of ints""" + return self.DataLengths + + def get_all_sample_numbers(self): + """return numpy array that contains the sample numbers of each antenna. Divide this by the sample frequency to get time + since the timestame of the first data point. Note that since these are, in general, different, they do NOT refer to sample + 0 of "get_data" in general """ + return self.SampleNumbers + + def get_nominal_sample_number(self): + """return the sample number of the 0th data sample returned by get_data. + Divide by sample_frequency to get time from timestamp of the 0th data sample""" + return self.nominal_sample_number + + def get_nominal_data_lengths(self): + """return the number of data samples that are usable for each antenna, accounting for different starting sample numbers. + returns array of ints""" + return self.nominal_DataLengths + + def get_ITRF_antenna_positions(self, copy=False): + """returns the ITRF positions of the antennas. Returns a 2D numpy array. If copy is False, then this just returns the internal array of values""" + if copy: + return np.array( self.ITRF_dipole_positions ) + else: + return self.ITRF_dipole_positions + + def get_LOFAR_centered_positions(self, out=None): + """returns the positions (as a 2D numpy array) of the antennas with respect to CS002. + if out is a numpy array, it is used to store the antenna positions, otherwise a new array is allocated""" + return md.convertITRFToLocal(self.ITRF_dipole_positions, out=out) + + def get_timing_callibration_phases(self): + """only a test function for the moment, do not use""" + fpath = os.path.dirname(self.filename) + '/'+self.StationName + phase_calibration = md.getStationPhaseCalibration(self.StationName, self.antennaSet, file_location=fpath ) + phase_calibration = phase_calibration[ self.antenna_filter ] + return phase_calibration + + def get_timing_callibration_delays(self): + """return the timing callibration of the anntennas, as a 1D np array. If not included in the metadata, will look + for a data file in the same directory as this file. Otherwise returns None""" + + if self.have_metadata and not self.forcemetadata_delays: + return self.calibrationDelays + else: + fpath = os.path.dirname(self.filename) + '/'+self.StationName + phase_calibration = md.getStationPhaseCalibration(self.StationName, self.antennaSet, file_location=fpath ) + phase_calibration = phase_calibration[ self.antenna_filter ] + return md.convertPhase_to_Timing(phase_calibration, 1.0/self.SampleFrequency) + + def get_data(self, start_index, num_points, antenna_index=None, antenna_ID=None): + """return the raw data for a specific antenna, as an 1D int16 numpy array, of length num_points. First point returned is + start_index past get_nominal_sample_number(). Specify the antenna by giving the antenna_ID (which is a string, same + as output from get_antenna_names(), or as an integer antenna_index. An antenna_index of 0 is the first antenna in + get_antenna_names().""" + + if antenna_index is None: + if antenna_ID is None: + raise LookupError("need either antenna_ID or antenna_index") + antenna_index = self.dipoleNames.index(antenna_ID) + else: + antenna_ID = self.dipoleNames[ antenna_index ] + + initial_point = self.sample_offsets[ antenna_index ] + start_index + final_point = initial_point+num_points + + return self.file[ self.stationKey ][ antenna_ID ][initial_point:final_point] + +class MultiFile_Dal1: + """A class for reading the data from one station from multiple files""" + def __init__(self, filename_list, force_metadata_ant_pos=False, polarization_flips=None, bad_antennas=[], additional_ant_delays=None, station_delay=0.0, only_complete_pairs=True): + """filename_list: list of filenames for this station for this event. + force_metadata_ant_pos -if True, then load antenna positions from a metadata file and not the raw data file. Default False + polarization_flips -list of even antennas where it is known that even and odd antenna names are flipped in file. This is assumed to apply both to data and timing calibration + bad_antennas -list of antennas that should not be used. Each item in the list is a tuple, first item of tuple is name of even antenna, second item is a 0 or 1 indicating if even or odd antenna is bad. + assumed to be BEFORE antenna flips are accounted for + additional_ant_delays -a dictionary. Each key is name of even antenna, each value is a tuple with additional even and odd antenna delays. This should rarely be needed. + assumed to be found BEFORE antenna flips are accounted for + station_delay -a single number that represents the clock offset of this station, as a delay + NOTE: polarization_flips, bad_antennas, additional_ant_delays, and station_delay can now be strings that are file names. If this is the case, they will be read automatically + only_complete_pairs -if True, discards antenna if the other in pair is not present or is bad. If False, keeps all good antennas with a 'none' value if other antenna in pair is missing""" + + self.files = [TBBData_Dal1(fname, force_metadata_ant_pos) for fname in filename_list] + + if isinstance(polarization_flips, str): + polarization_flips = read_antenna_pol_flips( polarization_flips ) + if isinstance(bad_antennas, str): + bad_antennas = read_bad_antennas( read_bad_antennas ) + if isinstance(additional_ant_delays, str): + additional_ant_delays = read_antenna_delays( additional_ant_delays ) + + #### get some data that should be constant #### TODO: change code to make use of getters + self.antennaSet = self.files[0].antennaSet + self.StationID = self.files[0].StationID + self.StationName = self.files[0].StationName + self.SampleFrequency = self.files[0].SampleFrequency + self.FilterSelection = self.files[0].FilterSelection + self.Time = self.files[0].Time + self.bad_antennas = bad_antennas + self.odd_pol_additional_timing_delay = 0.0 # anouther timing delay to add to all odd-polarized antennas + + if isinstance(station_delay, str): + station_delay = read_station_delays( station_delay )[ self.StationName ] + + self.station_delay = station_delay + + #### check consistancy of data #### + for TBB_file in self.files: + if TBB_file.antennaSet != self.antennaSet: + raise IOError("antenna set not the same between files for station: "+self.StationName) + if TBB_file.StationID != self.StationID: + raise IOError("station ID not the same between files for station: "+self.StationName) + if TBB_file.StationName != self.StationName: + raise IOError("station name not the same between files for station: "+self.StationName) + if TBB_file.FilterSelection != self.FilterSelection: + raise IOError("filter selection not the same between files for station: "+self.StationName) + if TBB_file.Time != self.Time: + raise IOError("antenna set not the same between files for station: "+self.StationName) + + ## check LBA outer antenna set + if self.antennaSet != "LBA_OUTER": + print("WARNING: antenna set on station", self.StationName, "is not LBA_OUTER") + + + #### find best files to get antennas from #### + ## require each antenna shows up once, and even pol is followed by odd pol + + self.dipoleNames = [] + self.antenna_to_file = [] ##each item is a tuple. First item is file object, second is antenna index in file + + unused_antenna_names = [] + unused_antenna_to_file = [] + bad_PolE_antennas = [ant for ant,pol in bad_antennas if pol==0] + bad_PolO_antennas = [ant for ant,pol in bad_antennas if pol==1] ## note that this is still the name of the even antenna, although it is the ODD antenna that is bad!!! + for TBB_file in self.files: + file_ant_names = TBB_file.get_antenna_names() + + for ant_i,ant_name in enumerate(file_ant_names): + if (ant_name in self.dipoleNames): + continue + ant_ID = int(ant_name[-3:]) + + if ant_ID%2 == 0: ##check if antenna is even + if ant_name in bad_PolE_antennas: + continue + + odd_ant_name = ant_name[:-3] + str(ant_ID+1).zfill(3) + if odd_ant_name in unused_antenna_names: ## we have the odd antenna + self.dipoleNames.append(ant_name) + self.dipoleNames.append(odd_ant_name) + + self.antenna_to_file.append( (TBB_file, ant_i) ) + odd_unused_index = unused_antenna_names.index( odd_ant_name ) + self.antenna_to_file.append( unused_antenna_to_file[ odd_unused_index ] ) + + unused_antenna_names.pop( odd_unused_index ) + unused_antenna_to_file.pop( odd_unused_index ) + else: ## we haven't found the odd antenna, so store info for now + unused_antenna_names.append(ant_name) + unused_antenna_to_file.append( (TBB_file, ant_i) ) + + else: ## antenna is odd + even_ant_name = ant_name[:-3] + str(ant_ID-1).zfill(3) + if even_ant_name in bad_PolO_antennas: ## note that have to check if EVEN antenna is in bad antenna names... + continue + + if even_ant_name in unused_antenna_names: ## we have the odd antenna + self.dipoleNames.append(even_ant_name) + self.dipoleNames.append(ant_name) + + even_unused_index = unused_antenna_names.index( even_ant_name ) + self.antenna_to_file.append( unused_antenna_to_file[ even_unused_index ] ) + + unused_antenna_names.pop( even_unused_index ) + unused_antenna_to_file.pop( even_unused_index ) + + self.antenna_to_file.append( (TBB_file, ant_i) ) + + else: ## we haven't found the odd antenna, so store info for now + unused_antenna_names.append(ant_name) + unused_antenna_to_file.append( (TBB_file, ant_i) ) + + if not only_complete_pairs: + for ant_name, to_file in zip(unused_antenna_names, unused_antenna_to_file): + ant_ID = int(ant_name[-3:]) + if ant_ID%2 == 0: ##check if antenna is even + + self.dipoleNames.append( ant_name ) + self.antenna_to_file.append( to_file ) + + self.dipoleNames.append( ant_name[:-3] + str(ant_ID+1).zfill(3) ) ## add the odd antenna + self.antenna_to_file.append( None ) ## doesn't exist in a file + + else: + + self.dipoleNames.append( ant_name[:-3] + str(ant_ID-1).zfill(3) ) ## add the even antenna + self.antenna_to_file.append( None ) ## doesn't exist in a file + + self.dipoleNames.append( ant_name ) + self.antenna_to_file.append( to_file ) + + + self.index_adjusts = np.arange( len(self.antenna_to_file) ) ##used to compensate for polarization flips + ### when given an antnna index to open data, use this index instead to open the correct data location + + + #### get sample numbers and offsets and lengths and other related stuff #### + self.SampleNumbers = [] + self.DataLengths = [] + for TBB_file, file_ant_i in self.antenna_to_file: + self.SampleNumbers.append( TBB_file.SampleNumbers[file_ant_i] ) + self.DataLengths.append( TBB_file.DataLengths[file_ant_i] ) + + self.SampleNumbers = np.array( self.SampleNumbers, dtype=int ) + self.DataLengths = np.array(self.DataLengths, dtype=int) + + self.nominal_sample_number = np.max( self.SampleNumbers ) + self.sample_offsets = self.nominal_sample_number - self.SampleNumbers + self.nominal_DataLengths = self.DataLengths - self.sample_offsets + + self.even_ant_pol_flips = None + if polarization_flips is not None: + self.set_polarization_flips( polarization_flips ) + self.additional_ant_delays = additional_ant_delays + + def set_polarization_flips(self, even_antenna_names): + """given a set of names(IDs) of even antennas, flip the data between the even and odd antennas""" + self.even_ant_pol_flips = even_antenna_names + for ant_name in even_antenna_names: + if ant_name in self.dipoleNames: + even_antenna_index = self.dipoleNames.index(ant_name) + + self.index_adjusts[even_antenna_index] += 1 + self.index_adjusts[even_antenna_index+1] -= 1 + + def set_odd_polarization_delay(self, new_delay): + self.odd_pol_additional_timing_delay = new_delay + + def set_station_delay(self, station_delay): + """ set the station delay, should be a number""" + self.station_delay = station_delay + + def find_and_set_polarization_delay(self, verbose=False, tolerance=1e-9): + fpath = os.path.dirname(self.files[0].filename) + '/'+self.StationName + phase_calibration = md.getStationPhaseCalibration(self.StationName, self.antennaSet, file_location=fpath ) + all_antenna_calibrations = md.convertPhase_to_Timing(phase_calibration, 1.0/self.SampleFrequency) + + even_delays = all_antenna_calibrations[::2] + odd_delays = all_antenna_calibrations[1::2] + odd_offset = odd_delays-even_delays + median_odd_offset = np.median( odd_offset ) + if verbose: + print("median offset is:", median_odd_offset) + below_tolerance = np.abs( odd_offset-median_odd_offset ) < tolerance + if verbose: + print(np.sum(below_tolerance), "antennas below tolerance.", len(below_tolerance)-np.sum(below_tolerance), "above.") + ave_best_offset = np.average( odd_offset[below_tolerance] ) + if verbose: + print("average of below-tolerance offset is:", ave_best_offset) + self.set_odd_polarization_delay( -ave_best_offset ) + + above_tolerance = np.zeros( len(all_antenna_calibrations), dtype=bool ) + above_tolerance[::2] = np.logical_not( below_tolerance ) + above_tolerance[1::2] = above_tolerance[::2] + above_tolerance = above_tolerance[ md.make_antennaID_filter(self.get_antenna_names()) ] + return [AN for AN, AT in zip(self.get_antenna_names(),above_tolerance) if AT] + + + #### GETTERS #### + def needs_metadata(self): + for TBB_file in self.files: + if TBB_file.needs_metadata(): + return True + return False + + def get_station_name(self): + """returns the name of the station, as a string""" + return self.StationName + + def get_station_ID(self): + """returns the ID of the station, as an integer. This is not the same as StationName. Mapping is givin in utilities""" + return self.StationID + + def get_antenna_names(self): + """return name of antenna as a list of strings. This is really the RCU id, and the physical antenna depends on the antennaSet""" + return self.dipoleNames + + def has_antenna(self, antenna_name): + """if only_complete_pairs is False, then we could have antenna names without the data. Return True if we actually have the antenna, False otherwise. Account for polarization flips.""" + if antenna_name in self.dipoleNames: + index = self.index_adjusts( self.dipoleNames.index(antenna_name) ) + if self.antenna_to_file[index] is None: + return False + else: + return True + else: + return False + + def get_antenna_set(self): + """return the antenna set as a string. Typically "LBA_OUTER" """ + return self.antennaSet + + def get_sample_frequency(self): + """gets samples per second. Typically 200 MHz.""" + return self.SampleFrequency + + def get_filter_selection(self): + """return a string that represents the frequency filter used. Typically "LBA_10_90""" + return self.FilterSelection + + def get_timestamp(self): + """return the POSIX timestamp of the first data point""" + return self.Time + + def get_full_data_lengths(self): + """get the number of samples stored for each antenna. Note that due to the fact that the antennas do not start recording + at the exact same instant (in general), this full data length is not all usable + returns array of ints""" + return self.DataLengths + + def get_all_sample_numbers(self): + """return numpy array that contains the sample numbers of each antenna. Divide this by the sample frequency to get time + since the timestame of the first data point. Note that since these are, in general, different, they do NOT refer to sample + 0 of "get_data" """ + return self.SampleNumbers + + def get_nominal_sample_number(self): + """return the sample number of the 0th data sample returned by get_data. + Divide by sample_frequency to get time from timestamp of the 0th data sample""" + return self.nominal_sample_number + + def get_nominal_data_lengths(self): + """return the number of data samples that are usable for each antenna, accounting for different starting sample numbers. + returns array of ints""" + return self.nominal_DataLengths + + def get_ITRF_antenna_positions(self, out=None): + """returns the ITRF positions of the antennas. Returns a 2D numpy array. + if out is a numpy array, it is used to store the antenna positions, otherwise a new array is allocated. + Does not account for polarization flips, but shouldn't need too.""" + if out is None: + out = np.empty( (len(self.dipoleNames), 3) ) + + for ant_i, (TBB_file,station_ant_i) in enumerate(self.antenna_to_file): + out[ant_i] = TBB_file.ITRF_dipole_positions[station_ant_i] + + return out + + def get_LOFAR_centered_positions(self, out=None): + """returns the positions (as a 2D numpy array) of the antennas with respect to CS002. + if out is a numpy array, it is used to store the antenna positions, otherwise a new array is allocated. + Does not account for polarization flips, but shouldn't need too.""" + if out is None: + out = np.empty( (len(self.dipoleNames), 3) ) + + md.convertITRFToLocal( self.get_ITRF_antenna_positions(), out=out ) + + return out + + def get_timing_callibration_phases(self): + """only a test function for the moment, do not use""" + + out = [None for i in range(len(self.dipoleNames))] + + for TBB_file in self.files: + ret = TBB_file.get_timing_callibration_phases() + if ret is None: + return None + + for ant_i, (TBB_fileA,station_ant_i) in enumerate(self.antenna_to_file): + if TBB_fileA is TBB_file: + out[ant_i] = ret[station_ant_i] + + return np.array(out) + + def get_timing_callibration_delays(self, out=None): + """return the timing callibration of the anntennas, as a 1D np array. If not included in the metadata, will look + for a data file in the same directory as this file. Otherwise returns None. + if out is a numpy array, it is used to store the antenna delays, otherwise a new array is allocated. + This takes polarization flips, and additional_ant_delays into account (assuming that both were found BEFORE the pol flip was found). + Also can account for a timing difference between even and odd antennas, if it is set.""" + + if out is None: + out = np.zeros( len(self.dipoleNames) ) + + for TBB_file in self.files: + ret = TBB_file.get_timing_callibration_delays() + if ret is None: + return None + + + for ant_i, adjust_i in enumerate(self.index_adjusts): + TBB_fileA,station_ant_i = self.antenna_to_file[adjust_i] + + if TBB_fileA is TBB_file: + out[ant_i] = ret[station_ant_i] + + if self.additional_ant_delays is not None: + ## additional_ant_delays stores only even antenna names for historical reasons. so we need to be clever here + antenna_polarization = 0 if (ant_i%2==0) else 1 + even_ant_name = self.dipoleNames[ ant_i-antenna_polarization ] + if even_ant_name in self.additional_ant_delays: + if even_ant_name in self.even_ant_pol_flips: + antenna_polarization = int(not antenna_polarization) + out[ant_i] += self.additional_ant_delays[ even_ant_name ][ antenna_polarization ] ## 90% sure this is correct sign + + out[1::2] += self.odd_pol_additional_timing_delay + + return out + + def get_total_delays(self, out=None): + """Return the total delay for each antenna, accounting for all antenna delays, polarization delay, station clock offsets, and trigger time offsets (nominal sample number). + This function should be prefered over 'get_timing_callibration_delays', but the offsets can have a large average. It is recomended to pick one antenna (on your referance station) + and use it as a referance antenna so that it has zero timing delay. Note: this creates two defintions of T=0. I will call 'uncorrected time' is when the result of this function is + used as-is, and a referance antenna is not choosen. (IE, the referance station can have a large total_delay offset), 'corrected time' will be otherwise.""" + + delays = self.get_timing_callibration_delays(out) + delays += self.station_delay - self.get_nominal_sample_number()*5.0E-9 + + return delays + + def get_geometric_delays(self, source_location, out=None, antenna_locations=None): + """Calculate travel time from a XYZ location to each antenna. out can be an array of length equal to number of antennas. + antenna_locations is the table of antenna locations, given by get_LOFAR_centered_positions(). If None, it is calculated. Note that antenna_locations CAN be modified in this function. + If antenna_locations is less then all antennas, then the returned array will be correspondingly shorter. + The output of this function plus get_total_delays plus emission time of the source is the time the source is seen on each antenna.""" + + if antenna_locations is None: + antenna_locations = self.get_LOFAR_centered_positions() + + if out is None: + out = np.empty( len(antenna_locations), dtype=np.double ) + + if len(out) != len(antenna_locations): + print("ERROR: arrays are not of same length in geometric_delays()") + return None + + antenna_locations -= source_location + antenna_locations *= antenna_locations + np.sum(antenna_locations, axis=1, out=out) + np.sqrt(out, out=out) + out /= util.v_air + return out + + def get_data(self, start_index, num_points, antenna_index=None, antenna_ID=None): + """return the raw data for a specific antenna, as an 1D int16 numpy array, of length num_points. First point returned is + start_index past get_nominal_sample_number(). Specify the antenna by giving the antenna_ID (which is a string, same + as output from get_antenna_names(), or as an integer antenna_index. An antenna_index of 0 is the first antenna in + get_antenna_names().""" + + if antenna_index is None: + if antenna_ID is None: + raise LookupError("need either antenna_ID or antenna_index") + antenna_index = self.dipoleNames.index(antenna_ID) + + antenna_index = self.index_adjusts[antenna_index] ##incase of polarization flips + + initial_point = self.sample_offsets[ antenna_index ] + start_index + final_point = initial_point+num_points + + to_file = self.antenna_to_file[antenna_index] + if to_file is None: + raise LookupError("do not have data for this antenna") + TBB_file,station_antenna_index = to_file + antenna_ID = self.dipoleNames[ antenna_index ] + + if final_point >= len( TBB_file.file[ TBB_file.stationKey ][ antenna_ID ] ): + print("WARNING! data point", final_point, "is off end of file", len( TBB_file.file[ TBB_file.stationKey ][ antenna_ID ] )) + + return TBB_file.file[ TBB_file.stationKey ][ antenna_ID ][initial_point:final_point] + + +if __name__ == "__main__": + import matplotlib.pyplot as plt + +# timeID = "D20170929T202255.000Z" +# station = "RS406" + antenna_id = 0 + + block_size = 2**16 + block_number = 2#3900 + +# raw_fpaths = filePaths_by_stationName(timeID) +# +# infile = MultiFile_Dal1(raw_fpaths[station]) + infile = MultiFile_Dal1(["/Users/anelles/Experiments/LOFAR/A_NuRadioReco_Test/L78862_D20121205T051644.039Z_CS002_R000_tbb.h5"]) + + print( infile.get_LOFAR_centered_positions() ) + + data = infile.get_data(block_number*block_size, block_size, antenna_index=antenna_id) + + plt.plot(data) + plt.show() +# + + + + + + + + + + + + diff --git a/NuRadioReco/modules/io/lofar/readLOFAR.py b/NuRadioReco/modules/io/lofar/readLOFAR.py new file mode 100644 index 00000000..2b33510a --- /dev/null +++ b/NuRadioReco/modules/io/lofar/readLOFAR.py @@ -0,0 +1,148 @@ +import sys +import numpy as np +import matplotlib.pyplot as plt + +#LOFAR library +import NuRadioReco.modules.io.lofar.raw_tbb_IO as IO +import NuRadioReco.modules.io.lofar.metadata as md + +import NuRadioReco.framework.event +import NuRadioReco.framework.station +import NuRadioReco.framework.channel + +import logging +logger = logging.getLogger("readLOFAR") +logging.basicConfig() + + +def to_lofar_id(station_id): + """ + Convert NuRadioReco event id (int) to LOFAR str event id + + Parameters: + ---------- + station_id: int + NuRadioReco ID + """ + id = str(station_id)[-3:] + if station_id > 3000: + logger.error("Station ID {} too high, not implemented".format(station_id)) + elif station_id > 2000: + LOFAR_id = 'RS'+id + else: + LOFAR_id = 'CS' + id + return LOFAR_id + + +def from_lofar_id(LOFAR_id): + """ + Convert LOFAR str event id to NuRadioReco event id + + Parameters: + ---------- + station_id: str + LOFAR id, in form of CS002 etc. + + """ + if 'RS' in LOFAR_id: + station_id = 2000 + station_id += int(LOFAR_id[-3]) + elif "CS" in LOFAR_id: + station_id = 1000 + station_id += int(LOFAR_id[-3:]) + else: + logger.error("Station type {} not implemented".format(LOFAR_id)) + return station_id + +def get_event_id_from_time(timestamp): + """ + Get LOFAR standard event ID + + Parameters: + ------------ + timestamp: int + unix timestamp of event + + """ + timestamp_0 = 1262304000 # Unix timestamp on Jan 1, 2010 (date -u --date "Jan 1, 2010 00:00:00" +"%s") + result = int(timestamp - timestamp_0) + + return result + + +class readLOFAR: + + """ + Reads LOFAR data + """ + + + def __init__(self,filenames = [], DAL=1): + """ + Init read LOFAR class + + ------------- + Parameters: + filenames: list + list of filenames belonging to one event + will be sorted by station + DAL: int + 1: data access library version of LOFAR + currently implemented 1 + default: 1 + + """ + self._DAL = DAL + if DAL not in [1]: + logger.error("Version of DAL {} not implemented".format(DAL)) + sys.exit() + + filenames_dict = {} + for filename in filenames: + if self._DAL == 1: + file = IO.TBBData_Dal1(filename) + LOFAR_id = (file.get_station_name()) + if LOFAR_id in filenames_dict.keys(): + filenames_dict[LOFAR_id].append(filename) + else: + filenames_dict[LOFAR_id] = [filename] + + + self.filenames = filenames_dict + + + def run(self): + for station in self.filenames.keys(): + if self._DAL == 1: + file_object = IO.MultiFile_Dal1(self.filenames[station]) + else: + logger.error("Version of DAL {} not implemented".format(self._DAL)) + + itrf_positions = file_object.get_ITRF_antenna_positions() + positions = md.convertITRFToLocal(itrf_positions) + + antenna_set = file_object.get_antenna_set() + if 'LBA' in antenna_set: + antenna_model = 'LBA' + elif 'HBA' in antenna_set: + antenna_model = 'HBA' + else: + logger.warning("Antenna set {} not recognized.".format(antenna_set)) + + evt_number = get_event_id_from_time(file_object.get_timestamp()) + + # LOFAR does not operate with run numbers + evt = NuRadioReco.framework.event.Event(1, evt_number) + + yield evt + + def end(self): + pass + + + +if __name__ == "__main__": + + f = readLOFAR(['/Users/anelles/Experiments/LOFAR/A_NuRadioReco_Test/L78862_D20121205T051644.039Z_CS002_R000_tbb.h5'],DAL=1) + + f.run() \ No newline at end of file diff --git a/NuRadioReco/modules/io/lofar/utilities.py b/NuRadioReco/modules/io/lofar/utilities.py new file mode 100644 index 00000000..0684b7d2 --- /dev/null +++ b/NuRadioReco/modules/io/lofar/utilities.py @@ -0,0 +1,317 @@ +#!/usr/bin/env python3 + +##ON APP MACHINE + +import sys + +from os import listdir, mkdir +from os.path import isdir, dirname, abspath + +import weakref + +from scipy import fftpack +import numpy as np + +## some global variables, this needs to be fixed at some point +default_raw_data_loc = None#"/exp_app2/appexp1/public/raw_data" +default_processed_data_loc = None#"/home/brian/processed_files" + +MetaData_directory = dirname(abspath(__file__)) + '/data' ## change this if antenna_response_model is in a folder different from this module + +#### constants +C = 299792458.0 +RTD = 180.0/3.1415926 ##radians to degrees +n_air = 1.000293 +v_air = C/n_air + +latlonCS002 = np.array([52.91512249, 6.869837540]) ## lattitude and longitude of CS002 in degrees + +#### log data to screen and to a file + +class logger(object): + class std_writer(object): + def __init__(self, logger): + self.logger_ref = weakref.ref(logger) + + def write(self, msg): + logger=self.logger_ref() + logger.out_file.write(msg) + if logger.to_screen: + logger.old_stdout.write(msg) + + def flush(self): + logger=self.logger_ref() + logger.out_file.flush() + + + def __init__(self): + + self.has_stderr = False + self.has_stdout = False + + self.old_stderr = sys.stderr + self.old_stdout = sys.stdout + + self.set("out_log") + + def set(self, fname, to_screen=True): + self.out_file = open(fname, 'w') + + self.set_to_screen( to_screen ) + + + def __call__(self, *args): + for a in args: + if self.to_screen: + self.old_stdout.write(str(a)) + self.old_stdout.write(" ") + + self.out_file.write(str(a)) + self.out_file.write(" ") + + self.out_file.write("\n") + if self.to_screen: + self.old_stdout.write("\n") + + self.out_file.flush() + self.old_stdout.flush() + + def set_to_screen(self, to_screen=True): + self.to_screen = to_screen + + def take_stdout(self): + + if not self.has_stdout: + sys.stdout = self.std_writer(self) + self.has_stdout = True + + def take_stderr(self): + + if not self.has_stderr: + sys.stderr = self.std_writer(self) + self.has_stderr = True + + def restore_stdout(self): + if self.has_stdout: + sys.stdout = self.old_stdout + self.has_stdout = False + + def restore_stderr(self): + if self.has_stderr: + sys.stderr = self.old_stderr + self.has_stderr = False + + def flush(self): + self.out_file.flush() + +# def __del__(self): +# self.restore_stderr() +# self.restore_stdout() + +#log = logger() + +def iterate_pairs(list_one, list_two, list_one_avoid=[], list_two_avoid=[]): + """returns an iterator that loops over all pairs of the two lists""" + for item_one in list_one: + if item_one in list_one_avoid: + continue + for item_two in list_two: + if item_two in list_two_avoid: + continue + yield (item_one, item_two) + + + +#### some file utils + +def Fname_data(Fpath): + """ takes both pulse data file names and h5 file names and returns UTC_time, station_name, Fpath""" + Fname = Fpath.split('/')[-1] + data = Fname.split('_') + timeID = data[1] + station_name = data[2] + + if len(data[3][1:])==0: + file_number = 0 + else: + file_number = int(data[3][1:]) + + return timeID, station_name, Fpath, file_number + + +##note that timeID is a string representing the datetime of a LOFAR trigger. such as: D20130619T094846.507Z +## the timeID is used to uniquely identify triggers + +def get_timeID(fname): + data=fname.split("_") + return data[1] + +def year_from_timeID(timeID): + return timeID[1:5] + +def raw_data_dir(timeID, data_loc=None): + """gives path to the raw data folder for a particular timeID, given location of data structure. Defaults to default_raw_data_loc""" + + if data_loc is None: + data_loc = default_raw_data_loc + + path = data_loc + '/' + year_from_timeID(timeID)+"/"+timeID + return path + +def processed_data_dir(timeID, data_loc=None): + """gives path to the analysis folders for a particular timeID, given location of data structure. Defaults to default_processed_data_loc + makes the directory if it doesn't exist""" + + if data_loc is None: + data_loc = default_processed_data_loc + + path=data_loc + "/" + year_from_timeID(timeID)+"/"+timeID + if not isdir(path): + mkdir(path) + return path + + +## a python list where the keys are the number of a station and the values are the station name +SId_to_Sname = [None]*209 #just to pre-initilize list, so syntax below is possible +SId_to_Sname[1] = "CS001" +SId_to_Sname[2] = "CS002" +SId_to_Sname[3] = "CS003" +SId_to_Sname[4] = "CS004" +SId_to_Sname[5] = "CS005" +SId_to_Sname[6] = "CS006" +SId_to_Sname[7] = "CS007" +#SId_to_Sname[8] = "CS008" +#SId_to_Sname[9] = "CS009" +#SId_to_Sname[10] = "CS010" +SId_to_Sname[11] = "CS011" +#SId_to_Sname[12] = "CS012" +SId_to_Sname[13] = "CS013" +#SId_to_Sname[14] = "CS014" +#SId_to_Sname[15] = "CS015" +#SId_to_Sname[16] = "CS016" +SId_to_Sname[17] = "CS017" +#SId_to_Sname[18] = "CS018" +#SId_to_Sname[19] = "CS019" +#SId_to_Sname[20] = "CS020" +SId_to_Sname[21] = "CS021" +#SId_to_Sname[22] = "CS022" +#SId_to_Sname[23] = "CS023" +SId_to_Sname[24] = "CS024" +#SId_to_Sname[25] = "CS025" +SId_to_Sname[26] = "CS026" +#SId_to_Sname[27] = "CS027" +SId_to_Sname[28] = "CS028" +#SId_to_Sname[29] = "CS029" +SId_to_Sname[30] = "CS030" +SId_to_Sname[31] = "CS031" +SId_to_Sname[32] = "CS032" +SId_to_Sname[101] = "CS101" +#SId_to_Sname[102] = "CS102" +SId_to_Sname[103] = "CS103" +SId_to_Sname[121] = "CS201" +SId_to_Sname[141] = "CS301" +SId_to_Sname[142] = "CS302" +SId_to_Sname[161] = "CS401" +SId_to_Sname[181] = "CS501" + +#SId_to_Sname[104] = "RS104" +#SId_to_Sname[105] = "RS105" +SId_to_Sname[106] = "RS106" +#SId_to_Sname[107] = "RS107" +#SId_to_Sname[108] = "RS108" +#SId_to_Sname[109] = "RS109" +#SId_to_Sname[122] = "RS202" +#SId_to_Sname[123] = "RS203" +#SId_to_Sname[124] = "RS204" +SId_to_Sname[125] = "RS205" +#SId_to_Sname[126] = "RS206" +#SId_to_Sname[127] = "RS207" +SId_to_Sname[128] = "RS208" +#SId_to_Sname[129] = "RS209" +SId_to_Sname[130] = "RS210" +#SId_to_Sname[143] = "RS303" +#SId_to_Sname[144] = "RS304" +SId_to_Sname[145] = "RS305" +SId_to_Sname[146] = "RS306" +SId_to_Sname[147] = "RS307" +#SId_to_Sname[148] = "RS308" +#SId_to_Sname[149] = "RS309" +SId_to_Sname[150] = "RS310" +SId_to_Sname[166] = "RS406" +SId_to_Sname[167] = "RS407" +SId_to_Sname[169] = "RS409" +SId_to_Sname[183] = "RS503" +SId_to_Sname[188] = "RS508" +SId_to_Sname[189] = "RS509" + +SId_to_Sname[201] = "DE601" +SId_to_Sname[202] = "DE602" +SId_to_Sname[203] = "DE603" +SId_to_Sname[204] = "DE604" +SId_to_Sname[205] = "DE605" +SId_to_Sname[206] = "FR606" +SId_to_Sname[207] = "SE607" +SId_to_Sname[208] = "UK608" + +## this just "inverts" the previous list, discarding unused values +Sname_to_SId_dict = {name:ID for ID,name in enumerate(SId_to_Sname) if name is not None} + +def even_antName_to_odd(even_ant_name): + even_num = int(even_ant_name) + odd_num = even_num + 1 + return str( odd_num ).zfill( 9 ) + +def antName_is_even(ant_name): + return not int(ant_name)%2 + +def odd_antName_to_even(odd_ant_name): + odd_num = int(odd_ant_name) + even_num = odd_num + 1 + return str( even_num ).zfill( 9 ) + + +#### plotting utilities #### +def set_axes_equal(ax): + '''Make axes of 3D plot have equal scale so that spheres appear as spheres, + cubes as cubes, etc.. This is one possible solution to Matplotlib's + ax.set_aspect('equal') and ax.axis('equal') not working for 3D. + + + + Input + ax: a matplotlib axis, e.g., as output from plt.gca(). + ''' + + x_limits = ax.get_xlim3d() + y_limits = ax.get_ylim3d() + z_limits = ax.get_zlim3d() + + x_range = abs(x_limits[1] - x_limits[0]) + x_middle = np.mean(x_limits) + y_range = abs(y_limits[1] - y_limits[0]) + y_middle = np.mean(y_limits) + z_range = abs(z_limits[1] - z_limits[0]) + z_middle = np.mean(z_limits) + + # The plot bounding box is a sphere in the sense of the infinity + # norm, hence I call half the max range the plot radius. + plot_radius = 0.5*max([x_range, y_range, z_range]) + + ax.set_xlim3d([x_middle - plot_radius, x_middle + plot_radius]) + ax.set_ylim3d([y_middle - plot_radius, y_middle + plot_radius]) + ax.set_zlim3d([z_middle - plot_radius, z_middle + plot_radius]) + + +### some math functions? ### + +def normalize_angle_radians( angle_radians ): + """For an angle in radians, return the equivalent angle that is garunteed be between -pi and pi""" + while angle_radians > np.pi: + angle_radians -= 2.0*np.pi + while angle_radians < -np.pi: + angle_radians += 2.0*np.pi + return angle_radians + + + + \ No newline at end of file