stardat = {{ 5191,   0.383433,   0.868904,118, 16 },
{  5054,   0.456958,   0.764102,112, 28 },
{  4905,   0.555319,   0.735841,104, 34 },
{  6536,  -0.430286,   0.812643,192, 40 },
{  6536,  -0.430286,   0.812643,192, 42 },
{  6705,  -0.529763,   0.817511,202, 50 },
{  4301,   0.918772,   0.580957, 86, 64 },
{  5563,   0.112046,   0.448420,140, 64 },
{  5563,   0.112046,   0.448420,140, 66 },
{  7796,  -1.138965,   0.930789,242,100 },
{   424,  -0.020781,   0.161240,138,102 },
{  7924,  -1.192617,   0.831810,232,108 },
{  7924,  -1.192617,   0.831810,234,108 },
{    21,  -1.980011,   0.446945,176,156 },
{  2088,   2.365176,   0.654424, 62,164 },
{  1708,   2.583464,   0.615932, 76,174 },
{  1017,  -3.100536,   0.527057,118,186 },
{   464,  -2.521721,   0.574505,158,188 },
{  1577,   2.717648,   0.831774, 68,202 },
{   603,  -2.677937,   0.673679,150,202 },
{   936,  -3.018307,   0.683470,122,206 },
{   921,  -3.004177,   0.720677,122,210 },
{   337,  -2.430564,   0.812391,176,210 }}

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
