Snooping, how to calculate object coordinates and how to get correct time

[scriptorron]

This is not an idea, it is a call for help!

I saw that the INDI CCD simulator writes object coordinates (Alt/Az and J2000 Ra/Dec) as metadata in the FITS. This is a cool feature for documentation, stacking and as starting point for field solving. Therefore I tried to implement the same in this driver.

An INDI feature "Snooping" is used to get the required site coordinates and the observation direction. With Snooping the camera driver can ask other drivers for information. I recorded the network traffic and compared it with the metadata in the FITS. For some of the metadata I already know how to calculate them. But for others I have no idea:

Metadata in FITS:

OBJCTAZ =         3.599877E+02 / Azimuth of center of image in Degrees          
OBJCTALT=         5.104544E+01 / Altitude of center of image in Degrees         
OBJCTRA = '23 59 46.01'        / Object J2000 RA in Hours                       
OBJCTDEC= '89 52 19.09'        / Object J2000 DEC in Degrees                    
RA      =         3.599417E+02 / Object J2000 RA in Degrees                     
DEC     =         8.987197E+01 / Object J2000 DEC in Degrees                    
EQUINOX =                 2000 / Equinox                                        
DATE-OBS= '2023-03-30T13:13:51.616' / UTC start date of observation 

When this FITS picture was taken the CCD simulator snooped the following information over network:

• observation site: Longitude = 13.75°, Latitude=51.05°, Elevation=117.8m
• observation direction (equatorial coordinates of the mount): RA=20.659287020031786852 hours, DEC=90°

Local time was 15:13. The UTC time can be found in FITS metadata.

My questions and issues are:

1. How can I calculate the object coordinates? I already searched the internet for it and all explanations were very complicated. I do not want to reinvent the wheel. The calculations are likely not more than 4 lines of code when using Astropy.
   Okay, with a little bit more effort I certainly can find an answer to that. But the second issue is a killer:
2. How do I get date and time of observation? The Raspberry Pi has no battery buffered realtime clock. It updates its system clock from an NTP server in internet. But when I stay at night on a meadow far from civilization I will not have internet! It is possible to set date and time from an other Linux device. But this means additional effort for the users.
   The driver has no way to find out if the Raspberry Pi clock is correct or not. I have concerns to implement something that writes wrong information silently.

[sajmons]

@scriptorron

1. I'm using this library in one of my project. Maybe you can find there what you need? https://github.com/cosinekitty/astronomy
   Python version here: https://github.com/cosinekitty/astronomy/tree/master/source/python
2. I guess your best hope for correct time would be GPS dongle, I'm using one with indi_gpsd for my setup

Maybe you can ask the avtor of that library itself?

[scriptorron]

Interesting link. I will have a look on it. The advantage of astropy is that it is anyway used to create the FITS. I am sure it can also do such basic things. I was just too lazy for the manual.

A GPS dongle is a great idea. Likely, I will bye one. But that will not help other people. But I remember that INDI has a GPS simulator. That could be a way to "snoop" the time from a computer/Pi which has a battery buffered clock.

[scriptorron]

I found a solution for my 1st issue. As already guessed it needs only a few lines of code with astropy:

from astropy.coordinates import SkyCoord,EarthLocation, AltAz, ICRS
from astropy import units as u
from astropy.time import Time

def calc_FITS(
    ra, dec, UTC,
    lon=13.75*u.deg, lat=51.05*u.deg, height=117.8*u.meter,
):
    ObsLoc = EarthLocation(lon=lon, lat=lat, height=height)
    #
    c = SkyCoord(ra=ra, dec=dec, frame="icrs")
    cAltAz = c.transform_to(AltAz(obstime = Time(UTC), location = ObsLoc))
    #
    print(f'AIRMASS = {cAltAz.secz}')
    print(f'OBJCTAZ = {cAltAz.az/u.deg}')
    print(f'OBJCTALT= {cAltAz.alt/u.deg}')
    #
    J2000 = cAltAz.transform_to(ICRS())
    print(f'OBJCTRA = {J2000.ra.to_string(unit=u.hour)}')
    print(f'OBJCTDEC= {J2000.dec.to_string(unit=u.deg)}')
    print(f'RA      = {J2000.ra.degree:f}')
    print(f'DEC     = {J2000.dec.degree:f}')

For my first test image the driver sooped:

ra=20.659287020031786852*u.hourangle
dec=90*u.deg
UTC='2023-03-30T13:13:51.616'

The function above outputs:

AIRMASS = 1.2877214713687333
OBJCTAZ = 0.12024057854305253
OBJCTALT= 50.947177762067746
OBJCTRA = 18h37m39.1535s
OBJCTDEC= 90d00m00s
RA      = 279.413140
DEC     = 90.000000

while the header data in the FITS are:

AIRMASS =         1.285375E+00 / Airmass                                        
OBJCTAZ =         3.599877E+02 / Azimuth of center of image in Degrees          
OBJCTALT=         5.104544E+01 / Altitude of center of image in Degrees         
OBJCTRA = '23 59 46.01'        / Object J2000 RA in Hours                       
OBJCTDEC= '89 52 19.09'        / Object J2000 DEC in Degrees                    
RA      =         3.599417E+02 / Object J2000 RA in Degrees                     
DEC     =         8.987197E+01 / Object J2000 DEC in Degrees  

The large deviation in the calculated J2000 RA caused me some headache. But then I realized that it is the same point on sky: with DEC=90° it is the north celestial pole, independent what RA says.

Then I made a second test image more distant from the pole. Snooped values are:

ra=4.6240856854862251168*u.hourangle
dec=16.688656799999993297*u.deg
UTC='2023-04-05T11:27:53.655'

Calculated from that:

AIRMASS = 1.6468410143452783
OBJCTAZ = 111.65204332744206
OBJCTALT= 37.388988840152656
OBJCTRA = 4h37m26.7085s
OBJCTDEC= 16d41m19.1645s
RA      = 69.361285
DEC     = 16.688657

Compared to the metadata in the FITS:

AIRMASS =         1.643007E+00 / Airmass                                        
OBJCTAZ =         1.121091E+02 / Azimuth of center of image in Degrees          
OBJCTALT=         3.744145E+01 / Altitude of center of image in Degrees         
OBJCTRA = ' 4 36 07.37'        / Object J2000 RA in Hours                       
OBJCTDEC= '16 30 26.02'        / Object J2000 DEC in Degrees                    
RA      =         6.903072E+01 / Object J2000 RA in Degrees                     
DEC     =         1.650723E+01 / Object J2000 DEC in Degrees                    

The calculated coordinates and airmass still do not perfectly fit. There is about 0.5° difference in RA / AZ and about 0.2° difference in DEC / ALT. Possible reasons could be:

• Astropy and the CCD simulator calculate in different way and make different errors.
• The snooped data for the calculations above is not what the CCD simulator used. I picked these data out of recorded network traffic. It is possible that I picked the telescope coordinates a minute before the FITS was taken.

It will be tedious to find the reasons for these small differences. I will implement this calculation in the driver.