stardat = {{4, 0.049047, 0.942966, 153, 12}, 
   {7, 0.4515900000000001, 0.952198, 106, 14}, 
   {1, -0.237294, 1.005112, 185, 19}, 
   {2, 0.4190170000000001, 0.6889870000000001, 115, 43}, 
   {7, 0.990026, 0.892706, 56, 47}, 
   {6, 0.957742, 0.7774950000000001, 68, 55}, 
   {2, -0.660084, 0.889447, 215, 58}, {1, -0.023801, 0.551493, 148, 60}, 
   {5, 1.00335, 0.7116410000000001, 70, 64}, 
   {1, -0.7703790000000001, 0.876744, 220, 69}, 
   {3, -0.737946, 0.798642, 211, 71}, 
   {3, 1.201075, 0.6099130000000001, 69, 84}, 
   {1, 1.047274, 0.4433880000000001, 91, 88}, 
   {23, 1.560093, 0.6083420000000001, 58, 108}, 
   {3, -1.205304, 0.904519, 239, 112}, {1, -1.310067, 0.636226, 209, 124}, 
   {2, -1.435982, 0.7776290000000001, 225, 135}, 
   {1, -1.829201, 0.424524, 177, 151}, {2, -2.122725, 0.345518, 161, 157}, 
   {3, 2.177909, 0.8453, 35, 164}, {9, -2.02051, 0.637389, 192, 176}, 
   {7, -1.83407, 0.916905, 228, 179}, {1, 2.307612, 0.949383, 30, 182}, 
   {1, -3.128213, 0.459826, 118, 182}, 
   {2, 3.068306, 0.5280160000000001, 111, 189}, 
   {2, -2.247108, 0.648191, 180, 189}}

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"

TableForm[ out] >> results
OutputForm[ mse] >>> results
ColumnForm[ parms] >>> results

!!results
