Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Preliminary modifications from Lech #9

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

[WIP] Preliminary modifications from Lech #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 58 additions & 41 deletions examples/simulation/shower-event.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,72 +2,89 @@
import astropy.units as u
from grand import ECEF, LTP
from grand.simulation import Antenna, ShowerEvent, TabulatedAntennaModel
#from grand.simulation import ShowerEvent, TabulatedAntennaModel
#from grand.simulation.antenna .generic import compute_voltage
import numpy as np
from matplotlib import pyplot as plt
from astropy.coordinates import BaseRepresentation, CartesianRepresentation

import os
grand_astropy = True
try:
if os.environ['GRAND_ASTROPY']=="0":
grand_astropy=False
except:
pass

# Note: This examples requires some test data to be downloaded localy. You can
# get these data from the grand store as:
#
# wget https://github.com/grand-mother/store/releases/download/101/HorizonAntenna_EWarm_leff_loaded.npy
#
# mkdir -p tests/simulation/data/zhaires
# cd tests/simulation/data/zhaires
# wget https://github.com/grand-mother/store/releases/download/101/zhaires-test.tar.gz
# tar -xzf zhaires-test.tar.gz
# cd -


# Load the radio shower simulation data
showerdir = '../../tests/simulation/data/zhaires'
shower = ShowerEvent.load(showerdir)
shower = ShowerEvent.load('tests/simulation/data/zhaires')
if shower.frame is None:
shower.localize(39.5 * u.deg, 90.5 * u.deg) # Coreas showers have no
# localization info. This must
# be set manually

print("Shower frame=",shower.frame)
print("Zenith (Zhaires?!) =",shower.zenith)
print("Azimuth (Zhaires?!) =",shower.azimuth)
print("Xmax=",shower.maximum)
#shower.core=CartesianRepresentation(0, 0, 2900, unit='m')
print("Core=",shower.core)


# Define an antenna model
#
# A tabulated model of the Butterfly antenna is used. Note that a single EW
# arm is assumed here for the sake of simplicity

antenna_model = TabulatedAntennaModel.load('./HorizonAntenna_EWarm_leff_loaded.npy')
antenna_model = TabulatedAntennaModel.load(
'HorizonAntenna_EWarm_leff_loaded.npy')

# Loop over electric fields and compute the corresponding voltages
for antenna_index, field in shower.fields.items():
# Compute the antenna local frame
#
# The antenna is placed within the shower frame. It is oriented along the
# local magnetic North by using an ENU/LTP frame (x: East, y: North, z: Upward)
antenna_location = shower.frame.realize_frame(field.electric.r)
print(antenna_index,"Antenna pos=",antenna_location)
# local magnetic North by using an ENU/LTP frame (x: East, y: North, z:
# Upward)

if grand_astropy:
print("field", field.electric.r, shower.frame)

antenna_frame = LTP(location=antenna_location, orientation='NWU',magnetic=True, obstime=shower.frame.obstime)
antenna = Antenna(model=antenna_model, frame=antenna_frame)
# LWP: joins field.electric.r and shower.frame
antenna_location = shower.frame.realize_frame(field.electric.r)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This wraps the Cartesian coordinates of the E field position, i.e. the antenna one, with an astropy.CoordinatesFrame allowing for frame transforms.

print(antenna_location)
antenna_frame = LTP(location=antenna_location, orientation='ENU',
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This creates an ENU local frame at the antenna position.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Even though we are now using a single GRAND frame, the antenna response is given in a local frame that might be rotated. In this example I am assuming that the antenna local frame is ENU.

magnetic=True, obstime=shower.frame.obstime)
# LWP: joins antenna_model and antenna_frame
antenna = Antenna(model=antenna_model, frame=antenna_frame)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The 'frame' parameter needs a replacement. E.g. it could be three angles: inclination, declination and twist indicating the antenna orientation w.r.t. the GRAND frame?

In the long run I would rather have those expressed in geographic coordinates that can be directly measured on the field and let the Antenna constructor properly compute the corresponding values in the GRAND frame. But this would contradict the credo that everything needs to be in a single frame. We might also leave to 'experts' the task of computing the antenna orientation in the GRAND frame.

print("antenna_frame", antenna_frame)
#print("antenna_model", antenna_model)
#print("antenna", antenna)
else:
# LWP: Move shower frame by field.electric.r, then somehow add this to Antenna without astropy
# This need deep moditication of the Antenna class, so for now do the astropy stuff
print("field", field.electric.r, shower.frame)

# LWP: joins field.electric.r and shower.frame
antenna_location = shower.frame.realize_frame(field.electric.r)
print(antenna_location)
antenna_frame = LTP(location=antenna_location, orientation='ENU',
magnetic=True, obstime=shower.frame.obstime)
# LWP: joins antenna_model and antenna_frame
antenna = Antenna(model=antenna_model, frame=antenna_frame)
print("antenna_frame", antenna_frame)
#print("antenna_model", antenna_model)
#print("antenna", antenna)

# Compute the voltage on the antenna
#
# The electric field is assumed to be a plane-wave originating from the
# shower axis at the depth of maximum development. Note that the direction
# of observation and the electric field components are provided in the
# shower frame. This is indicated by the `frame` named argument.
direction = shower.maximum - field.electric.r
print("Direction to Xmax = ",direction)
#print(antenna_frame.realize_frame(direction))
Exyz = field.electric.E.represent_as(CartesianRepresentation)

field.voltage = antenna.compute_voltage(direction, field.electric,frame=shower.frame)

plt.figure()
plt.subplot(211)
plt.plot(field.electric.t,Exyz.x,label='Ex')
plt.plot(field.electric.t,Exyz.y,label='Ey')
plt.plot(field.electric.t,Exyz.z,label='Ez')
plt.xlabel('Time ('+str(field.electric.t.unit)+')')
plt.ylabel('Efield ('+str(Exyz.x.unit)+')')
plt.legend(loc='best')
plt.subplot(212)
plt.plot(field.voltage.t,field.voltage.V,label='V$_{EW}$')
plt.xlabel('Time ('+str(field.voltage.t.unit)+')')
plt.ylabel('Voltage ('+str(field.voltage.V.unit)+')')
plt.legend(loc='best')
plt.show()
direction = shower.maximum - field.electric.r
print("computing voltage")
field.voltage = antenna.compute_voltage(direction, field.electric,
frame=shower.frame)
print("computed voltage", field.voltage)
exit()
55 changes: 45 additions & 10 deletions grand/simulation/antenna/generic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,15 @@

from ... import io, ECEF, LTP

import os
grand_astropy = True
try:
if os.environ['GRAND_ASTROPY']=="0":
grand_astropy=False
except:
pass


__all__ = ['Antenna', 'AntennaModel', 'ElectricField', 'MissingFrameError',
'Voltage']

Expand Down Expand Up @@ -92,6 +101,7 @@ class Antenna:
def compute_voltage(self, direction: Union[ECEF, LTP, BaseRepresentation],
field: ElectricField, frame: Union[ECEF, LTP, None]=None) \
-> Voltage:

# Uniformise the inputs
if self.frame is None:
antenna_frame = None
Expand All @@ -105,7 +115,6 @@ def compute_voltage(self, direction: Union[ECEF, LTP, BaseRepresentation],
else:
antenna_frame = cast(Union[ECEF, LTP], self.frame)
frame_required = False

if field.frame is None:
E_frame, frame_required = frame, True
else:
Expand All @@ -129,8 +138,15 @@ def irfft(q):
def fftfreq(n, t):
dt = (t[1] - t[0]).to_value('s')
return numpy.fft.fftfreq(n, dt) * u.Hz

E = field.E.represent_as(CartesianRepresentation)

# print("E not repr", field.E)
# LWP: Seems unnecessary at least in the tested case - field.E already in cartesian
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed. My understanding is that now coordinates would always be Cartesian in the local GRAND frame by default. I.e. it is the user's responsibility to provide proper coordinates.

if grand_astropy:
E = field.E.represent_as(CartesianRepresentation)
else:
E = field.E
# print("E cart", E)
# exit()
Ex = rfft(E.x)
Ey = rfft(E.y)
Ez = rfft(E.z)
Expand All @@ -139,22 +155,41 @@ def fftfreq(n, t):
if dir_frame is not None:
# Change the direction to the antenna frame
if isinstance(direction, BaseRepresentation):
direction = dir_frame.realize_frame(direction)
direction = direction.transform_to(antenna_frame) # Now compute direction vector data in antenna frame
direction = direction.data
print("dir pre", direction, dir_frame, antenna_frame)
# LWP: this joins direction and dir_frame - not needed if we get rid of astropy
if grand_astropy:
direction = dir_frame.realize_frame(direction)
else:
E = field.E
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this line. I guess it is a mistake? Maybe just pass as direction is left unchanged here.

Note that dir_frame should not be needed anymore now if the direction is always assumed to be expressed in the GRAND frame.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

E is undefined before that, and it is needed just below for calculating Ex, Ey, Ez


print("dir post", direction)
#exit()

# LWP: This adds the difference of the bases (shower_frame-antenna_frame) to the direction, which can be done manually
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The fact that direction gets offset was a bug. The direction should only be rotated actually since it is a vector, not a point.

if grand_astropy:
direction = direction.transform_to(antenna_frame).data
print("dire postpost", direction)
else:
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as comment above. In this case the direction should be rotated, not translated.

# Bases are in meters
base_diff = numpy.array([dir_frame.location.x.value-antenna_frame.location.x.value, dir_frame.location.y.value-antenna_frame.location.y.value, dir_frame.location.z.value-antenna_frame.location.z.value])/1000
# LWP: The frames are NWU and ENU orientation, yet their values are as similar as they had the same orientation...
print(antenna_frame, dir_frame)
print(base_diff)
# LWP: ToDo: Need to check the orientations of frames and perhaps reposition the x,y,z
direction = numpy.array([direction.x.value,direction.y.value,direction.z.value])+base_diff
print(direction)
#exit()

Leff:CartesianRepresentation
Leff = self.model.effective_length(direction, f)

if antenna_frame is not None:
# Change the effective length to the E-field frame
tmp = antenna_frame.realize_frame(Leff)
tmp = tmp.transform_to(E_frame)
# Transorm from antenna_frame to E_frame is a translation.
# The Leff vector norm should therefeore not be modified, but it is because it is defined as a point in astropy coordinates
Leff = tmp.cartesian

# Here we have to do an ugly patch for Leff values to be correct
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was Olivier's hack in order to patch my bug. If astropy is removed this hack should not be needed anymore.

V = irfft(Ex * (Leff.x - Leff.x[0]) + Ey * (Leff.y - Leff.y[0]) + Ez * (Leff.z - Leff.z[0]))
V = irfft(Ex * Leff.x + Ey * Leff.y + Ez * Leff.z)
t = field.t
t = t[:V.size]

Expand Down
87 changes: 78 additions & 9 deletions grand/simulation/antenna/tabulated.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,20 @@
PhysicsSphericalRepresentation
import astropy.units as u
import numpy
import h5py

from .generic import AntennaModel
from ... import io

import os
grand_astropy = True
try:
if os.environ['GRAND_ASTROPY']=="0":
grand_astropy=False
except:
pass


__all__ = ['DataTable', 'TabulatedAntennaModel']


Expand Down Expand Up @@ -70,6 +80,8 @@ def load(cls, source: Union[str, Path, io.DataNode]) \
source = Path(source)
if source.suffix == '.npy':
loader = '_load_from_numpy'
elif source.suffix == '.hdf5':
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed :)

loader = '_load_from_hdf5'
else:
loader = '_load_from_datafile'

Expand Down Expand Up @@ -120,17 +132,67 @@ def _load_from_numpy(cls, path: Union[Path, str]) -> TabulatedAntennaModel:
leff_phi = leffp, phase_phi = phasep)
return cls(table=t)

@classmethod
def _load_from_hdf5(cls, path: Union[Path, str], arm: str = "EW") -> TabulatedAntennaModel:
print("tuutu")
# Open the HDF5 file with antenna info
ant_file = h5py.File(path, 'r')
# Select the requested arm from the file
ant_arm = ant_file[arm]
# Load the antenna parameters (to numpy arrays, thus [:])
f = ant_arm["frequency"][:]
R = ant_arm["resistance"][:]
X = ant_arm["reactance"][:]
theta = ant_arm["theta"][:]
phi = ant_arm["phi"][:]
lefft = ant_arm["leff_theta"][:]
leffp = ant_arm["leff_phi"][:]
phaset = ant_arm["phase_theta"][:]
phasep = ant_arm["phase_phi"][:]
print("all shape", f.shape, R.shape, X.shape, theta.shape, phi.shape, lefft.shape, leffp.shape, phaset.shape, phasep.shape)

n_f = f.shape[0]
n_theta = len(numpy.unique(theta[0,:]))
n_phi = int(R.shape[1] / n_theta)
shape = (n_f, n_phi, n_theta)

dtype = 'f4'
f = f[:,0].astype(dtype) * u.MHz
theta = theta[0, :n_theta].astype(dtype) * u.deg
phi = phi[0, ::n_theta].astype(dtype) * u.deg
R = R.reshape(shape).astype(dtype) * u.Ohm
X = X.reshape(shape).astype(dtype) * u.Ohm
lefft = lefft.reshape(shape).astype(dtype) * u.m
leffp = leffp.reshape(shape).astype(dtype) * u.m
phaset = numpy.deg2rad(phaset).reshape(shape).astype(dtype) * u.rad
phasep = numpy.deg2rad(phasep).reshape(shape).astype(dtype) * u.rad

t = DataTable(frequency = f, theta = theta, phi = phi, resistance = R,
reactance = X, leff_theta = lefft, phase_theta = phaset,
leff_phi = leffp, phase_phi = phasep)
return cls(table=t)


def effective_length(self, direction: BaseRepresentation,
frequency: u.Quantity) -> CartesianRepresentation:

direction = direction.represent_as(PhysicsSphericalRepresentation)
theta, phi = direction.theta, direction.phi

# LWP: Replace with manual conversion to spherical, assuming direction is a numpy array
if grand_astropy:
direction = direction.represent_as(PhysicsSphericalRepresentation)
theta, phi = direction.theta, direction.phi
else:
theta = numpy.arctan2(numpy.sqrt(direction[0]**2+direction[1]**2),direction[2])
phi = numpy.arctan2(direction[1],direction[0])
Comment on lines +184 to +185
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indee, this would be enough for now. With the proposed coordinates type extending numpy.ndarray it could be something like:

spherical = direction.parametrize_as(SphericalVector)
theta, phi = spherical.theta, spherical.phi


# Interpolate using a tri-linear interpolation in (f, phi, theta)
t = self.table

dtheta = t.theta[1] - t.theta[0]
rt1 = ((theta - t.theta[0]) / dtheta).to_value(u.one)
# LWP: subtracting values from numpy array
if grand_astropy:
rt1 = ((theta - t.theta[0]) / dtheta).to_value(u.one)
else:
rt1 = ((theta - t.theta[0].to_value('rad')) / dtheta.to_value('rad'))
Comment on lines +191 to +195
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In principle t.theta and theta should both have the same unit since we decided to use a single system of units. Hence, this would simply be:

rt1 = (theta - t.theta[0]) / dtheta

it0 = int(numpy.floor(rt1) % t.theta.size)
it1 = it0 + 1
if it1 == t.theta.size: # Prevent overflow
Expand All @@ -140,7 +202,12 @@ def effective_length(self, direction: BaseRepresentation,
rt0 = 1 - rt1

dphi = t.phi[1] - t.phi[0]
rp1 = ((phi - t.phi[0]) / dphi).to_value(u.one)
# LWP: subtracting values from numpy array
if grand_astropy:
rp1 = ((phi - t.phi[0]) / dphi).to_value(u.one)
else:
rp1 = ((phi - t.phi[0].to_value('rad')) / dphi.to_value('rad'))

Comment on lines +205 to +210
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as previous comment

ip0 = int(numpy.floor(rp1) % t.phi.size)
ip1 = ip0 + 1
if ip1 == t.phi.size: # Results are periodic along phi
Expand All @@ -156,8 +223,6 @@ def interp(v):
rp0 * rt1 * v[:, ip0, it1] + rp1 * rt1 * v[:, ip1, it1]
return numpy.interp(x, xp, fp, left=0, right=0)



ltr = interp(t.leff_theta.to_value('m'))
lta = interp(t.phase_theta.to_value('rad'))
lpr = interp(t.leff_phi.to_value('m'))
Expand All @@ -167,8 +232,12 @@ def interp(v):
lt = ltr * numpy.exp(1j * lta)
lp = lpr * numpy.exp(1j * lpa)

t, p = theta.to_value('rad'), phi.to_value('rad')

# LWP: to_value not needed anymore
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes :)

if grand_astropy:
t, p = theta.to_value('rad'), phi.to_value('rad')
else:
t, p = theta, phi

ct, st = numpy.cos(t), numpy.sin(t)
cp, sp = numpy.cos(p), numpy.sin(p)
lx = lt * ct * cp - sp * lp
Expand Down