(ANSWERED) Daily aspectarian - unexpected error :(

[gipsea]

Dear all,
following the latest discussion I had here regarding zodiac entry and find_minima(), I'm now studing the daily aspects between 2 planets during a period of time. This is my function I though should work:

def sin_30(ephemeris, base ,t1,t2):
    """aspectarian modulus 30degree or pi/6 (7 occurrencies in total
     t1 =ephemeris(target1)
     t2 =ephemeris(target2)
    """
    if base == ephemeris['sun']:
        sun_at = base.at
        def sin_30_at(t):
            s = sun_at(t)
            _, t1lon, _ = s.observe(t1).ecliptic_latlon()
            _, t2lon, _ = s.observe(t2).ecliptic_latlon()

 #delta is the long difference which range is 0 - pi
            if abs(t1lon.radians-t2lon.radians)<=pi:
                delta = abs(t1lon.radians-t2lon.radians)
            else:
                delta = abs (min(t1lon.radians,t2lon.radians)+360-max(t1lon.radians,t2lon.radians))
            return (delta// (pi/ 7) % 7).astype(int)
    else:
         earth_at = base.at
         def sin_30_at(t):
             e = earth_at(t)
             _, t1lon, _ = e.observe(t1).apparent().ecliptic_latlon
             _, t2lon, _ = e.observe(t2).apparent().ecliptic_latlon
             if abs(t1lon.radians-t2lon.radians)<=pi:
                delta = abs(t1lon.radians-t2lon.radians)
             else:
                delta = abs (min(t1lon.radians,t2lon.radians)+360-max(t1lon.radians,t2lon.radians))
             return (delta// (pi / 7) % 7).astype(int)
    sin_30_at.step_days = 1.0
    return sin_30_at

At the moment I receive this error:

File "/home/mauro/anaconda3/lib/python3.8/site-packages/skyfield/almanac.py", line 373, in sin_30_at
if abs(t1lon.radians-t2lon.radians)<=pi:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I was expecting _, t1lon, _ = e.observe(t1).apparent().ecliptic_latlon would return one number but it is returning an array whose size depend on time span.

Being for some reason an array rather a single value I receive the above problem.

Can someone suggest why I receive an array insted of a single number and how to resolve the error?

Cheers

[brandon-rhodes]

I was expecting _, t1lon, _ = e.observe(t1).apparent().ecliptic_latlon would return one number but it is returning an array whose size depend on time span.

The array t1lon should have exactly as many items in it as t1 has times. Maybe you could print(len(t1)) and print(t1lon.shape) and verify whether you are getting back the name number of longitudes as the number of times you are asking about?

[gipsea]

Hi @brandon-rhodes , thanks for your time.

Since last time I posted the above code I moved on trying to study other features and I always got stuck somehow which ring a bell I might have misinterpretate some of skyfiled code.
The feature I'm studying now are all based more or less on apparent position and using find_minima or find_maxima function to find something like Moon apogee/perigee, when a planet change from direct to retrograde and the study of synodic aspects. The above function was actually a 'difference' of longitudes found by instance by the usual function
ra, decl, dist = earth.at(t).observe(moon).apparent.radec()

Below there is the simplest code I can think of to do some tests:

from skyfield.api import load
import datetime
from skyfield.searchlib import find_discrete, find_minima, find_maxima

ts = load.timescale()
planets = load('de421.bsp')

sun = planets['sun']
earth = planets['earth']
moon = planets['moon']

data_inizio = datetime.date(year=2021,month=5, day=1)
data_fine = datetime.date(year=2021,month=5,day=10)
t1 = ts.utc(data_inizio)
t2 = ts.utc(data_fine)

def apogee(t):
    e= earth.at(t)
    apparent = e.observe(moon).apparent()
    _, _,dist = apparent.radec()
    print("Distance inside function")
    print(dist)
    return dist
apogee.step_days=1
print("Distance outside function")
print(apogee(t1))
t, values = find_maxima(t1, t2,apogee)
for ti,vi in zip(t,values):
    print(ti.utc_strftime('%Y-%m-%d %H:%M:%S '), '%.2f' % vi,
          'Distance moon - Earth')

this is the output on spider console:

Distance outside function
Distance inside function
0.00245349 au
0.00245349 au
Distance inside function
[0.00242968 0.00245547 0.00248475 0.00251583 0.00254718 0.00257748
 0.00260572 0.0026312  0.00265346 0.00267232 0.00268772 0.00269975
 0.00270857 0.0027143 ] au
Traceback (most recent call last):

  File "/home/mauro/Documents/ephem/untitled0.py", line 35, in <module>
    t, values = find_maxima(t1, t2,apogee)

  File "/home/mauro/anaconda3/lib/python3.8/site-packages/skyfield/searchlib.py", line 163, in find_maxima
    left, right = _choose_brackets(y)

  File "/home/mauro/anaconda3/lib/python3.8/site-packages/skyfield/searchlib.py", line 183, in _choose_brackets
    dsd = diff(sign(diff(y)))

  File "<__array_function__ internals>", line 5, in diff

  File "/home/mauro/anaconda3/lib/python3.8/site-packages/numpy/lib/function_base.py", line 1247, in diff
    raise ValueError("diff requires input that is at least one dimensional")

ValueError: diff requires input that is at least one dimensional

I definetively use something incorrectly and I can't get my head around. How would any of you calculate by instance the max and/or min distance of the moon from the Earth? Would you use find_maxima() and find_minima()?

thank you all

[gipsea]

As usual...as soon I I post something the answer come along.

for the above function I think the error was:
return dist
where the distance unit was not specified...

At the momen with this function working I'm going ahead with most of my missing feature.

Soon I'll be back to the first function posted and return possibly with a solution:)

Cheers

[brandon-rhodes]

I'm glad you worked out what the problem was, even ahead of a Skyfield contributor having time to try out your script! Good luck as you proceed.

[gipsea]

Sorry everybody but I'm back to the same post...
I'm developing what I think being similar function and the finding about last function I thought was valid for the sin_30 function posted at the beginning of this thread but I was wrong.
I prepared the below script for anyone who has the change to play with it. At least now I'm down to were 'apparently' this gets stuck..

from skyfield.api import load
import datetime
from skyfield.searchlib import find_discrete
from skyfield.framelib import ecliptic_frame as ef

ts = load.timescale()
planets = load('de421.bsp')

sun = planets['sun']
earth = planets['earth']
mer = planets['mercury']
ven = planets['venus']

data_inizio = datetime.date(year=2021,month=5, day=1)
data_fine = datetime.date(year=2021,month=5,day=30)
t1 = ts.utc(data_inizio)
t2 = ts.utc(data_fine)

sun_at = sun.at
def sinodici(t):
    t1long=t2long= None
    s = sun_at(t)
    t1lat, t1long, t1dist = s.observe(mer).frame_latlon(ef)
    t1l = t1long.degrees
    print("t1long within function:" ,t1l)
    t2lat, t2long,t2dist = s.observe(ven).frame_latlon(ef)
    t2l =t2long.degrees
    print("t2long within function: ",t2l)
    delta = abs(t1l-t2l)   
    #if abs(t1l-t2l)<=180:
    #   delta = abs(t1l-t2l)
    #else:
    #    delta = abs (min(t1l,t2l)+360-max(t1l,t2l))

    return (delta// (30) % 7).astype(int)
sinodici.step_days=1.0
    
time,value= find_discrete(t1, t2, sinodici)

for ti,vi in zip(time,value):
    print('%s %.4f' % (ti.utc_strftime('%b %d'), vi))   

This function as it is works but I actualy need to use the commented code:

    if abs(t1l-t2l)<=180:
         delta = abs(t1l-t2l)
    else:
        delta = abs (min(t1l,t2l)+360-max(t1l,t2l))

If I uncomment it I receive this error:

File "/home/mauro/Documents/ephem/delete_moon.py", line 57, in
time,value= find_discrete(t1, t2, sinodici)

File "/home/mauro/anaconda3/lib/python3.8/site-packages/skyfield/searchlib.py", line 40, in find_discrete
return _find_discrete(ts, jd, f, epsilon, num)

File "/home/mauro/anaconda3/lib/python3.8/site-packages/skyfield/searchlib.py", line 52, in _find_discrete
y = f(t)

File "/home/mauro/Documents/ephem/delete_moon.py", line 41, in sinodici
if abs(t1l-t2l)<=180:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

showing that everything stops somehow with the if statement.
Does anybody has some consideration about this issue? I'm open to change approach to find the sinodic angles if this help.

I know this post is a little confused and I can split/re-post it possibly to reduce the confusion

Thank you all

[gipsea]

little update:
with the last code printing at screen different info I realized the find_discrete function pass an array of time to the called function whilst I thought itwould pass one specific time at the time. This was the origin of my confusion.
After spendig few hours I manipulate the the function sinodici

def sinodici(t)
    delta=np.zeros(len(t))
    s = sun_at(t)
    _, t1long, _ = s.observe(mer).frame_latlon(ef)
    t1l = t1long.degrees
    _, t2long,_ = s.observe(ven).frame_latlon(ef)
    t2l =t2long.degrees
    for i in range(len(t)):
        d=abs(t1l[i]-t2l[i])
        if d>180:
            delta[i]=abs(min(t1l[i],t2l[i])+360-max(t1l[i],t2l[i]))
        else:
            delta[i]=d
    return ((delta//30) % 7).astype(int)

The above function in conjunction with find_discrete function has 2 limitations:

1. doesn't find conjunction and opposition aspects
2. when the aspect moves from opposition back to conjunction it fails to display the correct aspect because degrees are moving backward from 180 to 0

For the time being I would say I found a way to bend the code the way I need with few workaround but if anybody has a better idea how to find aspect with specific modulus (30º in this case) is welcome to jump in

[brandon-rhodes]

You might find the NumPy unwrap() function useful — while it requires your units to be radians, not degrees, which is a little annoying, it edits out any 360° jumps in an angle. I often have to use it in almanac functions so that it doesn't look like one planet is passing another every time one of them goes 359°→0°.

[gipsea]

Hi @brandon-rhodes ,thanks for your valuable help. I'll take a look at the function and see if I can use in the future when I'll revise my code. I believe your suggestion was to help finding conjunction and opposition together with the other aspets...At the only for these 2 apsect I integrated my code starting from moon phase conj and opposition and adjusted for my need.