<jupyter_start><jupyter_text>Fink case study: Ephemerides for Solar System ScienceThis notebook shows how to use the Miriade IMCCE service with Fink Solar System Object data.Authors: <NAME>, <NAME>, <NAME>.<jupyter_code>import requests
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from sbpy.photometry import HG, HG1G2, HG12, HG12_Pen16
from astropy.modeling.fitting import LevMarLSQFitter
from sbpy.data import Obs
import astropy.units as u

from scipy.optimize import curve_fit

import seaborn as sns
sns.set_context('talk')

import time

from astropy.coordinates import SkyCoord
import astropy.units as u
# Define a simple ephemerides query
def miriade(ident, jd, observer='I41', rplane='1', tcoor=5):
    '''Gets asteroid ephemerides from IMCCE Miriade for a suite of JD for a single SSO
    Original function by <NAME>

    :ident: int, float, str - asteroid identifier
    :jd: array - dates to query
    :observer: str - IAU Obs code - default to ZTF: https://minorplanetcenter.net//iau/lists/ObsCodesF.html
    :returns: pd.DataFrame - Input dataframe with ephemerides columns appended
              False - If query failed somehow

    '''
    
    # ------
    # Miriade URL 
    url = 'https://ssp.imcce.fr/webservices/miriade/api/ephemcc.php'
    
    if rplane=='2':
        tcoor='1'
        
    # Query parameters
    params = {
        '-name': f'a:{ident}',
        '-mime': 'json',
        '-rplane': rplane,
        '-tcoor': tcoor,
        '-output': '--jd,--colors(SDSS:r,SDSS:g)',
        '-observer': observer, 
        '-tscale': 'UTC'
    }
    
    # Pass sorted list of epochs to speed up query
    files = {'epochs': ('epochs', '\n'.join(['%.6f' % epoch
                                             for epoch in jd]))}
    # Execute query
    try:
        r = requests.post(url, params=params, files=files, timeout=2000)
    except requests.exceptions.ReadTimeout:
        return False

    j = r.json()

    # Read JSON response
    try:
        ephem = pd.DataFrame.from_dict(j['data'])
    except KeyError:
        return False

    return ephem
def lin_fit(x, a, b):
    """ Linear function ax + b
    """
    return a*x + b

def sine_fit(x, a, b):
    """ Sinusoidal function a*sin( 2*(x-b) )
    :x: float - in degrees
    :a: float - Amplitude
    :b: float - Phase offset

    """
    return a*np.sin( 2*np.radians(x-b) )<jupyter_output><empty_output><jupyter_text>Let's define some constants:<jupyter_code># SSO id
ssname = '8467'
#'1041'
# Disctionary for filters
filters = {1:'g', 2:'R', 3:'i'}<jupyter_output><empty_output><jupyter_text>Fink data<jupyter_code># Get data from FINK 
r = requests.post(
  'http://134.158.75.151:24000/api/v1/sso',
  json={
    'n_or_d': ssname,
    'output-format': 'json'
  }
)

# Sort by epoch - critical for ephemerides speed up
sso = pd.read_json(r.content)
sso = sso.sort_values(by='i:jd').reset_index(drop=True)
# Select only a few columns
cols = [
    'd:cdsxmatch', 'd:rf_kn_vs_nonkn', 'd:mulens',
    'd:rf_snia_vs_nonia', 'd:roid', 'd:snn_sn_vs_all', 'd:snn_snia_vs_nonia',
    'd:tracklet', 'i:jd', 'i:ra', 'i:dec', 'i:fid', 'i:magpsf', 'i:sigmapsf',
    'i:magnr', 'i:sigmagnr', 'i:isdiffpos', 'i:magzpsci'
]

fink = sso[cols]<jupyter_output><empty_output><jupyter_text>Ephemerides<jupyter_code># Query ephemerides
start_t = time.time()
eph = miriade(ssname, fink['i:jd'], observer='I41' )
end_t = time.time()
print(f"This took {end_t - start_t:.3} seconds.")

sc = SkyCoord(eph['RA'], eph['DEC'], unit=(u.deg, u.deg) )

eph = eph.drop(columns=['RA','DEC'])
eph['RA'] = sc.ra.value*15
eph['Dec'] = sc.dec.value

# Add Ecliptic coordinates
eph_ec = miriade(ssname, fink['i:jd'], rplane='2' )

sc = SkyCoord(eph_ec['Longitude'], eph_ec['Latitude'], unit=(u.deg, u.deg) )
eph['Longitude'] = sc.ra.value
eph['Latitude'] = sc.dec.value


# Merge fink & Eph
info = pd.concat( [eph, fink], axis=1 )

# Compute magnitude reduced to unit distance
info['i:magpsf_red'] = info['i:magpsf'] - 5*np.log10( info['Dobs']*info['Dhelio'] )<jupyter_output>This took 1.78 seconds.<jupyter_text>Magnitude comparison<jupyter_code># Simple plot
filts = sso['i:fid'].unique()

fig, ax = plt.subplots(
    3, 1, 
    figsize=(12, 10), 
#     sharex=True, 
    gridspec_kw={'height_ratios': [3, 1.5, 1.5]}
)

colors = ['#15284F', '#F5622E']
for i,f in enumerate(filts):
    cond = sso['i:fid']==f
    
    ax[0].scatter(
        sso.loc[cond,'i:jd'], 
        sso.loc[cond,'i:magpsf'], 
        label='ZTF {:}'.format(filters[f]),
        color=colors[i],
        alpha=1.0
    )

    ax[0].plot(
        eph.Date[cond], 
        eph['SDSS:{}'.format(filters[f].lower())][cond], 
        ls='', 
        marker='s', 
        label='SDSS:{} from ephem'.format(filters[f].lower()), alpha=0.2)

ax[0].set_ylim(ax[0].get_ylim()[::-1])
ax[0].legend(loc='best')
ax[0].set_ylabel('Magnitude')
ax[0].set_xlabel('Epoch (JD)');

for i,f in enumerate(filts):
    cond = sso['i:fid']==f
    
    # Magnitude
    mag = sso.loc[cond,'i:magpsf'][cond]
    err= sso.loc[cond,'i:sigmapsf'][cond]
    
    ephemData = eph['SDSS:{}'.format(filters[f].lower())][cond]

    diff = mag - ephemData
    
    popt, pcov = curve_fit(sine_fit, eph.Longitude[cond], diff)
    perr = np.sqrt(np.diag(pcov))
    print(popt, perr)
    
    ax[i+1].errorbar(
        eph.Longitude[cond], 
        diff, 
        yerr=err, 
        ls='', marker='.', label=r'$\Delta$ {}'.format(filters[f].lower()), color=colors[i])
    
    longitude_arr = np.linspace(0,360,num=100)
    ax[i+1].plot(
        longitude_arr, 
        sine_fit(longitude_arr, *popt),
        label='Sine fit: a={:.0E}, b={:.0E}'.format(*popt),
        color=colors[i],
        ls='-',
        alpha=0.5,
        lw=3
    )
    
    ax[i+1].legend(loc='best', ncol=2)

    ax[i+1].set_ylim(-0.5, 0.5)
    ax[i+1].set_xlim(0,360)
    ax[i+1].axhline(0, ls='--', color='grey')
    ax[i+1].set_ylabel('residuals')
ax[2].set_xlabel('Ecliptic longitude ($^o$)');
# Simple plot
filts = info['i:fid'].unique()

fig, ax = plt.subplots(
    3, 1, 
    figsize=(12, 10), 
#     sharex=True, 
    gridspec_kw={'height_ratios': [3, 1.5, 1.5]}
)

colors = ['#15284F', '#F5622E']
for i,f in enumerate(filts):
    cond = info['i:fid']==f
    
    ax[0].scatter(
        info.loc[cond,'i:jd'], 
        info.loc[cond,'i:magpsf'], 
        label='ZTF {:}'.format(filters[f]),
        color=colors[i],
        alpha=1.0
    )

    ax[0].plot(
        info.Date[cond], 
        info['SDSS:{}'.format(filters[f].lower())][cond], 
        ls='', 
        marker='s', 
        label='SDSS:{} from ephem'.format(filters[f].lower()), alpha=0.2)

ax[0].set_ylim(ax[0].get_ylim()[::-1])
ax[0].legend(loc='best')
ax[0].set_ylabel('Magnitude')
ax[0].set_xlabel('Epoch (JD)');

for i,f in enumerate(filts):
    cond = info['i:fid']==f
    
    # Magnitude
    mag = info.loc[cond,'i:magpsf'][cond]
    err= info.loc[cond,'i:sigmapsf'][cond]
    
    ephemData = info['SDSS:{}'.format(filters[f].lower())][cond]

    diff = mag - ephemData
    
    popt, pcov = curve_fit(sine_fit, info.Longitude[cond], diff)
    perr = np.sqrt(np.diag(pcov))
    print(popt, perr)
    
    ax[i+1].errorbar(
        info.Longitude[cond], 
        diff, 
        yerr=err, 
        ls='', marker='.', label=r'$\Delta$ {}'.format(filters[f].lower()), color=colors[i])
    
    longitude_arr = np.linspace(0,360,num=100)
    ax[i+1].plot(
        longitude_arr, 
        sine_fit(longitude_arr, *popt),
        label='Sine fit: a={:.0E}, b={:.0E}'.format(*popt),
        color=colors[i],
        ls='-',
        alpha=0.5,
        lw=3
    )
    
    ax[i+1].legend(loc='best', ncol=2)

    ax[i+1].set_ylim(-0.5, 0.5)
    ax[i+1].set_xlim(0,360)
    ax[i+1].axhline(0, ls='--', color='grey')
    ax[i+1].set_ylabel('residuals')
ax[2].set_xlabel('Ecliptic longitude ($^o$)');<jupyter_output>[-0.06743838 14.64085576] [ 0.02284717 22.09911907]
[-0.09013547 -9.715322  ] [0.03908945 5.02009608]<jupyter_text>The residuals between observed and predicted magnitude are most-likely due to the difference of aspect angle: the objets is not a perfect sphere. We are seeing its oblateness here. Position comparison<jupyter_code>fig, ax = plt.subplots(
    figsize=(10, 10), 
    sharex=True,
)

ax.scatter(sso['i:ra'], sso['i:dec'], label='ZTF', alpha=0.2, color=colors[1])
    
ax.plot(eph.RA, eph.Dec, ls='', color='black', marker='x', alpha=0.2, label='Ephemerides')
ax.legend(loc='best')
ax.set_xlabel('RA ($^o$)')
ax.set_ylabel('DEC ($^o$)')

axins = ax.inset_axes([0.2, 0.2, 0.45, 0.45])

deltaRAcosDEC = (sso['i:ra'] - eph.RA)*np.cos(np.radians(sso['i:dec']))*3600
deltaDEC = (sso['i:dec']- eph.Dec)*3600
axins.plot( deltaRAcosDEC, deltaDEC, ls='', color=colors[0], marker='x', alpha=0.5)
axins.errorbar( np.mean(deltaRAcosDEC), np.mean(deltaDEC), xerr=np.std(deltaRAcosDEC), yerr=np.std(deltaDEC) )
axins.axhline(0, ls='--', color='black')
axins.axvline(0, ls='--', color='black')
axins.set_xlabel(r'$\Delta$RA ($^{\prime\prime}$)')
axins.set_ylabel(r'$\Delta$DEC ($^{\prime\prime}$)')<jupyter_output><empty_output><jupyter_text>Phase curve<jupyter_code># SSO Color index
# BCa: I will prepare LSST/ZTF colors to have them at hand
V_minus_g = -0.32
V_minus_r = 0.13
print(V_minus_g, V_minus_r)
fitter = LevMarLSQFitter(calc_uncertainties=True)
fig, ax = plt.subplots(figsize=(12, 5))

kind = 'HG'

for i,f in enumerate(filts):
    cond = info['i:fid']==f
    
    # Color conversion
    if filters[f]=='g':
        color_sso = V_minus_g
    else:
        color_sso = V_minus_r
    
    obs = Obs.from_dict(
        {
            'alpha': info.loc[cond,'Phase'].values * u.deg, 
            'mag': (info.loc[cond,'i:magpsf_red'] + color_sso).values * u.mag,
        }
    )
    
    if kind == 'HG1G2':
        fitfunc = HG1G2
        label = 'H={:.2f}, G1={:.2f}, G2={:.2f}'
    elif kind == 'HG12':
        fitfunc = HG12
        label = 'H={:.2f}, G12={:.2f}'
    elif kind == 'HG':
        fitfunc = HG
        label = 'H={:.2f}, G={:.2f}'
        
    model2 = fitfunc.from_obs(obs, fitter, 'mag', weights=info.loc[cond, 'i:sigmapsf'] * u.mag)
        
    if model2.cov_matrix is not None:
        perr = np.sqrt(np.diag(model2.cov_matrix.cov_matrix))
    else:
        perr = [np.nan, np.nan, np.nan]
    
    ax.errorbar(
        info.loc[cond,'Phase'], 
        info.loc[cond,'i:magpsf_red'] + color_sso, 
        yerr=info.loc[cond,'i:sigmapsf'], 
        label='{:} -> V_obs'.format(filters[f]), 
        linestyle='None', 
        marker='o', 
        color=colors[i], 
        alpha=0.7
    )
    
    ax.plot(
        info.loc[cond,'Phase'], 
        model2.to_mag(info.loc[cond, 'Phase'].values * u.deg, unit=u.deg).to_value(),
        label=label.format(*model2.parameters),
        color=colors[i]
    )
    
ax.legend(loc='best')
ax.set_ylim(ax.get_ylim()[::-1])
ax.set_xlabel('Phase angle ($^o$)')
ax.set_ylabel('V_obs [mag]')
ax.set_title('Per-band fit');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18861047835990888, 'fraction_numerical': 0.029703872437357632, 'mean_word_length': 3.1325301204819276, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 32, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '37687', 'license': 'Apache-2.0', 'n_tokens_mistral': 4628, 'n_tokens_neox': 4374, 'n_words': 802, 'path': 'sso/fink_sso_imcce.ipynb', 'repo_name': 'astrolabsoftware/fink-tutorials', 'size': '10996'}