Index: /trunk/Ohana/src/opihi/dvo/pmeasure.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/pmeasure.c	(revision 10793)
+++ /trunk/Ohana/src/opihi/dvo/pmeasure.c	(revision 10794)
@@ -1,11 +1,9 @@
 # include "dvoshell.h"
-// XXX this function is inefficient if it needs to plot stars for many catalog files
-// it would be better to accumulate the vectors and occasionally send them to the 
-// screen if the number of points exceeds some limit (say 10k - 100k)
+# define NCHUNK 10000
 
 int pmeasure (int argc, char **argv) {
   
   FILE *f;
-  int i, j, k, m, N, status, VERBOSE;
+  int i, j, k, m, Narg, Npts, NPTS, status, VERBOSE;
   double Mz, Mr, mag;
   double Radius, Rmin, Rmax, R, D;
@@ -13,4 +11,5 @@
   unsigned FlagChoice, FlagClip, flags;
   int PhotcodeClip;
+  float *Xvec, *Yvec, *Zvec;
 
   SkyTable *sky;
@@ -18,5 +17,4 @@
   Catalog catalog;
   Graphdata graphmode;
-  Vector Xvec, Yvec, Zvec;
 
   if (!InitPhotcodes ()) return (FALSE);
@@ -30,45 +28,45 @@
 
   PhotcodeClip = -1;
-  if ((N = get_argument (argc, argv, "-p"))) {
-    remove_argument (N, &argc, argv);
-    PhotcodeClip = GetPhotcodeCodebyName (argv[N]);
-    remove_argument (N, &argc, argv);
+  if ((Narg = get_argument (argc, argv, "-p"))) {
+    remove_argument (Narg, &argc, argv);
+    PhotcodeClip = GetPhotcodeCodebyName (argv[Narg]);
+    remove_argument (Narg, &argc, argv);
   }
   IDchoice = 0;
   IDclip = FALSE;
-  if ((N = get_argument (argc, argv, "-ID"))) {
-    remove_argument (N, &argc, argv);
-    IDchoice  = atoi(argv[N]);
-    remove_argument (N, &argc, argv);
+  if ((Narg = get_argument (argc, argv, "-ID"))) {
+    remove_argument (Narg, &argc, argv);
+    IDchoice  = atoi(argv[Narg]);
+    remove_argument (Narg, &argc, argv);
     IDclip = TRUE;
   }
   FlagChoice = 0;
   FlagClip = FALSE;
-  if ((N = get_argument (argc, argv, "-flag"))) {
-    remove_argument (N, &argc, argv);
-    FlagChoice  = atoi(argv[N]);
-    remove_argument (N, &argc, argv);
+  if ((Narg = get_argument (argc, argv, "-flag"))) {
+    remove_argument (Narg, &argc, argv);
+    FlagChoice  = atoi(argv[Narg]);
+    remove_argument (Narg, &argc, argv);
     FlagClip = TRUE;
   }
 
   LimExclude = FALSE;
-  if ((N = get_argument (argc, argv, "-x"))) {
-    remove_argument (N, &argc, argv);
+  if ((Narg = get_argument (argc, argv, "-x"))) {
+    remove_argument (Narg, &argc, argv);
     LimExclude = TRUE;
   }
 
   VERBOSE = FALSE;
-  if ((N = get_argument (argc, argv, "-v"))) {
-    remove_argument (N, &argc, argv);
+  if ((Narg = get_argument (argc, argv, "-v"))) {
+    remove_argument (Narg, &argc, argv);
     VERBOSE = TRUE;
   }
 
-  if ((N = get_argument (argc, argv, "-m"))) {
-    remove_argument (N, &argc, argv);
-    Mr  = atof(argv[N]);
-    remove_argument (N, &argc, argv);
-    Mz = atof(argv[N]);
+  if ((Narg = get_argument (argc, argv, "-m"))) {
+    remove_argument (Narg, &argc, argv);
+    Mr  = atof(argv[Narg]);
+    remove_argument (Narg, &argc, argv);
+    Mz = atof(argv[Narg]);
     Mr = Mr - Mz;
-    remove_argument (N, &argc, argv);
+    remove_argument (Narg, &argc, argv);
   }
 
@@ -78,4 +76,7 @@
   }
   SetGraph (graphmode);
+  graphmode.style = 2; /* set style to points */
+  graphmode.size = -1; /* point size determined by Zvec */
+  graphmode.etype = 0; /* no errorbars */
 
   Radius = MAX (fabs(graphmode.xmax), fabs(graphmode.ymax));
@@ -85,4 +86,11 @@
   skylist = SkyListByRadius (sky, -1, graphmode.coords.crval1, graphmode.coords.crval2, Radius);
   
+  /* storage for plotting the points */
+  Npts = 0;
+  NPTS = 1000;
+  ALLOCATE (Xvec, float, NPTS);
+  ALLOCATE (Yvec, float, NPTS);
+  ALLOCATE (Zvec, float, NPTS);
+
   for (j = 0; j < skylist[0].Nregions; j++) {
     catalog.filename = skylist[0].filename[j];
@@ -97,15 +105,6 @@
     dvo_catalog_unlock (&catalog);
 
-    /* data has been loaded, get ready to plot it */
-    Xvec.Nelements = catalog.Nmeasure;
-    Yvec.Nelements = catalog.Nmeasure;
-    Zvec.Nelements = catalog.Nmeasure;
-    ALLOCATE (Xvec.elements, float, Xvec.Nelements);
-    ALLOCATE (Yvec.elements, float, Yvec.Nelements);
-    ALLOCATE (Zvec.elements, float, Zvec.Nelements);
-
     /* project stars to screen display coords */
-    Xvec.Nelements = 0;
-    for (N = i = 0; i < catalog.Naverage; i++) {
+    for (i = 0; i < catalog.Naverage; i++) {
       if (IDclip && (catalog.average[i].code != IDchoice)) continue;
       while (catalog.average[i].R < Rmin) catalog.average[i].R += 360.0;
@@ -119,7 +118,7 @@
 	if ((PhotcodeClip != -1) && (catalog.measure[m+k].source != PhotcodeClip)) continue;
 	mag = PhotCat (&catalog.measure[m+k]);
-	Zvec.elements[N] = MIN (1.0, MAX (0.01, (mag - Mz) / Mr));
-	if (LimExclude && (Zvec.elements[N] > 0.99)) continue;
-	if (Zvec.elements[N] < 0.011) continue;
+	Zvec[Npts] = MIN (1.0, MAX (0.01, (mag - Mz) / Mr));
+	if (LimExclude && (Zvec[Npts] > 0.99)) continue;
+	if (Zvec[Npts] < 0.011) continue;
 	R = catalog.average[i].R - catalog.measure[m+k].dR_PS/3600.0;
 	D = catalog.average[i].D - catalog.measure[m+k].dD_PS/3600.0;
@@ -131,29 +130,31 @@
 	  free (date);
 	}
-	status = fRD_to_XY (&Xvec.elements[N], &Yvec.elements[N], R, D, &graphmode.coords);
+	status = fRD_to_XY (&Xvec[Npts], &Yvec[Npts], R, D, &graphmode.coords);
 	if (!status) {
 	  fprintf (stderr, ".");
 	  continue;
 	}
-	N ++;
+	Npts ++;
+
+	if (Npts == NPTS - 1) {
+	    NPTS += 1000;
+	    REALLOCATE (Xvec, float, NPTS);
+	    REALLOCATE (Yvec, float, NPTS);
+	    REALLOCATE (Zvec, float, NPTS);
+	}
+	if (Npts > NCHUNK) {
+	    PlotVectorTriplet (Npts, Xvec, Yvec, Zvec, &graphmode);
+	    Npts = 0;
+	}
       }
     }
-    Zvec.Nelements = Yvec.Nelements = Xvec.Nelements = N;
-
-    graphmode.style = 2; /* set style to points */
-    graphmode.size = -1; /* point size determined by Zvec */
-    graphmode.etype = 0; /* no errorbars */
-    PrepPlotting (N, &graphmode);
-    
-    PlotVector (N, Xvec.elements);
-    PlotVector (N, Yvec.elements);
-    PlotVector (N, Zvec.elements);
-
-    free (Xvec.elements);
-    free (Yvec.elements);
-    free (Zvec.elements);
-
     dvo_catalog_free (&catalog);
   }
+  if (Npts > 0) PlotVectorTriplet (Npts, Xvec, Yvec, Zvec, &graphmode);
+
+  free (Xvec);
+  free (Yvec);
+  free (Zvec);
+
   return (TRUE);