#!/usr/bin/env python

from numpy import *

# data jsou cteny ze souboru se strukturou:
# a d x y
# ..
# a,d je rektascence a deklinace ve stupnich
# x,y jsou kartezske souradnice hvezd

file = loadtxt("astrometrie.dat")

a = file[:,0]
d = file[:,1]
xx = file[:,2]
yy = file[:,3]

# pocet hvezd
n = size(a)

# definice stredu snimku
a0 = mean(a)
d0 = mean(d)

# stred CCD snimku, upravit dle skutecneho
# rozmeru snimku
x0 = 765/2
y0 = 510/2

# radian
rad = 180/pi

print("#      X0         Y0       A=cos      B=sin ") 

for j in range(1,10):
  # range (prvni cislo, kolikrat)

  # prepocet na gnomonicke souradnice
  u = - (a - a0)*cos(d0/rad)
  v = d - d0

  # prepocet na stred snimku
  x = xx - x0
  y = yy - y0

  # meritko
  cc = empty(n-2)
  for i in range(1,n-2):
     # i = 1 ... n-2

     # vzdalenosti hvezd
     d1 = sqrt((x[i+1] - x[i])**2 + (y[i+1] - y[i])**2)
     d2 = sqrt((u[i+1] - u[i])**2 + (v[i+1] - v[i])**2)
     cc[i] = d1/d2

     # v pripade hledani spatne identifikace odkomentovat:
#     print(i,cc[i])

  # meritko [pixely na stupen]
  c = mean(cc[1:n-2])
  # stat. chyba meritka 
  dc = std(cc[1:n-2])/sqrt(n-3)

  # prepocet u,v na stejne meritko
  u = c*u
  v = c*v

  # vypocty sum v matici nejmensich ctvercu
  sx = sum(u)
  sy = sum(v)
  sx2 = sum(u**2) + sum(v**2)
  sa = sum(x)
  sb = sum(y)
  sa2 = sum(u*x + v*y)
  sab = sum(v*x - u*y) 

  mnc = array([ [n,0,sx,sy], [0,n,sy,-sx], [sx,sy,sx2,0], [sy,-sx,0,sx2] ])
  b = array([sa,sb,sa2,sab])

  # vypocty parametru
  # t[ X0, Y0, A, B ]
  t = linalg.solve(mnc,b)

  # prubezny vypis
  print("{0} {1:10.5f} {2:10.5f} {3:8.5f} {4:8.5f}".format(j,t[0],t[1],t[2],t[3]))

  # korekce posunu stredu
  a0 = a0 + t[0]/c/cos(d0/rad)
  d0 = d0 - t[1]/c

# fin

print(" rektascenze deklinace    x      y    dx[arcsec]    dy[arcsec]")
s = 0
for i in range(0,n):
  dx = t[0] + t[2]*u[i] + t[3]*v[i] - x[i]
  dy = t[1] - t[3]*u[i] + t[2]*v[i] - y[i]
  s = s + dx**2 + dy**2
  print("{0:10.5f} {1:10.5f}   {2:5.1f} {3:5.1f} {4:10.1f} {5:10.1f}"\
         .format(a[i],d[i],xx[i],yy[i],3600*dx/c,3600*dy/c))

# inversni matice
mnc1 = linalg.inv(mnc)
dt = sqrt(s*diag(mnc1)/(n-4))

print("Meritko: {0:5.3f} +- {1:5.3f} [arcsec/pix]".format(3600/c,3600*dc/(c*c)))

print("Reseni:")
for i in range(0,size(t)):
  print("{0} {1:10.5f} +- {2:10.5}".format(i,t[i],dt[i]))