Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 41665)
+++ /trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 41666)
@@ -42,4 +42,5 @@
 $(SRC)/fixcols.$(ARCH).o	   \
 $(SRC)/fiximage.$(ARCH).o	   \
+$(SRC)/forcedphot.$(ARCH).o	   \
 $(SRC)/gauss.$(ARCH).o		   \
 $(SRC)/getvel.$(ARCH).o	           \
@@ -90,7 +91,8 @@
 $(SRC)/imfit-qgauss-psf.$(ARCH).o	   \
 $(SRC)/imfit-sgauss.$(ARCH).o	   \
-$(SRC)/imfit-sgauss-psf.$(ARCH).o	   \
+$(SRC)/imfit-sgauss-psf.$(ARCH).o  \
 $(SRC)/imfit-qfgauss.$(ARCH).o	   \
-$(SRC)/imfit-qrgauss.$(ARCH).o	   
+$(SRC)/imfit-qrgauss.$(ARCH).o	   \
+$(SRC)/imfit-trail.$(ARCH).o
 
 # dependancy rules for include files ########################
Index: /trunk/Ohana/src/opihi/cmd.astro/forcedphot.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/forcedphot.c	(revision 41666)
+++ /trunk/Ohana/src/opihi/cmd.astro/forcedphot.c	(revision 41666)
@@ -0,0 +1,146 @@
+# include "astro.h"
+
+float Nsigma = 4.0;
+float sigma = 4.0 / 2.35; // FWHM = 4 pixels, 1 arcsec on PS1 GPC1
+float calc_forcedphot (Buffer *buf, float x, float y);
+float norm = 2.0;
+
+int forcedphot (int argc, char **argv) {
+
+  int N = 0;
+  Buffer *buf = NULL;
+  Vector *output = NULL;
+  Vector *xvec = NULL;
+  Vector *yvec = NULL;
+
+  int VERBOSE = TRUE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    VERBOSE = FALSE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-quiet"))) {
+    VERBOSE = FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-nsigma"))) {
+    remove_argument (N, &argc, argv);
+    Nsigma  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  if ((N = get_argument (argc, argv, "-sigma"))) {
+    remove_argument (N, &argc, argv);
+    sigma  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  int FLUXNORMAL = TRUE;
+  if ((N = get_argument (argc, argv, "-unit-normal"))) {
+    FLUXNORMAL = FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
+  // XXX this should be the output name of a variable for 
+  char *outname = NULL;
+  if ((N = get_argument (argc, argv, "-output"))) {
+    remove_argument (N, &argc, argv);
+    outname = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  } else {
+    outname = strcreate ("fflux");
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: psfqf (buffer) x y [-sigma sigma] [-nsigma nsigma] [-output outname]\n");
+    gprint (GP_ERR, "  measures forced-photometry using defined PSF for supplied coordinate(s)\n");
+    gprint (GP_ERR, "  supplied coordinates may be scalar or vector(s)\n");
+    gprint (GP_ERR, "  output values are saved in $fflux or vector fflux, or supplied name\n");
+    goto escape;
+  }
+
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) goto escape;
+
+  if (FLUXNORMAL) {
+    norm = 2.0;
+  } else {
+    norm = 1.0 / (2.0 * M_PI * SQ(sigma));
+  }
+
+  // Scalar coordinate input:
+  double Xin = NAN;
+  double Yin = NAN;
+  if (SelectScalar (argv[2], &Xin)) {
+    if (!SelectScalar (argv[3], &Yin)) {
+      gprint (GP_ERR, "syntax error: mixed vector and scalars for (x y)\n");
+      goto escape;
+    }
+    float flux = calc_forcedphot (buf, Xin, Yin);
+    if (VERBOSE) gprint (GP_LOG, "forced flux at %10.6f %10.6f = %f\n", Xin, Yin, flux);
+    set_variable (outname, flux);
+    free (outname);
+    return TRUE;
+  }
+
+  // Vector coordinate input:
+  if ((xvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((yvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if (xvec[0].Nelements != yvec[0].Nelements) {
+    fprintf (stderr, "mis-matched vector lengths\n");
+    goto escape;
+  }
+  REQUIRE_VECTOR_FLT (xvec, FALSE); 
+  REQUIRE_VECTOR_FLT (yvec, FALSE); 
+  
+  if ((output = SelectVector (outname, ANYBUFFER, TRUE)) == NULL) goto escape;
+  
+  if (output) ResetVector (output, OPIHI_FLT, xvec->Nelements);
+
+  for (int i = 0; i < xvec->Nelements; i++) {
+    Xin = xvec->elements.Flt[i];
+    Yin = yvec->elements.Flt[i];
+    float flux = calc_forcedphot (buf, Xin, Yin);
+    output->elements.Flt[i] = flux;
+  }
+  if (VERBOSE) gprint (GP_LOG, "forced flux calculated for %d points (output is %s)\n", output->Nelements, output->name);
+  
+  free (outname);
+  return (TRUE);
+
+ escape:
+  free (outname);
+  return FALSE;
+}
+
+// XXX things to consider:
+// if source is totally off image, flux should be nan, not 0.0
+// if any pixel is nan, skip or nan the source?
+float calc_forcedphot (Buffer *buf, float x, float y) {
+
+  int Nx = buf->matrix.Naxis[0];
+  int Ny = buf->matrix.Naxis[1];
+
+  int Xmin = x - Nsigma*sigma;
+  int Xmax = x + Nsigma*sigma;
+  int Ymin = y - Nsigma*sigma;
+  int Ymax = y + Nsigma*sigma;
+
+  float *bvalue = (float *)buf->matrix.buffer;
+  float f1 = -0.5 / SQ(sigma);
+
+  float sum = 0.0;
+  for (int iy = Ymin; iy <= Ymax; iy++) {
+    if (iy < 0) continue;
+    if (iy >= Ny) continue;
+    for (int ix = Xmin; ix <= Xmax; ix++) {
+      if (ix < 0) continue;
+      if (ix >= Nx) continue;
+      float r = SQ(ix - x) + SQ(iy - y);
+      
+      float gvalue = bvalue[ix + iy*Nx]*exp(f1*r);
+      sum += gvalue;
+    }
+  }
+  sum *= norm;
+  return sum;
+}
Index: /trunk/Ohana/src/opihi/cmd.astro/imfit-trail.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/imfit-trail.c	(revision 41666)
+++ /trunk/Ohana/src/opihi/cmd.astro/imfit-trail.c	(revision 41666)
@@ -0,0 +1,133 @@
+# include "imfit.h"
+
+opihi_flt trailTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *);
+void  trailCL ();
+
+# define PAR_X      0
+# define PAR_Y      1
+# define PAR_THETA  2
+# define PAR_SIGMA  3
+# define PAR_LENGTH 4
+# define PAR_I0     5
+# define PAR_SKY    6
+
+void trail_setup (char *name) {
+
+  if (strcmp(name, "trail")) return;
+
+  fitfunc = trailTD;
+  imfit_cleanup = trailCL;
+  Npar = 7;
+  Nfpar = 0;
+
+  /* allocate free and fixed parameters */
+  ALLOCATE (par, opihi_flt, MAX (Npar, 1));
+  bzero (par, Npar*sizeof(opihi_flt));
+  ALLOCATE (fpar, opihi_flt, MAX (Nfpar, 1));
+  bzero (fpar, Nfpar*sizeof(opihi_flt));
+
+  par[PAR_X      ] = get_variable_default ("Xg",      0.0);
+  par[PAR_Y      ] = get_variable_default ("Yg",      0.0);
+  par[PAR_THETA  ] = get_variable_default ("Tg",      0.0);
+  par[PAR_SIGMA  ] = get_variable_default ("Wg",      2.0);
+  par[PAR_LENGTH ] = get_variable_default ("Lg",     10.0);
+  par[PAR_I0     ] = get_variable_default ("Zpk", 10000.0);
+  par[PAR_SKY    ] = get_variable_default ("Sg",      0.0);
+  sky = &par[PAR_SKY];
+}
+
+void trailCL () {
+  set_variable ("Xg",  par[PAR_X     ]);
+  set_variable ("Yg",  par[PAR_Y     ]);
+  set_variable ("Wg",  par[PAR_SIGMA ]);
+  set_variable ("Tg",  par[PAR_THETA ]);
+  set_variable ("Lg",  par[PAR_LENGTH]);
+  set_variable ("Zpk", par[PAR_I0    ]);
+  set_variable ("Sg",  par[PAR_SKY   ]);
+}
+
+/* trailed gaussian -- x, y, Sigma, Theta, Length, I, sky */
+opihi_flt trailTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  OHANA_UNUSED_PARAM(Npar);
+
+  opihi_flt X = x - par[PAR_X];
+  opihi_flt Y = y - par[PAR_Y];
+  
+  opihi_flt S2 = 2.0 * SQ(par[PAR_SIGMA]);
+
+  opihi_flt ST = sin(RAD_DEG*par[PAR_THETA]);
+  opihi_flt CT = cos(RAD_DEG*par[PAR_THETA]);
+
+  opihi_flt Zp = (X*CT + Y*ST + 0.5*par[PAR_LENGTH]) / sqrt(S2);
+  opihi_flt Zm = (X*CT + Y*ST - 0.5*par[PAR_LENGTH]) / sqrt(S2);
+
+  opihi_flt Ep = erf(Zp);
+  opihi_flt Em = erf(Zm);
+
+  opihi_flt Rxy = Y*CT - X*ST;
+  opihi_flt Gxy = exp(-Rxy*Rxy/S2);
+
+  opihi_flt Pxy = Gxy * (Ep - Em);
+  opihi_flt f = Pxy * par[PAR_I0] + par[PAR_SKY];
+
+  if (dpar != NULL) {
+    dpar[PAR_SKY]    = 1.0;
+    dpar[PAR_I0]     = Pxy;
+
+    float dGdR = -2.0 * Rxy * Gxy / S2; // -R Gxy / (2 Sigma^2)
+
+    // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo
+    float dRdX = +ST;
+    float dRdY = -CT;
+    float dRdT = -Y*ST - X*CT;
+
+    float dGdX = dGdR * dRdX;
+    float dGdY = dGdR * dRdY;
+    float dGdT = dGdR * dRdT;
+
+    // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo
+    float dZpdX = -CT / sqrt(S2);
+    float dZmdX = dZpdX; // float dZmdX = -CT / sqrt(S2); dZmdX = dZpdX
+
+    float dZpdY = -ST / sqrt(S2); // float dZmdY = -ST / sqrt(S2); dZmdY = dZpdY
+    float dZmdY = dZpdY;
+
+    float dZpdL = +0.5 / sqrt(S2);
+    float dZmdL = -0.5 / sqrt(S2);
+
+    float dZpdT = (-X*ST + Y*CT) / sqrt(S2);
+    float dZmdT = dZpdT; // dZpdT = dZmdT
+
+    float dEdZp = exp (-Zp*Zp) * M_2_SQRTPI;
+    float dEdZm = exp (-Zm*Zm) * M_2_SQRTPI;
+
+    float dEpdX = dEdZp * dZpdX;
+    float dEmdX = dEdZm * dZmdX;
+
+    float dEpdY = dEdZp * dZpdY;
+    float dEmdY = dEdZm * dZmdY;
+
+    float dEpdL = dEdZp * dZpdL;
+    float dEmdL = dEdZm * dZmdL;
+
+    float dEpdT = dEdZp * dZpdT;
+    float dEmdT = dEdZm * dZmdT;
+
+    float dPdX = dGdX * (Ep - Em) + Gxy * (dEpdX - dEmdX);
+    float dPdY = dGdY * (Ep - Em) + Gxy * (dEpdY - dEmdY);
+
+    // dGdL is 0.0 because dRdL is 0.0
+    float dPdL = Gxy * (dEpdL - dEmdL);
+
+    float dPdT = dGdT * (Ep - Em) + Gxy * (dEpdT - dEmdT);
+
+    dpar[PAR_X]      = par[PAR_I0] * dPdX;
+    dpar[PAR_Y]      = par[PAR_I0] * dPdY;
+
+    dpar[PAR_LENGTH] = par[PAR_I0] * dPdL;
+    dpar[PAR_THETA]  = par[PAR_I0] * dPdT;
+    dpar[PAR_SIGMA]  = 0;	// we don't actually allow this to vary, so we do not need to calculate it
+  }
+
+  return (f);
+}
Index: /trunk/Ohana/src/opihi/cmd.astro/imfit.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 41665)
+++ /trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 41666)
@@ -65,4 +65,5 @@
     qfgauss_setup (argv[N]);
     qrgauss_setup (argv[N]);
+    trail_setup (argv[N]);
     if (fitfunc == NULL) {
       gprint (GP_ERR, "unknown function %s\n", argv[N]);
@@ -124,8 +125,12 @@
   dchisq = ochisq;
   chisq  = ochisq;
-  for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
+
+//for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
+
+  for (i = 0; (i < 25); i++) {
     chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
     dchisq = ochisq - chisq;
     ochisq = chisq;
+    fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
   }  
   set_int_variable ("Niter",  i);
Index: /trunk/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 41665)
+++ /trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 41666)
@@ -27,4 +27,5 @@
 int fixcols                 PROTO((int, char **));
 int fixrows                 PROTO((int, char **));
+int forcedphot              PROTO((int, char **));
 int gauss                   PROTO((int, char **));
 int gaussfit                PROTO((int, char **));
@@ -100,4 +101,5 @@
   {1, "fixcols",     fixcols,      "fix bad columns by comparing with others"},
   {1, "fixrows",     fixrows,      "fix bad rows by comparing with others"},
+  {1, "forcedphot",  forcedphot,   "forced photometry on a star or set of stars, assuming gaussian profile"},
   {1, "gauss",       gauss,        "get statistics on a star, assuming gaussian profile"},
   {1, "getvel",      getvel,       "rotcurve to velocities"},