stardat = { { 5054,    0.174857,    0.782538, 110,  34},
{ 6705,   -0.824635,    0.780336, 198,  64},
{ 4301,    0.640339,    0.622669,  80,  66},
{ 4301,    0.640339,    0.622669,  80,  68},
{ 7001,   -1.037074,    0.973810, 230,  70},
{ 4434,    0.491103,    0.505226,  98,  74},
{ 5563,   -0.104488,    0.448645, 132,  74},
{ 5563,   -0.104488,    0.448645, 134,  74},
{  424,    0.006584,    0.160812, 126, 110},
{  424,    0.006584,    0.160812, 128, 110},
{ 7796,   -1.466421,    0.873605, 234, 118},
{ 2088,    1.983570,    0.703262,  48, 164},
{ 2088,    1.983570,    0.703262,  46, 166},
{ 1708,    2.185239,    0.656690,  60, 176},
{  168,   -2.597055,    0.429010, 156, 176},
{ 1708,    2.185239,    0.656690,  58, 178},
{ 1708,    2.185239,    0.656690,  60, 178},
{ 2286,    1.980009,    1.102883,   4, 182},
{ 2286,    1.980009,    1.102883,   6, 182},
{ 2216,    2.016280,    1.097424,   8, 186},
{ 1017,    2.743637,    0.537599, 100, 192},
{ 8775,   -2.233389,    0.968167, 220, 202},
{ 1577,    2.335125,    0.866012,  50, 204},
{ 8775,   -2.233389,    0.968167, 220, 204},
{  603,   -3.090811,    0.658580, 130, 212},
{  921,    2.871999,    0.724596, 102, 216},
{  921,    2.871999,    0.724596, 102, 218},
{  337,   -2.820777,    0.783121, 156, 222},
{  337,   -2.820777,    0.783121, 158, 222},
{  617,   -3.114902,    0.987747, 130, 250} }

nstars = Length[ stardat]
ix = Table[ i, {i, nstars}]
azimuths = Part[ stardat, ix, 2]
zenangs = Part[ stardat, ix, 3]
cols = Part[ stardat, ix, 4]
rows = Part[ stardat, ix, 5]

x0 = Cos[ -azimuths] Sin[ zenangs]
y0 = Sin[ -azimuths] Sin[ zenangs]

xp[ca_] := x0 Cos[ -ca] + y0 Sin[ -ca]
yp[ca_] := -x0 Sin[ -ca] + y0 Cos[ -ca]

rho[r_,k_] := r Sin[ k zenangs]
phi[ca_] := +ArcTan[ xp[ca], yp[ca]]

col[c0_,r_,k_,ca_] := c0 + rho[ r, k] Sin[ phi[ ca]]
row[r0_,r_,k_,ca_] := r0 - rho[ r, k] Cos[ phi[ ca]]

colerr[c0_,r_,k_,ca_] := cols - col[ c0, r, k, ca]
rowerr[r0_,r_,k_,ca_] := rows - row[ r0, r, k, ca]

sqerror[c0_,r0_,r_,k_,ca_] := colerr[ c0, r, k, ca]^2 +
			      rowerr[ r0, r, k, ca]^2

totsqerror[c0_,r0_,r_,k_,ca_] := Apply[ Plus,
			                sqerror[ c0, r0, r, k, ca]]

fit = FindMinimum[ totsqerror[c0,r0,r,k,ca], {c0,130.936},{r0,127.054},
                   {r,199.853},{k,0.642907},{ca,0.230475},
                   AccuracyGoal -> 6]

parms = Part[ fit, 2]

fc0 = Part[ parms, 1, 2]
fr0 = Part[ parms, 2, 2]
fr = Part[ parms, 3, 2]
fk = Part[ parms, 4, 2]
fca = Part[ parms, 5, 2]

error = Sqrt[ sqerror[fc0,fr0,fr,fk,fca]]
out=Transpose[ Partition[ Join[ Flatten[ Transpose[ stardat]], error], nstars]]

mse = "mean squared error = " <> ToString[ Part[ fit, 1] / nstars] <>
      " pixel^2"

OutputForm[ "Star  Azimuth    Zenith    Col  Row  PixError"] >> results
OutputForm[ TableForm[ out, TableSpacing -> {0, 2}]] >>> results
OutputForm[ mse] >>> results
ColumnForm[ parms] >>> results

!!results
