Study of pointing accuracy of 3.6m Devasthal Optical Telescope

Back

Study of pointing accuracy of 3.6m Devasthal Optical

Telescope

Mohan Agrawal

July 9, 2017

Contents

1 Introduction 2

1.1 Ways to point a telescope . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 About DOT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Corrections involved 4

3 Modeling Process 7

3.1 Kinematic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3.2 Geometrical errors in kinematic model . . . . . . . . . . . . . . . . . . . 10

4 Model Fitting 12

4.1 Multiple Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4.1.1 Criteria of estimate . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4.1.2 The Gauss-Markov assumptions . . . . . . . . . . . . . . . . . . . 14

4.1.3 The Gauss-Markov Theorem and related properties . . . . . . . . 14

4.1.4 ANOVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

4.1.5 Hypothesis Testing . . . . . . . . . . . . . . . . . . . . . . . . . . 18

4.2 Regressing the pointing data . . . . . . . . . . . . . . . . . . . . . . . . . 18

5 Obtaining statistics 20

5.1 Residual analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

5.2 Revisiting assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

5.3 Model validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

6 Practical considerations 29

6.1 Oﬀ axis imager . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

6.2 Star selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

6.3 Sources of uncertainity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

7 Conclusion 34

8 Acknowledgements 34

Note: All the data, code and graphs used in the report can be found at

https://github.com/mohanagr/pointing-model

Abstract

Pointing accuracy is the ﬁrst and foremost requirement for any astronomical observation.

With this report an eﬀort has been made to study and document various ways of modeling

an alt-az telescope’s pointing accuracy, obatining statistical inferences from the ﬁt models,

peering into the physical causes of the modeled errors and providing a comprehensive list

of shortcomings and limitations of the approach used. The study has been made using

the data obtained from 3.6m Devasthal Telescope. A kinematic model has been used to

model the pointing errors, vector form of least square method used to ﬁt the model and

various statistical parameters obtained that allow one to asses the goodness of ﬁt of the

model and provide solution to some common pointing diﬃculties that may arise during

an observation e.g. the imager being used is not aligned with the telescope’s rotator’s

axis. At the end several other sources/causes of error that may aﬀect the accuracy, have

been summarized.

1 Introduction

With amateur telescopes when one wants to look at a bright object, say a planet, all that

needs to be done is, swing the telescope to the location where the planet is approximately

located and manually move it until it is focused in the eyepiece. But what if the target

object is faint? i.e. a deep sky object. Then the usual method is to use a ﬁnding-chart.

This approach (although still being used at some places) is highly ineﬃcient simply

because of the sheer amount of wasted time.

Instead, large and professional telescopes have stepper motors & encoders (either

absolute or incremental) and have physical readouts. Here all an observer needs to do

is enter the desired location, usually as (RA, Dec), and the telescope (hopefully) points

to that direction. Although this approach points the telescope with minimum eﬀort, it

is riddled with errors - the celestial target may or may not appear in the FOV (Field

Of View) and even if it does, it is certain that it will not appear in the center of the

FOV/crosshair. Taking the simplest of examples, it is possible that the encoders have

some zero-errors or the bearings have runout. The direction where one wants to look

is diﬀerent from the direction the telescope’s looking. Therefore there arises a need

to introduce some correcting oﬀset in the values that are to be input. Or still better,

generate a model that takes into account various errors that can be present (both static

and dynamic), takes the desired coordinates as input and gives the corrected coordinates

as output - the coordinates which should actually be entered in order for the image to be

centered in the FOV.

1.1 Ways to point a telescope

The term ’pointing’ itself refers to the initial acquisition of the target. Depending on the

type of target (bright/faint) and the accuracy demands (large observatories often need

sub-arcsecond pointing accuracies) there are diﬀerent ways to point the telescopes.

1. Blind pointing : Also known as dead reckoning, this is the most basic way wherein

one enters the target coordinates and the telescope moves (usually to some diﬀerent

coordinates since certain errors always present). This is generally used when the

accuracy required is not too demanding (several arc seconds is acceptable) e.g. in

case of wide-ﬁeld imaging.

2. Oﬀset pointing : When the pointing demands are higher than what the system is

capable of, a common way is to introduce a constant oﬀset in the input values. This

oﬀset can be calculated by ﬁrst entering the actual coordinates of the (bright) object,

then centering the object in the FOV and noting down the resultant coordinates

shown by telescope readouts. The diﬀerence then gives the oﬀset. Using this oﬀset

for the rest of the observations often gives a decent performance with absolute

accuracy being less than 10

and often nears 1

− 2

if the observations are being

made only over a small patch of the sky.

3. Blind oﬀset pointing : When the target object itself is faint, an oﬀset is calculated

as above using a bright object in the vicinity of the target object.

4. Pointing models : This method involves constructing a mathematical model using

a set of calibration stars that can then predict correct encoder demands for any

future observations. Pointing models are used in all modern TCS (Telescope Control

Systems) and are capable of catering to sub-arcsecond pointing demands. The full-

blown pointing models (consisting of about 20 terms) generated for the 3.6m DOT

(Devasthal Optical Telescope) give a pointing accuracy of 0.5

− 0.6

This report is concerned with the development and usage of method no. 4. A pointing

model has been ﬁt manually and the subsequent results are compared with that of the

results obtained from TPOINT software, to assess the performance of TPOINT’s model

and to calculate certain statistics that TPOINT does not provide. A python script

has also been written to compare various TPOINT models and generate bar plots to

enable the user to have a quick look at the performance of various models and choose

appropriately.

1.2 About DOT

The 3.6m Devasthal Optical Telescope (DOT) is an alt-azimuth, Ritchey −Cr´etien type

telescope. The optical diagram is ﬁgure 1. The ’Main Port’ sits at FPA axial. There

are two FPA-side, holding the two ’Side Ports’. The pointing data used here has been

obtained from a CCD imager mounted at Main Port. Quick comparison can be done for

pointing models for diﬀerent ports (although not much variation is expected) as shown

in ﬁgure 19 later in the report.

Figure 1: Optical diagram for the telescope featuring hyperbolic primary and secondary

mirrors. The focal length with the corrector is 32795.8mm. Image adapted from telescope

manuals.

2 Corrections involved

Before proceeding to build a model, it is necessary to know various corrections have to

be considered by a pointing model. These corrections can be broadly categorized into:

1. Corrections for apparent position : These include eﬀects of gravitonal ﬁelds

of Sun and Moon on the earth - precession an nutation, proper motion of the

stars, annual parallax (signiﬁcant for nearby stars), light deﬂection (because of

gravitational lensing by the Sun) and Earth’s rotation. These are the geometrical

corrections in the positions of the stars on the celestial sphere i.e. not for someone

on Earth.

2. Corrections for topocentric position : Topocentric means the positions which

will be observed by an observer on earth. These include correcting for diurnal

abberation and accounting for atmospheric refraction eﬀect.

3. Telescope speciﬁc corrections : These corrections account for the imperfections

of the telescope (to be discussed in subsequent sections). Their aim is to mop up

any signiﬁcant diﬀerence left between actual coordinates and the encoder demands

after the aforementioned corrections have been applied.

The corrections are shown in the following ﬁgure

[1]

Figure 2: The list of corrections that lead from the actual star coordiantes obtained from

a catalogue to the ’instrumental place’ - coordinates that will be actually input in the

telescope.

The apparent [α, δ] is obtained after applying celestial corrections. RA is then con-

verted to HA (Hour angle) using LST (Local Sidereal Time). Apparent [HA, Dec] then

gets corrected for earth speciﬁc eﬀects including refraction and is called the Observed

Place [h, δ]. This is then converted to Horizontal Coordinate System i.e. [Az, El] repre-

senting Azimuth & Elevation. These alt-az coordinates then undergo corrections by the

telescope’s pointing model (e.g. the one generated by TPOINT) and the ﬁnal result is

used to orient the telescope, called the Instrument Place

[2]

Some deﬁnitions consider ’apparent place/position’ to be the one which is observed

from earth, here that is called ’observed place’. ’Apparent position’ is the one which is

obtained after considering the eﬀect of proper motion - from catalogue epoch to present

date; eﬀects of ’general precession’ which is the well-known precession of Earth’s pole

around the pole of the ecliptic with a 26, 000 year period (consequently the vernal equinox

moves around the celestial equator, advancing by 50.25

/year) which produces an eﬀect

on star positions of upto 50

(also being the reason as to why most catalogues specify

an epoch); and the nutation eﬀects (residual wobbles, generally because of moon) which

aﬀect the pointing by at most 10

. These eﬀects combined, give us the ’apparent place’.

Eﬀects such as light deﬂection (< 2

), Diurnal abberration(< 0.3

) and parallax (< 1

for nearest of the stars) can be neglected. Further most TCS’s allow the user to enter

eﬀects such as parallax etc., if user does want them to be considered.

While calculating observed place, the most signiﬁcant contribution is by that of re-

fraction. As shown below, continuous variation in the refractive index of the air raises

the apparent position of the star.

Figure 3: The eﬀect of refraction on pointing

Usually, this eﬀect scales with cot ζ (where ζ is the distance to zenith, 90

◦

−Elevation)

and can be as large as 1

for ζ = 45

◦

. This eﬀect raises the z coordiante of the aim vector

(position vector of the star in the sky) and the change is given by

[3]

∆z = K/z

where K(P, T ) ≈ 2.77 ×10

−7

P × (1.0419 − 0.00380T )

P is the pressure at the observatory in mbar and T is the ambient temperature in

◦

For DOT this calculation is done by the TCS.

Most catalogues give what is called here the ’apparent place’ which is then input into

the TCS and steps 2 & 3 are computed; even if apparent place is not available, the TCS

will give user the option to specify the epoch of the observation, parallax and proper

motion and then the TCS will compute all 3 steps.

This study is concerned speciﬁcally with step 3. The data that is used to ﬁt the

observations looks as shown in the following ﬁgure:

Figure 4: Format of data used to ﬁt pointing models. The two leftmost values represent

Apparent (Az, Alt) and the two rightmost represent Instrument (Az, Alt), the one shown

by encoders

It is important to note that the diﬀerence of the two sets of values shown above is

assumed to have come solely from telescope’s ineﬃcacies; which is to be corrected by step

3. In the sections that follow, this step will be discussed in detail, as to how these errors

are modeled. For calculating eﬀect of precession, nutation, aberration etc. explicitly one

may refer [2].

3 Modeling Process

The pointing error is represented as a pair (∆Az, ∆E). In order to build a model, we

need for formulate those errors i.e. we need to be able to represent the errors as:

∆A = f(A, E) + η (1)

∆E = g(A, E) + ξ (2)

η and ξ represent inherent irreducible random error in the model.

The goal of modeling is to generate a deterministic model that relates telescope mount

axis positions to the pointing error. Although no random variables are involved in the

above representation, the model is not truly deterministic. Some stochastic variables are

either excluded or too diﬃcult to integrate into the model. They contribute as a random

error e.g. wind eﬀects, local thermal variations etc. It is the choice of f and g that

diﬀerentiates the two common approaches.

One approach is to ﬁt a spherical harmonic function model for both the errors. Zonal

harmonics are taken upto 4

order and the tesseral terms only to the ﬁrst order. The

function may look something like:

f(A, E) = a

+ a

sin A + a

cos A cos E + a

sin A cos E + a

sin

+ a

cos A sin E cos E + a

sin A sin E cos E + a

sin

E...

It has been shown

[4]

that such a model lacks stability, often has high correlation

between the terms and performs badly when extrapolated. Further from the point of

view of a non-statiscian it provides absolutely no physical insights i.e. what caused those

errors or relating the errors with physical misalignment/errors of the telescope.

Other, more popular approach is to use a kinematic model. In such a model we derive

expressions for diﬀerent types of geometrical errors, misalignments that might actually

be present in the telescope and their eﬀect on pointing i.e. contribution to (∆A, ∆E).

3.1 Kinematic model

A kinematic model incorporates geometrical errors that can be induced in a telescope

because of mechanical misalignments, ﬂexures and low order bearing eﬀects. Former two

are by far the dominating eﬀects. There is nothing emprical in such a model, every error

has a physical basis, which is why as will be seen later, every term in the model comes

out to be highly signiﬁcant. (very low p-values in hypothesis testing).

A telescope employs a two degree-of-freedom motion system (rotations about Az and

El axes) which is why the basic kinematic model also consists of two rotation motions

only. A number of errors are possible with such a system.

(a) Non perpendicularity

(b) Collimation error

Figure 5: Figure representing diﬀerent geometrical errors. Blue vector in (a) represents

original axis and green shows the tilted one. Red vector shows the originally expected

pointing direction and the brown one shows the displaced direction. (Image adapted from

wildacard-innovations.com)

• The Azimuth and Elevation axes might not be perpendicular.

• There might be collimation error. The optical path might not be parallel to tele-

scope’s boresight e.g. the optical tube is not perpendicular to the elevation axis.

• The telescope pier might be tilted (NS or EW).

• Telescope tube might droop down because of gravitational loading. (ﬂexure)

In a kinematic model, the errors each of the above mentioned factors introduce in

(∆A, ∆E) can be calculated using 3-D geometry, which then form the required f(A, E)

and g(A, E).

There are various ways to represent a kinematic model. Here we will use a rotation-

only derivative of D-H (Denavit-Hartenberg) representation. This is because translational

errors are hardly present in a telescope & even if they are, they hardly aﬀect pointing of

the telescope. It is the rotational alignment which is error prone and has an on-sky eﬀect.

For all practical purposes the D-H convention is equivalent to using Euler angles or more

speciﬁcally Tait-Bryan angles for rotation. They are used to represent the orientation

of a rigid body with respect to a ﬁxed coordinate system. The rotations in Tait-Bryan

system are about x-y-z axes. These rotations can be of two type:

1. Intrinsic : About XYZ axes which is attached to the rigid body (i.e. it’s orienta-

tion changes w.r.t. ground).

2. Extrinsic : About xyz axes which remains ﬁxed (to the ground).

Here we will use extrinsic rotations. It is worthwhile to note that the XYZ axes

represent telescope’s mount axes (ones which will change orientation e.g. when the mount

is tilted). For our purposes the rigid body talked about earlier will be the telescope’s

boresight or the Optical Tube Assembly (OTA) (unless speciﬁed). Rotation about xyz

axes are represented as follows:





1 0 0

0 cos x −sin x

0 sin x cos x





≈





1 0 0

0 1 −x

0 x 1





(3)





cos y 0 sin y

0 1 0

−sin y 0 cos y





≈





1 0 y

0 1 0

−y 0 1





(4)





cos z −sin z 0

sin z cos z 0

0 0 1





≈





1 −z 0

z 1 0

0 0 1





(5)

For R

the signs have been reversed compared to the other two because the rotations

are considered positive in the anti-clock direction when looking towards the axis. As

shown in ﬁgure 6, image b. −θ for y axis implies that the rotation is into the page. As

per Euler’s rotation theorem, any orientation of a rigid body in space can be achieved by

3 elemental rotations about these axes.

Figure 6: Principal axes rotations

For an ideal telescope when it’s elevation and azimuth angles are set to zero it should

point at Northern horizon, which in our case represents the pointing vector



0 1 0



It is necessary to remember that all rotations are extrinsic (about xyz, attached to

ground). Now when the telescope is moved to a position (A, E), it is given two rotations

using R

and R

such that:





cos A sin A 0

−sin A cos A 0

0 0 1





z = −A has been substituted in R

since the rotation is clockwise (towards east or

x-axis). And for rotation around elevation axis, x = E in R





1 0 0

0 cos E −sin E

0 sin E cos E





Because rotating the telescope about azimuth changes the orientation of elevation

axis, R

should be applied before R

. Hence the total rotation matrix is:

= R





cos A sin A cos E −sin A sin E

−sin A cos A cos E −cos A sin E

0 sin E cos E





(6)

and produces the ideal, undeﬂected line of sight vector P

= R





sin A cos E

cos A cos E

sin E





(7)

From any pointing vector P, azimuth (A) and elevation angles can be found out by:

A = tan

−1

(

P 1

P 2

) (8)

and

E = sin

−1

(P 3) (9)

3.2 Geometrical errors in kinematic model

Each of the error discussed in previous section produces a small rotation of the line of

sight vector. The amount by which it changes (A, E) of the line of sight can be calculated.

The geometrical errors which have been used in this model have been summarized in table

1. As an example, derivation of the non-perpendicularity error has been shown.

To avoid any confusion it is necessary to understand that while visualizing these ro-

tations, one must keep the physical eﬀects that cause these errors (as seen by someone

standing on the ground) separate from the way the ﬁnal position is achieved mathemat-

ically. E.g. consider the non perpendicularity error. E.g. it is obvious for an observer

to see that the motions of Azimuth and Elevation axes will inﬂuence each other. When

the telescope tube is raised to some elevation, it will not rise straight up instead it will

rise at an angle, as shown in ﬁgure 5(a). In order to get to this position mathematically,

we ﬁrst rotate the vector by desired amount (the elevation input) around the x axis and

Eﬀect Symbol Matrix Extrinsic rotation axis

Optic axis misalignment CA R

Tube Flexure TF R

Elevation setting E R

Az-El non perpendicularity NPAE R

Azimuth setting A R

Mount tilt(NS) AN R

Mount tilt(EW) AW R

Table 1: Summary of geometrical error eﬀects

then rotate it about y axis by the amount of non-perpendicularity (NPAE ) to get the

desired tilt eﬀect. The ﬁnal orientation of tube will be exactly the same as one would

visualize when a single rotation is provided to an imperfect telescope. The xyz axes do

not represent the telescope.

On a telescope that suﬀers from NPAE, the ﬁnal telescope vector can be given as:

= R





NP AE sin A cos E + cos E sin A

−NP AE sin E sin A + cos A cos E

sin E





(10)

Using eqn. 6 and 7:

n sin E cos A + cos E sin A

−n sin E sin A + cos E cos A

(11)

= E (12)

Although the NP error will inﬂuence the vertical as well as azimuthal orientation,

given the fact that the errors are usually very small in magnitude, R

has been taken as

given below.





1 0 NP AE

0 1 0

−NP AE 0 1





(13)

This is the reason why the eﬀect of NP over elevation can be neglected and E

comes

out to be independent of E in eqn 12.

Using the diﬀerence formula for tangent

tan(A

− A) =

tan A

− tan A

1 + tan A

tan A

(14)

and again taking the assumption that the error in Azimuth is small, letting

tan(A

− A) ≈ A

− A (15)

On substituting eqn. 11 in 14 (using 15) and simplyfying, we obtain the result:

∆A = NP AE tan E (16)

Continuting the previous discussion, similar is the case with NS and EW mount tilt

errors. Although it is obvious for someone on the ground that the mount tilt will eﬀect

Eﬀect ∆A ∆E

A encoder index error IA 0

E encoder index error 0 IE

Optic axis misalignment −CA sec E 0

Tube Flexure 0 −T F cot E

Az-El non perpendicularity NP AE tan E 0

Mount tilt(NS) AN sin A tan E AN cos A

Mount tilt(EW) AW cos A tan E AW sin A

Table 2: Error terms

both the telescope’s mount axis and the line of sight (tilting the telescope plane), we are

not considered with the telescope’s plane. Our only concern is what happens to the line

of sight because of the error. Hence we rotate the LOS vector by the respective error

amounts around x and y axis to emulate the eﬀect of mount tilt. A natural question

then arises, what if some other rotation is given after mount tilt, now that in reality, the

telescope’s alt and az axes have tilted? The answer is, although we rotate the LOS vector

about the ﬁxed extrinsic axes, we also keep in mind, what physical changes the actual

errors are making in reality. So, no rotation is given about an axis has previously been

transformed. E.g. Elevation setting is done before including NPAE since NPAE tilts the

elevation axis. The only rotations a telescope is capable of having is around Alt and Az

axes. They have already been given before the tilt eﬀect is included.

The gist of this discussion is, there is a speciﬁc order to the rotations of the line of

sight vector if the eﬀect of all the errors are to be considered.

The ﬁnal telescope vector is given as:

= R

(17)

If all the error amounts listed in table 1 were known one could in principle calculate

the actual pointing direction using eqn. 17. Since they are not, a common strategy is to

evaluate the eﬀect of each error on altitude and azimuth (as was done for NPAE, eqn.

16), combine all the eﬀects and then do some sort of statistical ﬁt using available pointing

data (ﬁgure 4) to determine the error estimates. These eﬀects are listed in table 2.

We ﬁnally obtain the relations talked about in eqn. 1 and 2.

∆A = f(A, E) = IA−CA sec E+NP AE tan E +AN sin A tan E +AW cos A tan E (18)

and

∆E = g(A, E) = IE −T F cot E + AN cos A + AW sin A (19)

∆A and

∆E has been used instead of ∆A and ∆E because here we have excluded the

random error terms η and ξ. The ’hat’ symbol represents that the functions are our

estimates of the actual error. More on this will follow in subsequent sections.

4 Model Fitting

A derivative of multiple regression model was ﬁt to the pointing data using least-squares

minimization technique. Before getting to the results obtained after the ﬁtting, it is

necessary to examine the form of the model, the assumptions involved and the er-

rors/limitations of the approach.

4.1 Multiple Regression

In a multiple regression model we have the following components:

• X being an n ×k matrix. n represents the number of observations and k represents

the number of independent variables for each of the n observation. E.g. sec E,

tan E etc. in our case for each of the observation in the dataset. (Accurate example

to follow later).

• y being a n ×1 vector of observations on the dependent variable (n diﬀerent obser-

vations of ∆A or ∆E in our case).

•  being a n × 1 vector of random disturbances (which are inherent to the model,

similar to η or ξ).

• β being a k ×1 vector of unknown population parameters that we want to estimate

(e.g. NP AE or CA).

This is can be represented as:

y = Xβ +  (20)

This is assumed to be an accurate description of the properties of the popuation

(all the stars in the sky). Our goal is to obtain estimates of β, the population

parameter vector.

4.1.1 Criteria of estimate

Our estimates for the population parameter β are given by

β or b. Now when we predict

the output ’y’ using the same values that were used to ﬁt the model, we get an estimate

of the observations of the dependent variable - ˆy i.e.

ˆy = Xb (21)

We want to minimize the diﬀerence between ˆy and y. There are many notions of such

least-distance. One of them is to minimize the summation of the square of errors known

as RSS (Residual Sum of Squares) given by

i=1

, where e

is simply y

− ˆy

In matrix form:

RSS = e

· e = (y − Xb)

· (y − Xb)

= y

y − 2b

y + b

(22)

Using matrix diﬀerentiation, diﬀerentiating above expression w.r.t. b and equating it

to 0, we get our estimates for β.

b = (X

−1

y (23)

4.1.2 The Gauss-Markov assumptions

These are the assumptions involved. These assumptions are needed to justify a number

of properties of least-squared estimators which make them such a popular choice. These

properties are listed in section 4.1.3.

It is important for one to not confuse e (residual) and  (inherent random disturbance)

and not use the results about one for the other.

1. The zero conditional-mean assumption which states that the random disturbances

 average out to 0.

E(|X) = 0 (24)

This also implies from eqn. 20 that:

E(y) = Xβ (25)

2. Homoskedasticity of the disturbances and no autocorrelation.

E( · 

|X) = σ

I (26)

 ·

simply gives the variance-covariance matrix of the disturbances. Homoskedas-

ticity implies that the disturbances have constant variance. No autocorrelation

implies that the disturbances in each of sample’s observation have no correlation

with each other. Hence the variance-covariance matrix is a diagonal matrix with

constant entries.

This assumption is important when one wants to calculate conﬁdence intervals for

the predicted values.

3. Linearity between X and y which enables us to write eqn. 20. One has badly be-

haved residuals (ideally they should be randomly distributed) when this assumption

is violated. The cure then is to transform either LHS or RHS of the equation, say

replace y by log y or

√

4. Normal distribution of the disturbances.

|X ≈ N(0, σ

) (27)

This assumption is taken to make hypothesis testing easier (obtaining p-value using

t-distribution) and calculating conﬁdence intervals for parameter estimates. It can

be justiﬁed on the basis of Central Limit Theorem.

4.1.3 The Gauss-Markov Theorem and related properties

This theorem states that the OLS (oridnary least squared) estimator is the Best Linear,

Unbiased, Eﬃecient estimator. This is why it is used for almost every ﬁtting purposes. It

is important to state this theorem and it’s properties because it is very easy in statistics

to mistake one formula for the other or use two terms interchangibly when they actually

mean diﬀerent things.

1. b is an unbiased estimator of β. That is, using eqn. 20 and 23

b = β + (X

−1

 (28)

E(b) = E(β) + (X

−1

E(X

)

= β

(29)

since E(X

) = 0

This tells us that mean of estimated value of β approaches the true value over

several samples taken from the population.

2. b is a linear function of disturbances as shown by equation 28.

3. b has minimum variance among all linear, unbiased estimators.

4.1.4 ANOVA

ANOVA refers to ’Analysis of Variances’. When a model is ﬁt to a data, the total

variability in the data can be divided in two parts: the variance which is explained by

the ﬁtted model and the variance which is left unexplained.

Figure 7: Variability in the data

SSM =

i=1

( ˆy

− ¯y) (30)

SSE =

i=1

( ˆy

− y

) (31)

SST = SSM + SSE (32)

DF M (Degrees of freedom for the model) = p − 1 (33)

DF E (Degrees of freedom for error) = n − p (34)

DF T (Total degrees of freedom) = n − 1 (35)

Hence,

MSM (Explained variance) = SSM/DF M (36)

MSE (Unexplained variance) = SSE/DF E (37)

This concept is showed in ﬁgure 7.

This is particularly important for this enables us to calculate the F-statistic or perform

a F-test, which is used to check the following null hypothesis which checks whether or

not the y values bear any signiﬁcance to the terms on RHS (feasibiity of the ﬁt).

: β

= β

= . . . β

= 0

: β

6= 0 for at least one j

F-statistic is calculated using MSM/M SE. Generally a high value of F-statistic

indicates that null hypothesis should be rejected (more speciﬁcally, one has to calculate

a (at least) 95% conﬁdence interval using F-distribution and then check if F-statistic lies

in the interval or not. If not, H

is rejected).

In our case of regressing pointing data, a ridiculously high F-statistic is obtained

(> 10

). It is expected because the terms which we are using to relate the errors are

actual, physical, geometrical errors. Such a high value shows that there is absolutely

nothing empirical in the model.

The goodness of ﬁt of the model is calculated by R

statistic which is nothing by

the ratio of SSM/SST i.e. the fraction of total variance explained by the model. It is

worthwhile noting here that, when more terms are added to the model, irrespective of

whether or not they bear any signiﬁcance, R

will increase. This is why many statisticians

prefer a statistic called adjusted R

which penalizes the model for adding more terms and

is given as

= 1 − (1 − R

)

n−1

n−p

Although R

can be calculated in our case, it too is very high (> 0.99) and bears less

signiﬁcance, more on which will be clear in next section.

By far the most important parameter of ANOVA is MSE. It is used to estimate the

variance of random disturbances in the model

[5]

- σ

= MSE =

RSS

n − p

(38)

MSE or RMSE (

√

MSE) is one absolute ﬁt of criterion. If any changes to the model

decrease the RMSE, those changes are desirable, whether it is adding/removing terms or

transforming the values. RMSE for the ﬁt model in this study was around 0.4 (better

than TPOINT estimate) arcseconds (although the way it was obtained diﬀers a little).

The variance-covariance matrix of the residuals can also be found out using MSE.

Following proof (39 through 41) which is often skipped in a lot of literature, reiterates

the importance of aforementioned assumptions.

e = y − ˆy = Xβ +  − Xb

= Xβ +  − X(X

−1

= Xβ +  − X(X

−1

(Xβ + )

= (I − H)

(39)

Matrix X(X

−1

is called the hat-maker matrix since it puts a hat on whatever it

operates on e.g. ˆy = Hy. It’s important property is that it’s idempotent and symmetric

i.e. H · H

= H.

Now,

V ar(e) = e · e

= (I − H)(I − H)

 · 

= (I − H) · 

= (I − H)σ

(Using assumption 4)

(40)

Since by eqn. 39, s

is an estimate of σ

Est. V ar(e) = s

(I − H) (41)

If our assumption is perfectly valid (which is not in most cases) eqn. 41 is a pretty

good estimate, but one will see that explicitly calculating LHS of eqn. 40 and RHS of

eqn. 41 will not yield the same result. This check has been performed in the study and

hence mentioned here.

Lastly, we are also interested in checking the variances of our parameter estimates i.e.

V ar(b|X) = E[(b − E(b|X))(b − E(b|X))

|X] (42)

= E[(b − β)(b − β)

|X] (43)

= σ

−1

(44)

The conditional variance tells us, how much variance is left if we use E(Y |X) to predict

SE (Standard Error) = Est. V ar(b|X) = s

−1

(45)

There is often a confusion between standard error and standard deviation. b is not

only an estimate of β but it is also a random variable itself, since it depends on the choice

of the sample X from the population. This is why we have used ’conditoned’ mean and

’conditioned’ variance everywhere. b(X) is an estimate (of β, as function of X) and b is a

random variable. So, the standard error of an estimate (b(x)) equals standard deviation

of the random variable (b). Standard error simply shows how far away your estimate is

from the conditoned expected value (mean of b over all possible X’s), which here, is also

the true value of b (β) as given by eqn. 29.

4.1.5 Hypothesis Testing

An example of hypothesis testing was shown in the previous section when F-test was

done. It is natural to question what if one wants to test the signiﬁcance of one particular

coeﬃcient in b. Using the method described here one can assess the signiﬁcance of a

coeﬃcient and calculate a conﬁdence interval for it.

From the discussion in previous section (esp. eqn. 28, 29 & 44) i.e. if we are willing

to assume that the disturbances() are normally distibuted, we can write the following:

b = N(β, σ

−1

) (46)

We state the following null-hypothesis:

: ˆy doesn

t depend on β

i.e. β

= 0

: β

6= 0

In order to check the signiﬁcance of β

we look at it’s estimate b

. We want to determine

if b

is far enough from 0 so that we can be conﬁdent that it is not zero. How far is far

enough? In practice, we calculate a t-statistic given by:

t =

SE(b

)

(47)

where SE(b

) is given by diagonal element (i,i) of eqn. 45. It shows the no. of standard

deviations that b

is far away from 0. If it really is independent eqn. 47 will have a

t-distribution with n − p degrees of freedom. Using this distribution we calculate the

probability of of b

≥ |t| which is called as the p-value. This is essentially asking "what

are the chances of obtaining an even more extreme value of b

when it is assumed to be

0?”. A small p-value indicates a relation between response and predictor variables.

Using the concept of SE and p-values we can calculate 95% conﬁdence interval for b

values. which is given as:

(1 − α)% confidence interval for b

= b

± SE

× t

,n−p

(48)

if α = 0.5 it tells us the interval in which we have 95% probability of ﬁnding the true

value of b

, it simply gives us an idea of how far away b

is from β

4.2 Regressing the pointing data

The form of multile regression discussed earlier (eqn. 29) can not be used directly to

ﬁt a model to pointing data. It would be ideal if we wanted to minimize the sum of

squared error in either of eqn. 18 or 19 only i.e. (f −

or (g − ˆg)

but the coordinates

represent a vector and we want to minimze the error distance between two points (true

and estimated) in the sky. The distance between two nearby points in the sky is given

as:

∆D =

(∆A cos E)

+ (∆E)

(49)

The above formula can be obtained fairly simply from spherical trigonometry. It is (∆D)

that we want to minimize. Further as shown in table 2, there are some coeﬃcients that

are common for both ∆A and ∆E. Running two diﬀerent models for them does not allow

for a common coeﬃcients.

The pointing data has the true star coordinates and the coordinates shown by the

encoder. For constructing a regression model we will calculate the diﬀerence between

those two values, (∆A, ∆E). This represents the y. Using the formulations developed in

eqn. 18 & 19, we will estimate those true errors as (

∆A,

∆E) using regression method

discussed below. This represents our ˆy. It is then trivial to see that the pointing errors

as in eqn. 49 is the error in model error estimates - ∆(

∆A). For brevity we shall adopt

to using ∆f and ∆g instead (as shown in eqn. 29).

In order to get across this diﬃculty, we will develop a model that looks like the

following

[6]

(1)

(0)

(−sec E)

(0)

. . . (cos A tan E)

n×p

(0)

(1)

(0)

(−cot E)

. . . (sin A)

n×p

(50)

i represents the number of observation, going from 1 to n

p×1



IA IE CA T F NP AE AN AW



(51)

f = X

β + η

g = X

β + ξ

(52)

It is worthwhile to note that each of X

β and X

β ﬁt our usual notion of a model (with

η and ξ representing the familiar ). The term in β which is not present in a particular X

has the column corresponding to it, all zeroes. E.g. column 3 in X

is 0 since CA does

not aﬀect ∆E and column 4 corresponding to X

is 0 since TF does not aﬀect ∆A (please

refer table 2). It is also noteworthy that it is not necessary for the corrections in altitude

and azimuth to share the same coeﬃcient (e.g. they can be diﬀerent AN

and AN

alt

)

and depends on the type of telescope

[7]

. While constructing two separate models, it was

seen that the two diiﬀerent coeﬃcients have nearly the same value (varying only by few

milli arcseconds). Hence it is a same assumption to assume same coeﬃcients for both the

correction in single model approach.

Since we do not want to run two separate models, in order to capture the vector-like

notion and minimize the squared error as shown in eq. 49, we formulate the model as:





2n×1





2n×p

p×1





2n×1

(53)

Both the models have been stacked on top on each other, giving us the familiar form of:

y = Xβ + 

The only diﬀerence being that the dimensions have changed from n to 2n. Replacing that

in equations mentioned in section 4.1 will give us the correct estimates (for MSE etc.).

One can easily verify that sum of squared error is as expected by eqn. 49:

RSS = e

· e = (y − ˆy)

· (y − ˆy)



(f −

f) (g − ˆg)





(f −

(g − ˆg)



= (∆f)

+ (∆g)

(54)

The only thing that seems oﬀ is the missing cos E term. That can be accounted for if eqn.

2 is multiplied by cos E on both sides (replace f by f cos E). That will make absolutely

no diﬀerence in the b values.

Replacing X in eqn. 23 using eqn. 53, we get the estimated value of β:

b = (X

+ X

)

−1

∆f + X

∆g) (55)

The dimensions of entities in equation 20. has changed, we have done nothing new

here. All the other estimates have been calculated using the methods described in previous

section. E.g. MSE is now given as:

MSE =

RSS

2n − p

(56)

and residual pointing error is given as:

∆P =

(∆f cos E)

+ (∆E)

2n − p

MSE (57)

5 Obtaining statistics

All calculations were done using ipython and Jupyter notebooks, utilising standard plot-

ting, machine learning and numerical computation libraries - scikit-learn, numpy, scipy,

matplotlib, pandas. Such models are usually ﬁt using TPOINT but manual statistical

analysis was done because TPOINT:

• does not give any information about the distribution of estimated parameters (con-

ﬁdence intervals etc).

• does not do any kind of cross-validtion to estimate MSE for the population (all-sky).

• may (and does) miss some signiﬁcant terms.

• gives no idea about underﬁtting/overﬁtting.

• does not check for normality of residuals (and other assumptions).

Initially two separate models were also ﬁt since the correlation between f and g was

nearly negligible (< 0.02). Based on the assumption that minimizing ∆f and ∆g would

minimize their squared sum. Although this approach produced comparable results, it is

not technically accurate since Corr(X, Y ) only measures the linear dependence between

two vectors. If two sets of values are independent then they have Corr(X, Y ) = 0 but

vice versa is usually not true. Hence we shall stick to the model developed above.

On ﬁtting only the 7 geometrical terms statistics were obtained as shown in ﬁg.

8. RMSE is about 27 arcsecond and sky pointing accuracy is about 36 arcseconds.