stardat= {{5054, 0.027166, 0.7856990000000001, 127, 24}, 
  {4905, 0.12717, 0.7659930000000001, 117, 27}, 
  {4301, 0.4975170000000001, 0.639556, 91, 51}, 
  {5563, -0.218755, 0.4435990000000001, 142, 66}, 
  {424, 0.020863, 0.161244, 130, 101}, 
  {2990, 1.416618, 1.100629, 2, 103}, {1708, 1.986804, 0.683963, 52, 157}, 
  {264, -2.922509, 0.340293, 141, 164}, {21, -2.637136, 0.380205, 155, 164}, 
  {403, -3.105621, 0.345961, 133, 166}, {168, -2.839571, 0.41613, 147, 172}, 
  {915, 2.610035, 0.481665, 100, 175}, 
  {936, 2.639086, 0.7004590000000001, 89, 199}, 
  {7001, -0.565796, 1.037398, 197, 14}, {6705, -0.384281, 0.830671, 168, 25}, 
  {5054, 0.5999240000000001, 0.749066, 77, 45}, 
  {7924, -1.034052, 0.858563, 221, 66}, {5563, 0.220768, 0.443478, 118, 67}, 
  {4301, 1.06345, 0.5560070000000001, 69, 88}, 
  {424, -0.03359299999999999, 0.162054, 131, 101}, 
  {21, -1.791178, 0.4753000000000001, 190, 134}, 
  {264, -1.999967, 0.415363, 179, 143}, {168, -1.970509, 0.496065, 189, 145}, 
  {403, -2.152744, 0.404589, 174, 149}, {337, -2.241165, 0.833902, 213, 185}, 
  {1017, -2.879888, 0.531555, 149, 186}, {603, -2.476633, 0.68963, 185, 189}, 
  {2088, 2.566344, 0.6357680000000001, 88, 190}, 
  {7924, -0.419635, 0.936802, 176, 16}, {6705, 0.194403, 0.841846, 109, 19}, 
  {5563, 0.6653599999999999, 0.392654, 99, 83}, 
  {21, -1.128647, 0.5921760000000001, 198, 89}, 
  {168, -1.275975, 0.613495, 205, 98}, {424, -0.069298, 0.168523, 132, 100}, 
  {4905, 1.302552, 0.619254, 55, 102}, 
  {264, -1.297521, 0.5326010000000001, 196, 102}, 
  {337, -1.527184, 0.94579, 247, 115}, {4301, 1.708221, 0.440395, 75, 130}, 
  {915, -1.925049, 0.557848, 197, 145}, 
  {1017, -2.051418, 0.6053650000000001, 199, 156}, 
  {936, -2.028977, 0.76791, 217, 163}, {1708, -2.637433, 0.611985, 169, 189}, 
  {2088, -2.871424, 0.617822, 153, 196}, {424, -0.070329, 0.177381, 132, 99}, 
  {5563, 1.1561, 0.297929, 96, 107}, 
  {1708, -1.84687, 0.7054960000000001, 216, 144}}

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"

Format[LineBreak[_]] = "\0"
OutputForm[ TableForm[ out, TableSpacing -> {0, 2}]] >> results
OutputForm[ mse] >>> results
ColumnForm[ parms] >>> results

!!results
