stardat = { { 4295,   -0.013908,    0.760328, 126,  28},
{ 4905,   -0.416017,    0.749545, 160,  36},
{ 4905,   -0.416017,    0.749545, 160,  38},
{ 5191,   -0.632489,    0.844499, 184,  38},
{ 4301,   -0.020042,    0.666604, 126,  40},
{ 5054,   -0.521660,    0.757735, 170,  40},
{ 2088,    1.142453,    0.846394,  26,  80},
{ 2088,    1.142453,    0.846394,  30,  80},
{ 5563,   -0.648658,    0.395340, 154,  82},
{ 6132,   -1.133459,    0.547815, 186,  92},
{ 1708,    1.314294,    0.797563,  28,  98},
{ 1708,    1.314294,    0.797563,  30,  98},
{ 1708,    1.314294,    0.797563,  28, 100},
{ 1708,    1.314294,    0.797563,  30, 100},
{  424,    0.062496,    0.166308, 124, 102},
{ 6536,   -1.489857,    0.652271, 206, 114},
{ 6705,   -1.609464,    0.645980, 206, 124},
{ 8974,    2.343600,    0.055317, 120, 128},
{ 1017,    1.766825,    0.648825,  46, 140},
{  915,    1.828990,    0.572670,  56, 142},
{ 7001,   -1.865862,    0.831684, 222, 150},
{ 7001,   -1.865862,    0.831684, 224, 150},
{  264,    2.399878,    0.370095,  94, 158},
{   21,    2.695747,    0.376828, 104, 166},
{  168,    2.531282,    0.433783,  94, 168},
{  603,    2.193356,    0.720953,  52, 176},
{ 7924,   -2.428339,    0.642118, 178, 182},
{ 7417,   -2.156045,    0.980633, 224, 186},
{  337,    2.481476,    0.807386,  64, 202},
{ 7949,   -2.495335,    0.835120, 188, 202},
{ 7949,   -2.495335,    0.835120, 188, 204},
{   15,    2.793784,    0.898701,  88, 226},
{ 8775,    3.112266,    0.907118, 122, 234} }

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
