@@ -456,7 +456,12 @@ def julian_ephemeris_millennium(julian_ephemeris_century):
456456@jcompile (nopython = True )
457457def sum_mult_cos_add_mult (arr , x ):
458458 # shared calculation used for heliocentric longitude, latitude, and radius
459- s = (arr [:, [0 ]] * np .cos (arr [:, [1 ]] + arr [:, [2 ]] * np .expand_dims (x , axis = 0 ))).sum (axis = 0 )
459+ if USE_NUMBA :
460+ s = 0.
461+ for row in range (arr .shape [0 ]):
462+ s += arr [row , 0 ] * np .cos (arr [row , 1 ] + arr [row , 2 ] * x )
463+ else :
464+ s = (arr [:, [0 ]] * np .cos (arr [:, [1 ]] + arr [:, [2 ]] * np .expand_dims (x , axis = 0 ))).sum (axis = 0 )
460465 return s
461466
462467@jcompile ('float64(float64)' , nopython = True )
@@ -557,20 +562,41 @@ def moon_ascending_longitude(julian_ephemeris_century):
557562 nopython = True )
558563def longitude_obliquity_nutation (julian_ephemeris_century , x0 , x1 , x2 , x3 , x4 ,
559564 out ):
560- a = NUTATION_ABCD_ARRAY [:, [0 ]]
561- b = NUTATION_ABCD_ARRAY [:, [1 ]]
562- c = NUTATION_ABCD_ARRAY [:, [2 ]]
563- d = NUTATION_ABCD_ARRAY [:, [3 ]]
564- arg = np .radians (
565- NUTATION_YTERM_ARRAY [:, [0 ]]* np .expand_dims (x0 , axis = 0 ) +
566- NUTATION_YTERM_ARRAY [:, [1 ]]* np .expand_dims (x1 , axis = 0 ) +
567- NUTATION_YTERM_ARRAY [:, [2 ]]* np .expand_dims (x2 , axis = 0 ) +
568- NUTATION_YTERM_ARRAY [:, [3 ]]* np .expand_dims (x3 , axis = 0 ) +
569- NUTATION_YTERM_ARRAY [:, [4 ]]* np .expand_dims (x4 , axis = 0 )
570- )
571-
572- delta_psi_sum = ((a + b * julian_ephemeris_century ) * np .sin (arg )).sum (axis = 0 )
573- delta_eps_sum = ((c + d * julian_ephemeris_century ) * np .cos (arg )).sum (axis = 0 )
565+
566+ if USE_NUMBA :
567+ delta_psi_sum = 0.0
568+ delta_eps_sum = 0.0
569+ for row in range (NUTATION_YTERM_ARRAY .shape [0 ]):
570+ a = NUTATION_ABCD_ARRAY [row , 0 ]
571+ b = NUTATION_ABCD_ARRAY [row , 1 ]
572+ c = NUTATION_ABCD_ARRAY [row , 2 ]
573+ d = NUTATION_ABCD_ARRAY [row , 3 ]
574+ arg = np .radians (
575+ NUTATION_YTERM_ARRAY [row , 0 ]* x0 +
576+ NUTATION_YTERM_ARRAY [row , 1 ]* x1 +
577+ NUTATION_YTERM_ARRAY [row , 2 ]* x2 +
578+ NUTATION_YTERM_ARRAY [row , 3 ]* x3 +
579+ NUTATION_YTERM_ARRAY [row , 4 ]* x4
580+ )
581+ delta_psi_sum += (a + b * julian_ephemeris_century ) * np .sin (arg )
582+ delta_eps_sum += (c + d * julian_ephemeris_century ) * np .cos (arg )
583+
584+ else :
585+ a = NUTATION_ABCD_ARRAY [:, [0 ]]
586+ b = NUTATION_ABCD_ARRAY [:, [1 ]]
587+ c = NUTATION_ABCD_ARRAY [:, [2 ]]
588+ d = NUTATION_ABCD_ARRAY [:, [3 ]]
589+ arg = np .radians (
590+ NUTATION_YTERM_ARRAY [:, [0 ]]* np .expand_dims (x0 , axis = 0 ) +
591+ NUTATION_YTERM_ARRAY [:, [1 ]]* np .expand_dims (x1 , axis = 0 ) +
592+ NUTATION_YTERM_ARRAY [:, [2 ]]* np .expand_dims (x2 , axis = 0 ) +
593+ NUTATION_YTERM_ARRAY [:, [3 ]]* np .expand_dims (x3 , axis = 0 ) +
594+ NUTATION_YTERM_ARRAY [:, [4 ]]* np .expand_dims (x4 , axis = 0 )
595+ )
596+
597+ delta_psi_sum = ((a + b * julian_ephemeris_century ) * np .sin (arg )).sum (axis = 0 )
598+ delta_eps_sum = ((c + d * julian_ephemeris_century ) * np .cos (arg )).sum (axis = 0 )
599+
574600 delta_psi = delta_psi_sum * 1.0 / 36000000
575601 delta_eps = delta_eps_sum * 1.0 / 36000000
576602 # seems like we ought to be able to return a tuple here instead
0 commit comments