Index: /trunk/Ohana/src/gastro2/Makefile =================================================================== --- /trunk/Ohana/src/gastro2/Makefile (revision 8299) +++ /trunk/Ohana/src/gastro2/Makefile (revision 8300) @@ -25,5 +25,4 @@ COORDTEST = \ $(SRC)/coordtest.$(ARCH).o \ -$(SRC)/gaussj.$(ARCH).o \ $(SRC)/gpairs.$(ARCH).o \ $(SRC)/polyfit.$(ARCH).o @@ -44,5 +43,4 @@ $(SRC)/lumfunc.$(ARCH).o \ $(SRC)/gregions2.$(ARCH).o \ -$(SRC)/gaussj.$(ARCH).o \ $(SRC)/sort.$(ARCH).o \ $(SRC)/misc.$(ARCH).o \ Index: /trunk/Ohana/src/gastro2/src/coordtest.c =================================================================== --- /trunk/Ohana/src/gastro2/src/coordtest.c (revision 8299) +++ /trunk/Ohana/src/gastro2/src/coordtest.c (revision 8300) @@ -2,5 +2,4 @@ # define A00 +500.00 -// # define A00 +0.00 # define A10 +4.00 # define A01 +1.00 @@ -8,11 +7,10 @@ # define A11 +0.000 # define A02 -0.001 -# define A30 +0.002 -# define A21 +0.001 -# define A12 -0.001 -# define A03 -0.002 +# define A30 +0.00002 +# define A21 +0.00001 +# define A12 -0.00001 +# define A03 -0.00002 # define B00 400.00 -// # define B00 +0.00 # define B10 -1.00 # define B01 4.00 @@ -20,8 +18,8 @@ # define B11 0.000 # define B02 -0.001 -# define B30 -0.002 -# define B21 -0.001 -# define B12 +0.001 -# define B03 +0.002 +# define B30 -0.00002 +# define B21 -0.00001 +# define B12 +0.00001 +# define B03 +0.00002 int main (int argc, char **argv) { @@ -109,4 +107,5 @@ double Lo, Mo, dL, dL2, dM, dM2; double Xo, Yo, dX, dX2, dY, dY2; + double Lx, Mx; dL = dL2 = dM = dM2 = 0.0; Index: unk/Ohana/src/gastro2/src/gaussj.c =================================================================== --- /trunk/Ohana/src/gastro2/src/gaussj.c (revision 8299) +++ (revision ) @@ -1,77 +1,0 @@ -# include - -int gaussj (double **a, int n, double **b, int m) { - - int *indxc,*indxr,*ipiv; - int i, icol, irow, j, k, l, ll; - double big,dum,pivinv; - - irow = icol = 0; - - ALLOCATE (indxc, int, n); - ALLOCATE (indxr, int, n); - ALLOCATE (ipiv, int, n); - for (j = 0; j < n; j++) - ipiv[j] = 0; - - for (i = 0; i < n; i++) { - big = 0.0; - for (j = 0; j < n; j++) { - if (ipiv[j] != 1) { - for (k = 0; k < n; k++) { - if (ipiv[k] == 0) { - if (fabs (a[j][k]) >= big) { - big = fabs (a[j][k]); - irow = j; - icol = k; - } - } - else - if (ipiv[k] > 1) { - fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n"); - return (0); - } - } - } - } - ipiv[icol]++; - if (irow != icol) { - for (l = 0; l < n; l++) - SWAP (a[irow][l], a[icol][l]); - for (l = 0; l < m; l++) - SWAP (b[irow][l], b[icol][l]); - } - indxr[i] = irow; - indxc[i] = icol; - if (a[icol][icol] == 0.0) { - fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n"); - return (0); - } - pivinv = 1.0 / a[icol][icol]; - a[icol][icol] = 1.0; - for (l = 0; l < n; l++) - a[icol][l] *= pivinv; - for (l = 0; l < m; l++) - b[icol][l] *= pivinv; - for (ll = 0; ll < n; ll++) { - if (ll != icol) { - dum = a[ll][icol]; - a[ll][icol] = 0.0; - for (l = 0; l < n; l++) - a[ll][l] -= a[icol][l]*dum; - for (l = 0; l < m; l++) - b[ll][l] -= b[icol][l]*dum; - } - } - } - - for (l = n - 1; l >= 0; l--) { - if (indxr[l] != indxc[l]) - for (k = 0; k < n; k++) - SWAP (a[k][indxr[l]], a[k][indxc[l]]); - } - free (ipiv); - free (indxr); - free (indxc); - return (1); -} Index: /trunk/Ohana/src/gastro2/src/gheader2.c =================================================================== --- /trunk/Ohana/src/gastro2/src/gheader2.c (revision 8299) +++ /trunk/Ohana/src/gastro2/src/gheader2.c (revision 8300) @@ -15,8 +15,31 @@ oldsize = header.size; - if (Target[0].answer.N < 6) { - gfits_modify (&header, "NASTRO", "%d", 1, 0); - gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); - goto skipstuff; + /* check for insufficient number of stars */ + switch (Target[0].coords.Npolyterms) { + case 0: + case 1: + if (Target[0].answer.N < 6) { + gfits_modify (&header, "NASTRO", "%d", 1, 0); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + goto skipstuff; + } + break; + case 2: + if (Target[0].answer.N < 12) { + gfits_modify (&header, "NASTRO", "%d", 1, 0); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + goto skipstuff; + } + break; + case 3: + if (Target[0].answer.N < 20) { + gfits_modify (&header, "NASTRO", "%d", 1, 0); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + goto skipstuff; + } + break; + default: + fprintf (stderr, "invalid order\n"); + exit (2); } Index: /trunk/Ohana/src/gastro2/src/gproject2.c =================================================================== --- /trunk/Ohana/src/gastro2/src/gproject2.c (revision 8299) +++ /trunk/Ohana/src/gastro2/src/gproject2.c (revision 8300) @@ -34,5 +34,5 @@ /* create Focal Plane to Tangent Plane transformation from input coords */ - strcpy (FPtoTP.ctype, "FP---LIN"); + strcpy (FPtoTP.ctype, "FP---PLY"); FPtoTP.crval1 = FPtoTP.crval2 = 0; Index: /trunk/Ohana/src/gastro2/src/lumfunc.c =================================================================== --- /trunk/Ohana/src/gastro2/src/lumfunc.c (revision 8299) +++ /trunk/Ohana/src/gastro2/src/lumfunc.c (revision 8300) @@ -97,5 +97,5 @@ b[1][0] += y[i]*x[i]; } - gaussj (c, 2, b, 1); + dgaussj (c, 2, b, 1); *C0 = b[0][0]; *C1 = b[1][0]; Index: /trunk/Ohana/src/gastro2/src/plots.c =================================================================== --- /trunk/Ohana/src/gastro2/src/plots.c (revision 8299) +++ /trunk/Ohana/src/gastro2/src/plots.c (revision 8300) @@ -375,5 +375,5 @@ fprintf (stderr, "residuals using RD_to_XY\n"); plot_resid_plot (0, xvect0, yvect0, Npair); - plot_resid_plot (1, xvect1, yvect1, Npair); + // plot_resid_plot (1, xvect1, yvect1, Npair); plot_fullfield_pairs (xvect2, yvect2, 2*Npair); @@ -381,6 +381,6 @@ fit_apply (&x, &y, st[idx1[i]].X, st[idx1[i]].Y); - dx = x - sr[idx2[i]].L; - dy = y - sr[idx2[i]].M; + dx = 5.0*(x - sr[idx2[i]].P); + dy = 5.0*(y - sr[idx2[i]].Q); xvect0[i] = st[idx1[i]].X; @@ -389,14 +389,14 @@ yvect1[i] = dy; - xvect2[2*i+0] = sr[idx2[i]].L; - yvect2[2*i+0] = sr[idx2[i]].M; + xvect2[2*i+0] = sr[idx2[i]].P; + yvect2[2*i+0] = sr[idx2[i]].Q; xvect2[2*i+1] = x; yvect2[2*i+1] = y; } - fprintf (stderr, "residuals using RD_to_XY\n"); - plot_resid_plot (0, xvect0, yvect0, Npair); - plot_resid_plot (1, xvect1, yvect1, Npair); - plot_fullfield_pairs (xvect2, yvect2, 2*Npair); + fprintf (stderr, "residuals using fit_apply\n"); + plot_resid_plot (1, xvect0, yvect0, Npair); + // plot_resid_plot (1, xvect1, yvect1, Npair); + // plot_fullfield_pairs (xvect2, yvect2, 2*Npair); free (xvect2); Index: /trunk/Ohana/src/gastro2/src/polyfit.c =================================================================== --- /trunk/Ohana/src/gastro2/src/polyfit.c (revision 8299) +++ /trunk/Ohana/src/gastro2/src/polyfit.c (revision 8300) @@ -101,5 +101,5 @@ vector[i][1] /= max; } - gaussj (matrix, NPARS, vector, 2); + dgaussj (matrix, NPARS, vector, 2); i = 0; for (m = 0; m < NPOWR; m++) { @@ -233,4 +233,6 @@ /* convert fit terms to coords and polyterms */ +/**** what do we do with the value of ctype??? + can we leave it alone? ****/ int fit_adjust (Coords *coords) { @@ -240,11 +242,11 @@ double Xo, Yo, det; int Np, Nv; - double **A, **B; + double **A, **B, Fx, Fy; /* start with the linear solution for Xo,Yo */ coords[0].cdelt1 = hypot (xsum[1][0], ysum[1][0]); coords[0].cdelt2 = hypot (xsum[0][1], ysum[0][1]); - - coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + // coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + det = 1.0 / (xsum[1][0]*ysum[0][1] - xsum[0][1]*ysum[1][0]); Xo = det*(ysum[0][0]*xsum[0][1] - xsum[0][0]*ysum[0][1]); @@ -265,10 +267,10 @@ for (i = 0; i < 10; i++) { - fit_apply (&B[0][0], &B[1][0], Xo, Yo); + fit_apply (&Fx, &Fy, Xo, Yo); fit_apply_dx (&A[0][0], &A[0][1], Xo, Yo); fit_apply_dy (&A[1][0], &A[1][1], Xo, Yo); - B[0][0] *= -1; - B[1][0] *= -1; - gaussj (A, 2, B, 1); + B[0][0] = -Fx; + B[1][0] = -Fy; + dgaussj (A, 2, B, 1); Xo += B[0][0]; Yo += B[1][0]; @@ -295,5 +297,4 @@ case 2: - a10 = xsum[1][0] + 2.0*xsum[2][0]*Xo + xsum[1][1]*Yo; a01 = xsum[0][1] + 2.0*xsum[0][2]*Yo + xsum[1][1]*Xo; @@ -326,10 +327,9 @@ case 3: - /* the linear solution can be analytically inverted */ - a10 = xsum[1][0] + 2*xsum[2][0]*Xo + xsum[1][1]*Yo + 3*xsum[3][0]*Xo*Xo + 2*xsum[2][1]*Xo*Yo + xsum[1][2]*Yo*Yo; - a01 = xsum[0][1] + 2*xsum[0][2]*Yo + xsum[1][1]*Xo + 3*xsum[0][3]*Yo*Yo + 2*xsum[1][2]*Xo*Yo + xsum[2][1]*Xo*Xo; - a20 = xsum[2][0] + 3*xsum[3][0]*Xo + xsum[2][1]*Yo; - a11 = xsum[1][1] + xsum[2][1]*Xo + xsum[1][2]*Yo; - a02 = xsum[0][2] + 3*xsum[0][3]*Yo + xsum[1][2]*Xo; + a10 = xsum[1][0] + 2*xsum[2][0]*Xo + xsum[1][1]*Yo + 3*xsum[3][0]*Xo*Xo + 2*xsum[2][1]*Xo*Yo + xsum[1][2]*Yo*Yo; + a01 = xsum[0][1] + 2*xsum[0][2]*Yo + xsum[1][1]*Xo + 3*xsum[0][3]*Yo*Yo + 2*xsum[1][2]*Xo*Yo + xsum[2][1]*Xo*Xo; + a20 = xsum[2][0] + 3*xsum[3][0]*Xo + xsum[2][1]*Yo; + a11 = xsum[1][1] + 2*xsum[2][1]*Xo + 2*xsum[1][2]*Yo; + a02 = xsum[0][2] + 3*xsum[0][3]*Yo + xsum[1][2]*Xo; a30 = xsum[3][0]; a21 = xsum[2][1]; @@ -337,9 +337,9 @@ a03 = xsum[0][3]; - b10 = ysum[1][0] + 2*ysum[2][0]*Xo + ysum[1][1]*Yo + 3*ysum[3][0]*Xo*Xo + 2*ysum[2][1]*Xo*Yo + ysum[1][2]*Yo*Yo; - b01 = ysum[0][1] + 2*ysum[0][2]*Yo + ysum[1][1]*Xo + 3*ysum[0][3]*Yo*Yo + 2*ysum[1][2]*Xo*Yo + ysum[2][1]*Xo*Xo; - b20 = ysum[2][0] + 3*ysum[3][0]*Xo + ysum[2][1]*Yo; - b11 = ysum[1][1] + ysum[2][1]*Xo + ysum[1][2]*Yo; - b02 = ysum[0][2] + 3*ysum[0][3]*Yo + ysum[1][2]*Xo; + b10 = ysum[1][0] + 2*ysum[2][0]*Xo + ysum[1][1]*Yo + 3*ysum[3][0]*Xo*Xo + 2*ysum[2][1]*Xo*Yo + ysum[1][2]*Yo*Yo; + b01 = ysum[0][1] + 2*ysum[0][2]*Yo + ysum[1][1]*Xo + 3*ysum[0][3]*Yo*Yo + 2*ysum[1][2]*Xo*Yo + ysum[2][1]*Xo*Xo; + b20 = ysum[2][0] + 3*ysum[3][0]*Xo + ysum[2][1]*Yo; + b11 = ysum[1][1] + 2*ysum[2][1]*Xo + 2*ysum[1][2]*Yo; + b02 = ysum[0][2] + 3*ysum[0][3]*Yo + ysum[1][2]*Xo; b30 = ysum[3][0]; b21 = ysum[2][1];