[WIP] Preliminary modifications from Lech #9
Conversation
| 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 |
There was a problem hiding this comment.
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: | ||
| direction = dir_frame.realize_frame(direction) | ||
| else: | ||
| E = field.E |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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 |
There was a problem hiding this comment.
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: |
There was a problem hiding this comment.
Same as comment above. In this case the direction should be rotated, not translated.
| Leff = tmp.cartesian | ||
|
|
||
| # Here we have to do an ugly patch for Leff values to be correct |
There was a problem hiding this comment.
This was Olivier's hack in order to patch my bug. If astropy is removed this hack should not be needed anymore.
| @@ -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': | |||
| theta = numpy.arctan2(numpy.sqrt(direction[0]**2+direction[1]**2),direction[2]) | ||
| phi = numpy.arctan2(direction[1],direction[0]) |
There was a problem hiding this comment.
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
| # 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')) |
There was a problem hiding this comment.
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
| # 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')) | ||
|
|
| @@ -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 | |||
| 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) |
There was a problem hiding this comment.
This wraps the Cartesian coordinates of the E field position, i.e. the antenna one, with an astropy.CoordinatesFrame allowing for frame transforms.
| # 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', |
There was a problem hiding this comment.
This creates an ENU local frame at the antenna position.
There was a problem hiding this comment.
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.
| 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) |
There was a problem hiding this comment.
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.
Fake PR, just for discussing ...