stardat = { { 5191,    0.107685,    0.882538, 112,  16},
{ 5054,    0.190350,    0.781941, 106,  28},
{ 4905,    0.289412,    0.758834,  96,  32},
{ 4905,    0.289412,    0.758834,  96,  34},
{ 6132,   -0.424957,    0.650896, 158,  48},
{ 6705,   -0.808157,    0.782732, 194,  54},
{ 4295,    0.694645,    0.711325,  68,  56},
{ 7001,   -1.019478,    0.976625, 224,  58},
{ 4301,    0.655426,    0.620685,  76,  62},
{ 4301,    0.655426,    0.620685,  78,  62},
{ 5563,   -0.092535,    0.448968, 130,  66},
{ 7310,   -0.944044,    0.467445, 172,  88},
{  424,    0.005062,    0.160793, 124, 102},
{ 7924,   -1.501756,    0.777146, 220, 114},
{ 2890,    1.645992,    0.991597,   6, 132},
{   21,   -2.371823,    0.400809, 160, 158},
{ 2088,    2.003912,    0.700264,  46, 160},
{  264,   -2.635125,    0.352589, 148, 162},
{  403,   -2.814345,    0.352220, 140, 164},
{  168,   -2.572222,    0.430748, 154, 168},
{  168,   -2.572222,    0.430748, 156, 168},
{ 1708,    2.206488,    0.654024,  58, 172},
{ 1708,    2.206488,    0.654024,  60, 172},
{  915,    2.869719,    0.468124, 110, 180},
{ 1017,    2.767499,    0.536361, 100, 186},
{ 1017,    2.767499,    0.536361, 102, 186},
{ 8775,   -2.213273,    0.970779, 220, 192},
{ 1791,    2.232795,    0.958974,  34, 196},
{  603,   -3.067825,    0.658784, 132, 204},
{  936,    2.875363,    0.687477, 104, 206},
{  337,   -2.798911,    0.784191, 158, 212},
{   15,   -2.515277,    0.919052, 192, 212},
{  337,   -2.798911,    0.784191, 158, 214},
{  617,   -3.093989,    0.987869, 132, 242} }

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

