IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 8241


Ignore:
Timestamp:
Aug 8, 2006, 3:50:18 PM (20 years ago)
Author:
eugene
Message:

substantial work to correctly handle higher order polynomial.
conversions between the internal fitting model (R = sum (Aij xi yj))
and the header representation were incorrect: they did not correctly
handle the shift to the CRPIXi value (ie, x -> (x - xo)). i fixed this
by making the fit explicitly work on the fit to the tangent-plane
coordinates, instead of fitting corrections, and calculating the correct
change of polynomial coefficients. this is not as simple as it seems at
first glance: the header version of the coefficients does not have a
zero-order term, and instead has the CRPIX value. to convert from one to
the other requires solving the equation L(Xo,Yo),M(Xo,Yo) = (0,0). I do
this using the Newton-Raphson method.

Location:
trunk/Ohana/src/gastro2
Files:
3 added
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/gastro2/Makefile

    r6426 r8241  
    2323$(SRC)/gregions2.$(ARCH).o
    2424
     25COORDTEST = \
     26$(SRC)/coordtest.$(ARCH).o \
     27$(SRC)/gaussj.$(ARCH).o \
     28$(SRC)/gpairs.$(ARCH).o \
     29$(SRC)/polyfit.$(ARCH).o
     30
    2531GASTRO = \
    2632$(SRC)/getptolemy.$(ARCH).o \
    … …  
    2935$(SRC)/gheader2.$(ARCH).o \
    3036$(SRC)/gfit2.$(ARCH).o \
     37$(SRC)/gpairs.$(ARCH).o \
     38$(SRC)/polyfit.$(ARCH).o \
    3139$(SRC)/plotstuff.$(ARCH).o \
    3240$(SRC)/rotate2.$(ARCH).o \
    … …  
    6270        $(CC) $(2MASS) -o $(BIN)/extr2mass.linux $(CCFLAGS)
    6371
     72coordtest: $(COORDTEST)
     73        $(CC) $(COORDTEST) -o $(BIN)/coordtest.linux $(CCFLAGS)
     74
    6475# dependancy rules for binary code ##########################
    6576$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH)

  • trunk/Ohana/src/gastro2/include/gastro2.h

    r7997 r8241  
    33
    44typedef struct {
    5   double R;
    6   double D;
    7   double X;
    8   double Y;
    9   double M;
    10   double dM;
     5  double R, D;
     6  double P, Q;
     7  double X, Y;
     8  double M, dM;
    119  int type;
    1210} StarData;
    13 
    14 # if (0)
    15 typedef struct {
    16   float X;
    17   float Y;
    18   float M;
    19   float dM;
    20   float Mgal;
    21   float Map;
    22   float sky;
    23   float fx;
    24   float fy;
    25   float df;
    26   char  dophot;
    27   char  dummy[3];
    28 } Stars;
    29 # endif
    3011
    3112typedef struct {
    … …  
    129110double MAGLIM_MIN;
    130111double MAGLIM_MAX;
    131 
     112int NGRID_PIX;
     113
     114int    ASCA;
    132115int    FORCE;
    133116double F_RA;
    … …  
    173156void      fit_lum_bin (double *x, double *y, int N, double *C0, double *C1);
    174157void      fit_norm ();
    175 double    fit_scat (StarData *st, StarData *sr);
     158double    fit_scat (StarData *st, StarData *sr, Coords *coords);
    176159int       gaussj (double **a, int n, double **b, int m);
    177160void      gcenter (CmpCatalog *Target, RefCatalog *Ref);
    … …  
    209192void      plot_init_gridplot ();
    210193void      plot_lumfunc (CmpCatalog *Target, RefCatalog *Ref);
    211 void      plot_resid (StarData *st, StarData *sr);
     194void      plot_resid (StarData *st, StarData *sr, Coords *coords);
    212195void      plot_resid_init (int version, double xmax);
    213196void      plot_resid_plot (int version, float *xvect, float *yvect, int Nvect);

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

    r7996 r8241  
    2727  ScanConfig (config, "ASEC_PIX",          "%lf", 0, &ASEC_PIX);
    2828  ScanConfig (config, "NFIELD",            "%lf", 0, &NFIELD);
     29  ScanConfig (config, "NGRID_PIX",         "%d",  0, &NGRID_PIX);
    2930  ScanConfig (config, "NPOLYTERMS",        "%d",  0, &NPOLYTERMS);
    3031  ScanConfig (config, "ROT_ZERO",          "%lf", 0, &ROT_ZERO);
    … …  
    5859      exit (1);
    5960  }
     61  if (NGRID_PIX == 0) NGRID_PIX = 50.0;
    6062
    6163  if (strcasecmp (ROUGH_ASTROMETRY, "header") &&

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

    r7997 r8241  
    122122  }
    123123
     124  /** XXX temporary trick to deal with the very wide-field ASCA images
     125      this alters the definition of the reference field boundaries */
     126  ASCA = FALSE;
     127  if ((N = get_argument (*argc, argv, "-asca"))) {
     128    ASCA = TRUE;
     129    remove_argument (N, argc, argv);
     130  }
     131
    124132  CATDUMP = FALSE;
    125133  if ((N = get_argument (*argc, argv, "-dump"))) {

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

    r3390 r8241  
    3030  }
    3131
    32   for (i = 0; i < 3; i++) {
     32  for (i = 0; i < 5; i++) {
    3333    gfit (&Target, &Ref, MIN (MAX (1, NPOLYTERMS), 3));
    3434    fprintf (stderr, "precision: %f\n", Target.answer.dR / sqrt(Target.answer.N));

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

    r2442 r8241  
    11# include "gastro2.h"
    2 
    3 static int NPAIR, Npair;
    4 static int *idx1, *idx2;
    52
    63void gfit (CmpCatalog *Target, RefCatalog *Ref, int order) {
    74
    85  int i, j, j0;
     6  int Npair, *idx1, *idx2;
    97  double Radius, Radius2;
    108  double dX, dY, dR;
    … …  
    1210  StarData *st, *sr;
    1311
     12  /* XXX why is this hardwired here? */
    1413  NFIELD = 0.1;
    1514  gproject (Target, Ref, &Subset);
    … …  
    2524  sort_stars_X (Subset.stars, Subset.N);
    2625
    27   Radius = MAX (2.0 * Target[0].answer.dR, 0.5);
     26  Radius = MAX (2.5 * Target[0].answer.dR, 0.5);
    2827  Radius2 = Radius*Radius;
    2928
    … …  
    6362  }
    6463 
    65 # if (0)
    66   pair_init ();
    67   for (i = 0; i < Target[0].N; i++) {
    68     pair_add (i, i);
    69   }
    70 # endif
    71 
    72   fit_init (order);
    73   fit_norm ();
    74 
    75   plot_resid_init (0, (double) Target[0].header.Naxis[0]);
    76   plot_resid_init (1, (double) Target[0].header.Naxis[1]);
    77   if (PLOTSTUFF) plot_resid (st, sr);
    78 
     64  Npair = pair_lists (&idx1, &idx2);
    7965  /* find fit for matched pairs */
    8066  fit_init (order);
    8167  for (i = 0; i < Npair; i++) {
    82     fit_add (st[idx1[i]].X, st[idx1[i]].Y, sr[idx2[i]].X, sr[idx2[i]].Y, 1.0);
     68    fit_add (st[idx1[i]].X, st[idx1[i]].Y, sr[idx2[i]].P, sr[idx2[i]].Q, 1.0);
    8369  }
    8470  fit_eval ();
    8571
    86   Target[0].answer.dR = fit_scat (st, sr);
     72  /* XXX this is weak: the fit_scat call requires the coords from the fit_adjust call */
    8773  Target[0].answer.N  = fit_adjust (&Target[0].coords);
     74  Target[0].answer.dR = fit_scat (st, sr, &Target[0].coords);
    8875
     76  if (PLOTSTUFF) fprintf (stderr, "ploting resid (2)\n");
    8977  plot_resid_init (0, (double) Target[0].header.Naxis[0]);
    9078  plot_resid_init (1, (double) Target[0].header.Naxis[1]);
    91   if (PLOTSTUFF) plot_resid (st, sr);
     79  if (PLOTSTUFF) plot_resid (st, sr, &Target[0].coords);
    9280  free (idx1);
    9381  free (idx2);
    9482}
    95 
    96 static int NTERM, NPOWR, NPARS, NORDER, Npts;
    97 static double **sum, **xsum, **ysum;
    98 static double **matrix, **vector;
    99 
    100 void fit_init (int order) {
    101 
    102   int i;
    103 
    104   NORDER = order;
    105   NPOWR = NORDER + 1;
    106   NTERM = 2*NORDER + 1;
    107   NPARS = (NORDER + 1)*(NORDER + 2) / 2;
    108   Npts  = 0;
    109 
    110   /* allocate arrays for fit solution */
    111   ALLOCATE (sum, double *, NTERM);
    112   ALLOCATE (xsum, double *, NTERM);
    113   ALLOCATE (ysum, double *, NTERM);
    114   for (i = 0; i < NTERM; i++) {
    115     ALLOCATE (sum[i], double, NTERM);
    116     bzero (sum[i], NTERM*sizeof(double));
    117     ALLOCATE (xsum[i], double, NTERM);
    118     bzero (xsum[i], NTERM*sizeof(double));
    119     ALLOCATE (ysum[i], double, NTERM);
    120     bzero (ysum[i], NTERM*sizeof(double));
    121   }
    122   ALLOCATE (matrix, double *, NPARS);
    123   ALLOCATE (vector, double *, NPARS);
    124   for (i = 0; i < NPARS; i++) {
    125     ALLOCATE (matrix[i], double, NPARS);
    126     ALLOCATE (vector[i], double, 2);
    127     bzero (vector[i], 2*sizeof(double));
    128     bzero (matrix[i], NPARS*sizeof(double));
    129   }
    130 
    131 }
    132 
    133 void fit_add (double x1, double y1, double x2, double y2, double wt) {
    134 
    135   int n, m;
    136   double xterm, yterm, term;
    137 
    138   xterm = 1;
    139   for (n = 0; n < NTERM; n++) {
    140     yterm = 1;
    141     for (m = 0; m < NTERM; m++) {
    142       term = xterm*yterm;
    143       if (n+m < NTERM) {
    144         sum[n][m] += term;
    145       }
    146       if (n+m < NPOWR) {
    147         xsum[n][m] += x2*term;
    148         ysum[n][m] += y2*term;
    149       }
    150       yterm *= y1;
    151     }
    152     xterm *= x1;
    153   }
    154   Npts ++;
    155 
    156 }
    157 
    158 void fit_eval () {
    159 
    160   int i, j, n, m, M, N;
    161   double max;
    162 
    163   fprintf (stderr, "npts: %d\n", Npts);
    164 
    165   i = 0;
    166   for (m = 0; m < NPOWR; m++) {
    167     for (n = 0; n < NPOWR - m; n++, i++) {
    168       vector[i][0] = xsum[n][m];
    169       vector[i][1] = ysum[n][m];
    170     }   
    171   }
    172   j = 0;
    173   for (M = 0; M < NPOWR; M++) {
    174     for (N = 0; N < NPOWR - M; N++, j++) {
    175       i = 0;
    176       for (m = 0; m < NPOWR; m++) {
    177         for (n = 0; n < NPOWR - m; n++, i++) {
    178           matrix[i][j] = sum[n+N][m+M];
    179         }       
    180       }
    181     }
    182   }       
    183   max = 0.0;
    184   for (i = 0; i < NPARS; i++) {
    185     for (j = 0; j < NPARS; j++) {
    186       max = MAX (max, fabs(matrix[i][j]));
    187     }
    188     max = MAX (max, fabs(vector[i][0]));
    189     max = MAX (max, fabs(vector[i][1]));
    190   }
    191   for (i = 0; i < NPARS; i++) {
    192     for (j = 0; j < NPARS; j++) {
    193       matrix[i][j] /= max;
    194     }
    195     vector[i][0] /= max;
    196     vector[i][1] /= max;
    197   }
    198 # if (0)
    199   for (i = 0; i < NPARS; i++) {
    200     for (j = 0; j < NPARS; j++) {
    201       fprintf (stderr, "%10.4e  ", matrix[i][j]);
    202     }
    203     fprintf (stderr, "  %10.4e  ", vector[i][0]);
    204     fprintf (stderr, "  %10.4e  \n", vector[i][1]);
    205   }
    206 # endif
    207   gaussj (matrix, NPARS, vector, 2);
    208 # if (0)
    209   fprintf (stderr, "\n\n");
    210   for (i = 0; i < NPARS; i++) {
    211     for (j = 0; j < NPARS; j++) {
    212       fprintf (stderr, "%11.4e  ", matrix[i][j]);
    213     }
    214     fprintf (stderr, "  %11.4e  ", vector[i][0]);
    215     fprintf (stderr, "  %11.4e  \n", vector[i][1]);
    216   }
    217 # endif
    218   i = 0;
    219   for (m = 0; m < NPOWR; m++) {
    220     for (n = 0; n < NPOWR - m; n++, i++) {
    221       xsum[n][m] = vector[i][0];
    222       ysum[n][m] = vector[i][1];
    223     }   
    224   }
    225   i = 0;
    226   for (m = 0; m < NPOWR; m++) {
    227     for (n = 0; n < NPOWR - m; n++, i++) {
    228       fprintf (stderr, "RA x^%dy^%d: %10.4g    DEC x^%dy^%d: %10.4g \n",
    229                n, m, vector[i][0], n, m, vector[i][1]);
    230     }   
    231   }
    232 }
    233 
    234 void fit_norm () {
    235 
    236   xsum[0][0] = 0;
    237   xsum[1][0] = 1;
    238   xsum[0][1] = 0;
    239 
    240   ysum[0][0] = 0;
    241   ysum[1][0] = 0;
    242   ysum[0][1] = 1;
    243 }
    244 
    245 void fit_apply (double *x, double *y, double X, double Y) {
    246 
    247   int m, n;
    248   double xterm, yterm;
    249   double Xo, Yo;
    250 
    251   Xo = Yo = 0;
    252   yterm = 1;
    253   for (m = 0; m < NPOWR; m++) {
    254     xterm = 1;
    255     for (n = 0; n < NPOWR - m; n++) {
    256       Xo += xterm*yterm*xsum[n][m];
    257       Yo += xterm*yterm*ysum[n][m];
    258       xterm *= X;
    259     }   
    260     yterm *= Y;
    261   }
    262  
    263   *x = Xo;
    264   *y = Yo;
    265 }
    266 
    267 
    268 void pair_init () {
    269 
    270   Npair = 0;
    271   NPAIR = 100;
    272   ALLOCATE (idx1, int, NPAIR);
    273   ALLOCATE (idx2, int, NPAIR);
    274 
    275 }
    276 
    277 void pair_add (int i1, int i2) {
    278 
    279   idx1[Npair] = i1;
    280   idx2[Npair] = i2;
    281 
    282   Npair ++;
    283 
    284   if (Npair == NPAIR) {
    285     NPAIR += 100;
    286     REALLOCATE (idx1, int, NPAIR);
    287     REALLOCATE (idx2, int, NPAIR);
    288   }
    289 
    290 }
    291 
    292 double fit_scat (StarData *st, StarData *sr) {
    293 
    294   int i;
    295   double x, y, dx, dy, dX, dY, dX2, dY2, dR;
    296  
    297   dX = dY = dX2 = dY2 = 0;
    298   for (i = 0; i < Npair; i++) {
    299 
    300     fit_apply (&x, &y, sr[idx2[i]].X, sr[idx2[i]].Y);
    301    
    302     dx = x - st[idx1[i]].X;
    303     dy = y - st[idx1[i]].Y;
    304    
    305     dX += dx;
    306     dY += dy;
    307     dX2 += dx*dx;
    308     dY2 += dy*dy;
    309    
    310   }
    311 
    312   dX = dX / Npair;
    313   dY = dY / Npair;
    314   fprintf (stderr, "scatter: %f, %f\n", sqrt(dX2/Npair - dX*dX), sqrt(dY2/Npair - dY*dY));
    315   fprintf (stderr, "precise: %f, %f\n", sqrt((dX2/Npair - dX*dX) / Npair), sqrt((dY2/Npair - dY*dY)/Npair));
    316 
    317   dR = 0.5 * sqrt(fabs(dX2/Npair - dX*dX)) + 0.5 * sqrt (fabs(dY2/Npair - dY*dY));
    318   return (dR);
    319 }
    320 
    321 /* convert new terms to adjustments in coords and to polyterms */
    322 int fit_adjust (Coords *coords) {
    323 
    324   int i, j, N;
    325   double S1, S2, p11, p12, p21, p22;
    326   double a0, a1, a2, b0, b1, b2, det;
    327   double X, Y;
    328   int Np, Nv;
    329    
    330   S1 = coords[0].cdelt1;
    331   S2 = coords[0].cdelt2;
    332   p11 = coords[0].pc1_1;    p12 = coords[0].pc1_2;
    333   p21 = coords[0].pc2_1;    p22 = coords[0].pc2_2;
    334    
    335   /* get the correct vector entries for the linear terms */
    336   N = mk_vector (0, 0, NORDER);
    337   a0 = vector[N][0];  b0 = vector[N][1];
    338   N = mk_vector (1, 0, NORDER);
    339   a1 = vector[N][0];  b1 = vector[N][1];
    340   N = mk_vector (0, 1, NORDER);
    341   a2 = vector[N][0];  b2 = vector[N][1];
    342 
    343   det = 1.0 / (a1*b2 - a2*b1);
    344 
    345   coords[0].pc1_1 = p11*a1 + p12*b1*(S2/S1);
    346   coords[0].pc2_1 = p21*a1 + p22*b1*(S2/S1);
    347    
    348   coords[0].pc1_2 = p12*b2 + p11*a2*(S1/S2);
    349   coords[0].pc2_2 = p22*b2 + p21*a2*(S1/S2);
    350    
    351   X = (coords[0].crpix1 - a0);
    352   Y = (coords[0].crpix2 - b0);
    353   coords[0].crpix1 = det*(X*b2 - Y*a2);
    354   coords[0].crpix2 = det*(Y*a1 - X*b1);
    355 
    356   coords[0].Npolyterms = NORDER;
    357   if (NORDER > 1) strcpy (coords[0].ctype, "DEC--PLY");
    358 
    359   /* generate higher order terms from vector */
    360   for (i = 0; i < NORDER + 1; i++) {
    361     for (j = 0; j < (NORDER - i + 1); j++) {
    362       if (i + j < 2) continue;
    363       Np = mk_polyterm (i, j, NORDER);
    364       Nv = mk_vector (i, j, NORDER);
    365       coords[0].polyterms[Np][0] = det*(vector[Nv][0]*b2  - vector[Nv][1]*a2);  /* x2 y0 */
    366       coords[0].polyterms[Np][1] = det*(vector[Nv][1]*a1  - vector[Nv][0]*b1);  /* x2 y0 */
    367     }
    368   }
    369   while (coords[0].crval1 < 0) coords[0].crval1 += 360.0;
    370   while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0;
    371 
    372   /* test for valid solution */
    373   return (Npts);
    374 
    375 }
    376 
    377 int mk_polyterm (int n, int m, int norder) {
    378  
    379   int i, nt, N;
    380  
    381   N = 0;
    382   nt = n + m;
    383   for (i = 2; i < nt; i++) {
    384     N += i + 1;
    385   }
    386   N += m;
    387   return (N);
    388 }
    389 
    390 int mk_vector (int n, int m, int norder) {
    391  
    392   int i, N;
    393  
    394   N = 0;
    395   for (i = 0; i < m; i++) {
    396     N += (norder - i + 1);
    397   }
    398   N += n;
    399   return (N);
    400 }
    401 
    402 void plot_resid (StarData *st, StarData *sr) {
    403 
    404   int i;
    405   double x, y, dx, dy;
    406   float *xvect0, *yvect0, *xvect1, *yvect1, *xvect2, *yvect2;
    407   int Nvect;
    408 
    409   Nvect = Npair;
    410   ALLOCATE (xvect0, float, Nvect);
    411   ALLOCATE (yvect0, float, Nvect);
    412   ALLOCATE (xvect1, float, Nvect);
    413   ALLOCATE (yvect1, float, Nvect);
    414 
    415   ALLOCATE (xvect2, float, 2*Nvect);
    416   ALLOCATE (yvect2, float, 2*Nvect);
    417  
    418   for (i = 0; i < Npair; i++) {
    419     fit_apply (&x, &y, st[idx1[i]].X, st[idx1[i]].Y);
    420    
    421     dx = x - sr[idx2[i]].X;
    422     dy = y - sr[idx2[i]].Y;
    423    
    424     xvect0[i] = sr[idx2[i]].X;
    425     xvect1[i] = sr[idx2[i]].Y;
    426     yvect0[i] = dx;
    427     yvect1[i] = dy;
    428 
    429     xvect2[2*i+0] = sr[idx2[i]].X;
    430     yvect2[2*i+0] = sr[idx2[i]].Y;
    431     xvect2[2*i+1] = x;
    432     yvect2[2*i+1] = y;
    433   }
    434 
    435   plot_resid_plot (0, xvect0, yvect0, Nvect);
    436   plot_resid_plot (1, xvect1, yvect1, Nvect);
    437   plot_fullfield_pairs (xvect2, yvect2, 2*Nvect);
    438 
    439   free (xvect2);
    440   free (yvect2);
    441   free (xvect0);
    442   free (yvect0);
    443   free (xvect1);
    444   free (yvect1);
    445 
    446 }

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

    r7997 r8241  
    44
    55  int i, N;
    6   double X, Y, M, Moff;
     6  double X, Y, P, Q, M, Moff;
    77  double XMIN, XMAX, YMIN, YMAX, MMIN, MMAX;
    8   Coords *coords;
     8  Coords TPtoSky, FPtoTP, *coords;
    99  StarData *in, *out;
    1010
    … …  
    1515  in     = Ref[0].stars;
    1616  out    = Subset[0].stars;
     17
    1718  coords = &Target[0].coords;
     19
     20  /* create Tangent Plane to Sky transformation from input coords */
     21  strcpy (TPtoSky.ctype, coords[0].ctype);
     22  TPtoSky.crval1 = coords[0].crval1;
     23  TPtoSky.crval2 = coords[0].crval2;
     24
     25  TPtoSky.cdelt1 = TPtoSky.cdelt2 = 1;
     26  TPtoSky.crpix1 = TPtoSky.crpix2 = 0;
     27  TPtoSky.pc1_1 = TPtoSky.pc2_2 = 1;
     28  TPtoSky.pc1_2 = TPtoSky.pc2_1 = 0;
     29  TPtoSky.Npolyterms = 0;
     30  for (i = 0; i < 7; i++) {
     31    TPtoSky.polyterms[i][0] = 0;
     32    TPtoSky.polyterms[i][1] = 0;
     33  }
     34
     35  /* create Focal Plane to Tangent Plane transformation from input coords */
     36  strcpy (FPtoTP.ctype, "FP---LIN");
     37  FPtoTP.crval1 = FPtoTP.crval2 = 0;
     38
     39  FPtoTP.cdelt1 = coords[0].cdelt1;
     40  FPtoTP.cdelt2 = coords[0].cdelt2;
     41  FPtoTP.crpix1 = coords[0].crpix1;
     42  FPtoTP.crpix2 = coords[0].crpix2;
     43  FPtoTP.pc1_1  = coords[0].pc1_1;
     44  FPtoTP.pc1_2  = coords[0].pc1_2;
     45  FPtoTP.pc2_1  = coords[0].pc2_1;
     46  FPtoTP.pc2_2  = coords[0].pc2_2;
     47
     48  FPtoTP.Npolyterms = coords[0].Npolyterms;
     49  for (i = 0; i < 7; i++) {
     50    FPtoTP.polyterms[i][0] = coords[0].polyterms[i][0];
     51    FPtoTP.polyterms[i][1] = coords[0].polyterms[i][1];
     52  }
    1853
    1954  Moff = Ref[0].Moff;
    … …  
    5590    }
    5691
     92    /* get tangent-plane coordinates as well */
     93    RD_to_XY (&P, &Q, in[i].R, in[i].D, &TPtoSky);
     94
    5795    out[N] = in[i];
    5896    out[N].X = X;
    5997    out[N].Y = Y;
     98    out[N].P = P;
     99    out[N].Q = Q;
    60100    out[N].M = M;
    61101    N++;
    … …  
    77117  REALLOCATE (Subset[0].stars, StarData, MAX (1, Subset[0].N));
    78118
     119  /* get tangent-plane coords for target as well */
     120  for (i = 0; i < Target[0].N; i++) {
     121    XY_to_RD (&P, &Q, Target[0].stars[i].X, Target[0].stars[i].Y, &FPtoTP);
     122    Target[0].stars[i].P = P;
     123    Target[0].stars[i].Q = Q;
     124  }
     125
    79126}
    80127

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

    r7962 r8241  
    6868  Yo = -0.5*NFIELD*NY;
    6969
    70   catstats[0].RA[0] = +360.0;
    71   catstats[0].RA[1] = -360.0;
    72   catstats[0].DEC[0] = +90.0;
    73   catstats[0].DEC[1] = -90.0;
     70  if (ASCA) {
     71      XY_to_RD (&R, &D, 0.5*NX, 0.5*NY, &Target[0].coords);
     72      catstats[0].RA[0]  = R - 90.0;
     73      catstats[0].RA[1]  = R + 90.0;
     74      catstats[0].DEC[0] = MAX (-90.0, D - 90.0);
     75      catstats[0].DEC[1] = MIN (+90.0, D + 90.0);
     76      if (VERBOSE) fprintf (stderr, "asca region: %f - %f, %f - %f\n",
     77                            catstats[0].RA[0], catstats[0].RA[1], catstats[0].DEC[0], catstats[0].DEC[1]);
     78      return;
     79  }
    7480
    7581  for (x = 0; x <= 1.0; x += 0.5) {
    … …  
    8995    }
    9096  }
     97
     98  /* is a pole in the image?  if so, include it... */
     99  status = RD_to_XY (&X, &Y, 0.0, 90.0, &Target[0].coords);
     100  if (status) catstats[0].DEC[1] = 90.0;
     101
     102  status = RD_to_XY (&X, &Y, 0.0, -90.0, &Target[0].coords);
     103  if (status) catstats[0].DEC[0] = -90.0;
     104
    91105  if (VERBOSE) fprintf (stderr, "full region: %f - %f, %f - %f\n",
    92106           catstats[0].RA[0], catstats[0].RA[1], catstats[0].DEC[0], catstats[0].DEC[1]);

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

    r7997 r8241  
    11# include "gastro2.h"
    22
    3 # define NPIX 50
     3/* this needs to be user-configured */
    44static int NX, NY, Nbin;
    55static double C0x, C1x, C0y, C1y;
    … …  
    9898    if (PLOTSTUFF) plot_done_gridplot (dX, dY);
    9999
    100     XMIN -= 0.5*NPIX;
    101     XMAX -= 0.5*NPIX;
    102     YMIN -= 0.5*NPIX;
    103     YMAX -= 0.5*NPIX;
     100    XMIN -= 0.5*NGRID_PIX;
     101    XMAX -= 0.5*NGRID_PIX;
     102    YMIN -= 0.5*NGRID_PIX;
     103    YMAX -= 0.5*NGRID_PIX;
    104104    gridfree ();
    105105  }
    … …  
    111111void gridinit (double XMIN, double XMAX, double YMIN, double YMAX, int Nr, int Nt) {
    112112
    113   NX = (XMAX - XMIN) / NPIX;
    114   NY = (YMAX - YMIN) / NPIX;
     113  NX = (XMAX - XMIN) / NGRID_PIX;
     114  NY = (YMAX - YMIN) / NGRID_PIX;
    115115
    116116  C1x =          NX / (XMAX - XMIN);

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

    r2442 r8241  
    169169
    170170  int i, Nvect;
    171   float *xvect, *yvect, *zvect, M, dM;
     171  float *xvect, *yvect, *zvect, M, dM, mRefMin, mRefMax;
    172172  char c;
    173173 
    … …  
    216216  ALLOCATE (yvect, float, Nvect);
    217217  ALLOCATE (zvect, float, Nvect);
     218
     219  mRefMin = 30.0;
     220  mRefMax = -5.0;
     221  for (i = 0; i < Nvect; i++) {
     222      mRefMin = MIN (mRefMin, Ref[0].stars[i].M);
     223      mRefMax = MAX (mRefMax, Ref[0].stars[i].M);
     224  }
     225  dM = mRefMax - mRefMin;
     226
    218227  for (i = 0; i < Nvect; i++) {
    219228    xvect[i] = Ref[0].stars[i].X;
    220229    yvect[i] = Ref[0].stars[i].Y;
    221     M = (Target[0].lum.Mmax - Ref[0].stars[i].M) / dM;
     230    M = (mRefMax - Ref[0].stars[i].M) / dM;
    222231    zvect[i] = MIN (1.0, MAX (0.01, M));
    223232  }
    … …  
    330339}
    331340
     341
     342void plot_resid (StarData *st, StarData *sr, Coords *coords) {
     343
     344  int i;
     345  double x, y, dx, dy;
     346  float *xvect0, *yvect0, *xvect1, *yvect1, *xvect2, *yvect2;
     347  int Npair, *idx1, *idx2;
     348
     349  Npair = pair_lists (&idx1, &idx2);
     350  ALLOCATE (xvect0, float, Npair);
     351  ALLOCATE (yvect0, float, Npair);
     352  ALLOCATE (xvect1, float, Npair);
     353  ALLOCATE (yvect1, float, Npair);
     354
     355  ALLOCATE (xvect2, float, 2*Npair);
     356  ALLOCATE (yvect2, float, 2*Npair);
     357 
     358  for (i = 0; i < Npair; i++) {
     359    RD_to_XY (&x, &y, sr[idx2[i]].R, sr[idx2[i]].D, coords);
     360   
     361    dx = x - st[idx1[i]].X;
     362    dy = y - st[idx1[i]].Y;
     363   
     364    xvect0[i] = st[idx1[i]].X;
     365    xvect1[i] = st[idx1[i]].Y;
     366    yvect0[i] = dx;
     367    yvect1[i] = dy;
     368
     369    xvect2[2*i+0] = st[idx1[i]].X;
     370    yvect2[2*i+0] = st[idx1[i]].Y;
     371    xvect2[2*i+1] = x;
     372    yvect2[2*i+1] = y;
     373  }
     374
     375  plot_resid_plot (0, xvect0, yvect0, Npair);
     376  plot_resid_plot (1, xvect1, yvect1, Npair);
     377  plot_fullfield_pairs (xvect2, yvect2, 2*Npair);
     378
     379  free (xvect2);
     380  free (yvect2);
     381  free (xvect0);
     382  free (yvect0);
     383  free (xvect1);
     384  free (yvect1);
     385
     386}
     387
    332388void fill_lumfunc (StarData *stars, int N, float *lbin, float *bin, int *nb) {
    333389
Note: See TracChangeset for help on using the changeset viewer.