<jupyter_start><jupyter_text>Case8-challenge00_raeff_pix2wcs Modified version from Case7-0challengef by <NAME>. In this note, we estimate the field parameters and distortion parameters from the observed positions on the focal plane in the overlapped plates. We also use reference stars (Gaia stars) whose sky coordinates are known in a certain accuracy. The SIP-convention distortion is considered in this note.  Preparation First, we load the data from https://github.com/xr0038/jasmine_warpfield/tree/master/challenge/case8.<jupyter_code>import astropy.io.ascii as asc
import astropy.units as u

objects = asc.read('/Users/dkawata/work/obs/projs/JASMINE-Mission/analysis-testing-e2e/jasmine_warpfield/challenge/case8/case8_challenge_00.txt') #consists of x (um), y (um), catalog_id, ra (deg), dec (deg), and field.
pointings = asc.read('/Users/dkawata/work/obs/projs/JASMINE-Mission/analysis-testing-e2e/jasmine_warpfield/challenge/case8/case8_challenge_00_pointing.txt') # consists of field, ra (deg), dec (deg), and pa (deg).<jupyter_output><empty_output><jupyter_text>We can convert the units of x and y from um to pix with assuming the pixel size to be 15 um. However, we will use um unit for the detector coordinates. The input data are created with Sip definition of crpix=[0,0] and origin=1, which map the origin to [0 um, 0 um].<jupyter_code>pix_size = 15.*u.um
# objects['x'] = (objects['x']/pix_size).si
# objects['y'] = (objects['y']/pix_size).si
# objects: x (px), y (px), catalog_id, ra (deg), dec (deg), and field.
# pointings: field, ra (deg), dec (deg), and pa (deg).<jupyter_output><empty_output><jupyter_text>Then, we change the ids for easy handling.<jupyter_code>from astropy.table import unique
import numpy as np

ids = unique(objects, keys='catalog_id')['catalog_id']
objects.add_column(-1, name='id')
for i in range(0, np.size(ids)):
    pos = np.where(objects['catalog_id']==ids[i])
    objects['id'][pos] = i

objects.remove_column('catalog_id')
objects.rename_column('id', 'catalog_id')<jupyter_output><empty_output><jupyter_text>Here, we make some arrays for further analysis. One new array is true_radec which stores true ra/dec values. Duplicated information (rows for the same object) is removed, and the rows are sorted with object ids. Another new array is observed_xy. It contains field ids, observed x/y positions on the focal plane, catalog ids. We rename ra, dec to ra_est, dec_est to store the estimated sky positions.<jupyter_code>true_radec = objects['catalog_id', 'ra', 'dec'].copy()
true_radec.sort('catalog_id')
true_radec = unique(true_radec, keys='catalog_id') # consists of catalog_id, ra (deg), and dec (deg).

observed_xy = objects['field', 'x', 'y', 'catalog_id', 'ra', 'dec'].copy()
# observed_xy.rename_column('ra', 'ra_est')
# observed_xy.rename_column('dec', 'dec_est')
observed_xy.add_column(observed_xy['ra'], name='ra_est')
observed_xy.add_column(observed_xy['dec'],name='dec_est')
# observed_xy will have field, x (px), y (px), catalog_id, and estimated ra/dec (deg).

# initializing ra_est and dec_est
observed_xy['ra_est']  = 0.0
observed_xy['dec_est'] = 0.0<jupyter_output><empty_output><jupyter_text>In addition, we make another array which stores field parameters, ra and dec (deg) of the origin of the pointing and position angle, pa (deg). The plate scale, um pixel scale to deg in the sky, is assumed to be the same value for all plates. At this time, an approximated (initial guess) value is stored in a variable (plate_scale).<jupyter_code>field_params = pointings.copy() # field, ra (deg), dec (deg), and pa (deg).
true_field_params = field_params.copy()
# field_params['pa'] -= 240.0 # offset?
# plate_scale = 8.e-6*u.deg*(pix_size/u.um).si # in deg/pix
plate_scale = 8.e-6*u.deg/u.um
print(plate_scale)<jupyter_output>8e-06 deg / um<jupyter_text>Let's check the object distribution on sky.<jupyter_code>import matplotlib.pylab as plt
import numpy as np

color = ['red', 'blue', 'green', 'orange']

for i in range(0, np.max(field_params['field'])+1):
    pos = np.where(objects['field']==i)
    plt.scatter(objects['ra'][pos], objects['dec'][pos], marker='o', facecolor='None', edgecolor=color[i], s=10*i+10)
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')<jupyter_output><empty_output><jupyter_text>We can see that the data consists of four image plates (different colours indicating the objects observd by the different plantes) and that the overlapped region has a size of about a 1/4 FoV. We select the objects in the overlapped region for further analysis. Here, we select the regions all 4 plates overlaps, but we can use the region overlapped with at least 2 plates.<jupyter_code>true_radec_overlap = true_radec.copy()
observed_xy_overlap = observed_xy.copy()
for cid in true_radec['catalog_id']:
    if np.count_nonzero(observed_xy['catalog_id']==cid)!=4:
    # if np.count_nonzero(observed_xy['catalog_id']==cid)<=1:
        pos = np.where(true_radec_overlap['catalog_id']==cid)[0]
        true_radec_overlap.remove_rows(pos)
        pos = np.where(observed_xy_overlap['catalog_id']==cid)[0]
        observed_xy_overlap.remove_rows(pos)
print(' The number of overlapped unique stars =', len(true_radec_overlap))
print(' The total number of observations of these overlapped stars =', len(observed_xy_overlap))<jupyter_output>The number of overlapped unique stars = 222
 The total number of observations of these overlapped stars = 888<jupyter_text>Let's check the distribution of the selected unique objects.<jupyter_code>plt.scatter(true_radec_overlap['ra'], true_radec_overlap['dec'], marker='o', facecolor='None', edgecolor='orange')
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')
print(len(true_radec_overlap['ra']))<jupyter_output>222<jupyter_text>These objects will be used for the following analysis. We again modify the catalog id for easy handling.<jupyter_code>ids = unique(true_radec_overlap, keys='catalog_id')['catalog_id']
true_radec_overlap.add_column(-1, name='id')
observed_xy_overlap.add_column(-1, name='id')
for i in range(0, np.size(ids)):
    pos = np.where(true_radec_overlap['catalog_id']==ids[i])
    true_radec_overlap['id'][pos] = i
    pos = np.where(observed_xy_overlap['catalog_id']==ids[i])
    observed_xy_overlap['id'][pos] = i

true_radec_overlap.remove_column('catalog_id')
true_radec_overlap.rename_column('id', 'catalog_id')
observed_xy_overlap.remove_column('catalog_id')
observed_xy_overlap.rename_column('id', 'catalog_id')<jupyter_output><empty_output><jupyter_text>First guess of the positions At first, we define a wcs constructor, including SIP polynomial distortion convention, https://irsa.ipac.caltech.edu/data/SPITZER/docs/files/spitzer/shupeADASS.pdf and https://docs.astropy.org/en/stable/api/astropy.wcs.Sip.html.<jupyter_code>from astropy.wcs import WCS
from astropy.wcs import Sip
import astropy.units as u

def wcs(ra_ptg, dec_ptg, pa_ptg, scale, a=None, b=None, ap=None, bp=None):
    w = WCS(naxis=2)
    w.wcs.crpix=[0,0]
    w.wcs.cdelt=np.array([-scale, scale])
    w.wcs.crval=[ra_ptg, dec_ptg]
    w.wcs.ctype=["RA---TAN-SIP", "DEC--TAN-SIP"]
    w.wcs.pc=[[ np.cos(pa_ptg*u.deg), -np.sin(pa_ptg*u.deg)],
              [np.sin(pa_ptg*u.deg), np.cos(pa_ptg*u.deg)]]
#    if a is not None and b is not None:
    w.sip = Sip(a, b, ap, bp, [0, 0])
    
    return w<jupyter_output><empty_output><jupyter_text>Then, we estimate the sky coordinates from the observed focal-plane positions and (approximated) field parameters. Here, we do not add the distorption, but naively convert pixel coordinate (x, y) to sky coordinate, ($\alpha$, $\delta$) (ra_est, dec_est).<jupyter_code>for i in range(0, np.size(field_params)):
    fp = field_params[i]
    w = wcs(fp['ra'], fp['dec'], fp['pa'], plate_scale.value)
    pos = np.where(observed_xy_overlap['field']==fp['field'])
    ret = w.all_pix2world(np.concatenate(([observed_xy_overlap[pos]['x']], [observed_xy_overlap[pos]['y']])).T, 0)
    observed_xy_overlap['ra_est'][pos]  = ret[:, 0]
    observed_xy_overlap['dec_est'][pos] = ret[:, 1]<jupyter_output><empty_output><jupyter_text>Let's check the true positions and estimated positions.<jupyter_code>plt.scatter(true_radec_overlap['ra'], true_radec_overlap['dec'], marker='x', label='True')
plt.scatter(observed_xy_overlap['ra_est'], observed_xy_overlap['dec_est'], marker='+', label='Estimated')
print(' number of stars used =', len(observed_xy_overlap['ra_est']))
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')
plt.legend()<jupyter_output>number of stars used = 888<jupyter_text>Test for distortion with A/B, using the following c and d.<jupyter_code>c = np.zeros(shape=(3, 3))
d = np.zeros(shape=(3, 3))

c[0,2]=-2.34153374723336e-09
c[1,1]=1.5792128155073916e-08
c[1,2]=7.674347291529089e-15
c[2,0]=-4.694743859349522e-09
c[2,1]=5.4256004358596465e-15
c[2,2]=-4.6341769281246224e-21
d[0,2]=-1.913280244657798e-08
d[1,1]=-5.622875292409728e-09
d[1,2]=-1.0128311203344238e-14
d[2,0]=3.1424733259527392e-09
d[2,1]=-9.08024075521211e-15
d[2,2]=-1.4123037013352912e-20<jupyter_output><empty_output><jupyter_text>We check if all_pix2world takes into account SIP parameters of A and B by comparing ($\alpha$, $\delta$) converted from (x, y) pixel coordinate without distortion (above observed_xy_overlap['ra_est'] and observed_xy_overlap['dec_est']) and ($\alpha$, $\delta$) converted from (x, y) pixel coordinate with A and B, ra_dist, dec_dist below.<jupyter_code># print(observed_xy_overlap['ra_est'])
c *= 100.0
ra_dist = np.zeros_like(observed_xy_overlap['ra_est'])
dec_dist = np.zeros_like(observed_xy_overlap['dec_est'])

for i in range(0, np.size(field_params)):
    fp = field_params[i]
    w = wcs(fp['ra'], fp['dec'], fp['pa'], plate_scale.value, a=c, b=d)
    pos = np.where(observed_xy_overlap['field']==fp['field'])
    ret = w.all_pix2world(np.concatenate(([observed_xy_overlap[pos]['x']], [observed_xy_overlap[pos]['y']])).T, 0)
    ra_dist[pos] = ret[:,0]
    dec_dist[pos] = ret[:,1]
    
print(' diff ra=', ra_dist-observed_xy_overlap['ra_est'])
print(' diff dec=', dec_dist-observed_xy_overlap['dec_est'])

plt.scatter(ra_dist, dec_dist, marker='x', label='Distorted')
plt.scatter(observed_xy_overlap['ra_est'], observed_xy_overlap['dec_est'], marker='+', label='No distortion')
print(' number of stars used =', len(observed_xy_overlap['ra_est']))
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')
plt.legend()<jupyter_output>diff ra=          ra_est        
          deg          
-----------------------
 1.2170394995791867e-05
  5.816765053623385e-06
 -9.162415494756715e-06
-1.7497795170129393e-05
-1.4185463442117907e-05
 -1.119997426712871e-05
  6.512692038995738e-06
  5.819727562084154e-07
 1.0868943661535013e-05
  8.796283736955957e-07
                    ...
-2.5157424943245132e-06
 -8.514140176885121e-06
  2.359488917136332e-06
  4.017218486751517e-06
 7.2278966172234504e-06
  3.721560091207721e-05
 1.5106658224794955e-05
 3.1841070267546456e-05
  4.521866071627301e-05
  9.789179591734865e-05
Length = 888 rows
 diff dec=         dec_est        
          deg          
-----------------------
  1.631389871903366e-05
   7.99703507325944e-06
-1.3330118619592213e-05
-2.6079412609902874e-05
 -2.134303151279937e-05
-1.6830686099922332e-05
  7.842507688593514e-06
  8.140157774505496e-07
  1.477228683555154e-05
  1.222328322114663e-06
                    ...
 -3.754488030693892e-06
-1.2820045057537754e-05
 [...]<jupyter_text>Check if these stars cover the large enough detector region by looking at their (x, y) position in the detector coordinate.<jupyter_code>plt.scatter(objects['x'], objects['y'], marker='x', label='All', s=5)
print(' number of all stars=', len(objects['x']))
plt.scatter(observed_xy_overlap['x'], observed_xy_overlap['y'], marker='+', label='Overlap')
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')<jupyter_output>number of all stars= 2563<jupyter_text>Here, there are four estimated (ignoring distortion) positions, (observed_xy_overlap['ra_est'], observed_xy_overlap['dec_est']), in the sky coordinate for each unique object. We take their mean values as the first-guess positions and store them in radec_est array.<jupyter_code>from astropy.table import Table

radec_est = Table(names=['catalog_id', 'ra_est', 'dec_est'], \
            dtype=['int64', 'float64', 'float64'])
#           units=[None, u.deg, u.deg], \
#           dtype=['int64', 'float64', 'float64'])

radec_est['ra_est'].unit = u.deg
radec_est['dec_est'].unit = u.deg

cat_ids = unique(observed_xy_overlap, 'catalog_id')['catalog_id']

for i in cat_ids:
    pos = np.where(observed_xy_overlap['catalog_id'] == i)
    ra = np.mean(observed_xy_overlap[pos]['ra_est'])*u.deg
    dec = np.mean(observed_xy_overlap[pos]['dec_est'])*u.deg
    radec_est.add_row([i, ra, dec])
    
# print('radec_est=', radec_est)<jupyter_output><empty_output><jupyter_text>Let's check the estimated positions.<jupyter_code>plt.scatter(true_radec_overlap['ra'], true_radec_overlap['dec'], marker='x', label='True')
plt.scatter(radec_est['ra_est'], radec_est['dec_est'], marker='+', label='First guess')
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')
plt.legend()<jupyter_output><empty_output><jupyter_text>Parameter adjustment At first, we define a function which calculates x/y positions from the ra/dec values estimated above and the field/catalog ids.<jupyter_code>def xy_calculator(observed_xy, field_params, plate_scale, ap, bp, radec_info):
    # observed_xy:  consists of field, x (px), y (px), catalog_id, ra_est (deg), and dec_est(deg).
    # field_params: consists of field, ra (deg), dec (deg), and pa (deg).
    # radec_info:   consists of catalog_id, ra_est (deg), and dec_est (deg).

    observed_xy_cp = observed_xy.copy()
    observed_xy_cp.rename_column('x', 'x_est')
    observed_xy_cp.rename_column('y', 'y_est')
    observed_xy_cp['x_est']   = None
    observed_xy_cp['y_est']   = None
    observed_xy_cp['ra_est']  = None
    observed_xy_cp['dec_est'] = None

    for i in range(0, np.size(radec_info)):
        pos = np.where(observed_xy_cp['catalog_id']==radec_info[i]['catalog_id'])
        observed_xy_cp['ra_est'][pos]  = radec_info[i]['ra_est']
        observed_xy_cp['dec_est'][pos] = radec_info[i]['dec_est']
    
    for i in range(0, np.size(field_params)):
        fp = field_params[i]
        w = wcs(fp['ra'], fp['dec'], fp['pa'], plate_scale, ap=ap, bp=bp)
        pos = np.where(observed_xy_cp['field']==fp['field'])
        radec0 = np.concatenate(([observed_xy_cp[pos]['ra_est']], [observed_xy_cp[pos]['dec_est']])).T
        ret = w.sip_foc2pix(w.wcs_world2pix(radec0, 1)-w.wcs.crpix, 1)
        observed_xy_cp['x_est'][pos] = ret[:, 0]
        observed_xy_cp['y_est'][pos] = ret[:, 1]

    return observed_xy_cp['x_est', 'y_est']<jupyter_output><empty_output><jupyter_text>Next, we define a function to map from (x, y) pixel coordinate to ($\alpha$, $\delta$), using A/B Sip distortion parameters using wcs.all_pix2world, https://docs.astropy.org/en/stable/api/astropy.wcs.WCS.htmlastropy.wcs.WCS.all_pix2world with input field parameters of $\alpha_{\rm ptgs}$ (deg), $\delta_{\rm ptgs}$ (deg) and pa$_{\rm ptgs}$ (deg) of each field (plate) pointing. This conversion is described as follows. Here, we follow the description at https://www.stsci.edu/itt/review/DrizzlePac/HTML/ch33.html  DefinitionCRVAL1: $\alpha_{\rm ptgs}$ right assension at the pointing centre.CRVAL2: $\delta_{\rm ptgs}$ declination at the pointing centre.CRPIX1: the x reference location of the image plate, corresponding to the pointing centre. We set CRPIX1=0.CRPIX2: the y reference location of the image plate, corresponding to the pointing centre. We set CRPIX2=0.CDELT1: x-scale. We set -scale.CDELT2: y-scale. We set scale. wcs compute the sky coordidate, ($\alpha$, $\delta$) of star at (x, y) on the detector as follows (This needs to be confirmed...).According to https://irsa.ipac.caltech.edu/data/SPITZER/docs/files/spitzer/shupeADASS.pdf , intermediate world coordinate, $(\alpha', \delta')$ are computed with$ \begin{pmatrix}\alpha' \\\delta' \\\end{pmatrix}=\begin{pmatrix}CD_{11} & CD_{12} \\CD_{21} & CD_{22} \\\end{pmatrix}\begin{pmatrix}x+f(x,y) \\y+g(x,y) \\\end{pmatrix},$where $$f(x,y) = \sum_{p,q} A_{pq} x^p, y^q \ \ \ \  p+q\leq {\rm A_{ORDER}}, \\g(x,y) = \sum_{p,q} B_{pq} x^p, y^q \ \ \ \  p+q\leq {\rm B_{ORDER}}. \\$$${A_{ORDER}}$ and ${B_{ORDER}}$ is the order of polynomial. For example, with ${A_{ORDER}}=3$, $f(x,y)$ can be described as,$f(x,y) = A_{20} x^2 + A_{02} y^2 + A_{11} x y + A_{21} x^2 y + A_{12} x y^2 + A_{30} x^3 + A_{03} y^3.$Note that $p,q\leq1$ components, i.e. $A_{00}=A_{01}=A_{10}$, are ignored here. **Therefore, we set $A_{00}=A_{01}=A_{10}$.** From the intermediate world coordinate $(\alpha', \delta')$ to the world coordinate $(\alpha, \delta)$ is taken care by WCS ([Greisen & Calabretta 2002, A&A, 395, 1061](https://ui.adsabs.harvard.edu/abs/2002A%26A...395.1061G/abstract); [Clabretta & Greisen 2002, A&A, 395, 1077](https://ui.adsabs.harvard.edu/abs/2002A%26A...395.1077C/abstract)). CD matrix above is considered to be described (**need to check** ).$\begin{pmatrix}CD_{11} & CD_{12} \\CD_{21} & CD_{22} \\\end{pmatrix}=\begin{pmatrix}{\rm CDELT1} \cos({\rm pa_{ptgs}}) & -{\rm CDELT2} \sin({\rm pa_{ptgs}}) \\{\rm CDELT1} \sin({\rm pa_{ptgs}}) & {\rm CDELT2} \cos({\rm pa_{ptgs}}),  \\\end{pmatrix}$according to equation (189) in [Clabretta & Greisen 2002, A&A, 395, 1077](https://ui.adsabs.harvard.edu/abs/2002A%26A...395.1077C/abstract). Note that here, we do not have a term to describe a linear skewness.<jupyter_code>def radec_calculator_ab(observed_xy, field_params, plate_scale, a, b):
    # observed_xy:  consists of field, x (px), y (px), catalog_id, ra_est (deg), and dec_est(deg).
    # field_params: consists of field, ra (deg), dec (deg), and pa (deg).

    
    observed_xy_cp = observed_xy.copy()
    # observed_xy_cp.rename_column('x', 'x_est')
    # observed_xy_cp.rename_column('y', 'y_est')
    # observed_xy_cp['x_est']   = None
    # observed_xy_cp['y_est']   = None
    observed_xy_cp['ra_est']  = None
    observed_xy_cp['dec_est'] = None
    
    for i in range(0, np.size(field_params)):
        fp = field_params[i]
        w = wcs(fp['ra'], fp['dec'], fp['pa'], plate_scale, a=a, b=b)
        pos = np.where(observed_xy_cp['field']==fp['field'])
        pix0 = np.concatenate(([observed_xy_cp[pos]['x']], [observed_xy_cp[pos]['y']])).T
        ret = w.all_pix2world(pix0, 1)
        # ret = w.sip_pix2foc(w.wcs_pix2world(pix0, 1)-w.wcs.crval, 1)
        observed_xy_cp['ra_est'][pos]  = ret[:, 0]
        observed_xy_cp['dec_est'][pos] = ret[:, 1]
    
    return observed_xy_cp['ra_est', 'dec_est']<jupyter_output><empty_output><jupyter_text>Using scipy.optimize least_squares, assuming the pointing sky coordinate, RA, DEC are accurately known. Define model function to solve with Least Squares. Here, we consider that the position of the pointings, $(\alpha_{\rm ptgs}, \delta_{\rm ptgs})$ are accurately known, and they are fixed. Following equation (1) of [Bernstein et al. (2017)](https://ui.adsabs.harvard.edu/abs/2017PASP..129g4503B/abstract), we minimise the following $\chi^2$ to fit the parameters of position angle of ${\rm pa_{ptgs,p}}$ for each image plate, p, and the common "scale" parameter to convert pixel unit of um to world coordinate unit of deg, in addition to the components of ${\rm A_p}$, ${\rm B_p}$ distortion matrices. $\chi^2 = \sum_i (w_i | \alpha_{i,s} - \bar{\alpha_s}|)^2,$where $\alpha_{i,s}({\rm pa_{ptgs,p}, {\rm scale}, A_p, B_p})$ is the world coordinate position estimated from $(x_i, y_i)$ for star s on a plate p, using the parameters for plate p, ${\rm pa_{ptgs,p}, A_p, B_p}$. $\bar{\alpha_s}$ is the mean sky position of star s, and described as$$ \bar{\alpha_s} = \frac{\sum_{i\in s} w_i \alpha_{i,s}}{\sum_{i\in s} w_i}.$$Here, we do not take into account any error, and set $w_i = 1$.<jupyter_code>import copy

# def model_func(params, n_fields, dim_sip, observed_xy):
def model_func(params, ra_ptgs, dec_ptgs, n_fields, dim_sip, observed_xy):
    # params = (ra_ptgs, dec_ptgs, pa_ptg..., scale, a..., b...)
    pa_ptgs, scale, a, b =\
          np.split(params, [n_fields, n_fields+1,\
                          n_fields+1+(dim_sip+1)**2])
    # ra_ptgs, dec_ptgs, pa_ptgs, scale, a, b =\
    #    np.split(params, [n_fields, 2*n_fields, 3*n_fields, 3*n_fields+1,\
    #                      3*n_fields+1+(dim_sip+1)**2])

    field_params = Table(data=[ra_ptgs, dec_ptgs, pa_ptgs, -np.ones(shape=(np.size(ra_ptgs)))],\
                         names=['ra', 'dec', 'pa', 'field'],\
                         dtype=['float64', 'float64', 'float64', 'int64'])
#                         names=['ra', 'dec', 'pa', 'field'],\
#                         units=[u.deg, u.deg, u.deg, None],\
#                         dtype=['float64', 'float64', 'float64', 'int64'])
    field_params['ra'].unit = u.deg
    field_params['dec'].unit = u.deg
    field_params['pa'].unit = u.deg
    field_params['field'] = np.arange(0, np.size(field_params))
    
    # use copy of observed_xy
    observed_xy_cp = observed_xy.copy()
    
    a_matrix = np.reshape(a, (dim_sip+1, dim_sip+1))
    b_matrix = np.reshape(b, (dim_sip+1, dim_sip+1))
    
    # mns = np.concatenate(((0, 1), np.arange(dim_sip+1, 2*(dim_sip)+1))) 
    # for mn in mns:
    #    for m in range(np.max([0, mn-dim_sip]), np.min([mn+1, dim_sip+1])):
    #        n = mn - m
    #        ap_matrix[m, n] = 0
    #        bp_matrix[m, n] = 0
    a_matrix[0, 0] = 0.0
    a_matrix[0, 1] = 0.0
    a_matrix[1, 0] = 0.0
    b_matrix[0, 0] = 0.0
    b_matrix[0, 1] = 0.0
    b_matrix[1, 0] = 0.0

    m, n = np.indices((dim_sip+1, dim_sip+1))
    mn = m + n
    a_matrix = a_matrix * (1.e-3**mn)
    b_matrix = b_matrix * (1.e-3**mn)
    
    # compute ra/dec from x/y with the parameters.
    ret = radec_calculator_ab(observed_xy_cp, field_params, scale[0], \
                                                              a_matrix, b_matrix)
    
    observed_xy_cp['ra_est'] = ret['ra_est']
    observed_xy_cp['dec_est'] = ret['dec_est']
    
    # compute the mean ra/dec for unique stars
    cat_ids = unique(observed_xy_cp, 'catalog_id')['catalog_id']
    
    ra_mean = np.zeros(len(observed_xy['ra_est']))
    dec_mean = np.zeros(len(observed_xy['dec_est']))
    
    # effective RA
    observed_xy_cp['ra_est'] = observed_xy_cp['ra_est']*np.cos(observed_xy_cp['dec_est'].data*u.deg)

    for i in cat_ids:
        pos = np.where(observed_xy_cp['catalog_id'] == i)
        ra_mean[pos] = np.mean(observed_xy_cp[pos]['ra_est'])*u.deg
        dec_mean[pos] = np.mean(observed_xy_cp[pos]['dec_est'])*u.deg
    

    radec_est = np.concatenate((observed_xy_cp['ra_est'], observed_xy_cp['dec_est']))
    radec_est_mean = np.concatenate((ra_mean, dec_mean))
    

    residuals = radec_est - radec_est_mean

    return residuals<jupyter_output><empty_output><jupyter_text>Next, we execute the least-square calculation to derive the field parameters and sky positions of the objects in the overlapped region.<jupyter_code>from scipy.optimize import least_squares
import time

dim_sip = 3
a = np.zeros(shape=(dim_sip+1, dim_sip+1))
b = np.zeros(shape=(dim_sip+1, dim_sip+1))

# constructing a_init (initial parameter set).
# a_init = np.array(np.concatenate((field_params['ra'], field_params['dec'], field_params['pa'], \
#                                   [plate_scale.value], a.flatten(), b.flatten()))) # This must be an ndarray.
a_init = np.array(np.concatenate((field_params['pa'], \
                                  [plate_scale.value], a.flatten(), b.flatten()))) # This must be an ndarray.
print(' # of fitting parameters =', len(a_init))
# constraining ra/dec values in 'observed' between -180 and 180 deg.
# measured = np.concatenate((observed_xy_overlap['x'], observed_xy_overlap['y']))
# print(' # of data points =', len(measured))
#pos = np.where(measured>180.)
#measured[pos] -= 360.
#pos = np.where(measured<-180.)
#measured[pos] += 360.

start = time.time()
# result = least_squares(model_func, a_init, loss='linear', args=(np.size(field_params), \
#                         dim_sip, observed_xy_overlap), \
#                       verbose=2)
result = least_squares(model_func, a_init, loss='linear', args=(field_params['ra'], \
                      field_params['dec'], np.size(field_params), dim_sip, observed_xy_overlap), \
                      verbose=2)
print(' time=',time.time()-start)

## pa should be a positive value between 0 and 360.
#if result[3] < 0:
#    result[3] = -result[3]
#    result[2] = result[2] + 180.0
#
#if result[2] > 360.0 or result[2] < 0.0:
#    result[2] = result[2]%360.0<jupyter_output># of fitting parameters = 37
   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality   
       0              1         4.8745e-03                                    6.54e+04    
       1              2         2.9443e-09      4.87e-03       6.03e+02       3.18e+00    
       2              3         2.7798e-09      1.65e-10       1.21e+03       1.35e-03    
       3              4         2.7798e-09      1.37e-16       2.41e+03       1.33e-06    
       4             19         2.7798e-09      0.00e+00       0.00e+00       1.33e-06    
`xtol` termination condition is satisfied.
Function evaluations 19, initial cost 4.8745e-03, final cost 2.7798e-09, first-order optimality 1.33e-06.
 time= 11.082858800888062<jupyter_text>Checking results  PreparationExtracting the resaltnat fitting parameter values.<jupyter_code>n_fields = np.size(field_params)
n_objects = np.size(radec_est)
true_ra_ptgs = true_field_params['ra'].data
true_dec_ptgs = true_field_params['dec'].data
true_pa_ptgs = true_field_params['pa'].data

# ra_ptgs, dec_ptgs, pa_ptgs, scale, a, b =\
#       np.split(result.x, [n_fields, 2*n_fields, 3*n_fields, 3*n_fields+1,\
#                           3*n_fields+1+(dim_sip+1)**2])
pa_ptgs, scale, a, b =\
       np.split(result.x, [n_fields, n_fields+1,\
                           n_fields+1+(dim_sip+1)**2])
ra_ptgs = field_params['ra'].data
dec_ptgs = field_params['dec'].data

a_matrix = np.reshape(a, (dim_sip+1, dim_sip+1))
b_matrix = np.reshape(b, (dim_sip+1, dim_sip+1))
# A/B scaling
m, n = np.indices((dim_sip+1, dim_sip+1))
mn = m + n
a_matrix = a_matrix * (1.e-3**mn)
b_matrix = b_matrix * (1.e-3**mn)

fit_field_params = Table(data=[ra_ptgs, dec_ptgs, pa_ptgs, -np.ones(shape=(np.size(ra_ptgs)))],\
                         names=['ra', 'dec', 'pa', 'field'],\
                         dtype=['float64', 'float64', 'float64', 'int64'])
fit_field_params['ra'].unit = u.deg
fit_field_params['dec'].unit = u.deg
fit_field_params['pa'].unit = u.deg
fit_field_params['field'] = np.arange(0, np.size(field_params))<jupyter_output><empty_output><jupyter_text>PointingsThe central position of the pointings are not fitted, but fixed with the true values. Hence, the differences should be zero.<jupyter_code>print(' pointing centre (fit) ra, dec (deg) =', ra_ptgs, dec_ptgs)
print(' pointing centre (true) ra, dec (deg) =', true_ra_ptgs, true_dec_ptgs)
print(' difference ra, dec (deg) =', ra_ptgs-true_ra_ptgs, dec_ptgs-true_dec_ptgs)<jupyter_output>pointing centre (fit) ra, dec (deg) = [265.6202439  265.70081783 265.4894155  265.56770499] [-28.85396419 -28.74323799 -28.78375368 -28.67010405]
 pointing centre (true) ra, dec (deg) = [265.6202439  265.70081783 265.4894155  265.56770499] [-28.85396419 -28.74323799 -28.78375368 -28.67010405]
 difference ra, dec (deg) = [0. 0. 0. 0.] [0. 0. 0. 0.]<jupyter_text>Pointings position angles without  A00=A01=A10=0<jupyter_code>print(' position angle (fit) (deg) =', pa_ptgs)
print(' position angle (true) (deg) =', true_pa_ptgs)
print(' difference =', pa_ptgs-true_pa_ptgs)<jupyter_output>position angle (fit) (deg) = [302.01808056 301.17384262 300.93373878 301.57720235]
 position angle (true) (deg) = [302.02408829 301.17958541 300.93717604 301.58002573]
 difference = [-0.00600773 -0.00574279 -0.00343726 -0.00282338]<jupyter_text>Pointing position angle with A00=A01=A10=0<jupyter_code>print(' position angle (fit) (deg) =', pa_ptgs)
print(' position angle (true) (deg) =', true_pa_ptgs)
print(' difference =', pa_ptgs-true_pa_ptgs)<jupyter_output>position angle (fit) (deg) = [302.01808056 301.17384262 300.93373878 301.57720235]
 position angle (true) (deg) = [302.02408829 301.17958541 300.93717604 301.58002573]
 difference = [-0.00600773 -0.00574279 -0.00343726 -0.00282338]<jupyter_text>Scale (deg/pix) with  A00=A01=A10=0<jupyter_code>print(' scale (fit) deg/um=', scale)
print(' true scale =',(1e-6/7.3/np.pi*180.0)*u.deg/u.um)
# print(' true scale =',(1e-6/7.3/np.pi*180.0)*u.deg*(pix_size/u.um).si)
print(' difference =', scale-(1e-6/7.3/np.pi*180.0))<jupyter_output>scale (fit) deg/um= [7.85066081e-06]
 true scale = 7.848736919600318e-06 deg / um
 difference = [1.92389171e-09]<jupyter_text>A/B<jupyter_code>print(' derived A/B matrices = ', a_matrix, b_matrix)<jupyter_output>derived A/B matrices =  [[ 0.00000000e+00  0.00000000e+00 -4.71606279e-09 -6.70732967e-14]
 [ 0.00000000e+00 -1.44099032e-09  1.56894290e-13  1.62253635e-10]
 [-9.93306006e-09 -6.47094309e-14 -2.16749008e-10  2.59837384e-13]
 [ 2.00058595e-13 -1.04182162e-09 -5.49703540e-13 -6.29226746e-16]] [[ 0.00000000e+00  0.00000000e+00 -4.26325960e-09 -1.94242607e-13]
 [ 0.00000000e+00 -5.27060768e-09  1.49196168e-13  6.31955220e-11]
 [-7.06469819e-10  4.81358364e-14  4.16572121e-10  1.25413380e-12]
 [ 4.50680831e-13  1.45572251e-09  1.99277731e-12  2.23491378e-15]]<jupyter_text>Object positions<jupyter_code>print(' field params=', fit_field_params)
radec_objs = radec_calculator_ab(observed_xy_overlap, fit_field_params, scale[0], a_matrix, b_matrix)

plt.scatter(true_radec_overlap['ra'], true_radec_overlap['dec'], marker='x', label='True')
plt.scatter(radec_est['ra_est'], radec_est['dec_est'], marker='+', label='Initial guess')
plt.scatter(radec_objs['ra_est'], radec_objs['dec_est'], marker='.', label='Final estimation')
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')
plt.title('Object positions')
plt.legend()<jupyter_output>field params=         ra                dec                 pa         field
       deg                deg                deg              
----------------- ------------------- ------------------ -----
265.6202439021891 -28.853964194125034  302.0180805640577     0
265.7008178261919   -28.7432379906527 301.17384261837486     1
265.4894154993913  -28.78375368278103 300.93373877569456     2
265.5677049936395 -28.670104050957786  301.5772023481107     3<jupyter_text>Position difference<jupyter_code>from astropy.coordinates import SkyCoord

distlist = []
print(np.shape(radec_objs))
for i in range(0, np.size(radec_objs)):
    c1 = SkyCoord(radec_objs['ra_est'][i]*u.deg, radec_objs['dec_est'][i]*u.deg)
    c2 = SkyCoord(observed_xy_overlap['ra'][i]*u.deg, observed_xy_overlap['dec'][i]*u.deg)
    distlist.append(c1.separation(c2).arcsec)
distlist = np.array(distlist)
#plt.hist(np.log10(distlist))
plt.hist(distlist)
plt.xlabel("Residual (arcsec)")
plt.ylabel("Number")
# use effective RA
dra = radec_objs['ra_est']*np.cos(radec_objs['dec_est'].data*u.deg) \
        -observed_xy_overlap['ra']*np.cos(observed_xy_overlap['dec'].data*u.deg)
ddec = ((radec_objs['dec_est']-observed_xy_overlap['dec']))
# print(' cos test=', np.cos(radec_objs['dec_est'].data*u.deg)-np.cos((radec_objs['dec_est'].data*u.deg).to_value(u.radian)))
dra_arcsec = (dra.data*u.deg).to_value(u.arcsec)
ddec_arcsec = (ddec.data*u.deg).to_value(u.arcsec)
plt.scatter(dra_arcsec, ddec_arcsec, marker='x')
plt.xlabel('dRA (arcsec)')
plt.ylabel('dDEC (arcsec)')
#plt.xlim([-0,8, 0.0])
#plt.ylim([-0.8, 0.0])<jupyter_output><empty_output><jupyter_text>Not fixing pointing RA, DEC, but use the reference stars.  First, we define the model function to evaluate the difference in the sky coordinate of i-th stars, (ra, dec)i, from the individual plate, j-th plate, coordinates, (x, y)ij and the residual between (ra, dec)i and (ra, dec)k, for k-th reference stars whose (ra, dec)k is known from the other observation, e.g. Gaia. We use effecitve RA $\alpha_* = \alpha \cos(\delta)$.<jupyter_code>def model_wrefs_func(params, n_fields, dim_sip, observed_xy, radec_refstars):
    # params = (ra_ptgs, dec_ptgs, pa_ptg..., scale, a..., b...)
    ra_ptgs, dec_ptgs, pa_ptgs, scale, a, b =\
        np.split(params, [n_fields, 2*n_fields, 3*n_fields, 3*n_fields+1,\
                          3*n_fields+1+(dim_sip+1)**2])

    field_params = Table(data=[ra_ptgs, dec_ptgs, pa_ptgs, -np.ones(shape=(np.size(ra_ptgs)))],\
                         names=['ra', 'dec', 'pa', 'field'],\
                         dtype=['float64', 'float64', 'float64', 'int64'])
#                         names=['ra', 'dec', 'pa', 'field'],\
#                         units=[u.deg, u.deg, u.deg, None],\
#                         dtype=['float64', 'float64', 'float64', 'int64'])
    field_params['ra'].unit = u.deg
    field_params['dec'].unit = u.deg
    field_params['pa'].unit = u.deg
    field_params['field'] = np.arange(0, np.size(field_params))
    
    # use copy of observed_xy
    observed_xy_cp = observed_xy.copy()
    
    a_matrix = np.reshape(a, (dim_sip+1, dim_sip+1))
    b_matrix = np.reshape(b, (dim_sip+1, dim_sip+1))
    
    # mns = np.concatenate(((0, 1), np.arange(dim_sip+1, 2*(dim_sip)+1))) 
    # for mn in mns:
    #    for m in range(np.max([0, mn-dim_sip]), np.min([mn+1, dim_sip+1])):
    #        n = mn - m
    #        ap_matrix[m, n] = 0
    #        bp_matrix[m, n] = 0
    a_matrix[0, 0] = 0.0
    a_matrix[0, 1] = 0.0
    a_matrix[1, 0] = 0.0
    b_matrix[0, 0] = 0.0
    b_matrix[0, 1] = 0.0
    b_matrix[1, 0] = 0.0

    # normalisation. 
    m, n = np.indices((dim_sip+1, dim_sip+1))
    mn = m + n
    a_matrix = a_matrix * (1.e-3**mn)
    b_matrix = b_matrix * (1.e-3**mn)
    
    # compute ra/dec from x/y with the parameters.
    ret = radec_calculator_ab(observed_xy_cp, field_params, scale[0], \
                                                              a_matrix, b_matrix)
    
    observed_xy_cp['ra_est'] = ret['ra_est']
    observed_xy_cp['dec_est'] = ret['dec_est']
    
    # compute the mean ra/dec for unique stars
    cat_ids = unique(observed_xy_cp, 'catalog_id')['catalog_id']
    
    ra_mean = np.zeros(len(observed_xy['ra_est']))
    dec_mean = np.zeros(len(observed_xy['dec_est']))

    # effective RA
    observed_xy_cp['ra_est'] = observed_xy_cp['ra_est']*np.cos(observed_xy_cp['dec_est'].data*u.deg)

    for i in cat_ids:
        pos = np.where(observed_xy_cp['catalog_id'] == i)
        ra_mean[pos] = np.mean(observed_xy_cp[pos]['ra_est'])*u.deg
        dec_mean[pos] = np.mean(observed_xy_cp[pos]['dec_est'])*u.deg

    # effective RA for reference stars
    radec_refstars_cp = radec_refstars.copy()
    radec_refstars_cp.rename_column('ra', 'ra*')
    radec_refstars_cp['ra*'] = radec_refstars_cp['ra*']*np.cos(radec_refstars_cp['dec'].data*u.deg)
    # print(radec_refstars['ra'], radec_refstars_cp['ra*'])
        
    # reference stars' measured mean ra, dec to be compared 
    # with the ra, dec of reference stars. 
    radec_est_refstars = radec_refstars.copy()
    radec_est_refstars.rename_column('ra', 'ra_est')
    radec_est_refstars.rename_column('dec', 'dec_est')
    
    for i,id in enumerate(radec_refstars['catalog_id']):
      # print('i, id=', i, id)
      # print(ra_mean[observed_xy_cp['catalog_id'] == id][0])
      radec_est_refstars[i]['ra_est'] = ra_mean[observed_xy_cp['catalog_id'] == id][0]
      radec_est_refstars[i]['dec_est'] = dec_mean[observed_xy_cp['catalog_id'] == id][0]

    
    radec_est = np.concatenate((observed_xy_cp['ra_est'], observed_xy_cp['dec_est'], \
                              radec_refstars_cp['ra*'], radec_refstars_cp['dec']))
    radec_est_mean = np.concatenate((ra_mean, dec_mean, radec_est_refstars['ra_est'], \
                                    radec_est_refstars['dec_est']))
    
    # print(' ref stars al diff=', np.sum(radec_refstars_cp['ra*']-radec_est_refstars['ra_est']))
    
    residuals = radec_est - radec_est_mean

    return residuals<jupyter_output><empty_output><jupyter_text>Pick the reference stars from true_radec_overlap of overlap stars.<jupyter_code># print(' true_radec_overlap =', true_radec_overlap)
print(' len =', len(true_radec_overlap))

# number of reference stars
n_refstars = 10
pos = np.random.choice(len(true_radec_overlap), size=n_refstars, replace=False)
radec_refstars = true_radec_overlap[pos]
print(radec_refstars)<jupyter_output>len = 222
        ra                 dec         catalog_id
       deg                 deg                   
------------------ ------------------- ----------
265.63919602229487 -28.700155552688123        168
 265.5751287163082 -28.748045938319695        152
265.58403117517577  -28.76026745169417        148
265.50674203925803 -28.813848543503703         79
265.57752654272946 -28.653709577335015        218
265.57217871776766 -28.648209008247456        217
 265.6104588526529 -28.720814374071875        163
265.65743669231915 -28.822639963970435         31
265.64537142662545 -28.784050026954677         49
 265.5976812951617 -28.829715640389534         23<jupyter_text>Now, let's run least_squares and get the distortion parameters with the reference stars' constraints.<jupyter_code>from scipy.optimize import least_squares
import time

dim_sip = 4
a = np.zeros(shape=(dim_sip+1, dim_sip+1))
b = np.zeros(shape=(dim_sip+1, dim_sip+1))

# constructing a_init (initial parameter set).
a_init = np.array(np.concatenate((field_params['ra'], field_params['dec'], \
                                  field_params['pa'], \
                                 [plate_scale.value], a.flatten(), b.flatten()))) # This must be an ndarray.
# a_init = np.array(np.concatenate((field_params['pa'], \
#                                  [plate_scale.value], a.flatten(), b.flatten()))) # This must be an ndarray.
print(' # of fitting parameters =', len(a_init))

print(' size of reference stars =', np.size(radec_refstars['catalog_id']))

start = time.time()
result = least_squares(model_wrefs_func, a_init, loss='linear', args= \
                       (np.size(field_params), dim_sip, observed_xy_overlap, \
                        radec_refstars), verbose=2)
print(' time=',time.time()-start)<jupyter_output># of fitting parameters = 63
 size of reference stars = 10
   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality   
       0              1         4.9202e-03                                    6.61e+04    
       1              2         3.3253e-10      4.92e-03       8.06e+02       3.25e+00    
       2              3         1.6627e-10      1.66e-10       1.61e+03       4.17e-04    
       3             18         1.6627e-10      0.00e+00       0.00e+00       4.17e-04    
`xtol` termination condition is satisfied.
Function evaluations 18, initial cost 4.9202e-03, final cost 1.6627e-10, first-order optimality 4.17e-04.
 time= 13.792067289352417<jupyter_text>Checking restuls Preparation<jupyter_code>n_fields = np.size(field_params)
n_objects = np.size(radec_est)
true_ra_ptgs = true_field_params['ra'].data
true_dec_ptgs = true_field_params['dec'].data
true_pa_ptgs = true_field_params['pa'].data

ra_ptgs, dec_ptgs, pa_ptgs, scale, a, b =\
       np.split(result.x, [n_fields, 2*n_fields, 3*n_fields, 3*n_fields+1,\
                           3*n_fields+1+(dim_sip+1)**2])
# pa_ptgs, scale, a, b =\
#       np.split(result.x, [n_fields, n_fields+1,\
#                           n_fields+1+(dim_sip+1)**2])
#ra_ptgs = field_params['ra'].data
# dec_ptgs = field_params['dec'].data

a_matrix = np.reshape(a, (dim_sip+1, dim_sip+1))
b_matrix = np.reshape(b, (dim_sip+1, dim_sip+1))
# A/B scaling
m, n = np.indices((dim_sip+1, dim_sip+1))
mn = m + n
a_matrix = a_matrix * (1.e-3**mn)
b_matrix = b_matrix * (1.e-3**mn)

fit_field_params = Table(data=[ra_ptgs, dec_ptgs, pa_ptgs, -np.ones(shape=(np.size(ra_ptgs)))],\
                         names=['ra', 'dec', 'pa', 'field'],\
                         dtype=['float64', 'float64', 'float64', 'int64'])
fit_field_params['ra'].unit = u.deg
fit_field_params['dec'].unit = u.deg
fit_field_params['pa'].unit = u.deg
fit_field_params['field'] = np.arange(0, np.size(field_params))<jupyter_output><empty_output><jupyter_text>Pointings RA, DEC, position angle and scale<jupyter_code>print(' pointing centre (fit) ra, dec (deg) =', ra_ptgs, dec_ptgs)
print(' pointing centre (true) ra, dec (deg) =', true_ra_ptgs, true_dec_ptgs)
print(' difference ra, dec (deg) =', ra_ptgs-true_ra_ptgs, dec_ptgs-true_dec_ptgs)
print(' position angle (fit) (deg) =', pa_ptgs)
print(' position angle (true) (deg) =', true_pa_ptgs)
print(' difference =', pa_ptgs-true_pa_ptgs)
print(' scale (fit, true) =', scale, (1e-6/7.3/np.pi*180.0)*u.deg/u.um)
print(' difference =', scale-(1e-6/7.3/np.pi*180.0))<jupyter_output>pointing centre (fit) ra, dec (deg) = [265.62028571 265.7008365  265.48944065 265.56770947] [-28.85395845 -28.74322645 -28.78374376 -28.67008675]
 pointing centre (true) ra, dec (deg) = [265.6202439  265.70081783 265.4894155  265.56770499] [-28.85396419 -28.74323799 -28.78375368 -28.67010405]
 difference ra, dec (deg) = [4.18069611e-05 1.86722713e-05 2.51542362e-05 4.47158408e-06] [5.74486724e-06 1.15429180e-05 9.92206167e-06 1.72978067e-05]
 position angle (fit) (deg) = [302.01858626 301.17544275 300.93198899 301.57583732]
 position angle (true) (deg) = [302.02408829 301.17958541 300.93717604 301.58002573]
 difference = [-0.00550203 -0.00414265 -0.00518704 -0.0041884 ]
 scale (fit, true) = [7.85063513e-06] 7.848736919600318e-06 deg / um
 difference = [1.89820786e-09]<jupyter_text>Object positions<jupyter_code>radec_objs = radec_calculator_ab(observed_xy_overlap, fit_field_params, scale[0], a_matrix, b_matrix)

plt.scatter(true_radec_overlap['ra'], true_radec_overlap['dec'], marker='x', label='True')
plt.scatter(radec_est['ra_est'], radec_est['dec_est'], marker='+', label='Initial guess')
plt.scatter(radec_objs['ra_est'], radec_objs['dec_est'], marker='.', label='Final estimation')
plt.scatter(radec_refstars['ra'], radec_refstars['dec'], marker='o', \
            label='Reference stars')
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')
plt.title('Object positions')
plt.legend()<jupyter_output><empty_output><jupyter_text>Position differences<jupyter_code>from astropy.coordinates import SkyCoord

distlist = []
print(np.shape(radec_objs))
for i in range(0, np.size(radec_objs)):
    c1 = SkyCoord(radec_objs['ra_est'][i]*u.deg, radec_objs['dec_est'][i]*u.deg)
    c2 = SkyCoord(observed_xy_overlap['ra'][i]*u.deg, observed_xy_overlap['dec'][i]*u.deg)
    distlist.append(c1.separation(c2).arcsec)
distlist = np.array(distlist)

#plt.hist(np.log10(distlist))
plt.hist(distlist)
plt.xlabel("Residual (arcsec)")
plt.ylabel("Number")
# dra = ((radec_objs['ra_est']-observed_xy_overlap['ra']).data)*u.deg
# effective RA
dra = radec_objs['ra_est']*np.cos(radec_objs['dec_est'].data*u.deg) \
        -observed_xy_overlap['ra']*np.cos(observed_xy_overlap['dec'].data*u.deg)
ddec = ((radec_objs['dec_est']-observed_xy_overlap['dec']).data)*u.deg
dra_arcsec = dra.to_value(u.arcsec)
ddec_arcsec = ddec.to_value(u.arcsec)
plt.scatter(dra_arcsec, ddec_arcsec, marker='x')
plt.xlabel('dRA (arcsec)')
plt.ylabel('dDEC (arcsec)')
#plt.xlim([-0,8, 0.0])
#plt.ylim([-0.8, 0.0])<jupyter_output><empty_output><jupyter_text>Apply the field parameters to all the objects.<jupyter_code>print(' total # of stars =', len(observed_xy))
radec_allobjs = radec_calculator_ab(observed_xy, fit_field_params, \
                                    scale[0], a_matrix, b_matrix)

plt.scatter(observed_xy['ra'], observed_xy['dec'], marker='x', label='True')
plt.scatter(radec_allobjs['ra_est'], radec_allobjs['dec_est'], marker='.', label='Final estimation')
plt.scatter(radec_refstars['ra'], radec_refstars['dec'], marker='o', \
            label='Reference stars')
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')
plt.title('Object positions')
plt.legend()
distlist = []
print(np.shape(radec_allobjs))
for i in range(0, np.size(radec_allobjs)):
    c1 = SkyCoord(radec_allobjs['ra_est'][i]*u.deg, radec_allobjs['dec_est'][i]*u.deg)
    c2 = SkyCoord(observed_xy['ra'][i]*u.deg, observed_xy['dec'][i]*u.deg)
    distlist.append(c1.separation(c2).arcsec)
distlist = np.array(distlist)

#plt.hist(np.log10(distlist))
plt.hist(distlist)
plt.xlabel("Residual (arcsec)")
plt.ylabel("Number")
# dra = ((radec_allobjs['ra_est']-observed_xy['ra']).data)*u.deg
# effective RA
dra = radec_allobjs['ra_est']*np.cos(radec_allobjs['dec_est'].data*u.deg) \
        -observed_xy['ra']*np.cos(observed_xy['dec'].data*u.deg)
ddec = ((radec_allobjs['dec_est']-observed_xy['dec']).data)*u.deg
dra_arcsec = dra.to_value(u.arcsec)
ddec_arcsec = ddec.to_value(u.arcsec)
plt.scatter(dra_arcsec, ddec_arcsec, marker='x')
plt.xlabel('dRA (arcsec)')
plt.ylabel('dDEC (arcsec)')
#plt.xlim([-0,8, 0.0])
#plt.ylim([-0.8, 0.0])<jupyter_output><empty_output><jupyter_text>With observatioal errors. We add the observational errors for both JASMINE observations and reference stars, Gaia stars. We first add the position error + displacement for observed_xy_overlap. Then, later we will add the noise to observed_xy (all observations). The displacement for the same observation of stars should be the same between observed_xy and observed_xy_overlap. However, for now for the simplification of set up, we use independent one.<jupyter_code># JASMINE pixel position uncertainty, let's set to 1/300 pix
pix_size = 15.*u.um
xy_error_jasmine = (1.0/300)*pix_size
print(' JASMINE pix error (um) =', xy_error_jasmine)
# Reference stars ra, dec error, let's set to 0.02 mas
radec_error_refstars = (0.2*u.mas).to(u.deg)
print(' Reference stars error (deg) =', radec_error_refstars)
# add errors to JASMINE pix position
# for overlap stars
observed_xy_overlap.rename_column('x', 'x0')
observed_xy_overlap.rename_column('y', 'y0')
observed_xy_overlap.add_column(observed_xy_overlap['x0'], name='x')
observed_xy_overlap.add_column(observed_xy_overlap['x0'], name='y')

observed_xy_overlap['x'] = np.random.normal(observed_xy_overlap['x0'], xy_error_jasmine)
observed_xy_overlap['y'] = np.random.normal(observed_xy_overlap['y0'], xy_error_jasmine)
# store the noise
observed_xy_overlap.add_column(observed_xy_overlap['x'], name='xy_err')
observed_xy_overlap['xy_err'] = xy_error_jasmine
# for all stars
observed_xy.rename_column('x', 'x0')
observed_xy.rename_column('y', 'y0')
observed_xy.add_column(observed_xy['x0'], name='x')
observed_xy.add_column(observed_xy['x0'], name='yt')

observed_xy['x'] = np.random.normal(observed_xy['x0'], xy_error_jasmine)
observed_xy['y'] = np.random.normal(observed_xy['y0'], xy_error_jasmine)
observed_xy.add_column(observed_xy['x'], name='xy_err')
observed_xy['xy_err'] = xy_error_jasmine
# add errors to reference stars
radec_refstars.rename_column('ra', 'ra0')
radec_refstars.rename_column('dec', 'dec0')
radec_refstars.add_column(radec_refstars['ra0'], name='ra')
radec_refstars.add_column(radec_refstars['dec0'], name='dec')


# print(' ra before noise =', radec_refstars['ra'])
radec_refstars['ra'] = np.random.normal(radec_refstars['ra0'], radec_error_refstars)
radec_refstars['dec'] = np.random.normal(radec_refstars['dec0'], radec_error_refstars)
# print(' ra w/added noise =', radec_refstars['ra'].to_value(u.mas))
# store the noise
radec_refstars.add_column(radec_refstars['ra'], name='radec_err')
radec_refstars['radec_err'] = radec_error_refstars<jupyter_output><empty_output><jupyter_text>We define the model function. We use effecitve RA $\alpha_* = \alpha \cos(\delta)$, when we compute residuals. Note that the errors are not needed to be adjusted with effective RA.<jupyter_code>def model_wrefs_werr_func(params, n_fields, dim_sip, observed_xy, radec_refstars):
    # params = (ra_ptgs, dec_ptgs, pa_ptg..., scale, a..., b...)
    ra_ptgs, dec_ptgs, pa_ptgs, scale, a, b =\
        np.split(params, [n_fields, 2*n_fields, 3*n_fields, 3*n_fields+1,\
                          3*n_fields+1+(dim_sip+1)**2])

    field_params = Table(data=[ra_ptgs, dec_ptgs, pa_ptgs, -np.ones(shape=(np.size(ra_ptgs)))],\
                         names=['ra', 'dec', 'pa', 'field'],\
                         dtype=['float64', 'float64', 'float64', 'int64'])
#                         names=['ra', 'dec', 'pa', 'field'],\
#                         units=[u.deg, u.deg, u.deg, None],\
#                         dtype=['float64', 'float64', 'float64', 'int64'])
    field_params['ra'].unit = u.deg
    field_params['dec'].unit = u.deg
    field_params['pa'].unit = u.deg
    field_params['field'] = np.arange(0, np.size(field_params))
    
    # use copy of observed_xy
    observed_xy_cp = observed_xy.copy()
    
    a_matrix = np.reshape(a, (dim_sip+1, dim_sip+1))
    b_matrix = np.reshape(b, (dim_sip+1, dim_sip+1))
    
    # mns = np.concatenate(((0, 1), np.arange(dim_sip+1, 2*(dim_sip)+1))) 
    # for mn in mns:
    #    for m in range(np.max([0, mn-dim_sip]), np.min([mn+1, dim_sip+1])):
    #        n = mn - m
    #        ap_matrix[m, n] = 0
    #        bp_matrix[m, n] = 0
    a_matrix[0, 0] = 0.0
    a_matrix[0, 1] = 0.0
    a_matrix[1, 0] = 0.0
    b_matrix[0, 0] = 0.0
    b_matrix[0, 1] = 0.0
    b_matrix[1, 0] = 0.0

    # normalisation. 
    m, n = np.indices((dim_sip+1, dim_sip+1))
    mn = m + n
    a_matrix = a_matrix * (1.e-3**mn)
    b_matrix = b_matrix * (1.e-3**mn)
    
    # compute ra/dec from x/y with the parameters.
    ret = radec_calculator_ab(observed_xy_cp, field_params, scale[0], \
                                                              a_matrix, b_matrix)
    
    observed_xy_cp['ra_est'] = ret['ra_est']
    observed_xy_cp['dec_est'] = ret['dec_est']
    
    # compute the mean ra/dec for unique stars
    cat_ids = unique(observed_xy_cp, 'catalog_id')['catalog_id']

    ra_mean = np.zeros(len(observed_xy['ra_est']))
    dec_mean = np.zeros(len(observed_xy['dec_est']))

    # compute weights from error in xy (um) -> radec (deg)
    w_observed_xy = 1.0/(observed_xy_cp['xy_err']*scale[0])
    
    # effective RA: error does not need to adjusted.
    observed_xy_cp['ra_est'] = observed_xy_cp['ra_est']*np.cos(observed_xy_cp['dec_est'].data*u.deg)

    for i in cat_ids:
        pos = np.where(observed_xy_cp['catalog_id'] == i)
        ra_mean[pos] = np.average(observed_xy_cp[pos]['ra_est'], \
                                 weights=w_observed_xy[pos])*u.deg
        dec_mean[pos] = np.average(observed_xy_cp[pos]['dec_est'], \
                                 weights=w_observed_xy[pos])*u.deg

    # effective RA for reference stars
    radec_refstars_cp = radec_refstars.copy()
    radec_refstars_cp.rename_column('ra', 'ra*')
    radec_refstars_cp['ra*'] = radec_refstars_cp['ra*']*np.cos(radec_refstars_cp['dec'].data*u.deg)
    # print(radec_refstars['ra'], radec_refstars_cp['ra*'])

    # reference stars' measured mean ra, dec to be compared 
    # with the ra, dec of reference stars. 
    radec_est_refstars = radec_refstars.copy()
    radec_est_refstars.rename_column('ra', 'ra*_est')
    radec_est_refstars.rename_column('dec', 'dec_est')
    # compute weights for reference stars
    w_refstars = 1.0/(radec_refstars['radec_err'])
    
    for i,id in enumerate(radec_refstars['catalog_id']):
      # print('i, id=', i, id)
      # print(ra_mean[observed_xy_cp['catalog_id'] == id][0])
      radec_est_refstars[i]['ra*_est'] = ra_mean[observed_xy_cp['catalog_id'] == id][0]
      radec_est_refstars[i]['dec_est'] = dec_mean[observed_xy_cp['catalog_id'] == id][0]

    radec_est = np.concatenate((observed_xy_cp['ra_est'], observed_xy_cp['dec_est'], \
                              radec_refstars_cp['ra*'], radec_refstars_cp['dec']))
    radec_est_mean = np.concatenate((ra_mean, dec_mean, radec_est_refstars['ra*_est'], \
                                    radec_est_refstars['dec_est']))
    w_all = np.concatenate((w_observed_xy, w_observed_xy, w_refstars, w_refstars))
    
    # print(w_all)
    
    residuals = w_all*(radec_est - radec_est_mean)

    return residuals<jupyter_output><empty_output><jupyter_text>Let's run least squares.<jupyter_code>from scipy.optimize import least_squares
from scipy.optimize import leastsq
import time

dim_sip = 4
a = np.zeros(shape=(dim_sip+1, dim_sip+1))
b = np.zeros(shape=(dim_sip+1, dim_sip+1))

# constructing a_init (initial parameter set).
a_init = np.array(np.concatenate((field_params['ra'], field_params['dec'], \
                                  field_params['pa'], \
                                 [plate_scale.value], a.flatten(), b.flatten()))) # This must be an ndarray.
# a_init = np.array(np.concatenate((field_params['pa'], \
#                                  [plate_scale.value], a.flatten(), b.flatten()))) # This must be an ndarray.
print(' # of fitting parameters =', len(a_init))

print(' size of reference stars =', np.size(radec_refstars['catalog_id']))

start = time.time()
result = least_squares(model_wrefs_werr_func, a_init, loss='linear', args= \
                        (np.size(field_params), dim_sip, observed_xy_overlap, \
                        radec_refstars), verbose=2)
# result = least_squares(model_wrefs_werr_func, a_init, args= \
#                       (np.size(field_params), dim_sip, observed_xy_overlap, \
#                        radec_refstars))
print(' time=',time.time()-start)<jupyter_output># of fitting parameters = 63
 size of reference stars = 10
   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality   
       0              1         4.5266e+10                                    6.01e+17    
       1              2         1.3364e+07      4.53e+10       8.06e+02       8.74e+15    
       2              3         1.2903e+04      1.34e+07       1.61e+03       4.40e+10    
       3              4         1.2903e+04      1.83e-01       3.22e+03       8.94e+09    
       4             18         1.2903e+04      2.39e-02       4.80e-05       9.47e+08    
       5             19         1.2903e+04      0.00e+00       0.00e+00       9.47e+08    
`xtol` termination condition is satisfied.
Function evaluations 19, initial cost 4.5266e+10, final cost 1.2903e+04, first-order optimality 9.47e+08.
 time= 34.5712468624115<jupyter_text>Checking results Extract the results.<jupyter_code>n_fields = np.size(field_params)
n_objects = np.size(radec_est)
true_ra_ptgs = true_field_params['ra'].data
true_dec_ptgs = true_field_params['dec'].data
true_pa_ptgs = true_field_params['pa'].data

ra_ptgs, dec_ptgs, pa_ptgs, scale, a, b =\
       np.split(result.x, [n_fields, 2*n_fields, 3*n_fields, 3*n_fields+1,\
                           3*n_fields+1+(dim_sip+1)**2])
# pa_ptgs, scale, a, b =\
#       np.split(result.x, [n_fields, n_fields+1,\
#                           n_fields+1+(dim_sip+1)**2])
#ra_ptgs = field_params['ra'].data
# dec_ptgs = field_params['dec'].data

print(' a and b matrices before scaling=', a, b)

a_matrix = np.reshape(a, (dim_sip+1, dim_sip+1))
b_matrix = np.reshape(b, (dim_sip+1, dim_sip+1))
# A/B scaling
m, n = np.indices((dim_sip+1, dim_sip+1))
mn = m + n
a_matrix = a_matrix * (1.e-3**mn)
b_matrix = b_matrix * (1.e-3**mn)

fit_field_params = Table(data=[ra_ptgs, dec_ptgs, pa_ptgs, -np.ones(shape=(np.size(ra_ptgs)))],\
                         names=['ra', 'dec', 'pa', 'field'],\
                         dtype=['float64', 'float64', 'float64', 'int64'])
fit_field_params['ra'].unit = u.deg
fit_field_params['dec'].unit = u.deg
fit_field_params['pa'].unit = u.deg
fit_field_params['field'] = np.arange(0, np.size(field_params))<jupyter_output>a and b matrices before scaling= [ 0.00000000e+00  0.00000000e+00 -3.97395033e-03 -4.90471591e-05
  1.59703913e-06  0.00000000e+00 -1.64052413e-03  1.58223449e-04
  6.67437986e-07 -1.58068873e+01 -1.04891186e-02 -5.51877616e-05
  1.77575220e-06  8.24310771e+02 -6.40165945e+02  1.83186442e-04
  5.20766875e-07  2.56289441e+02  1.42868197e+03 -2.97605857e+01
  7.55774468e-07  7.52113934e+02 -1.10783978e+03 -4.59397333e+01
  1.00740222e+03] [ 0.00000000e+00  0.00000000e+00 -4.67719155e-03 -9.31279501e-05
  7.55167787e-08  0.00000000e+00 -7.99102030e-03  7.25312767e-05
  1.51684733e-06  9.53929242e+02 -1.51568538e-03  9.96867287e-05
  2.33145432e-06 -4.37074033e+02  8.92234996e+01  3.52660007e-04
  1.83594733e-06 -7.04127471e+02  9.76867472e+01  1.54897834e+03
 -9.51777298e-07  6.91948839e+02  1.89079607e+02  1.83323321e+03
  3.48971122e+02]<jupyter_text>Evaluate fitting. We follow https://www.fixes.pub/program/444521.html.<jupyter_code>from scipy import linalg, optimize

chi2dof= np.sum(result.fun**2)/(result.fun.size -result.x.size)
print(' Xi^2/dof =', chi2dof)

J= result.jac
print(' shape of J =', np.shape(J))
# this does not work.
# cov= np.linalg.inv(J.T.dot(J))
# var= np.sqrt(np.diagonal(cov))
# print(' parameter variances =', var)

U, s, Vh= linalg.svd(result.jac, full_matrices=False)
tol= np.finfo(float).eps*s[0]*max(result.jac.shape)
w= s > tol
cov= (Vh[w].T/s[w]**2) @ Vh[w]  # robust covariance matrix
cov *= chi2dof
perr= np.sqrt(np.diag(cov))     # 1sigma uncertainty on fitted parameters

# extract errors

ra_ptgs_err, dec_ptgs_err, pa_ptgs_err, scale_err, a_err, b_err =\
       np.split(perr, [n_fields, 2*n_fields, 3*n_fields, 3*n_fields+1,\
                           3*n_fields+1+(dim_sip+1)**2])

# A/B scaling
a_err_matrix = np.reshape(a_err, (dim_sip+1, dim_sip+1))
b_err_matrix = np.reshape(b_err, (dim_sip+1, dim_sip+1))
# A/B scaling
m, n = np.indices((dim_sip+1, dim_sip+1))
mn = m + n
a_err_matrix = a_err_matrix * (1.e-3**mn)
b_err_matrix = b_err_matrix * (1.e-3**mn)

print(' parameter values =', ra_ptgs, dec_ptgs, pa_ptgs, scale, a_matrix, b_matrix)
print(' parameter variances =', ra_ptgs_err, dec_ptgs_err, pa_ptgs_err, scale_err, \
     a_err_matrix, b_err_matrix)<jupyter_output>Xi^2/dof = 14.891181458481045
 shape of J = (1796, 63)
 parameter values = [265.62027541 265.70083679 265.48943621 265.56771628] [-28.85395491 -28.74322559 -28.78374369 -28.67008935] [302.01809948 301.17509424 300.93140956 301.57520988] [7.85067393e-06] [[ 0.00000000e+00  0.00000000e+00 -3.97395033e-09 -4.90471591e-14
   1.59703913e-18]
 [ 0.00000000e+00 -1.64052413e-09  1.58223449e-13  6.67437986e-19
  -1.58068873e-14]
 [-1.04891186e-08 -5.51877616e-14  1.77575220e-18  8.24310771e-13
  -6.40165945e-16]
 [ 1.83186442e-13  5.20766875e-19  2.56289441e-13  1.42868197e-15
  -2.97605857e-20]
 [ 7.55774468e-19  7.52113934e-13 -1.10783978e-15 -4.59397333e-20
   1.00740222e-21]] [[ 0.00000000e+00  0.00000000e+00 -4.67719155e-09 -9.31279501e-14
   7.55167787e-20]
 [ 0.00000000e+00 -7.99102030e-09  7.25312767e-14  1.51684733e-18
   9.53929242e-13]
 [-1.51568538e-09  9.96867287e-14  2.33145432e-18 -4.37074033e-13
   8.92234996e-17]
 [ 3.52660007e-13  1.83594733e-18 -7.04127471e-13  9.76867472e-1[...]<jupyter_text>Pointings RA, DEC, position angle and scale<jupyter_code>print(' pointing centre (fit) ra, dec (deg) =', ra_ptgs, dec_ptgs)
print(' pointing centre (true) ra, dec (deg) =', true_ra_ptgs, true_dec_ptgs)
print(' difference ra, dec (deg) =', ra_ptgs-true_ra_ptgs, dec_ptgs-true_dec_ptgs)
print(' uncertainty ra, dec pointings =', ra_ptgs_err, dec_ptgs_err)
print(' position angle (fit) (deg) =', pa_ptgs)
print(' position angle (true) (deg) =', true_pa_ptgs)
print(' difference =', pa_ptgs-true_pa_ptgs)
print(' uncertainty =', pa_ptgs_err)
print(' scale (fit, true) =', scale, (1e-6/7.3/np.pi*180.0)*u.deg/u.um)
print(' difference =', scale-(1e-6/7.3/np.pi*180.0))
print(' uncertainty =', scale_err)<jupyter_output>scale (fit, true) = [7.85067393e-06] 7.848736919600318e-06 deg / um
 difference = [1.93700585e-09]
 uncertainty = [7.35420055e-12]<jupyter_text>Objects positions<jupyter_code>radec_objs = radec_calculator_ab(observed_xy_overlap, fit_field_params, scale[0], a_matrix, b_matrix)

plt.scatter(true_radec_overlap['ra'], true_radec_overlap['dec'], marker='x', label='True')
plt.scatter(radec_est['ra_est'], radec_est['dec_est'], marker='+', label='Initial guess')
plt.scatter(radec_objs['ra_est'], radec_objs['dec_est'], marker='.', label='Final estimation')
plt.scatter(radec_refstars['ra0'], radec_refstars['dec0'], marker='o', \
            label='Reference stars')
plt.xlabel('ra (deg)')
plt.ylabel('dec (deg)')
plt.title('Object positions')
plt.legend()
distlist = []
print(np.shape(radec_objs))
for i in range(0, np.size(radec_objs)):
    c1 = SkyCoord(radec_objs['ra_est'][i]*u.deg, radec_objs['dec_est'][i]*u.deg)
    c2 = SkyCoord(observed_xy_overlap['ra'][i]*u.deg, observed_xy_overlap['dec'][i]*u.deg)
    distlist.append(c1.separation(c2).arcsec)
distlist = np.array(distlist)

#plt.hist(np.log10(distlist))
plt.hist(distlist)
plt.xlabel("Residual (arcsec)")
plt.ylabel("Number")
# dra = ((radec_objs['ra_est']-observed_xy_overlap['ra']).data)*u.deg
# effective RA
dra = radec_objs['ra_est']*np.cos(radec_objs['dec_est'].data*u.deg) \
        -observed_xy_overlap['ra']*np.cos(observed_xy_overlap['dec'].data*u.deg)
ddec = ((radec_objs['dec_est']-observed_xy_overlap['dec']).data)*u.deg
dra_arcsec = dra.to_value(u.arcsec)
ddec_arcsec = ddec.to_value(u.arcsec)
plt.scatter(dra_arcsec, ddec_arcsec, marker='x')
plt.xlabel('dRA (arcsec)')
plt.ylabel('dDEC (arcsec)')<jupyter_output><empty_output><jupyter_text>Apply to all the data, taking into account uncertainties of their position and parameter uncertainties. We shall run Monte Carlo by randomly displacing the position of stars and distortion parameters.<jupyter_code>n_mc = 100

n_stars = len(observed_xy)
print(' total # of stars =', n_stars)
ra_allobjs_samp = np.empty((n_stars, n_mc))
dec_allobjs_samp = np.empty((n_stars, n_mc))

observed_xy_try = observed_xy.copy()
# flattened uncertainties of a, b matrix
a_flat = a_matrix.flatten()
b_flat = b_matrix.flatten()
a_err = a_err_matrix.flatten()
b_err = b_err_matrix.flatten()

for i in range(n_mc):
    # displace observed_xy positions
    observed_xy_try['x'] = np.random.normal(observed_xy['x'], observed_xy['xy_err'])
    observed_xy_try['y'] = np.random.normal(observed_xy['y'], observed_xy['xy_err'])

    # displace the parameters
    ra_ptgs_try = np.random.normal(ra_ptgs, ra_ptgs_err)
    dec_ptgs_try = np.random.normal(dec_ptgs, dec_ptgs_err)
    pa_ptgs_try = np.random.normal(pa_ptgs, pa_ptgs_err)
    scale_try = np.random.normal(scale, scale_err) 
    a_try = np.random.normal(a_flat, a_err)
    b_try = np.random.normal(b_flat, b_err)
    a_matrix_try = np.reshape(a_try, (dim_sip+1, dim_sip+1))
    b_matrix_try = np.reshape(b_try, (dim_sip+1, dim_sip+1))
    
    fit_field_params_try = Table(data=[ra_ptgs_try, dec_ptgs_try, pa_ptgs_try, \
                                   -np.ones(shape=(np.size(ra_ptgs)))],\
                         names=['ra', 'dec', 'pa', 'field'],\
                         dtype=['float64', 'float64', 'float64', 'int64'])
    fit_field_params_try['ra'].unit = u.deg
    fit_field_params_try['dec'].unit = u.deg
    fit_field_params_try['pa'].unit = u.deg
    fit_field_params_try['field'] = np.arange(0, np.size(field_params))
    
    radec_allobjs_try = radec_calculator_ab(observed_xy_try, fit_field_params_try, \
                                    scale_try[0], a_matrix_try, b_matrix_try)
    # use effective RA
    ra_allobjs_samp[:, i] = radec_allobjs_try['ra_est']*np.cos(radec_allobjs_try['dec_est'].data*u.deg)
    # ra_allobjs_samp[:, i] = radec_allobjs_try['ra_est']
    dec_allobjs_samp[:, i] = radec_allobjs_try['dec_est']
    

ra_allobjs_mean = np.mean(ra_allobjs_samp, axis=1)*u.deg
ra_allobjs_std = np.std(ra_allobjs_samp, axis=1)
dec_allobjs_mean = np.mean(dec_allobjs_samp, axis=1)*u.deg
dec_allobjs_std = np.std(dec_allobjs_samp, axis=1)

# we convert true obserfved ra to effective ra
raeff_observed_xy = observed_xy['ra']*np.cos(observed_xy['dec'].data*u.deg)
# raeff_observed_xy = observed_xy['ra']
dec_observed_xy = observed_xy['dec']

print(' mean dec correction for effective RA =', np.mean(np.cos(observed_xy['dec'].data*u.deg)))

# print(raeff_observed_xy, ra_allobjs_mean)

# error from the true value
ra_allobjs_err = ra_allobjs_mean-raeff_observed_xy
dec_allobjs_err = dec_allobjs_mean-dec_observed_xy

plt.scatter(ra_allobjs_err.to_value(u.arcsec), (ra_allobjs_std*u.deg).to_value(u.arcsec), marker='x', label='RA')
plt.scatter(dec_allobjs_err.to_value(u.arcsec), (dec_allobjs_std*u.deg).to_value(u.arcsec), marker='.', label='DEC')

print(' RA mean standard deviation of measurements (arcsec) =', \
      (np.mean(ra_allobjs_std)*u.deg).to_value(u.arcsec))
print(' RA standard deviation from the true values (arcsec) =', 
      (np.std(ra_allobjs_err)).to_value(u.arcsec))
print(' DEC mean standard deviation of measurements (arcsec) =', 
      (np.mean(dec_allobjs_std)*u.deg).to_value(u.arcsec))
print(' DEC standard deviation from the true values (arcsec)=', 
      (np.std(dec_allobjs_err)).to_value(u.arcsec))

plt.xlabel('deviatoin from the true radec (arcsec)')
plt.ylabel('standar deviation of measurement (arcsec)')
plt.title('Object positions')
plt.legend()<jupyter_output>total # of stars = 2563
 mean dec correction for effective RA = 0.876660634963193
 RA mean standard deviation of measurements (arcsec) = 0.006280788912085892
 RA standard deviation from the true values (arcsec) = 0.02474307609116223
 DEC mean standard deviation of measurements (arcsec) = 0.003751085599114781
 DEC standard deviation from the true values (arcsec)= 0.014589443344355368
StarcoderdataJupyterStructuredCleanDedup