<blockquote>
<p>Discrepancies of 0.000147842 degrees in AZ, and 0.00002705 degrees in EL.</p>
</blockquote>
<p>Those discrepancies are very small; they are essentially the same numbers. There is one difference: your CSPICE calculations were made using DE440 while Horizons used DE441. From <a href="https://iopscience.iop.org/article/10.3847/1538-3881/abd414" rel="nofollow noreferrer">The JPL Planetary and Lunar Ephemerides DE440 and DE441</a>,</p>
<blockquote>
<p>The ephemerides DE440 and DE441 are fit to the same data set, but DE441 assumes no damping between the lunar liquid core and the solid mantle, which avoids a divergence when integrated backward in time. Therefore, DE441 is less accurate than DE440 for the current century, but covers a much longer duration of years −13,200 to +17,191, compared to DE440 covering years 1550–2650.</p>
</blockquote>
<p>Try using DE441 rather than DE440 in your CSPICE calculations to see if that reduces the discrepancy.</p>

CC BY-SA 4.0
2023-03-11T14:01:15.620Z
DE441 and DE440 yield the exact same results, down to at least 8 decimals (AZ/EL).
phil5
CC BY-SA 4.0
2023-03-11T15:40:23.420Z
"Very small" is a relative term.  The accuracy of both DE's are less than 1 meter.  The errors mentioned above are about .5 arcsec and .1 arcsec.  It should be possible to get several more decimals worth of accuracy.
Greg Miller
CC BY-SA 4.0
2023-03-11T23:08:23.060Z
@GregMiller I wonder if a very slightly different epoch, something like a few tens of milliseconds could realign the two positions? Is there any chance they interpret the input time differently, perhaps assuming different timescales? Is `obstim = "2003 OCT 13 06:00:00.000000 UTC";` what you used in both cases?
uhoh
CC BY-SA 4.0
2023-03-11T23:09:42.030Z
@phil5 ditto...
uhoh
CC BY-SA 4.0
2023-03-11T10:14:54.590Z
2023-03-11T10:14:54.590Z
David Hammen
<p>This CSPICE <code>spkcpo_c</code> routine <a href="https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/spkcpo_c.html#Examples" rel="nofollow noreferrer">example</a> outputs:</p>
<pre><code>AZ/EL computed without frame kernel:

       Solar azimuth (deg):           316.67141786
       Solar elevation (deg):         -54.85253216
</code></pre>
<p>On my machine, I've computed the same position with the new <code>azlcpo_c</code> routine:</p>
<pre><code>furnsh_c ( &quot;de440.bsp&quot; );
furnsh_c ( &quot;pck00010.tpc&quot; );
furnsh_c ( &quot;naif0012.tls&quot; );
furnsh_c ( &quot;earth_000101_230601_230308.bpc&quot; );

str2et_c ( &quot;2003-Oct-13 06:00:00 UTC&quot;, &amp;et );

obslat = 35.425901111 * rpd_c();
obslon = -116.889537582 * rpd_c();
obsalt = 1.00179621332;

bodvrd_c ( &quot;EARTH&quot;, &quot;RADII&quot;, 3, &amp;n, radii );

re = radii[0];
rp = radii[2];
f  = ( re - rp ) / re;

georec_c ( obslon, obslat, obsalt, re, f, obspos );

azlcpo_c ( &quot;ELLIPSOID&quot;, &quot;SUN&quot;, et, &quot;CN+S&quot;, SPICEFALSE, SPICETRUE,
           obspos, &quot;EARTH&quot;, &quot;ITRF93&quot;, azlsta, &amp;lt );

printf ( &quot;\n&quot;
         &quot; Solar azimuth (deg):   %20.8f\n&quot;
         &quot; Solar elevation (deg): %20.8f\n&quot;
         &quot;\n&quot;,
         azlsta[1] * dpr_c(), azlsta[2] * dpr_c() );
</code></pre>
<p>which outputs:</p>
<pre><code>Solar azimuth (deg):           316.67141785
Solar elevation (deg):         -54.85253206
</code></pre>
<p>So far, so good.</p>
<p>Here's Horizons System's output (query <a href="https://ssd.jpl.nasa.gov/api/horizons.api?format=text&amp;COMMAND=%2710%27&amp;OBJ_DATA=%27NO%27&amp;MAKE_EPHEM=%27YES%27&amp;EPHEM_TYPE=%27OBSERVER%27&amp;CENTER=%27coord@399%27&amp;COORD_TYPE=%27GEODETIC%27&amp;SITE_COORD=%27-116.889537582,35.425901111,1.00179621332%27&amp;START_TIME=%272003-10-13%2006:00:00%20UT%27&amp;STOP_TIME=%272003-10-13%2006:00:01%27&amp;STEP_SIZE=%271%20d%27&amp;QUANTITIES=%274%27&amp;ANG_FORMAT=%27DEG%27&amp;EXTRA_PREC=%27YES%27" rel="nofollow noreferrer">URL</a>):</p>
<pre><code>*******************************************************************************
Ephemeris / WWW_USER Fri Mar 10 08:31:14 2023 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Sun (10)                        {source: DE441}
Center body name: Earth (399)                     {source: DE441}
Center-site name: (user defined site below)
*******************************************************************************
Start time      : A.D. 2003-Oct-13 06:00:00.0000 UT      
Stop  time      : A.D. 2003-Oct-13 06:00:01.0000 UT      
Step-size       : 1440 minutes
*******************************************************************************
Target pole/equ : IAU_SUN                         {East-longitude positive}
Target radii    : 696000.0, 696000.0, 696000.0 km {Equator_a, b, pole_c}       
Center geodetic : 243.1104624, 35.4259011, 1.0018 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 243.1104624,5203.99728,3677.05257 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : ITRF93                          {East-longitude positive}
Center radii    : 6378.137, 6378.137, 6356.752 km {Equator_a, b, pole_c}       
Target primary  : Sun
Vis. interferer : MOON (R_eq= 1737.400) km        {source: DE441}
Rel. light bend : Sun                             {source: DE441}
Rel. lght bnd GM: 1.3271E+11 km^3/s^2                                          
Atmos refraction: NO (AIRLESS)
RA format       : DEG
Time format     : BOTH
Calendar mode   : Mixed Julian/Gregorian
EOP file        : eop.230308.p230601                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2023-MAR-08. PREDICTS-&gt; 2023-MAY-31
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass (&gt;38.000=NO), Daylight (NO )
Table cut-offs 2: Solar elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table cut-offs 3: RA/DEC angular rate (     0.0=NO )                           
*******************************************************************************
 Date__(UT)__HR:MN Date_________JDUT     Azimuth__(a-app)__Elevation
********************************************************************
$$SOE
 2003-Oct-13 06:00 2452925.750000000  m  316.671270018 -54.852559210
$$EOE
*******************************************************************************
</code></pre>
<p>Discrepancies of 0.000147842 degrees in AZ, and 0.00002705 degrees in EL.</p>
<p>That seems like a significant difference. As a comparison, a 7~8 seconds difference in delta T (TT - UT) will result in similar discrepancies; equivalent to the Earth moving over 200 km around its orbit.</p>
<p>What am I missing?</p>

CC BY-SA 4.0
2023-03-10T23:50:46.333Z
It's been a while since I looked at SPICE.  Am I reading it right that you are not using the built in Earth model?  As there are warnings about its accuracy.  But I agree, that is a much larger discrepancy than I'd expect.
Greg Miller
CC BY-SA 4.0
2023-03-11T00:49:19.653Z
@GregMiller AFAIK, I'm using the ITRF93 high accuracy Earth rotation model (binary file). The model to avoid is IAU_EARTH, provided in text a file.
phil5
CC BY-SA 4.0
2023-03-11T10:24:49.360Z
@phil5 ITRF93 is not exactly "high accuracy", and the JPL reimplementation of it in terms of Chebyshev polynomials is even less accurate. People who do sub-arcsecond astronomy use more recent models from the IAU.
David Hammen
CC BY-SA 4.0
2023-03-11T14:31:18.940Z
@DavidHammen I was just using the official wording by [NAIF](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/Tutorials/pdf/individual_docs/23_lunar-earth_pck-fk.pdf). I'm more interested in knowing the source of error, given that Horizons System and SPICE supposedly use the same EOP files and planetary ephemeris.
phil5
CC BY-SA 4.0
2023-03-11T15:47:36.653Z
Regarding "a 7~8 seconds difference in delta T (TT - UT) will result in similar discrepancies; equivalent to the Earth moving over 200 km around its orbit": You are forgetting the Earth's one revolution per sidereal day rotation. Even a small time discrepancy in DUT1 can result in a significant az-el discrepancy. That said, I would not consider a half an arcsecond discrepancy between Horizons and CSPICE "significant". Lesson learned: Don't use CSPICE or Horizons for milliarcsecond astronomy.
David Hammen
CC BY-SA 4.0
2023-03-11T15:54:35.610Z
You can compare other data that Horizons can generate to see where the discrepincy creeps in.  Eg j2000 ra/dec, (of date) ra/dec, sidereal time, etc.
Greg Miller
CC BY-SA 4.0
2023-03-13T20:29:47.940Z
@GregMiller The rectangular coordinates of the target (Sun) with respect to the observer (DSS-14) are accurate to within 1 meter, which means that the rotation model and timing are OK. The error seems to either come from the transformation from geocentric (ITRF93) to horizontal, or the conversion from rectangular to spherical coordinates. I've implemented my own matrix transformations and the result coincides with SPICE exactly.
phil5
CC BY-SA 4.0
2023-03-13T20:34:56.327Z
@GregMiller Another thing making me suspicious of the Horizons results, is that the higher the observer's latitude on the surface of the Earth (near ±90 degrees), the bigger the discrepancy. Perhaps they are using some buggy implementation of arctan, but that's only a wild guess.
phil5
CC BY-SA 4.0
2023-03-13T21:32:08.630Z
@phil5 or possibly a different shape for the Earth's surface and local vertical. If the discrepancy is in alt/az but not RA/Dec, that might be indicative of something along those lines
uhoh
CC BY-SA 4.0
2023-03-13T22:10:17.247Z
@uhoh SPICE uses a unit vector in the direction of the observer's planetographic latitude and longtiude on the Earth's reference triaxial ellipsoid as a local normal. I guess it's possible that Horizons uses a different shape and slightly adjusts the normal vector.
phil5
CC BY-SA 4.0
2023-03-13T22:16:54.240Z
@uhoh However I don't see how you would get the same rectangular coordinates. I admit I don't have a strongest intuition for this problem.
phil5
CC BY-SA 4.0
2023-03-14T17:28:34.180Z
@uhoh Converting from [Earth-centered, Earth-fixed coordinates](https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_geodetic_coordinates) (ECEF) to geodetic coordinates is actually a non-trivial problem. It can be sensitive to small errors and in general requires more complicated numerical methods. This seems a more likely scenario for error to creep in, rather than Horizons aligning to the Earth geoid, or buggy math functions. (I'm slightly embarrassed of having considered it.)
phil5
CC BY-SA 4.0
2023-03-14T18:42:26.913Z
...and after some testing, I don't think that's it either. Even the most rudimentary iteration finds the correct geodetic coordinates.
phil5
<coordinate><positional-astronomy><ephemeris><algorithm><azimuth>
CSPICE vs JPL Horizons discrepancy in AZ/EL
phil5
# CSPICE vs JPL Horizons discrepancy in AZ/EL

(Asked by: username_0 on 2023-03-10)

This CSPICE spkcpo_c routine example outputs:
AZ/EL computed without frame kernel:

       Solar azimuth (deg):           316.67141786
       Solar elevation (deg):         -54.85253216

On my machine, I've computed the same position with the new azlcpo_c routine:
furnsh_c ( "de440.bsp" );
furnsh_c ( "pck00010.tpc" );
furnsh_c ( "naif0012.tls" );
furnsh_c ( "earth_000101_230601_230308.bpc" );

str2et_c ( "2003-Oct-13 06:00:00 UTC", &et );

obslat = 35.425901111 * rpd_c();
obslon = -116.889537582 * rpd_c();
obsalt = 1.00179621332;

bodvrd_c ( "EARTH", "RADII", 3, &n, radii );

re = radii[0];
rp = radii[2];
f  = ( re - rp ) / re;

georec_c ( obslon, obslat, obsalt, re, f, obspos );

azlcpo_c ( "ELLIPSOID", "SUN", et, "CN+S", SPICEFALSE, SPICETRUE,
           obspos, "EARTH", "ITRF93", azlsta, < );

printf ( "\n"
         " Solar azimuth (deg):   %20.8f\n"
         " Solar elevation (deg): %20.8f\n"
         "\n",
         azlsta[1] * dpr_c(), azlsta[2] * dpr_c() );

which outputs:
Solar azimuth (deg):           316.67141785
Solar elevation (deg):         -54.85253206

So far, so good.
Here's Horizons System's output (query URL):
*******************************************************************************
Ephemeris / WWW_USER Fri Mar 10 08:31:14 2023 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Sun (10)                        {source: DE441}
Center body name: Earth (399)                     {source: DE441}
Center-site name: (user defined site below)
*******************************************************************************
Start time      : A.D. 2003-Oct-13 06:00:00.0000 UT      
Stop  time      : A.D. 2003-Oct-13 06:00:01.0000 UT      
Step-size       : 1440 minutes
*******************************************************************************
Target pole/equ : IAU_SUN                         {East-longitude positive}
Target radii    : 696000.0, 696000.0, 696000.0 km {Equator_a, b, pole_c}       
Center geodetic : 243.1104624, 35.4259011, 1.0018 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 243.1104624,5203.99728,3677.05257 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : ITRF93                          {East-longitude positive}
Center radii    : 6378.137, 6378.137, 6356.752 km {Equator_a, b, pole_c}       
Target primary  : Sun
Vis. interferer : MOON (R_eq= 1737.400) km        {source: DE441}
Rel. light bend : Sun                             {source: DE441}
Rel. lght bnd GM: 1.3271E+11 km^3/s^2                                          
Atmos refraction: NO (AIRLESS)
RA format       : DEG
Time format     : BOTH
Calendar mode   : Mixed Julian/Gregorian
EOP file        : eop.230308.p230601                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2023-MAR-08. PREDICTS-> 2023-MAY-31
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass (>38.000=NO), Daylight (NO )
Table cut-offs 2: Solar elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table cut-offs 3: RA/DEC angular rate (     0.0=NO )                           
*******************************************************************************
 Date__(UT)__HR:MN Date_________JDUT     Azimuth__(a-app)__Elevation
********************************************************************
$$SOE
 2003-Oct-13 06:00 2452925.750000000  m  316.671270018 -54.852559210
$$EOE
*******************************************************************************

Discrepancies of 0.000147842 degrees in AZ, and 0.00002705 degrees in EL.
That seems like a significant difference. As a comparison, a 7~8 seconds difference in delta T (TT - UT) will result in similar discrepancies; equivalent to the Earth moving over 200 km around its orbit.
What am I missing?

**Tags:** coordinate, positional-astronomy, ephemeris, algorithm, azimuth

> **Comments:**
> **username_1**: It's been a while since I looked at SPICE.  Am I reading it right that you are not using the built in Earth model?  As there are warnings about its accuracy.  But I agree, that is a much larger discrepancy than I'd expect.  
> **username_0**: @GregMiller AFAIK, I'm using the ITRF93 high accuracy Earth rotation model (binary file). The model to avoid is IAU_EARTH, provided in text a file.  
> **username_2**: @username_0 ITRF93 is not exactly "high accuracy", and the JPL reimplementation of it in terms of Chebyshev polynomials is even less accurate. People who do sub-arcsecond astronomy use more recent models from the IAU.  
> **username_0**: @DavidHammen I was just using the official wording by [NAIF](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/Tutorials/pdf/individual_docs/23_lunar-earth_pck-fk.pdf). I'm more interested in knowing the source of error, given that Horizons System and SPICE supposedly use the same EOP files and planetary ephemeris.  
> **username_2**: Regarding "a 7~8 seconds difference in delta T (TT - UT) will result in similar discrepancies; equivalent to the Earth moving over 200 km around its orbit": You are forgetting the Earth's one revolution per sidereal day rotation. Even a small time discrepancy in DUT1 can result in a significant az-el discrepancy. That said, I would not consider a half an arcsecond discrepancy between Horizons and CSPICE "significant". Lesson learned: Don't use CSPICE or Horizons for milliarcsecond astronomy.  
> **username_1**: You can compare other data that Horizons can generate to see where the discrepincy creeps in.  Eg j2000 ra/dec, (of date) ra/dec, sidereal time, etc.  
> **username_0**: @GregMiller The rectangular coordinates of the target (Sun) with respect to the observer (DSS-14) are accurate to within 1 meter, which means that the rotation model and timing are OK. The error seems to either come from the transformation from geocentric (ITRF93) to horizontal, or the conversion from rectangular to spherical coordinates. I've implemented my own matrix transformations and the result coincides with SPICE exactly.  
> **username_0**: @GregMiller Another thing making me suspicious of the Horizons results, is that the higher the observer's latitude on the surface of the Earth (near ±90 degrees), the bigger the discrepancy. Perhaps they are using some buggy implementation of arctan, but that's only a wild guess.  
> **username_3**: @username_0 or possibly a different shape for the Earth's surface and local vertical. If the discrepancy is in alt/az but not RA/Dec, that might be indicative of something along those lines  
> **username_0**: @username_3 SPICE uses a unit vector in the direction of the observer's planetographic latitude and longtiude on the Earth's reference triaxial ellipsoid as a local normal. I guess it's possible that Horizons uses a different shape and slightly adjusts the normal vector.  
> **username_0**: @username_3 However I don't see how you would get the same rectangular coordinates. I admit I don't have a strongest intuition for this problem.  
> **username_0**: @username_3 Converting from [Earth-centered, Earth-fixed coordinates](https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_geodetic_coordinates) (ECEF) to geodetic coordinates is actually a non-trivial problem. It can be sensitive to small errors and in general requires more complicated numerical methods. This seems a more likely scenario for error to creep in, rather than Horizons aligning to the Earth geoid, or buggy math functions. (I'm slightly embarrassed of having considered it.)  
> **username_0**: ...and after some testing, I don't think that's it either. Even the most rudimentary iteration finds the correct geodetic coordinates.  

---

## 1 answers

### Answer by username_2 on 2023-03-11. Score: 1

> Discrepancies of 0.000147842 degrees in AZ, and 0.00002705 degrees in EL.

Those discrepancies are very small; they are essentially the same numbers. There is one difference: your CSPICE calculations were made using DE440 while Horizons used DE441. From [The JPL Planetary and Lunar Ephemerides DE440 and DE441](https://iopscience.iop.org/article/10.3847/1538-3881/abd414),

> The ephemerides DE440 and DE441 are fit to the same data set, but DE441 assumes no damping between the lunar liquid core and the solid mantle, which avoids a divergence when integrated backward in time. Therefore, DE441 is less accurate than DE440 for the current century, but covers a much longer duration of years −13,200 to +17,191, compared to DE440 covering years 1550–2650.

Try using DE441 rather than DE440 in your CSPICE calculations to see if that reduces the discrepancy.

> **Comments:**
> **username_0**: DE441 and DE440 yield the exact same results, down to at least 8 decimals (AZ/EL).  
> **username_1**: "Very small" is a relative term.  The accuracy of both DE's are less than 1 meter.  The errors mentioned above are about .5 arcsec and .1 arcsec.  It should be possible to get several more decimals worth of accuracy.  
> **username_3**: @GregMiller I wonder if a very slightly different epoch, something like a few tens of milliseconds could realign the two positions? Is there any chance they interpret the input time differently, perhaps assuming different timescales? Is `obstim = "2003 OCT 13 06:00:00.000000 UTC";` what you used in both cases?  
> **username_3**: @username_0 ditto...  

---

