hDirectionTriangulationCartesian(direction, Error, positions_1, positions_2, positions_3, time1, time2, time3, sign_of_solution)
Calculates the direction of one source using the arrival times of the signal in three antennas returning the result in cartesian coordinates
Parameters
direction 3-Vector of X,Y,Z coordindates of the direction vector towards the source. Error Closure error - if non-zero then there was no valid solution. positions_1 3-Vector (x,y,z) with cartesian coordinates for the 1st antenna in meters positions_2 3-Vector (x,y,z) with cartesian coordinates for the 2nd antenna in meters positions_3 3-Vector (x,y,z) with cartesian coordinates for the 3rd antenna in meters time1 Arrival times of the signal at 1st antenna in seconds time2 Arrival times of the signal at 2nd antenna in seconds time3 Arrival times of the signal at 3rd antenna in seconds sign_of_solution There can be one or two solutions returned: use +/-1 for choosing one of the two, or 0 for returning both (in AzElEr)
Description
Given three antenna positions, and (pulse) arrival times for each antenna, get a direction of arrival (az, el) assuming a source at infinity (plane wave). From three antennas we get in general two solutions for the direction, unless the arrival times are out of bounds, i.e. larger than the light-time between two antennas. In the latter case a non-zero closure error is returned, which contains the magntiude of the error made.
Usually, when the 3 antennas are in a more or less horizontal plane, one of the solutions will appear to come from below the horizon (el < 0) and can be discarded. With the parameter sign_of_solution one can choose either the positive or the negative solution.
Output is returned in the scalar variables az, el, er. Typically this does not work in python, so this function should be called from within c++. In Python use hDirectionTriangulations, which returs the values in vectors.
Usage
``hDirectionTriangulation(az,el,err,[x1,y1,z1],[x1,y1,z1],[x1,y1,z1],t1,t2,t3,+/-1) -> az,el, err``
See also
hDirectionTriangulation(Azimuth, Elevation, Error, positions_1, positions_2, positions_3, time1, time2, time3, sign_of_solution)
Calculates the direction of one source using the arrival times of the signal in three antennas
Parameters
Azimuth Azimuth of the source. Elevation Elevation of the source. Error Closure error - if non-zero then there was no valid solution. positions_1 3-Vector (x,y,z) with cartesian coordinates for the 1st antenna in meters positions_2 3-Vector (x,y,z) with cartesian coordinates for the 2nd antenna in meters positions_3 3-Vector (x,y,z) with cartesian coordinates for the 3rd antenna in meters time1 Arrival times of the signal at 1st antenna in seconds time2 Arrival times of the signal at 2nd antenna in seconds time3 Arrival times of the signal at 3rd antenna in seconds sign_of_solution There can be one or two solutions returned: use +/-1 for choosing one of the two, or 0 for returning both (in AzElEr)
Description
Given three antenna positions, and (pulse) arrival times for each antenna, get a direction of arrival (az, el) assuming a source at infinity (plane wave). From three antennas we get in general two solutions for the direction, unless the arrival times are out of bounds, i.e. larger than the light-time between two antennas. In the latter case a non-zero closure error is returned, which contains the magntiude of the error made.
Usually, when the 3 antennas are in a more or less horizontal plane, one of the solutions will appear to come from below the horizon (el < 0) and can be discarded. With the parameter sign_of_solution one can choose either the positive or the negative solution.
Output is returned in the scalar variables az, el, er. Typically this does not work in python, so this function should be called from within c++. In Python use hDirectionTriangulations, which returs the values in vectors.
Usage
``hDirectionTriangulation(az,el,err,[x1,y1,z1],[x1,y1,z1],[x1,y1,z1],t1,t2,t3,+/-1) -> az,el, err``
See also
hDirectionTriangulations(Azimuth, Elevation, errors, origins, weights, positions, times, sign_of_solution)
Calculates all the directions of one source using the arrival times of the signal in all possible antenna triangles from a list of antennas as input
Parameters
Azimuth Vector containing the calculated azimuth directions towards the source for all triangles Elevation Vector containing the calculated elevation directions towards the source for all triangles errors Vector with closure errors for each triangle origins Vector of 3-tuples (x,y,z) with cartesian coordinates of the center point of each triangle weights Vector with weights indicating the estimated position accuracy of each triangle positions Vector of 3-tuples with cartesian coordinates for all antennas in meters times Vector with arrival times of the signal at each antenna in seconds sign_of_solution There can be one or two solutions returned: use +/-1 for choosing one of the two
Description
This function calculates a list of arrival directions (azimuth and elevation) of a pulse/source from all possible triangles of (groups of three) antennas and the measured arrival times at these antennas.
For the direction finding it is assumed that the signal arrives as a plane wave. If one triangle of antennas does not give a valid solution a solution at the horizion is returned with the complex part of the solution (the “closure error”) returned in err.
To cope with near-field sources for each triangle the coordinates of its center are also calculated. This allows one to later look for points where the dircetions from all triangles intercept.
The inverse of the average lengths of the triangle baselines are returned in weights giving some indication how accurate the measurement in principle can be.
From three antennas we get in general two solutions for the direction, unless the arrival times are out of bounds, i.e. larger than the light-time between two antennas. Usually, when the 3 antennas are in a more or less horizontal plane, one of the solutions will appear to come from below the horizon (el < 0) and can be discarded. With sign_of_solution one can pick either the positive of the negative elevation.
As input one only needs to provide arrival times and a vector with the corresponding antenna positions x,y,z (Cartesian coordinates) as a 3-tuple.
If N triangles are provided then the azimuth, elevation, origins, and weights vectors have N*(N-1)*(N-2)/6 results (times number of components of each result, i.e. times three for origins and times one for the rest.)
Output: [az1,az2,...], [el1,el2,...], [err1,err2,...],[(x1,y1,z1),...],[weight1, weights2, ....]) in radians for azimuth and elevation.
Usage
``hDirectionTriangulations([az],[el],[err],[origins],[weights],[positions],[times],+/-1) -> az,el,err,origins,weights``
See also
Example
file=open('/Users/falcke/LOFAR/usg/data/lofar/oneshot_level4_CS017_19okt_no-9.h5')
file['ANTENNA_SET']='LBA_OUTER'
file['BLOCKSIZE']=2^17
oddantennas=['017000001', '017000002', '017000005', '017000007', '017001009', '017001010', '017001012', '017001015', '017002017', '017002019', '017002020', '017002023', '017003025', '017003026', '017003029', '017003031', '017004033', '017004035', '017004037', '017004039', '017005041', '017005043', '017005045', '017005047', '017006049', '017006051', '017006053', '017006055', '017007057', '017007059', '017007061', '017007063', '017008065', '017008066', '017008069', '017008071', '017009073', '017009075', '017009077', '017009079', '017010081', '017010083', '017010085', '017010087', '017011089', '017011091', '017011093', '017011095']
file['SELECTED_DIPOLES']=oddantennas
positions=file['ANTENNA_POSITIONS']
timeseries_data=file['TIMESERIES_DATA']
(start,end)=trun('LocatePulseTrain',rf.TimeBeamIncoherent(timeseries_data),nsigma=7,maxgap=3)
start-=16; end=start+int(2^ceil(log(end-start,2)));
timeseries_data_cut=hArray(float,[timeseries_data.shape()[-2],end-start])
timeseries_data_cut[...].copy(timeseries_data[...,start:end])
timeseries_data_cut.abs()
timeseries_data_cut[...].runningaverage(7,hWEIGHTS.GAUSSIAN)
mx=task('FitMaxima')(timeseries_data_cut,doplot=True,refant=0,sampleinterval=5,plotend=48,peak_width=10,legend=file['SELECTED_DIPOLES'])
mx.lags *= 10^-9
N=48
NT=N*(N-1)*(N-2)/6 # Number of Triangles ....
az=hArray(float,[NT])
el=hArray(float,[NT])
er=hArray(float,[NT])
weights=hArray(float,[NT])
origins=hArray(float,[NT,3])
hDirectionTriangulations(az,el,er,origins,weights,positions[-N:],hArray(mx.lags[-N:]),+1);
az/=deg; el/=deg;
indx=hArray(int,[NT])
ngood=indx.findlessthan(er,1e-10)
scrt=hArray(float,[NT])
scrt.copy(az,indx,ngood)
scrt[0:ngood.val()].mean()
# -> Vector(float, 1, fill=[147.618])
scrt.copy(el,indx,ngood)
mean_el=scrt[0:ngood.val()].mean()
rms_el=scrt[0:ngood.val()].stddev()
ngood=indx.findbetween(scrt,mean_el-rms_el,mean_el+rms_el)
scrt.copy(scrt,indx,ngood)
az=hArray(float,[NT]);az[...].copy(azel[...,0])
el=hArray(float,[NT]);el[...].copy(azel[...,1])
plt.clf(); plt.plot(az,el,marker='o',linestyle='None')
print 'az =',az.vec().mean(),'+/-',az.vec().stddev(),'el =',el.vec().mean(),'+/-',el.vec().stddev()
el. 86 -> NAN
make sure +elev is always pos.
closure time: dt1 + dt2 -dt3
azel -> hArray(float, [4L], fill=[-211.972,-83.0455,-211.822,83.0441]) # len=4 slice=[0:4])
positions[-3] -> hArray(float, [48L, 3L], fill=[23.568,-27.952,-0.0004]) # len=144 slice=[135:138])
positions[-2] -> hArray(float, [48L, 3L], fill=[12.106,-39.642,-0.003]) # len=144 slice=[138:141])
positions[-1] -> hArray(float, [48L, 3L], fill=[-28.833,-25.937,-0.007]) # len=144 slice=[141:144])
mx.lags[-3:] -> Vector(float, 3, fill=[-2.25e-08,-2.40625e-08,-1.0625e-08])
import pycrtools.srcfind as srcfind
p=np.array(positions[-3:].vec())
t=np.array(mx.lags[-3:].vec())
(az1,el1,az2,el2)=srcfind.directionFromThreeAntennas(p,t)
hDirectionTriangulationsCartesian(directions, errors, origins, weights, positions, times, sign_of_solution)
Calculates all the directions of one source using the arrival times of the signal in all possible antenna triangles from a list of antennas as input
Parameters
directions Vector containing the calculated direction vectors in Cartesian coordinates (3 values: x,y,z) towards the source for all triangles errors Vector with closure errors for each triangle origins Vector of 3-tuples (x,y,z) with cartesian coordinates of the center point of each triangle weights Vector with weights indicating the estimated position accuracy of each triangle positions Vector of 3-tuples with cartesian coordinates for all antennas in meters times Vector with arrival times of the signal at each antenna in seconds sign_of_solution There can be one or two solutions returned: use +/-1 for choosing one of the two
Description
This function calculates a list of arrival directions in Cartesian coordinates of a pulse/source from all possible triangles of (groups of three) antennas and the measured arrival times at these antennas.
For the direction finding it is assumed that the signal arrives as a plane wave. If one triangle of antennas does not give a valid solution a solution at the horizion is returned with the complex part of the solution (the “closure error”) returned in err.
To cope with near-field sources for each triangle the coordinates of its center are also calculated. This allows one to later look for points where the dircetions from all triangles intercept.
The inverse of the average lengths of the triangle baselines are returned in weights giving some indication how accurate the measurement in principle can be.
From three antennas we get in general two solutions for the direction, unless the arrival times are out of bounds, i.e. larger than the light-time between two antennas. Usually, when the 3 antennas are in a more or less horizontal plane, one of the solutions will appear to come from below the horizon (el < 0) and can be discarded. With sign_of_solution one can pick either the positive of the negative elevation.
As input one only needs to provide arrival times and a vector with the corresponding antenna positions x,y,z (Cartesian coordinates) as a 3-tuple.
If N triangles are provided then the Directions, Origins, and weights vectors have N*(N-1)*(N-2)/6 results (times number of components of each result, i.e. times three for origins and times one for the rest.)
Output: [az1,az2,...], [el1,el2,...], [err1,err2,...],[(x1,y1,z1),...],[weight1, weights2, ....]) in radians for azimuth and elevation.
Usage
``hDirectionTriangulationsCartesian([directions],[errors],[origins],[weights],[positions],[times],+/-1) -> cartesian directions,err,origins,weights``
See also
hDirectionTriangulation(), hDirectionTriangulationsCartesian()
Example
file=open('/Users/falcke/LOFAR/usg/data/lofar/oneshot_level4_CS017_19okt_no-9.h5')
file['ANTENNA_SET']='LBA_OUTER'
file['BLOCKSIZE']=2^17
oddantennas=['017000001', '017000002', '017000005', '017000007', '017001009', '017001010', '017001012', '017001015', '017002017', '017002019', '017002020', '017002023', '017003025', '017003026', '017003029', '017003031', '017004033', '017004035', '017004037', '017004039', '017005041', '017005043', '017005045', '017005047', '017006049', '017006051', '017006053', '017006055', '017007057', '017007059', '017007061', '017007063', '017008065', '017008066', '017008069', '017008071', '017009073', '017009075', '017009077', '017009079', '017010081', '017010083', '017010085', '017010087', '017011089', '017011091', '017011093', '017011095']
file['SELECTED_DIPOLES']=oddantennas
positions=file['ANTENNA_POSITIONS']
timeseries_data=file['TIMESERIES_DATA']
(start,end)=trun('LocatePulseTrain',rf.TimeBeamIncoherent(timeseries_data),nsigma=7,maxgap=3)
start-=16; end=start+int(2^ceil(log(end-start,2)));
timeseries_data_cut=hArray(float,[timeseries_data.shape()[-2],end-start])
timeseries_data_cut[...].copy(timeseries_data[...,start:end])
timeseries_data_cut.abs()
timeseries_data_cut[...].runningaverage(7,hWEIGHTS.GAUSSIAN)
mx=task('FitMaxima')(timeseries_data_cut,doplot=True,refant=0,sampleinterval=5,plotend=48,peak_width=10,legend=file['SELECTED_DIPOLES'])
mx.lags *= 10^-9
N=48
NT=N*(N-1)*(N-2)/6 # Number of Triangles ....
directions=hArray(float,[NT,3])
er=hArray(float,[NT])
weights=hArray(float,[NT])
origins=hArray(float,[NT,3])
hDirectionTriangulationsCartesian(directions,er,origins,weights,positions[-N:],hArray(mx.lags[-N:]),+1);
indx=hArray(int,[NT])
ngood=indx.findlessthan(er,1e-10).val()
scrt=hArray(float,[ngood,3])
scrt.copyvec(directions,indx,ngood,3)
meandirection=hArray(float,[3])
meandirection.mean(scrt)
meandirection /= meandirection.vectorlength().val()
#meandirection - > hArray(float, [3L], fill=[0.0824118,-0.111287,0.990365]) # len=3 slice=[0:3])
#pytmf.spherical2cartesian(1, pi2-81.7932*deg,pi2-143.4092*deg) -> (0.085090553169460376, -0.11461297030462134, 0.98975929639446536)
sp=pytmf.cartesian2spherical(meandirection[0],meandirection[1],meandirection[2]); azel=(180-(sp[2]+pi2)/deg,90-(sp[1])/deg)
# azel -> (143.47870717075114, 82.040161833968313)
hSkewLinesDistanceToClosestApproach(distance, baseline, origin1, direction1, origin2, direction2, max_distance)
Calculates the point of closest approach of two lines in 3D as distance along the first input vector. Multiple lines can be given.
Parameters
distance Distance of point of closest approach along first line from origin of first line (in Meters). baseline Baseline (separation) between the two points of origin (in Meters) - output value. origin1 3-Vector of Cartesian (X,Y,Z) coordindates of the point of origin of the first lines(normalized) direction vector of the first line. direction1 3-Vector of Cartesian (X,Y,Z) coordindates of the (normalized) direction vector of the first line. origin2 3-Vector of Cartesian (X,Y,Z) coordindates of the point of origin of the first lines(normalized) direction vector of the second line. direction2 3-Vector of Cartesian (X,Y,Z) coordindates of the (normalized) direction vector of the second line. max_distance Maximum allowed distance to return - is also used to voids divide by zero.
Description
Typically three-dimensional lines do not intersect (skew lines), however, there is a point of closest approach. In this function four vectors are provided as input, which provide the points of origin and the normalized (length=1) direction vector of the two lines. distance``will then contain the distance from the origin of the first vector (``origin1) along its direction vector (direction1) until the point of closest approach is reached. Per vector 3 Cartesian coordinates have to be provided, hence the length of the input vector has to be a mutliple of three.
Looping behaviour: The four vectors can contain multiple vectors. Calculations proceed until the end of origin2 or direction2 are reached. If direction1 and origin1 are shorter they will be wrapped und taken repeatedly. For example, by providing a single vector in origins1, directions1 and multiple ones in origin2, direction2 one obtains a list of closest approach distances along the first vector realtive to all the other vectors.
The Math is taken from http://www.mathworks.com/matlabcentral/newsreader/view_thread/246420
In general you can find the two points, A0 and B0, on the respective lines A1A2 and B1B2 which are closest together and determine if they coincide or else see how far apart they lie.
The following is a specific formula using matlab for determining these closest points A0 and B0 in terms of vector cross products and dot products. Assume that all four points are represented by three-element column vectors or by three-element row vectors. Do the following:
nA = dot(cross(B2-B1,A1-B1),cross(A2-A1,B2-B1)); nB = dot(cross(A2-A1,A1-B1),cross(A2-A1,B2-B1)); d = dot(cross(A2-A1,B2-B1),cross(A2-A1,B2-B1)); A0 = A1 + (nA/d)*(A2-A1); B0 = B1 + (nB/d)*(B2-B1);
Here we actually use a point of origin and a direction vector, which is direction1 = A2-A1 and direction2 = B2-B1. Also, we do not calculate nB.
Usage
``hSkewLinesDistanceToClosestApproach(distance,baseline,origin1,direction1,origin2,direction2,max_distance) -> distance, baseline``
See also
hDirectionTriangulations(), hDirectionTriangulationCartesian()
Example
#simple example first
distance=hArray(float,[1])
error_distance=hArray(float,[1])
origin1=hArray([0.,0.,0.])
direction1=hArray([0.,0.,1.])
origin2=hArray([1.,1.,0.])
direction2=hArray([-2.^-0.5,0.,2.^-0.5])
hSkewLinesDistanceToClosestApproach(distance,error_distance,origin1,direction1,origin2,direction2,10.,1.);
#more invovled example
file=open('/Users/falcke/LOFAR/usg/data/lofar/oneshot_level4_CS017_19okt_no-9.h5')
file['ANTENNA_SET']='LBA_OUTER'
file['BLOCKSIZE']=2^17
oddantennas=['017000001', '017000002', '017000005', '017000007', '017001009', '017001010', '017001012', '017001015', '017002017', '017002019', '017002020', '017002023', '017003025', '017003026', '017003029', '017003031', '017004033', '017004035', '017004037', '017004039', '017005041', '017005043', '017005045', '017005047', '017006049', '017006051', '017006053', '017006055', '017007057', '017007059', '017007061', '017007063', '017008065', '017008066', '017008069', '017008071', '017009073', '017009075', '017009077', '017009079', '017010081', '017010083', '017010085', '017010087', '017011089', '017011091', '017011093', '017011095']
file['SELECTED_DIPOLES']=oddantennas
positions=file['ANTENNA_POSITIONS']
timeseries_data=file['TIMESERIES_DATA']
(start,end)=trun('LocatePulseTrain',rf.TimeBeamIncoherent(timeseries_data),nsigma=7,maxgap=3)
start-=16; end=start+int(2^ceil(log(end-start,2)));
timeseries_data_cut=hArray(float,[timeseries_data.shape()[-2],end-start])
timeseries_data_cut[...].copy(timeseries_data[...,start:end])
timeseries_data_cut.abs()
timeseries_data_cut[...].runningaverage(7,hWEIGHTS.GAUSSIAN)
mx=task('FitMaxima')(timeseries_data_cut,doplot=True,refant=0,sampleinterval=5,plotend=48,peak_width=10,legend=file['SELECTED_DIPOLES'])
mx.lags *= 10^-9
N=48
NT=N*(N-1)*(N-2)/6 # Number of Triangles ....
directions=hArray(float,[NT,3])
er=hArray(float,[NT])
weights=hArray(float,[NT])
origins=hArray(float,[NT,3])
hDirectionTriangulationsCartesian(directions,er,origins,weights,positions[-N:],hArray(mx.lags[-N:]),+1);
indx=hArray(int,[NT])
ngood=indx.findlessthan(er,1e-10).val()
scrt=hArray(float,[ngood,3])
scrt.copyvec(directions,indx,ngood,3)
meandirection=hArray(float,[3])
meandirection.mean(scrt)
meandirection /= meandirection.vectorlength().val()
#meandirection - > hArray(float, [3L], fill=[0.0824118,-0.111287,0.990365]) # len=3 slice=[0:3])
#pytmf.spherical2cartesian(1, pi2-81.7932*deg,pi2-143.4092*deg) -> (0.085090553169460376, -0.11461297030462134, 0.98975929639446536)
sp=pytmf.cartesian2spherical(meandirection[0],meandirection[1],meandirection[2]); azel=(180-(sp[2]+pi2)/deg,90-(sp[1])/deg)
# azel -> (143.47870717075114, 82.040161833968313)
scrt.copyvec(origins,indx,ngood,3)
meanorigin=hArray(float,[3])
meanorigin.mean(scrt)
distance=hArray(float,[NT])
error_distance=hArray(float,[NT])
hSkewLinesDistanceToClosestApproach(distance,error_distance,meanorigin,meandirection,origins,directions,10000.,100.);
scrt2=hArray(float,[ngood])
scrt2.copyvec(distance,indx,ngood,1)
hSpectralPower(outvec, vecin)
Calculates the abs (square of square root) of a complex vector and adds it to the output vector.
Parameters
outvec Output vector containing the abs (square root of square) of the input vector added to itself. vecin Input vector containing the complex spectrum. (Looping parameter)
Description
The fact that the result is added to the output vector allows one to call the function multiple times and get a summed spectrum. If you need it only once, just fill the vector with zeros.
The number of loops (if used with an hArray) is here determined by the second parameter!
Usage
vecout.spectralpower(vecin) -> vecout_i=sqrt(real((vecin_i)*conj(vecin_i)))
See also
hSquareAdd(), hSpectralPower2()
Example
>>> spectrum=hArray(float,[1,128])
>>> cplxfft=hArray(complex,[10,128],fill=1+0j)
>>> spectrum[...].spectralpower(cplxfft[...])
hSpectralPower2(outvec, vecin)
Calculates the power of a complex spectrum and adds it to an output vector.
Parameters
outvec Output vector containing the square of the vector added to itself. vecin Input vector containing the complex spectrum. (Looping parameter)
Description
The fact that the result is added to the output vector allows one to call the function multiple times and get a summed spectrum. If you need it only once, just fill the vector with zeros.
The number of loops (if used with an hArray) is here determined by the second parameter!
Usage
vecout.spectralpower2(vecin) -> vecout_i=real((vecin_i)*conj(vecin_i))
See also
hSquareAdd(), hSpectralPower()
Example
>>> spectrum=hArray(float,[1,128])
>>> cplxfft=hArray(complex,[10,128],fill=1+0j)
>>> spectrum[...].spectralpower2(cplxfft[...])