IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 8300


Ignore:
Timestamp:
Aug 11, 2006, 4:56:32 PM (20 years ago)
Author:
eugene
Message:

moved gaussj to libohana, fixed polyterm to coords for higher order, added test program

Location:
trunk/Ohana/src/gastro2
Files:
1 deleted
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/gastro2/Makefile

    r8241 r8300  
    2525COORDTEST = \
    2626$(SRC)/coordtest.$(ARCH).o \
    27 $(SRC)/gaussj.$(ARCH).o \
    2827$(SRC)/gpairs.$(ARCH).o \
    2928$(SRC)/polyfit.$(ARCH).o
    … …  
    4443$(SRC)/lumfunc.$(ARCH).o \
    4544$(SRC)/gregions2.$(ARCH).o \
    46 $(SRC)/gaussj.$(ARCH).o \
    4745$(SRC)/sort.$(ARCH).o \
    4846$(SRC)/misc.$(ARCH).o \

  • trunk/Ohana/src/gastro2/src/coordtest.c

    r8241 r8300  
    22
    33# define A00 +500.00
    4 // # define A00 +0.00
    54# define A10 +4.00
    65# define A01 +1.00
    … …  
    87# define A11 +0.000
    98# define A02 -0.001
    10 # define A30 +0.002
    11 # define A21 +0.001
    12 # define A12 -0.001
    13 # define A03 -0.002
     9# define A30 +0.00002
     10# define A21 +0.00001
     11# define A12 -0.00001
     12# define A03 -0.00002
    1413
    1514# define B00 400.00
    16 // # define B00 +0.00
    1715# define B10  -1.00
    1816# define B01  4.00
    … …  
    2018# define B11  0.000
    2119# define B02 -0.001
    22 # define B30 -0.002
    23 # define B21 -0.001
    24 # define B12 +0.001
    25 # define B03 +0.002
     20# define B30 -0.00002
     21# define B21 -0.00001
     22# define B12 +0.00001
     23# define B03 +0.00002
    2624
    2725int main (int argc, char **argv) {
    … …  
    109107    double Lo, Mo, dL, dL2, dM, dM2;
    110108    double Xo, Yo, dX, dX2, dY, dY2;
     109    double Lx, Mx;
    111110
    112111    dL = dL2 = dM = dM2 = 0.0;

  • trunk/Ohana/src/gastro2/src/gheader2.c

    r7995 r8300  
    1515  oldsize = header.size;
    1616
    17   if (Target[0].answer.N < 6) {
    18     gfits_modify (&header, "NASTRO", "%d", 1, 0);
    19     gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry");
    20     goto skipstuff;
     17  /* check for insufficient number of stars */
     18  switch (Target[0].coords.Npolyterms) {
     19    case 0:
     20    case 1:
     21      if (Target[0].answer.N < 6) {
     22        gfits_modify (&header, "NASTRO", "%d", 1, 0);
     23        gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry");
     24        goto skipstuff;
     25      }
     26      break;
     27    case 2:
     28      if (Target[0].answer.N < 12) {
     29        gfits_modify (&header, "NASTRO", "%d", 1, 0);
     30        gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry");
     31        goto skipstuff;
     32      }
     33      break;
     34    case 3:
     35      if (Target[0].answer.N < 20) {
     36        gfits_modify (&header, "NASTRO", "%d", 1, 0);
     37        gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry");
     38        goto skipstuff;
     39      }
     40      break;
     41    default:
     42      fprintf (stderr, "invalid order\n");
     43      exit (2);
    2144  }
    2245 

  • trunk/Ohana/src/gastro2/src/gproject2.c

    r8241 r8300  
    3434
    3535  /* create Focal Plane to Tangent Plane transformation from input coords */
    36   strcpy (FPtoTP.ctype, "FP---LIN");
     36  strcpy (FPtoTP.ctype, "FP---PLY");
    3737  FPtoTP.crval1 = FPtoTP.crval2 = 0;
    3838

  • trunk/Ohana/src/gastro2/src/lumfunc.c

    r2442 r8300  
    9797    b[1][0] += y[i]*x[i];
    9898  }
    99   gaussj (c, 2, b, 1);
     99  dgaussj (c, 2, b, 1);
    100100  *C0 = b[0][0];
    101101  *C1 = b[1][0];

  • trunk/Ohana/src/gastro2/src/plots.c

    r8249 r8300  
    375375  fprintf (stderr, "residuals using RD_to_XY\n");
    376376  plot_resid_plot (0, xvect0, yvect0, Npair);
    377   plot_resid_plot (1, xvect1, yvect1, Npair);
     377  // plot_resid_plot (1, xvect1, yvect1, Npair);
    378378  plot_fullfield_pairs (xvect2, yvect2, 2*Npair);
    379379
    … …  
    381381    fit_apply (&x, &y, st[idx1[i]].X, st[idx1[i]].Y);
    382382   
    383     dx = x - sr[idx2[i]].L;
    384     dy = y - sr[idx2[i]].M;
     383    dx = 5.0*(x - sr[idx2[i]].P);
     384    dy = 5.0*(y - sr[idx2[i]].Q);
    385385   
    386386    xvect0[i] = st[idx1[i]].X;
    … …  
    389389    yvect1[i] = dy;
    390390
    391     xvect2[2*i+0] = sr[idx2[i]].L;
    392     yvect2[2*i+0] = sr[idx2[i]].M;
     391    xvect2[2*i+0] = sr[idx2[i]].P;
     392    yvect2[2*i+0] = sr[idx2[i]].Q;
    393393    xvect2[2*i+1] = x;
    394394    yvect2[2*i+1] = y;
    395395  }
    396396
    397   fprintf (stderr, "residuals using RD_to_XY\n");
    398   plot_resid_plot (0, xvect0, yvect0, Npair);
    399   plot_resid_plot (1, xvect1, yvect1, Npair);
    400   plot_fullfield_pairs (xvect2, yvect2, 2*Npair);
     397  fprintf (stderr, "residuals using fit_apply\n");
     398  plot_resid_plot (1, xvect0, yvect0, Npair);
     399  // plot_resid_plot (1, xvect1, yvect1, Npair);
     400  // plot_fullfield_pairs (xvect2, yvect2, 2*Npair);
    401401
    402402  free (xvect2);

  • trunk/Ohana/src/gastro2/src/polyfit.c

    r8241 r8300  
    101101    vector[i][1] /= max;
    102102  }
    103   gaussj (matrix, NPARS, vector, 2);
     103  dgaussj (matrix, NPARS, vector, 2);
    104104  i = 0;
    105105  for (m = 0; m < NPOWR; m++) {
    … …  
    233233
    234234/* convert fit terms to coords and polyterms */
     235/**** what do we do with the value of ctype???
     236      can we leave it alone? ****/
    235237int fit_adjust (Coords *coords) {
    236238
    … …  
    240242  double Xo, Yo, det;
    241243  int Np, Nv;
    242   double **A, **B;
     244  double **A, **B, Fx, Fy;
    243245   
    244246  /* start with the linear solution for Xo,Yo */
    245247  coords[0].cdelt1 = hypot (xsum[1][0], ysum[1][0]);
    246248  coords[0].cdelt2 = hypot (xsum[0][1], ysum[0][1]);
    247 
    248   coords[0].cdelt1 = coords[0].cdelt2 = 1.0;
     249  // coords[0].cdelt1 = coords[0].cdelt2 = 1.0;
     250
    249251  det = 1.0 / (xsum[1][0]*ysum[0][1] - xsum[0][1]*ysum[1][0]);
    250252  Xo = det*(ysum[0][0]*xsum[0][1] - xsum[0][0]*ysum[0][1]);
    … …  
    265267
    266268    for (i = 0; i < 10; i++) {
    267       fit_apply (&B[0][0], &B[1][0], Xo, Yo);
     269      fit_apply (&Fx, &Fy, Xo, Yo);
    268270      fit_apply_dx (&A[0][0], &A[0][1], Xo, Yo);
    269271      fit_apply_dy (&A[1][0], &A[1][1], Xo, Yo);
    270       B[0][0] *= -1;
    271       B[1][0] *= -1;
    272       gaussj (A, 2, B, 1);
     272      B[0][0] = -Fx;
     273      B[1][0] = -Fy;
     274      dgaussj (A, 2, B, 1);
    273275      Xo += B[0][0];
    274276      Yo += B[1][0];
    … …  
    295297
    296298    case 2:
    297      
    298299      a10 = xsum[1][0] + 2.0*xsum[2][0]*Xo + xsum[1][1]*Yo;
    299300      a01 = xsum[0][1] + 2.0*xsum[0][2]*Yo + xsum[1][1]*Xo;
    … …  
    326327     
    327328    case 3:
    328       /* the linear solution can be analytically inverted */
    329       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;
    330       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;
    331       a20 = xsum[2][0] + 3*xsum[3][0]*Xo + xsum[2][1]*Yo;
    332       a11 = xsum[1][1] +   xsum[2][1]*Xo + xsum[1][2]*Yo;
    333       a02 = xsum[0][2] + 3*xsum[0][3]*Yo + xsum[1][2]*Xo;
     329      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;
     330      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;
     331      a20 = xsum[2][0] + 3*xsum[3][0]*Xo +   xsum[2][1]*Yo;
     332      a11 = xsum[1][1] + 2*xsum[2][1]*Xo + 2*xsum[1][2]*Yo;
     333      a02 = xsum[0][2] + 3*xsum[0][3]*Yo +   xsum[1][2]*Xo;
    334334      a30 = xsum[3][0];
    335335      a21 = xsum[2][1];
    … …  
    337337      a03 = xsum[0][3];
    338338
    339       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;
    340       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;
    341       b20 = ysum[2][0] + 3*ysum[3][0]*Xo + ysum[2][1]*Yo;
    342       b11 = ysum[1][1] +   ysum[2][1]*Xo + ysum[1][2]*Yo;
    343       b02 = ysum[0][2] + 3*ysum[0][3]*Yo + ysum[1][2]*Xo;
     339      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;
     340      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;
     341      b20 = ysum[2][0] + 3*ysum[3][0]*Xo +   ysum[2][1]*Yo;
     342      b11 = ysum[1][1] + 2*ysum[2][1]*Xo + 2*ysum[1][2]*Yo;
     343      b02 = ysum[0][2] + 3*ysum[0][3]*Yo +   ysum[1][2]*Xo;
    344344      b30 = ysum[3][0];
    345345      b21 = ysum[2][1];
Note: See TracChangeset for help on using the changeset viewer.