Changeset 12220

Timestamp:

Mar 4, 2007, 3:45:18 PM (19 years ago)

Author:

eugene

Message:

finishing up the basic code; almost compiles

Location:

branches/dvo-mods-2007-02/Ohana/src/relastro

Files:

• Makefile (modified) (1 diff)
• include/relastro.h (modified) (8 diffs)
• src/FitMosaic.c (modified) (1 diff)
• src/FitPM.c (modified) (3 diffs)
• src/FitPMandPar.c (modified) (3 diffs)
• src/FitSimple.c (modified) (1 diff)
• src/ImageOps.c (modified) (5 diffs)
• src/ParFactor.c (modified) (2 diffs)
• src/StarOps.c (modified) (9 diffs)
• src/UpdateChips.c (modified) (2 diffs)
• src/UpdateMosaic.c (modified) (2 diffs)
• src/UpdateObjects.c (modified) (8 diffs)
• src/UpdateSimple.c (modified) (2 diffs)
• src/args.c (modified) (1 diff)
• src/bcatalog.c (modified) (5 diffs)
• src/fitpoly.c (modified) (8 diffs)
• src/mkpolyterm.c (modified) (1 diff)
• src/plot_scatter.c (modified) (1 diff)
• src/plotstuff.c (modified) (1 diff)
• src/relastro.c (modified) (2 diffs)

branches/dvo-mods-2007-02/Ohana/src/relastro/Makefile

@@ -21 +21 @@
 
 RELASTRO = \
-$(SRC)/relastro.$(ARCH).o        \
-$(SRC)/initialize.$(ARCH).o      \
 $(SRC)/ConfigInit.$(ARCH).o      \
-$(SRC)/args.$(ARCH).o            \
-$(SRC)/liststats.$(ARCH).o       \
-$(SRC)/load_images.$(ARCH).o     \
-$(SRC)/name_region.$(ARCH).o     \
-$(SRC)/find_images.$(ARCH).o     \
-$(SRC)/find_regions.$(ARCH).o    \
-$(SRC)/get_regions.$(ARCH).o     \
-$(SRC)/load_catalogs.$(ARCH).o   \
-$(SRC)/gcatalog.$(ARCH).o        \
-$(SRC)/bcatalog.$(ARCH).o        \
-$(SRC)/GridOps.v2.$(ARCH).o      \
-$(SRC)/StarOps.$(ARCH).o         \
+$(SRC)/FitChip.$(ARCH).o         \
+$(SRC)/FitMosaic.$(ARCH).o       \
+$(SRC)/FitPM.$(ARCH).o           \
+$(SRC)/FitPMandPar.$(ARCH).o     \
+$(SRC)/FitSimple.$(ARCH).o       \
 $(SRC)/ImageOps.$(ARCH).o        \
 $(SRC)/MosaicOps.$(ARCH).o       \
+$(SRC)/ParFactor.$(ARCH).o       \
+$(SRC)/SetSignals.$(ARCH).o      \
+$(SRC)/Shutdown.$(ARCH).o        \
+$(SRC)/StarOps.$(ARCH).o         \
+$(SRC)/UpdateChips.$(ARCH).o     \
+$(SRC)/UpdateMosaic.$(ARCH).o    \
+$(SRC)/UpdateObjects.$(ARCH).o   \
+$(SRC)/UpdateSimple.$(ARCH).o    \
+$(SRC)/args.$(ARCH).o            \
+$(SRC)/bcatalog.$(ARCH).o        \
+$(SRC)/dvo_astrom_ops.$(ARCH).o  \
+$(SRC)/fitpoly.$(ARCH).o         \
 $(SRC)/global_stats.$(ARCH).o    \
+$(SRC)/initialize.$(ARCH).o      \
+$(SRC)/liststats.$(ARCH).o       \
+$(SRC)/load_catalogs.$(ARCH).o   \
+$(SRC)/load_images.$(ARCH).o     \
+$(SRC)/mkpolyterm.$(ARCH).o      \
 $(SRC)/plot_scatter.$(ARCH).o    \
 $(SRC)/plotstuff.$(ARCH).o       \
-$(SRC)/misc.$(ARCH).o            \
+$(SRC)/relastro.$(ARCH).o        \
 $(SRC)/reload_catalogs.$(ARCH).o \
-$(SRC)/image-db.$(ARCH).o \
-$(SRC)/free_catalogs.$(ARCH).o   \
-$(SRC)/check_permissions.$(ARCH).o   \
+$(SRC)/setExclusions.$(ARCH).o   \
+$(SRC)/write_coords.$(ARCH).o
+
+DROP = \
 $(SRC)/setMrelFinal.$(ARCH).o    \
-$(SRC)/write_coords.$(ARCH).o    \
-$(SRC)/setExclusions.$(ARCH).o   \
-$(SRC)/Shutdown.$(ARCH).o        \
-$(SRC)/SetSignals.$(ARCH).o      \
-$(SRC)/wcatalog.$(ARCH).o
+
 
 $(RELASTRO): $(INC)/relastro.h

branches/dvo-mods-2007-02/Ohana/src/relastro/include/relastro.h

@@ -4 +4 @@
 # include <signal.h>
 
+typedef struct {
+  double R, D;  /* Sky Coords    - degrees */
+  double P, Q;  /* Tangent Plane - pixels  */
+  double L, M;  /* Focal Plane   - pixels  */
+  double X, Y;  /* Chip Coords   - pixels  */
+  double Mag, dMag;
+  int mask;
+} StarData;
+
 // structure to hold coordinate fitting terms
 typedef struct {
     int Npts;
     int Nterms;
-    int Npower;
     int Norder;
+    int Nsums;
+    int Nelems;
+    double **sum;
     double **xsum;
     double **ysum;
-
-    double **vector;
-    double **matrix;
+    double **xfit;
+    double **yfit;
 } CoordFit;
 
@@ -28 +38 @@
   double chisq;
   int Nfit;
-} PMfit;
+} PMFit;
 
 typedef struct {
@@ -62 +72 @@
 char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
 char   CameraConfig[256];
+char   SKY_TABLE[256];
+int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */
 
 double MAG_LIM;
@@ -108 +120 @@
 enum {black, white, red, orange, yellow, green, blue, indigo, violet};
 
+int DophotSelect, DophotValue;
+
 int TimeSelect;
 time_t TSTART, TSTOP;
 
-# ifdef GRID_V1
-int setGridMeasure (int meas, int cat, double X, double Y);
-# endif
-
-# ifdef GRID_V2
-int setGridMeasure (int meas, int cat, double X, double Y, int ccdnum);
-# endif
+int FITTING_MODE;
+enum {FIT_AVERAGE, FIT_PM_ONLY, FIT_PM_AND_PAR};
+
+SkyRegion UserPatch;
+int UserPatchSelect;
+
+int DoUpdateObjects;
+int DoUpdateSimple;
+int DoUpdateChips;
+int DoUpdateMosaics;
 
 /*** relphot prototypes ***/
@@ -126 +143 @@
 int          *SelectRefMosaic     PROTO((Mosaic **refmosaic, int *Nimage));
 int           args                PROTO((int argc, char **argv));
-int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog, GSCRegion *fullregion));
+int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
 void          clean_images        PROTO(());
 void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final));
@@ -171 +188 @@
 void          initstats           PROTO((char *mode));
 int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
-Catalog      *load_catalogs       PROTO((GSCRegion *region, int Nregion, GSCRegion *fullregion));
-GSCRegion    *load_images         PROTO((FITS_DB *db, char *seed, int *nregion, GSCRegion *fullregion));
+Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog));
+SkyList      *load_images         PROTO((FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect));
+Image        *select_images       PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage));
 
 void check_permissions (char *basefile);
@@ -195 +213 @@
 void          plot_star_coords    PROTO((Catalog *catalog, int Ncatalog));
 void          plot_stars          PROTO((Catalog *catalog, int Ncatalog));
-void          reload_catalogs     PROTO((GSCRegion *region, int Nregion, GSCRegion *fullregion));
+void          reload_catalogs     PROTO((SkyList *skylist));
 int           setExclusions       PROTO((Catalog *catalog, int Ncatalog));
 void          setMcal             PROTO((Catalog *catalog, int Poor));
@@ -226 +244 @@
 void          wimages             PROTO(());
 void          write_coords        PROTO((Header *header, Coords *coords));
+
+
+double **array_init (int Nx, int Ny);
+void array_free (double **array, int Nx);
+CoordFit *fit_init (int order);
+void fit_free (CoordFit *fit);
+void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt);
+void fit_eval (CoordFit *fit);
+double **poly2d_dx (double **poly, int Nx, int Ny);
+double **poly2d_dy (double **poly, int Nx, int Ny);
+double **poly2d_copy (double **poly, int Nx, int Ny);
+double poly2d_eval (double **poly, int Nx, int Ny, double x, double y);
+void fit_apply_coords (CoordFit *fit, Coords *coords);
+int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo);
+CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo);
+void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void initObjectData (Catalog *catalog, int Ncatalog);
+int UpdateObjects (Catalog *catalog, int Ncatalog);
+int UpdateSimple (Catalog *catalog, int Ncatalog);
+int UpdateChips (Catalog *catalog, int Ncatalog);
+int UpdateMosaic (Catalog *catalog, int Ncatalog);
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
+int ParFactor (double *pR, double *pD, double R, double D, time_t T);
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts);
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts);
+
+Coords *getMosaicForImage (int N);
+
+StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic);
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic);
+
+double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
+double getMeanD (Measure *measure, Average *average, SecFilt *secfilt);

branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitMosaic.c

@@ -4 +4 @@
 
   int i;
+  CoordFit *fit;
 
-  fit_init (coords[0].Npolyterms);
+  fit = fit_init (coords[0].Npolyterms);
   for (i = 0; i < Nmatch; i++) {
     if (raw[i].mask) continue;
-    fit_add (raw[i].L, raw[i].M, ref[i].P, ref[i].Q);
+    fit_add (fit, raw[i].L, raw[i].M, ref[i].P, ref[i].Q, 1.0);
   }
-  fit_eval ();
-  fit_apply_coords (coords);
+  fit_eval (fit);
+  fit_apply_coords (fit, coords);
+  fit_free (fit);
   /* FitChip and FitSimple update the coords in different ways? maybe not... */
 }

branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c

@@ -2 +2 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) {
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts) {
+
+  int i;
 
   double **A, **B;
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
 
   /* do I need to do this as 2 2x2 matrix equations? */
-  ALLOCATE (A, double *, 4);
-  ALLOCATE (B, double *, 4);
-  for (i = 0; i < 4; i++) {
-    ALLOCATE (A[i], double, 4);
-    ALLOCATE (B[i], double, 1);
-    bzero (A[i], 4*sizeof(double));
-    bzero (B[i], 1*sizeof(double));
-  }
+  A = array_init (4, 4);
+  B = array_init (4, 1);
 
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
   for (i = 0; i < Npts; i++) {
     /* handle case where dX or dY = 0.0 */
@@ -56 +54 @@
   dgaussj (A, 4, B, 1);
 
-  fit[0].Ro = B[0];
-  fit[0].uR = B[1];
-  fit[0].Do = B[2];
-  fit[0].uD = B[3];
+  fit[0].Ro = B[0][0];
+  fit[0].uR = B[1][0];
+  fit[0].Do = B[2][0];
+  fit[0].uD = B[3][0];
   fit[0].p  = 0.0;
   
@@ -68 +66 @@
   fit[0].dp  = 0.0;
   
-  for (i = 0; i < 4; i++) {
-    free (A[i]);
-    free (B[i]);
-  }
-  free (A);
-  free (B);
+  array_free (A, 4);
+  array_free (B, 4);
 
   return (TRUE);

branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPMandPar.c

@@ -2 +2 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts) {
+
+  int i;
 
   double **A, **B;
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+  double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
 
-  ALLOCATE (A, double *, 5);
-  ALLOCATE (B, double *, 5);
-  for (i = 0; i < 5; i++) {
-    ALLOCATE (A[i], double, 5);
-    ALLOCATE (B[i], double, 1);
-    bzero (A[i], 5*sizeof(double));
-    bzero (B[i], 1*sizeof(double));
-  }
+  A = array_init (5, 5);
+  B = array_init (5, 1);
 
+  PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0;
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
   for (i = 0; i < Npts; i++) {
     /* handle case where dX or dY = 0.0 */
@@ -78 +78 @@
   dgaussj (A, 5, B, 1);
 
-  fit[0].Ro = B[0];
-  fit[0].uR = B[1];
-  fit[0].Do = B[2];
-  fit[0].uD = B[3];
-  fit[0].p  = B[4];
+  fit[0].Ro = B[0][0];
+  fit[0].uR = B[1][0];
+  fit[0].Do = B[2][0];
+  fit[0].uD = B[3][0];
+  fit[0].p  = B[4][0];
   
   fit[0].dRo = sqrt(A[0][0]);
@@ -90 +90 @@
   fit[0].dp  = sqrt(A[4][4]);
   
-  for (i = 0; i < 5; i++) {
-    free (A[i]);
-    free (B[i]);
-  }
-  free (A);
-  free (B);
+  array_free (A, 5);
+  array_free (B, 5);
 
   /* get the chisq from the matrix values */

branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c

@@ -4 +4 @@
 
   int i;
+  CoordFit *fit;
 
-  fit_init (coords[0].Npolyterms);
+  fit = fit_init (coords[0].Npolyterms);
   for (i = 0; i < Nmatch; i++) {
     if (raw[i].mask) continue;
-    fit_add (raw[i].X, raw[i].Y, ref[i].L, ref[i].M);
+    fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0);
   }
-  fit_eval ();
-  fit_apply_coords (coords);
+  fit_eval (fit);
+  fit_apply_coords (fit, coords);
+  fit_free (fit);
 }
 

branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c

@@ -94 +94 @@
 void matchImage (Catalog *catalog, int meas, int cat) {
 
-  int i, ave, ccdnum;
-  double ra, dec, X, Y;
-  char *pname, *filter, *p, base[256];
+  int i;
   Measure *measure;
   
@@ -178 +176 @@
 }
 
-StarData *getImageRaw (int im, int *Nstars, int isMosaic) {
-
+/* XXX re-write to use image->nstars if that is possible */
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) {
+
+  int i, m, c;
+
+  Coords *mosaic;
   StarData *raw;
   
   ALLOCATE (raw, StarData, Nlist[im]);
 
+  mosaic = NULL;
   if (isMosaic) {
     mosaic = getMosaicForImage (im);
@@ -202 +205 @@
     XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     if (isMosaic) {
-      XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, &mosaic);
-      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &mosaic);
+      XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, mosaic);
+      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, mosaic);
     } else {
       raw[i].P = raw[i].L;
@@ -215 +218 @@
 }
 
-StarData *getImageRef (int im, int *Nstars, int isMosaic) {
-
+StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) {
+
+  int i, m, c, n;
+
+  Coords *mosaic;
   StarData *ref;
   
   ALLOCATE (ref, StarData, Nlist[im]);
 
+  mosaic = NULL;
   if (isMosaic) {
     mosaic = getMosaicForImage (im);
@@ -237 +244 @@
     RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
     if (isMosaic) {
-      LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, &mosaic);
-      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &mosaic);
+      LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, mosaic);
+      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, mosaic);
     } else {
       ref[i].L = ref[i].P;

branches/dvo-mods-2007-02/Ohana/src/relastro/src/ParFactor.c

@@ -28 +28 @@
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {
 
-  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
+  double n, L, g;
+
+# define J2000 2451545.       /* Julian date at standard epoch */
 
   n = jd - J2000;
@@ -43 +45 @@
 
   double jd;
-  double L, B;
+  double L, B, E, e, s, r, d;
 
   /* given a time T in UNIX seconds, determine the solar longitude S */

branches/dvo-mods-2007-02/Ohana/src/relastro/src/StarOps.c

@@ -28 +28 @@
   if (catalog[cat].average[ave].code & STAR_BAD) return (NO_MAG);  
  
-  value = PhotPrimary ? catalog[cat].average[ave].M : catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
+  value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
   return (value);
 }
@@ -65 +65 @@
       liststats (list, dlist, N, &stats);
 
-      if (PhotPrimary) {
-        catalog[i].average[j].M = stats.mean;
-        catalog[i].average[j].dM = stats.sigma;
-        catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
-      } else {
-        catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
-        catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma;
-        catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
-      }      
+      catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
+      catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma;
+      catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
     }
   }
@@ -121 +115 @@
 
       /* use sigma or error in dM for output? */
-      if (PhotPrimary) {
-        catalog[i].average[j].M = stats.mean;
-        catalog[i].average[j].dM = MAX (stats.sigma, stats.error);
-        catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
-      } else {
-        catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
-        catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma);
-        catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
-      }      
+      catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
+      catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma);
+      catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
     }
   }
@@ -183 +171 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
       if (catalog[i].average[j].code & STAR_BAD) continue;
-      Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
       if (Xm == -1) continue;
       Chisq = pow (10.0, 0.01*Xm);
       xlist[Ntot] = Chisq;
-      slist[Ntot] = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+      slist[Ntot] = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
       dlist[Ntot] = 1;
       Ntot ++;
@@ -203 +191 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
-      Xm = PhotPrimary ? catalog[i].average[j].Xm    : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
       Chisq = pow (10.0, 0.01*Xm);
       mark = (dM > MaxScatter) || (Xm == NO_MAG) || (Chisq > MaxChisq);
@@ -387 +375 @@
       if (catalog[i].average[j].code & STAR_BAD) continue;  
 
-      Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
       if (Xm == NO_MAG) continue;
       list[n] = pow (10.0, 0.01*Xm);
@@ -423 +411 @@
       if (catalog[i].average[j].code & STAR_BAD) continue;  
 
-      dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
       list[n] = dM;
       dlist[n] = 1;
@@ -452 +440 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
       if (catalog[i].average[j].code & STAR_BAD) continue;  
-      dMrel = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+      dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
       bin = dMrel / 0.00025;
       bin = MAX (0, MIN (NBIN-1, bin));
@@ -483 +471 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
       if (catalog[i].average[j].code & STAR_BAD) continue;
-      xlist[N] = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M;
-      value    = PhotPrimary ? catalog[i].average[j].Xm   : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      xlist[N] = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
+      value    = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
       if (value == NO_MAG) continue;
       ylist[N] = 0.01*value;

branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c

@@ -1 +1 @@
 # include "relastro.h"
 
-int UpdateChips () {
+int UpdateChips (Catalog *catalog, int Ncatalog) {
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-
+  int i, Nimage, Nstars;
   Image *image;
-  int Nimage;
+  StarData *raw, *ref;
 
   image = getimages (&Nimage);
@@ -18 +18 @@
 
     /* convert measure coordinates to raw entries */
-    raw = getImageRaw (&image[i], &Nstars, TRUE);
+    raw = getImageRaw (catalog, Ncatalog, i, &Nstars, TRUE);
 
     /* convert average coordinates to ref entries */
-    ref = getImageRef (&image[i], &Nstars, TRUE);
+    ref = getImageRef (catalog, Ncatalog, i, &Nstars, TRUE);
 
     /* XXX grab this code from mosastro */

branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c

@@ -1 +1 @@
 # include "relastro.h"
 
-int UpdateMosaic () {
+int UpdateMosaic (Catalog *catalog, int Ncatalog) {
 
   /* we can measure new image parameters for each mosaic independently */
-
+  int i, Nimage, Nstars;
   Image *image;
-  int Nimage;
+  StarData *raw, *ref;
 
   image = getimages (&Nimage);
@@ -18 +18 @@
 
     /* convert measure coordinates to raw entries */
-    raw = getImageRaw (&image[i], &Nstars, TRUE);
+    raw = getImageRaw (catalog, Ncatalog, i, &Nstars, TRUE);
 
     /* convert average coordinates to ref entries */
-    ref = getImageRef (&image[i], &Nstars, TRUE);
+    ref = getImageRef (catalog, Ncatalog, i, &Nstars, TRUE);
 
     /* XXX grab this code from mosastro */

branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c

@@ -6 +6 @@
 static double *R, *dR;
 static double *D, *dD;
+static double *pX;
+static double *pY;
 static time_t *T;
 static double *dT;
@@ -25 +27 @@
   ALLOCATE (X, double, MAX (1, Nmax));
   ALLOCATE (Y, double, MAX (1, Nmax));
-  ALLOCATE (pR, double, MAX (1, Nmax));
-  ALLOCATE (pD, double, MAX (1, Nmax));
 
   ALLOCATE (dR, double, MAX (1, Nmax));
@@ -33 +33 @@
   ALLOCATE (dX, double, MAX (1, Nmax));
   ALLOCATE (dY, double, MAX (1, Nmax));
+
+  ALLOCATE (pX, double, MAX (1, Nmax));
+  ALLOCATE (pY, double, MAX (1, Nmax));
 }  
 
@@ -38 +41 @@
 
   int i, j, k, m, N;
-  float chisq;
   StatType statsR, statsD;
   Coords coords;
@@ -88 +90 @@
       
       /* project all of the R,D coordinates to a plane centered on this coordinate */
-      /* calculate pR[i], pD[i] for each point */
       for (k = 0; k < N; k++) {
           RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
@@ -100 +101 @@
 
           fit.Ro = statsR.mean;
-          fit.dR = statsR.sigma;
+          fit.dRo = statsR.sigma;
 
           fit.Do = statsD.mean;
-          fit.dD = statsD.sigma;
+          fit.dDo = statsD.sigma;
 
           fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
@@ -117 +118 @@
             ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
           }
-          FitPMandPar (&fit, X, dX, Y, dY, T, pR, pD, N);
+          FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N);
           break;
         default:
@@ -137 +138 @@
       catalog[i].average[j].dP  = fit.dp;
 
-      catalog[i].average[j].Xp  = (fit.Npts > 1) ? 100.0*log10(fit.chisq) : NO_MAG;
+      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NO_MAG;
     }
   }

branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c

@@ -1 +1 @@
 # include "relastro.h"
 
-int UpdateSimple () {
+int UpdateSimple (Catalog *catalog, int Ncatalog) {
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-
+  int i, Nimage, Nstars;
   Image *image;
-  int Nimage;
+  StarData *raw, *ref;
 
   image = getimages (&Nimage);
@@ -19 +19 @@
 
     /* convert measure coordinates to raw entries */
-    raw = getImageRaw (&image[i], &Nstars, FALSE);
+    raw = getImageRaw (catalog, Ncatalog, i, &Nstars, FALSE);
 
     /* convert average coordinates to ref entries */
-    ref = getImageRef (&image[i], &Nstars, FALSE);
+    ref = getImageRef (catalog, Ncatalog, i, &Nstars, FALSE);
 
     /* XXX grab this code from mosastro */

branches/dvo-mods-2007-02/Ohana/src/relastro/src/args.c

@@ -174 +174 @@
 
   /* possible operations */
-  UpdateObjects = FALSE;
+  DoUpdateObjects = FALSE;
   if ((N = get_argument (argc, argv, "-update-objects"))) {
     remove_argument (N, &argc, argv);
-    UpdateObjects = TRUE;
-  }
-  UpdateImages = FALSE;
-  if ((N = get_argument (argc, argv, "-update-images"))) {
-    remove_argument (N, &argc, argv);
-    UpdateImages = TRUE;
+    DoUpdateObjects = TRUE;
+  }
+  DoUpdateSimple = FALSE;
+  if ((N = get_argument (argc, argv, "-update-simple"))) {
+    remove_argument (N, &argc, argv);
+    DoUpdateSimple = TRUE;
+  }
+  DoUpdateChips = FALSE;
+  if ((N = get_argument (argc, argv, "-update-chips"))) {
+    remove_argument (N, &argc, argv);
+    DoUpdateChips = TRUE;
+  }
+  DoUpdateMosaics = FALSE;
+  if ((N = get_argument (argc, argv, "-update-mosaics"))) {
+    remove_argument (N, &argc, argv);
+    DoUpdateMosaics = TRUE;
+  }
+
+  DophotSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-dophot"))) {
+    remove_argument (N, &argc, argv);
+    DophotValue = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    DophotSelect = TRUE;
   }
 

branches/dvo-mods-2007-02/Ohana/src/relastro/src/bcatalog.c

@@ -19 +19 @@
     if (catalog[0].average[i].Nm < 2) continue; 
 
-    /* XXX this limitation is absurd
-    if (catalog[0].average[i].R < fullregion[0].RA[0]) continue;
-    if (catalog[0].average[i].R > fullregion[0].RA[1]) continue;
-    if (catalog[0].average[i].D < fullregion[0].DEC[0]) continue;
-    if (catalog[0].average[i].D > fullregion[0].DEC[1]) continue;
-    */
-
     /* start with all stars good */
     subcatalog[0].average[Naverage] = catalog[0].average[i];
@@ -34 +27 @@
 
     if (RESET) {
-      float *p;
-      p = (PhotPrimary) ? &subcatalog[0].average[Naverage].M : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M;
-      *p = NO_MAG;
-      p = (PhotPrimary) ? &subcatalog[0].average[Naverage].dM : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM;
-      *p = NO_MAG;
+      // XXX reset the ra,dec coords?
       subcatalog[0].average[Naverage].code &= ~ID_STAR_FEW;
       subcatalog[0].average[Naverage].code &= ~ID_STAR_POOR;
@@ -59 +48 @@
 
       /* select measurements by quality */
-      if (catalog[0].measure[offset].dophot != 1) continue;
+      if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue;
 
       /* select measurements by mag limit */
@@ -66 +55 @@
 
       /* select measurements by measurement error */
-      if (catalog[0].measure[offset].dM > SIGMA_LIM) continue;
+      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) continue;
 
       /* select measurements by mag limit */
@@ -111 +100 @@
   subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
 
+  // XXX really do this or not???
   if (catalog[0].Naverage > 0) free (catalog[0].average);
   if (catalog[0].Nmeasure > 0) free (catalog[0].measure);

branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c

@@ -12 +12 @@
 */
 
+double **array_init (int Nx, int Ny) {
+
+  int i;
+  double **array;
+
+  ALLOCATE (array, double *, Nx);
+  for (i = 0; i < Nx; i++) {
+    ALLOCATE (array[i], double, Ny);
+    memset (array[i], 0, Ny*sizeof(double));
+  }    
+  return (array);
+}
+
+void array_free (double **array, int Nx) {
+
+  int i;
+
+  for (i = 0; i < Nx; i++) {
+    free (array[i]);
+  }    
+  free (array);
+}
+
 // XXX define a fit structure and drop the file static variables?
 CoordFit *fit_init (int order) {
 
-  int i;
+  CoordFit *fit;
 
   ALLOCATE (fit, CoordFit, 1);
@@ -29 +52 @@
   // xsum[i][j] holds \sum (x2 wt x1^i y1^j)
   // ysum[i][j] holds \sum (y2 wt x1^i y1^j)
-  ALLOCATE (fit[0].xsum, double *, fit[0].Nterms);
-  ALLOCATE (fit[0].ysum, double *, fit[0].Nterms);
-  for (i = 0; i < fit[0].Nterms; i++) {
-    ALLOCATE (fit[0].xsum[i], double, fit[0].Nterms);
-    ALLOCATE (fit[0].ysum[i], double, fit[0].Nterms);
-    bzero (fit[0].xsum[i], fit[0].Nterms*sizeof(double));
-    bzero (fit[0].ysum[i], fit[0].Nterms*sizeof(double));
-  }      
+  fit[0].xsum = array_init (fit[0].Nterms, fit[0].Nterms);
+  fit[0].ysum = array_init (fit[0].Nterms, fit[0].Nterms);
+
+  // xfit[i][j] holds x2 coeff for x1^i y1^j
+  // yfit[i][j] holds y2 coeff for x1^i y1^j
+  fit[0].xfit = array_init (fit[0].Nterms, fit[0].Nterms);
+  fit[0].yfit = array_init (fit[0].Nterms, fit[0].Nterms);
 
   // sum[i][j] holds \sum (wt x1^i y1^j)
-  ALLOCATE (fit[0].sum, double *, fit[0].Nsums);
-  for (i = 0; i < fit[0].Nsums; i++) {
-    ALLOCATE (fit[0].sum[i], double, fit[0].Nsums);
-    bzero (fit[0].sum[i], fit[0].Nsums*sizeof(double));
-  }
-
-  // matrix, vector hold the final linear system
-  ALLOCATE (fit[0].matrix, double *, fit[0].Nelems);
-  ALLOCATE (fit[0].vector, double *, fit[0].Nelems);
-  for (i = 0; i < fit[0].Nelems; i++) {
-    ALLOCATE (fit[0].matrix[i], double, fit[0].Nelems);
-    ALLOCATE (fit[0].vector[i], double, fit[0].Nelems);
-    bzero (fit[0].matrix[i], fit[0].Nelems*sizeof(double));
-    bzero (fit[0].vector[i], 2*sizeof(double));
-  }
+  fit[0].sum = array_init (fit[0].Nsums, fit[0].Nsums);
+
   return (fit);
 }
 
+void fit_free (CoordFit *fit) {
+
+  array_free (fit[0].xfit, fit[0].Nterms);
+  array_free (fit[0].yfit, fit[0].Nterms);
+
+  array_free (fit[0].sum, fit[0].Nsums);
+  array_free (fit[0].xsum, fit[0].Nterms);
+  array_free (fit[0].ysum, fit[0].Nterms);
+  
+  free (fit);
+}
+  
 // XXX use implicit masks as below or explicit masks (with function to set?)
 // XXX eg, add a global mask to this file and 
@@ -84 +105 @@
 void fit_eval (CoordFit *fit) {
 
-  int i, j, n, m, M, N;
+  int i, j, ix, iy, jx, jy;
+  double **matrix, **vector;
 
   if (fit[0].Npts == 0) {
     fprintf (stderr, "warning: no valid pts\n");
   }
+
+  // matrix, vector hold the final linear system
+  matrix = array_init (fit[0].Nelems, fit[0].Nelems);
+  vector = array_init (fit[0].Nelems, 2);
 
   /* remap the xsum,ysum terms into the vector */
@@ -94 +120 @@
     ix = i % fit[0].Nterms;
     iy = i / fit[0].Nterms;
-    vector[i][0] = xsum[ix][iy];
-    vector[i][1] = ysum[ix][iy];
+    vector[i][0] = fit[0].xsum[ix][iy];
+    vector[i][1] = fit[0].ysum[ix][iy];
 
     for (j = 0; j < fit[0].Nelems; j++) {
       jx = j % fit[0].Nterms;
       jy = j / fit[0].Nterms;
-      matrix[i][j] = sum[ix+jx][iy+jy];
-    }
-  }
-
-  dgaussj (matrix, NPARS, vector, 2); 
+      matrix[i][j] = fit[0].sum[ix+jx][iy+jy];
+    }
+  }
+
+  dgaussj (matrix, fit[0].Nelems, vector, 2); 
 
   for (i = 0; i < fit[0].Nelems; i++) {
@@ -112 +138 @@
              ix, iy, vector[i][0], ix, iy, vector[i][1]);
   }     
+
+  /* remap the vector terms into xfit,yfit */
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    fit[0].xfit[ix][iy] = vector[i][0];
+    fit[0].yfit[ix][iy] = vector[i][1];
+  }
+
+  array_free (matrix, fit[0].Nelems);
+  array_free (vector, fit[0].Nelems);
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_dx (double **poly, int Nx, int Ny) {
+
+  int i, j, Nxout, Nyout;
+  double **dpoly;
+
+  Nxout = Nx - 1;
+  Nyout = Ny;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  dpoly = array_init (Nxout, Nyout);
+
+  for (i = 0; i < Nxout; i++) {
+    for (j = 0; j < Nyout; j++) {
+      dpoly[i][j] = poly[i+1][j] * (i+1);
+    }
+  }
+  return dpoly;
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_dy (double **poly, int Nx, int Ny) {
+
+  int i, j, Nxout, Nyout;
+  double **dpoly;
+
+  Nxout = Nx;
+  Nyout = Ny - 1;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  dpoly = array_init (Nxout, Nyout);
+
+  for (i = 0; i < Nxout; i++) {
+    for (j = 0; j < Nyout; j++) {
+      dpoly[i][j] = poly[i][j+1] * (j+1);
+    }
+  }
+  return dpoly;
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_copy (double **poly, int Nx, int Ny) {
+
+  int i, j;
+  double **out;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  out = array_init (Nx, Ny);
+
+  for (i = 0; i < Nx; i++) {
+    for (j = 0; j < Ny; j++) {
+      out[i][j] = poly[i][j+1];
+    }
+  }
+  return out;
+}
+
+double poly2d_eval (double **poly, int Nx, int Ny, double x, double y) {
+
+  int i, j;
+  double xsum, ysum, sum;
+
+  sum = 0;
+  xsum = ysum = 1.0;
+
+  for (i = 0; i < Nx; i++) {
+    ysum = xsum;
+    for (j = 0; j < Ny; j++) {
+      sum += poly[i][j] * ysum;
+      ysum *= y;
+    }
+    xsum *= x;
+  }
+  return (sum);
 }
 
@@ -119 +232 @@
 void fit_apply_coords (CoordFit *fit, Coords *coords) {
 
-  int i, j, Np, Nv, N;
   double c11, c12;
   double c21, c22;
-  double R;
+  double Xo, Yo, R;
+  CoordFit *modfit;
 
   /* I have L,M = fit(X,Y). set corresponding terms for coords */
@@ -129 +242 @@
   // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ...
 
-  coords[0].crpix1,crpix2 = CoordsGetCenter (fit);
-  modfit = CoordsSetCenter (fit, xo, yo);
-
-  // set crpix1,2
-  fit_to_coordterms (coords, modfit, 0, 0);
+  CoordsGetCenter (fit, 0.001, &Xo, &Yo);
+  coords[0].crpix1 = Xo;
+  coords[0].crpix2 = Yo;
+  
+  // resulting fit should have zero constant terms
+  modfit = CoordsSetCenter (fit, Xo, Yo);
 
   // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0)
-  fit_to_coordterms (coords, modfit, 0, 1);
-  fit_to_coordterms (coords, modfit, 1, 0);
+  coords[0].pc1_1 = modfit[0].xfit[1][0];
+  coords[0].pc1_2 = modfit[0].xfit[0][1];
+  coords[0].pc2_1 = modfit[0].yfit[1][0];
+  coords[0].pc2_2 = modfit[0].yfit[0][1];
 
   // set the polyterm elements 
-  fit_to_coordterms (coords, modfit, 0, 2);
-  fit_to_coordterms (coords, modfit, 1, 1);
-  fit_to_coordterms (coords, modfit, 2, 0);
-
-  fit_to_coordterms (coords, modfit, 0, 3);
-  fit_to_coordterms (coords, modfit, 1, 2);
-  fit_to_coordterms (coords, modfit, 2, 1);
-  fit_to_coordterms (coords, modfit, 3, 0);
+  coords[0].polyterms[0][0] = modfit[0].xfit[2][0];
+  coords[0].polyterms[1][0] = modfit[0].xfit[1][1];
+  coords[0].polyterms[2][0] = modfit[0].xfit[0][2];
+
+  coords[0].polyterms[0][1] = modfit[0].yfit[2][0];
+  coords[0].polyterms[1][1] = modfit[0].yfit[1][1];
+  coords[0].polyterms[2][1] = modfit[0].yfit[0][2];
+
+  // I need to validate Norder
+  coords[0].polyterms[3][0] = modfit[0].xfit[3][0];
+  coords[0].polyterms[4][0] = modfit[0].xfit[2][1];
+  coords[0].polyterms[5][0] = modfit[0].xfit[1][2];
+  coords[0].polyterms[6][0] = modfit[0].xfit[0][3];
+
+  coords[0].polyterms[3][1] = modfit[0].yfit[3][0];
+  coords[0].polyterms[4][1] = modfit[0].yfit[2][1];
+  coords[0].polyterms[5][1] = modfit[0].yfit[1][2];
+  coords[0].polyterms[6][1] = modfit[0].yfit[0][3];
 
   /* we do not modify crval1,2: these are kept at the default values */
@@ -164 +290 @@
   coords[0].pc2_2  = c22*R;
 
+  fit_free (modfit);
   /* keep the order and type from initial values */
 }
 
-void fit_free (CoordFit *fit) {
-
-  int i;
-
-  for (i = 0; i < fit[0].Nterms; i++) {
-    free (fit[0].xsum[i]);
-    free (fit[0].ysum[i]);
-  }      
-  free (fit[0].xsum);
-  free (fit[0].ysum);
-
-  for (i = 0; i < fit[0].Nsums; i++) {
-    free (fit[0].sum[i]);
-  }
-  free (fit[0].sum);
-
-  for (i = 0; i < fit[0].Nelems; i++) {
-    free (fit[0].matrix[i]);
-    free (fit[0].vector[i]);
-  }
-  free (fit[0].matrix);
-  free (fit[0].vector);
-
-  free (fit);
-}
-  

branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c

@@ -1 @@
-# include "mosastro.h"
+# include "relastro.h"
 
-fit_to_coordterms (Coords *coords, CoordFit *fit, int nx, int ny) {
+double factorial (int N) {
 
-  int nsum;
+  int i;
+  double F;
 
-  nsum = nx + ny;
-
-  switch (nsum) {
-    case 0:
-      xterm = &coords[0].crpix1;
-      yterm = &coords[0].crpix2;
-      return;
-    case 1:
-      if (ny == 0) {
-        xterm = &coords[0].pc1_1;
-        yterm = &coords[0].pc2_1;
-      } else {
-        xterm = &coords[0].pc1_2;
-        yterm = &coords[0].pc2_2;
-      }
-      return;
-    case 2:
-      xterm = &coords[0].polyterm[ny][0];
-      yterm = &coords[0].polyterm[ny][1];
-      return;
-    case 3:
-      xterm = &coords[0].polyterm[ny+3][0];
-      yterm = &coords[0].polyterm[ny+3][1];
-      return;
-    default:
-      fprintf (stderr, "programming error %s:%d\n", __func__, __line__);
-      exit (2);
+  F = 1;
+  for (i = N; i > 1; i--) {
+    F *= i;
   }
-  fprintf (stderr, "programming error %s:%d\n", __func__, __line__);
-  exit (2);
+  return F;
 }
 
-
-// XXXX use a separate structure for the map (x2,y2 = f,g(x1,y1)) ?
-
-CoordsGetCenter () {
+// XXX use a separate structure for the map (x2,y2 = f,g(x1,y1)) ?
 
 // given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y
 // for which L,M = 0,0. tol is the allowed error on x,y.
-psPlane *psPlaneTransformGetCenter (psPlaneTransform *trans, double tol)
-{
+int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo) {
 
-    // crpix1,2 = X,Y(crval1,2)
-    // start with linear solution for Xo,Yo:
-    double R  = (trans->x->coeff[1][0]*trans->y->coeff[0][1] - trans->x->coeff[0][1]*trans->y->coeff[1][0]);
-    double Xo = (trans->y->coeff[0][0]*trans->x->coeff[0][1] - trans->x->coeff[0][0]*trans->y->coeff[0][1])/R;
-    double Yo = (trans->x->coeff[0][0]*trans->y->coeff[1][0] - trans->y->coeff[0][0]*trans->x->coeff[1][0])/R;
+  int i, Nx, Ny;
+  double R, Xo, Yo, dPos;
+  double **XdX, **XdY, **YdX, **YdY, **alpha, **beta;
+  double **xfit, **yfit;
 
-    // iterate to actual solution: requires small non-linear terms
-    if (trans->x->nX > 1) {
-        psPolynomial2D *XdX = psPolynomial2D_dX(NULL, trans->x);
-        psPolynomial2D *XdY = psPolynomial2D_dY(NULL, trans->x);
+  xfit = fit[0].xfit;
+  yfit = fit[0].yfit;
 
-        psPolynomial2D *YdX = psPolynomial2D_dX(NULL, trans->y);
-        psPolynomial2D *YdY = psPolynomial2D_dY(NULL, trans->y);
+  // solve for Xo,Yo where L,M(Xo,Yo) = 0,0
+  // start with linear solution for Xo,Yo:
+  R  = (xfit[1][0]*yfit[0][1] - xfit[0][1]*yfit[1][0]);
+  Xo = (yfit[0][0]*xfit[0][1] - xfit[0][0]*yfit[0][1])/R;
+  Yo = (xfit[0][0]*yfit[1][0] - yfit[0][0]*xfit[1][0])/R;
 
-        psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32);
-        psVector *Beta = psVectorAlloc (2, PS_DATA_F32);
+  Nx = fit[0].Nterms;
+  Ny = fit[0].Nterms;
 
-        /* this loop uses the Newton-Raphson method to solve for Xo,Yo
-        * it needs the high order terms to be small 
-        * Xo,Yo are in pixels;
-        */
-        double dPos = tol + 1;
-        for (int i = 0; (dPos > tol) && (i < 20); i++) {
-            // NOTE: order for Alpha is: [y][x]
-            Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo);
-            Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo);
-            Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo);
-            Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);
+  // iterate to actual solution: requires small non-linear terms
+  if (fit[0].Norder > 1) {
+    XdX = poly2d_dx (xfit, Nx, Ny);
+    XdY = poly2d_dy (xfit, Nx, Ny);
+    YdX = poly2d_dx (yfit, Nx, Ny);
+    YdY = poly2d_dy (yfit, Nx, Ny);
 
-            Beta->data.F32[0] = psPolynomial2DEval (trans->x, Xo, Yo);
-            Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo);
+    alpha = array_init (2, 2);
+    beta  = array_init (2, 1);
 
-            psMatrixGJSolveF32 (Alpha, Beta);
+    /* this loop uses the Newton-Raphson method to solve for Xo,Yo
+     * it needs the high order terms to be small 
+     * Xo,Yo are in pixels;
+     */
+    dPos = tol + 1;
+    for (i = 0; (dPos > tol) && (i < 20); i++) {
+      // NOTE: order for alpha is: [y][x]
+      alpha[0][0] = poly2d_eval (XdX, Nx-1, Ny,   Xo, Yo);
+      alpha[1][0] = poly2d_eval (XdY, Nx,   Ny-1, Xo, Yo);
+      alpha[0][1] = poly2d_eval (YdX, Nx-1, Ny,   Xo, Yo);
+      alpha[1][1] = poly2d_eval (YdY, Nx,   Ny-1, Xo, Yo);
 
-            Xo -= Beta->data.F32[0];
-            Yo -= Beta->data.F32[1];
-            dPos = hypot(Beta->data.F32[0], Beta->data.F32[1]);
-        }
-        psFree (Alpha);
-        psFree (Beta);
-        psFree (XdX);
-        psFree (XdY);
-        psFree (YdX);
-        psFree (YdY);
+      beta[0][0] = poly2d_eval (xfit, Nx, Ny, Xo, Yo);
+      beta[1][0] = poly2d_eval (yfit, Nx, Ny, Xo, Yo);
+
+      dgaussj (alpha, 2, beta, 1);
+
+      Xo -= beta[0][0];
+      Yo -= beta[1][0];
+      dPos = hypot(beta[0][0], beta[1][0]);
     }
-    psPlane *center = psPlaneAlloc ();
-    center->x = Xo;
-    center->y = Yo;
-
-    return center;
+    array_free (alpha, 2);
+    array_free (beta, 2);
+    array_free (XdX, Nx - 1);
+    array_free (XdY, Nx);
+    array_free (YdX, Nx - 1);
+    array_free (YdY, Nx);
+  }
+  *xo = Xo;
+  *yo = Yo;
+  return TRUE;
 }
 
 // convert a transformation L(x,y) to L'(x-xo,y-yo)
-psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, double Xo, double Yo)
-{
+CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo) {
 
-    // validate fit order
+  int i, j, Nx, Ny;
+  double **xPx, **yPx, **xPy, **yPy, **tmp;
 
-    if (output == NULL) {
-        output = psPlaneTransformAlloc(input->x->nX, input->x->nY);
+  CoordFit *output;
+  output = fit_init (input->Norder);
+
+  /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
+   * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and 
+   * noting that the constant term in each is the coefficient in the case of L(x,y) and is the 
+   * value of L'(-xo,-yo) in the second case.
+   */
+
+  Nx = input->Nterms;
+  Ny = input->Nterms;
+
+  xPx = poly2d_copy (input->xfit, Nx, Ny);
+  yPx = poly2d_copy (input->yfit, Nx, Ny);
+
+  for (i = 0; i <= input->Nterms; i++) {
+    xPy = poly2d_copy (xPx, Nx, Ny);
+    yPy = poly2d_copy (yPx, Nx, Ny);
+    for (j = 0; j <= input->Nterms; j++) {
+      output->xfit[i][j] = poly2d_eval (xPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j);
+      output->yfit[i][j] = poly2d_eval (yPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j);
+
+      // take the next derivative wrt y, catch output (is NULL on last pass)
+      if (Ny > 0) {
+        tmp = poly2d_dy(xPy, Nx, Ny);
+        array_free (xPy, Nx);
+        xPy = tmp;
+
+        tmp = poly2d_dy(yPy, Nx, Ny);
+        array_free (yPy, Nx);
+        yPy = tmp;
+      } else {
+        array_free (xPy, Nx);
+        array_free (yPy, Nx);
+      }
+
+      Ny --;
     }
-
-    /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
-     * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and 
-     * noting that the constant term in each is the coefficient in the case of L(x,y) and is the 
-     * value of L'(-xo,-yo) in the second case.
-     */
-
-    psPolynomial2D *tmp;
-
-    psPolynomial2D *xPx = psPolynomial2DCopy (NULL, input->x);
-    psPolynomial2D *yPx = psPolynomial2DCopy (NULL, input->y);
-
-    for (int i = 0; i <= input->x->nX; i++) {
-        psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx);
-        psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
-        for (int j = 0; j <= input->x->nY; j++) {
-            output->x->mask[i][j] = input->x->mask[i][j];
-            output->y->mask[i][j] = input->y->mask[i][j];
-            output->x->coeff[i][j] = (output->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
-            output->y->coeff[i][j] = (output->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
-
-            // take the next derivative wrt y, catch output (is NULL on last pass)
-            tmp = psPolynomial2D_dY(NULL, xPy);
-            psFree (xPy);
-            xPy = tmp;
-            tmp = psPolynomial2D_dY(NULL, yPy);
-            psFree (yPy);
-            yPy = tmp;
-        }
-        // take the next derivative wrt x, catch output (is NULL on last pass)
-        tmp = psPolynomial2D_dX(NULL, xPx);
-        psFree (xPx);
-        xPx = tmp;
-        tmp = psPolynomial2D_dX(NULL, yPx);
-        psFree (yPx);
-        yPx = tmp;
+    Ny = input->Nterms;
+    // take the next derivative wrt x, catch output (is NULL on last pass)
+    if (Nx > 0) {
+      tmp = poly2d_dx(xPx, Nx, Ny);
+      array_free (xPx, Nx);
+      xPx = tmp;
+      tmp = poly2d_dx(yPx, Nx, Ny);
+      array_free (yPx, Nx);
+      yPx = tmp;
+      Nx --;
+    } else {
+      array_free (xPx, Nx);
+      array_free (yPx, Nx);
     }
-    return output;
-}
-
-int mkpolyterm (int n, int m) {
-  
-  int i, nt, N;
-  
-  N = 0;
-  nt = n + m;
-  for (i = 2; i < nt; i++) {
-    N += i + 1;
   }
-  N += m;
-  return (N);
+  return output;
 }
 
 /*
-  coords uses an odd, rigid sequence for the coefficients:
+  Coords uses a rigid sequence for the coefficients:
 
   x^0 y^0 : crpix1,2 (but note this is applied before higher order terms)

branches/dvo-mods-2007-02/Ohana/src/relastro/src/plot_scatter.c

@@ -32 +32 @@
         if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
 
-        Mrel = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M;
+        Mrel = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
         xlist[N] = Mrel;
         ylist[N] = PhotSys  (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel;

branches/dvo-mods-2007-02/Ohana/src/relastro/src/plotstuff.c

@@ -49 +49 @@
   if (Xgraph[N] == 0) return;
 
-  KiiPS (Xgraph[N], TRUE, KAPA_NEWPLOT, filename);
+  KiiPS (Xgraph[N], TRUE, KAPA_PS_NEWPLOT, filename, "default");
   return;
 }

branches/dvo-mods-2007-02/Ohana/src/relastro/src/relastro.c

@@ -3 +3 @@
 int main (int argc, char **argv) {
 
-  int i, status, Ncatalog;
+  int status, Ncatalog;
   Catalog *catalog;
   FITS_DB db;
@@ -43 +43 @@
 
   /* major modes */
-  if (UpdateObjects) {
+  if (DoUpdateObjects) {
     UpdateObjects (catalog, Ncatalog);
   }
-  if (UpdateSimple) {
-    UpdateSimple ();
+  if (DoUpdateSimple) {
+    UpdateSimple (catalog, Ncatalog);
   }
-  if (UpdateChips) {
-    UpdateChips ();
+  if (DoUpdateChips) {
+    UpdateChips (catalog, Ncatalog);
   }
-  if (UpdateMosaics) {
-    UpdateMosaics ();
+  if (DoUpdateMosaics) {
+    UpdateMosaic (catalog, Ncatalog);
   }