rise time not found (too close from reference time ?)

[tstolarczyk]

Hello,

astroplan 0.7 dev1236

I found several cases where obs.target_rise_time miss the closest next rise time and return the next one.
It seems it happens when the reference time is too close from the rise time to be found.
Below an example (jupyter notebook) that reproduces this anomaly.
Thanks in advance for your help,
Thierry S.

from   astropy.coordinates   import SkyCoord
from   astropy.coordinates   import EarthLocation
from   astropy.coordinates   import AltAz
import astropy.units as u
from astroplan import Observer, FixedTarget
from astropy.time import Time
xyz = EarthLocation.from_geocentric( 5327285.09211954, 
                                          -1718777.11250295, 
                                          3051786.7327476,
                                          unit="m")

radec = SkyCoord(ra=15.44379741*u.degree, dec=-26.1439*u.degree, frame='icrs') 
tref = Time('2046-01-30 13:30:50.505',format="iso",scale="utc")
horizon = 10*u.deg

obs = Observer(location  = xyz, name = "test", timezone ="utc")   
target = FixedTarget(coord=radec, name="mysource")

high = obs.target_is_up(tref, target, horizon=10*u.deg )
if (not high): print("It is below the horizon, compute next rise")
t_rise = obs.target_rise_time(tref,target,which="next",horizon=10*u.deg)
print(t_rise.iso)   

It is below the horizon, compute next rise
2046-01-31 13:28:17.143

# Previous rise
t_rise_bef = obs.target_rise_time(tref,target,which="previous",horizon=10*u.deg)
print(t_rise_bef.iso) 

2046-01-29 13:36:09.289

# Find rise time "manually"
import numpy as np
trange = 500
dt = np.linspace(-trange/2,trange/2,10)
tsamp = tref + dt*u.s
altaz = radec.transform_to(AltAz(obstime  = tsamp,location = xyz))
t_rise_guess = tsamp[np.where(altaz.alt.value > 10)][0]
t_rise_guess

<Time object: scale='utc' format='iso' value=2046-01-30 13:32:13.838>

[bmorris3]

Hi, you're correct – if the reference time is very close to the rise/set time you're trying to find, you can "miss" the rise/set time that you're looking for with our simple grid search technique. This is sort of mentioned in the docs here but I would agree it could be made clearer.

Have you tried adjusting the keyword argument n_grid_points to a larger value, like n_grid_points=1000? This will slow down the algorithm but increase the precision in the search for the rise/set time.

[tstolarczyk]

Thanks for having a look into this problem.

Yes I have tried with n_grid_points=10000 and it does not change anything.
In fact the problem is not in the bin size of the first search, but simply in the algorithm in the function def _horiz_cross(self, t, alt, rise_set, horizon=0*u.degree) that intrinsically omits the first bin, and therefore skips any crossing in that bin.

It is clearly a bug, or at least an omitted feature.

My feeling is that the algorithm has been build under the assumption that only one crossing is possible in the investigated period, and it is not true !

With the example mentioned earlier I have checked that indeed in condition, the first of the n_grid_point=150 bins is True, as well as the last one (True means a crossing was detected).

From this, the indices of the intervals in the time array t is stored in before_indices and after_indices.
In this part the weird things start since the after_indices are computed only for the second True interval of the list

       before_indices = np.array(np.nonzero(condition))
        # we want to add an vector like (0, 1, ...) to get after indices
        after_indices = before_indices.copy()
        after_indices[1, :] += 1 # my comment : only the first element, whereas 2 are valid

one get :

before_indices
array([[  0,   0],
       [  0, 148],
       [  0,   0]], dtype=int64)
after_indices
array([[  0,   0],
       [  1, 149],
       [  0,   0]], dtype=int64)

whereas one would expect [0,1] for the first element if I understand the spirit of this action.
Anyhow, the altitude limits are obtained correctly (did not think too much why), this is not where we lose the crossing.
We have for the altitude (I have a chosen a 10 deg horizon), and this is fine :

• start :
  al1 <Quantity [9.76204762, 8.7489172 ] deg>
• stop :
  al2 <Quantity [11.44538854, 10.45077748] deg>

And we also get the twoarrays for the 2 time bins :

- tl1 <Time object: scale='utc' format='jd' value=[2468376.06308455 2468377.05637314]>
- tl2 <Time object: scale='utc' format='jd' value=[2468376.06979596 2468377.06308455]>

The drama happens at the following lines where the first crossing is dropped!
The function indeed requests only one crossing to be passed.
and the first one is omitted.

        alt_lims1[tuple(np.delete(before_indices, 1, 0))] = al1
        alt_lims2[tuple(np.delete(before_indices, 1, 0))] = al2
        jd_lims1[tuple(np.delete(before_indices, 1, 0))] = tl1.utc.jd
        jd_lims2[tuple(np.delete(before_indices, 1, 0))] = tl2.utc.jd

Sincerely, this code looks very complicated for what it does.
My guess (hope) is that it is due to find crossings in multiple cases (sources or whatever).
For this reason, I am not sure how to correct the function in the most general case.

What do you think?

[tstolarczyk]

Any progress on this issue? Thanks !

[tstolarczyk]

@bmorris3 kind reminder on this issue and my last post...