--- grads-2.0.a2/src/gagx.c	2008-05-02 15:21:50.000000000 +0200
+++ work/src/gagx.c	2008-08-06 17:14:07.000000000 +0200
@@ -2914,20 +2914,24 @@ char *umask, *vmask, *cmask=NULL;
     flag = 1;
     pgrc = pcm->result[2].pgr;
     if (pgrc->idim!=pgru->idim || pgrc->jdim!=pgru->jdim ||
         gagchk(pgrc,pgru,pgru->idim) || gagchk(pgrc,pgru,pgru->jdim)) {
       flag = 0;
-      gaprnt (0,"Error plotting streamlines:  Invalid color grid");
-      gaprnt (0,"   Color grid ignored -- has different scaling");
+      gaprnt (0,"Error plotting streamlines:  Invalid color grid\n");
+      gaprnt (0,"   Color grid ignored -- has different scaling\n");
     }
   }
 
   if ( (pcm->rotate && (pgru->idim!=2 || pgru->jdim!=3)) ||
       (!pcm->rotate && pgru->idim==2 && pgru->jdim==3)) {
-    pgru = gaflip(pgru,pcm);
-    pgrv = gaflip(pgrv,pcm);
-    if (flag) pgrc = gaflip(pgrc,pcm);
+    if (pgru->idim==0 && pgru->jdim==1) {
+      gaprnt (0,"Warning:  xyrev not supported with streamlines\n");
+    } else {
+      pgru = gaflip(pgru,pcm);
+      pgrv = gaflip(pgrv,pcm);
+      if (flag) pgrc = gaflip(pgrc,pcm);
+    }
   }
 
   gxstyl (1);
   gxwide (1);
   gas2d (pcm, pgru, 1);     /* Set up scaling, draw map if apprprt */
@@ -2936,11 +2940,11 @@ char *umask, *vmask, *cmask=NULL;
   idiv = 1.0; jdiv = 1.0;
 
   if (flag) {
     gamnmx (pgrc);
     if (pgrc->umin==0) {
-      gaprnt (0,"Connot color vectors -- Color grid all undefined\n");
+      gaprnt (0,"Cannot color vectors -- Color grid all undefined\n");
       flag = 0;
     } else gaselc(pcm,pgrc->rmin,pgrc->rmax);
   }
 
   u = pgru->grid;
@@ -2957,17 +2961,13 @@ char *umask, *vmask, *cmask=NULL;
   gxcolr (lcol);
   gxstyl(pcm->cstyle);
   gxwide(pcm->cthick);
   gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
 
-  if (flag) {
-    gxstrm (u,v,c,pgru->isiz,pgru->jsiz,umask,vmask,cmask,
-       flag,pcm->shdlvs,pcm->shdcls,pcm->shdcnt,pcm->strmden);
-  } else {
-    gxstrm (u,v,NULL,pgru->isiz,pgru->jsiz,umask,vmask,0,
-       flag,pcm->shdlvs,pcm->shdcls,pcm->shdcnt,pcm->strmden);
-  }
+  gxstrm (u,v,c,pgru->isiz,pgru->jsiz,umask,vmask,cmask,
+     flag,pcm->shdlvs,pcm->shdcls,pcm->shdcnt,pcm->strmden,
+     pcm->ahdsiz,pgru,pcm->mproj,&mpj);
 
   gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
   gxwide (pcm->annthk-3);
   gxcolr (pcm->anncol);
   if (pcm->pass==0 && pcm->grdsflg)
--- grads-2.0.a2/src/gx.h	2008-01-10 17:42:05.000000000 +0100
+++ work/src/gx.h	2008-08-06 21:25:50.000000000 +0200
@@ -158,10 +158,23 @@ struct xinfo {
                                          -999 does double plot */
   gaint ccol;                  /* Contour color */
   char *label;                 /* Contour label */
 };
 
+/* Structure for passing information to gxstrm subroutines */
+
+struct strmparam {
+  gadouble (*igrab) (gadouble *, gadouble), (*jgrab) (gadouble *, gadouble);
+  gadouble *u, *v, *c, *shdlvs, *ivals, *jvals, *dxt, *dyt, *xfact;
+  gaint *shdcls;
+  char *umask, *vmask, *cmask, *it;
+  gadouble dis, rscl, asiz;
+  gaint shdcnt, is, js, iss, jss, imin, jmin, cflag, isgeo, wrap;
+};
+
+struct gagrid;
+
 /* Function prototypes for GX library routines  */
 
 /* Functions in gxdev:
    gxdbgn: Initialize hardware
    gxdend: Terminate hardware
@@ -415,12 +428,14 @@ void  shdcmp (void);
 gaint shdmax (void);
 
 /* routines in gxstrm:  gxstrm (do streamlines) */
 
 void gxstrm (gadouble *, gadouble *, gadouble *, gaint, gaint, char *, char *, char *,
-	     gaint, gadouble *, gaint *, gaint, gaint);
-void strmar (gadouble, gadouble, gadouble, gadouble);
+	     gaint, gadouble *, gaint *, gaint, gaint,
+	     gadouble, struct gagrid *, gaint, struct mapprj *);
+gaint strmhalf (gadouble, gadouble, gaint, gaint, const struct strmparam *);
+void strmar (gadouble, gadouble, gadouble, gadouble, gadouble);
 gaint gxshdc (gadouble *, gaint *, gaint, gadouble);
 
 /* Routines in gxwmap:
    gxwmap: Draw world map
    gxnmap: Draw medium res n.am. map
--- grads-2.0.a2/src/gxstrm.c	2008-01-10 17:42:05.000000000 +0100
+++ work/src/gxstrm.c	2008-08-06 22:31:16.000000000 +0200
@@ -1,10 +1,10 @@
 /*  Copyright (C) 1988-2008 by Brian Doty and the 
     Institute of Global Environment and Society (IGES).  
     See file COPYRIGHT for more information.   */
 
-/* Authored by B. Doty */
+/* Authored by B. Doty, rewritten by Andreas Schneider */
 
 #ifdef HAVE_CONFIG_H
 #include "config.h"
 
 /* If autoconfed, only include malloc.h when it's presen */
@@ -18,253 +18,632 @@
 
 #endif /* HAVE_CONFIG_H */
 
 #include <stdio.h>
 #include <math.h>
-#include "gatypes.h"
+#include <string.h>
+#include "grads.h"
 #include "gx.h"
 
-void gxstrm (gadouble *u, gadouble *v, gadouble *c, gaint is, gaint js,
-   char *umask, char *vmask, char *cmask, gaint flag, gadouble *shdlvs,
-   gaint *shdcls, gaint shdcnt, gaint den) {
-gadouble x,y,xx,yy;
-gadouble *up, *vp, *cp, cv1,cv2,cv;
-gadouble uv1,uv2,uv,vv1,vv2,vv,auv,avv,xsav,ysav,xold=0.0,yold=0.0;
-gadouble fact,rscl,xxsv,yysv;
-gaint i,ii,jj,ii1,ij1,i2,j2,ipt,acnt,icol,scol,dis;
-gaint *it,siz,iacc,iisav,iscl,imn,imx,jmn,jmx,iz,jz,iss,jss,bflg;
-char *upmask,*vpmask,*cpmask;
-  scol = -9;
-  icol = 1;
+#define DEFDEN 5      /* input strmden value for default density */
+#define DENSCL 1.25992105 /* density factor for strmden+1: pow(2.0,1.0/3.0) */
+#define STRMDEN 270   /* desired number of flag boxes at default density along \
+                         X plus Y axis (270 means 2deg grid for a global map) */
+#define DEFDIS 1.5    /* desired minimum distance in flag box units */
+#define STARTFACT 1.6 /* factor for required minimum distance to start a line */
+#define STEPSPERBOX 4 /* number of steps per flag box width or height */
+#define UNDEFCOL 15   /* color for undefined values */
+#define ARRDIST 0.8   /* distance between arrows in output units (inches) */
+#define ARRMINLEN 0.2 /* approximate minimum length to plot an arrow */
+#define DEG2RAD (pi/180.0)
+#define MASK ((galint)0x0101010101010101LL)
+
+#define DEBUG 1
+
+/* We set up a grid of flag boxes to mark where we already have a streamline.
+   Grid points lie at the corners of flag boxes.  The area between 4 adjacent
+   grid points may be divided into more than one flag box to have always a
+   similar number of boxes within the visible area of the grid; horizontal and
+   vertical scaling are the same.
+
+   Streamlines can continue until another streamline has crossed a flag box
+   whose center is less than a certain distance away.  We don't take different
+   horizontal and vertical resolution into account for that purpose, but we
+   account for denser meridians at higher latitudes for map projections where
+   this is visible.
+
+   A streamline can get very close to itself for simpler calculation and to
+   allow for nearly closed streamlines.  It just can't touch the same flag box
+   again once it has left it.
+
+   Streamlines start from the center of flag boxes in both directions.  The
+   distance from other streamlines needs to be greater than that required to
+   continue, to avoid too much very short streamlines.  */
 
-  /* Figure out the interval for the flag grid */
-
-  i = is;
-  if (js>i) i = js;
-  iscl = 200/i;
-  iscl = iscl + den - 5;
-  if (iscl<1) iscl=1;
-  if (iscl>10) iscl=10;
-  fact = 0.5/((gadouble)iscl);
-  rscl = (gadouble)iscl;
- 
-  /* Allocate memory for the flag grid */
-
-  iss = is*iscl; jss = js*iscl;
-  siz = iss*jss;
-  it = (gaint *)malloc(sizeof(gaint) * siz);
-  if (it==NULL) {
+void gxstrm (gadouble *u, gadouble *v, gadouble *c, gaint is, gaint js,
+   char *umask, char *vmask, char *cmask, gaint cflag, gadouble *shdlvs,
+   gaint *shdcls, gaint shdcnt, gaint den, gadouble asiz, struct gagrid *pgr,
+   gaint mproj, struct mapprj *mpj) {
+gadouble x,y,xx,yy,wscl,sdis,dis2,min,xf;
+gaint iz,jz,i,j,ilo,ihi,jlo,jhi,imn,imx,jmn,jmx,di,dj,idis,jdis,ii;
+gaint iscl,steps,ilim,siz,lsiz,denscl,cnt,**pdis,*pd;
+struct strmparam param;
+char *jt0,*jt;
+galint *pt,*et;
+
+  param.u = u;
+  param.v = v;
+  param.c = c;
+  param.is = is;
+  param.js = js;
+  param.umask = umask;
+  param.vmask = vmask;
+  param.cmask = cmask;
+  param.cflag = cflag;
+  param.shdlvs = shdlvs;
+  param.shdcls = shdcls;
+  param.shdcnt = shdcnt;
+  param.asiz = asiz;
+  param.isgeo = 0;
+  param.wrap = 0;
+  denscl = 0;
+
+  if (mproj!=0 && pgr->idim==0 && pgr->jdim==1) { /* X/Y grid */
+    param.isgeo = 1;
+    param.imin = pgr->dimmin[0];
+    param.jmin = pgr->dimmin[1];
+    param.igrab = pgr->igrab;
+    param.jgrab = pgr->jgrab;
+    param.ivals = pgr->ivals;
+    param.jvals = pgr->jvals;
+    x = pgr->iabgr(pgr->iavals,mpj->lnmx) - pgr->iabgr(pgr->iavals,mpj->lnmn);
+    y = pgr->iabgr(pgr->iavals,mpj->ltmx) - pgr->iabgr(pgr->iavals,mpj->ltmn);
+    if (mproj>2) {
+      /* Account for denser grid near the poles */
+      denscl = 1;
+      if (mpj->ltmn>0.0) x *= cos(mpj->ltmn*DEG2RAD);
+      else if (mpj->ltmx<0.0) x *= cos(mpj->ltmx*DEG2RAD);
+      if (mproj==3 || mproj==4) { /* nps or sps */
+        if (param.igrab(param.ivals,pgr->dimmax[0])>=
+            param.igrab(param.ivals,param.imin)+359.999999) param.wrap = 1;
+        xx = mpj->lnmx - mpj->lnmn;
+        yy = mpj->ltmx - mpj->ltmn;
+        if (xx>=180.0) {
+          x *= (360.0/pi) / xx;
+          y *= (mproj==3 ? 90.0-mpj->ltmn : mpj->ltmx+90.0) *
+               (1.0-cos(xx*DEG2RAD/2.0)) / yy;
+        } else {
+          x *= sin(xx*DEG2RAD/2.0)*2.0 / (xx*DEG2RAD);
+          if (mproj==3)
+            y *= ((90.0-mpj->ltmn) - (90.0-mpj->ltmx)*cos(xx*DEG2RAD/2.0)) / yy;
+          else
+            y *= ((mpj->ltmx+90.0) - (mpj->ltmn+90.0)*cos(xx*DEG2RAD/2.0)) / yy;
+        }
+      }
+      /* XXX Tweek density for other other projection types too */
+    }
+  } else {
+    x = (gadouble)(is-1);
+    y = (gadouble)(js-1);
+  }
+  wscl = (STRMDEN * pow(DENSCL,den-DEFDEN)) / (x+y);
+  iscl = (gaint)(wscl + 0.5);
+  if (iscl<1) iscl = 1;
+  param.dis = DEFDIS * (gadouble)iscl/wscl;
+  sdis = param.dis*STARTFACT;
+  jdis = (gaint)sdis;
+  param.rscl = (gadouble)iscl;
+  steps = iscl*STEPSPERBOX;
+  param.iss = (is-1)*iscl;
+  param.jss = (js-1)*iscl;
+  ilim = param.iss-1;
+  siz = param.iss*param.jss;
+  lsiz = (siz+sizeof(galint)-1) / sizeof(galint);
+  param.it = (char *)malloc(sizeof(galint) * lsiz +
+                            sizeof(gadouble) * (is+js-2+param.jss) +
+                            sizeof(gaint *) * param.jss +
+                            sizeof(gaint) * param.jss*(jdis*2+1));
+  if (param.it==NULL) {
     printf ("Cannot allocate memory for streamline function\n");
     return;
   }
-  for (i=0; i<siz; i++) *(it+i) = 0;
+#ifdef DEBUG
+  printf("Streamline grid size: %dx%d, visible extent %.3fx%.3f\n",is,js,x,y);
+  printf("Streamline flag grid: scaling %d (wanted %.3f),",iscl,wscl);
+  printf(" size %dx%d (%.3f MiB)\n",param.iss,param.jss,
+         (gadouble)(lsiz*sizeof(galint))/1048576.0);
+  printf("Streamline distance:  start %.3f, continue %.3f (flag grid units)\n",
+         sdis,param.dis);
+#endif
+  memset(param.it,0,sizeof(galint)*lsiz);
+  param.dxt = (gadouble *)((galint *)param.it + lsiz);
+  param.dyt = param.dxt + is-1;
+  param.xfact = param.dyt + js-1;
+  pdis = (gaint **)(param.xfact + param.jss);
+  pd = (gaint *)(pdis + param.jss);
+  if (param.isgeo) {
+    xx = param.igrab(param.ivals,param.imin);
+    for (iz=0; iz<is-1; iz++) {
+      x = param.igrab(param.ivals,param.imin+iz+1);
+      param.dxt[iz] = x-xx;
+      xx = x;
+    }
+    yy = param.jgrab(param.jvals,param.jmin);
+    for (jz=0; jz<js-1; jz++) {
+      y = param.jgrab(param.jvals,param.jmin+jz+1);
+      param.dyt[jz] = y-yy;
+      yy = y;
+    }
+  } else {
+    for (iz=0; iz<is-1; iz++)
+      param.dxt[iz] = 1.0;
+    for (jz=0; jz<js-1; jz++)
+      param.dyt[jz] = 1.0;
+  }
+  if (param.isgeo && denscl) {
+    min = sdis/param.iss;
+    for (jz=0; jz<param.jss; jz++) {
+      xf = cos(param.jgrab(param.jvals,param.jmin+(jz+0.5)/param.rscl)*DEG2RAD);
+      if (xf<min) xf = min;
+      param.xfact[jz] = 1.0/xf;
+    }
+  } else {
+    for (jz=0; jz<param.jss; jz++)
+      param.xfact[jz] = 1.0;
+  }
+  dis2 = sdis*sdis;
+  for (jz=0; jz<param.jss; jz++) {
+    pdis[jz] = pd + jdis;
+    xf = param.xfact[jz];
+    for (dj=-jdis; dj<=jdis; dj++)
+      *pd++ = (gaint)(sqrt(dis2 - (gadouble)(dj*dj))*xf);
+  }
 
-  /* Loop through flag grid to look for start of streamlines.  
-     To start requires no streams drawn within surrounding 
-     flag boxes.  */
-
-  i2 = 0;
-  j2 = 0;
-  for (i=0; i<siz; i++) {
-    dis = 2;
-    if (den<5) dis = 3;
-    if (den>5) dis = 1;
-    imn = i2-dis; imx = i2+dis+1; 
-    jmn = j2-dis; jmx = j2+dis+1;
-    if (imn<0) imn = 0;
-    if (imx>iss) imx = iss;
+  /* Loop through flag grid to look for start of streamlines.
+     We start in the middle of the grid to avoid spurious short lines at the
+     edges that can just start there because the previous streamline is cut off
+     at the edge.  */
+
+  dj = 0;
+  j = jlo = jhi = param.jss>>1;
+  if (param.jss&1) goto start2;
+  else goto start1;
+  do {
+    j += ++dj;
+  start1:
+    pd = pdis[j];
+    jmn = j-jdis;
+    jmx = j+jdis;
     if (jmn<0) jmn = 0;
-    if (jmx>jss) jmx = jss;
-    iacc = 0;
-    for (jz=jmn; jz<jmx; jz++) {
-      ipt = jz*iss+imn;
-      for (iz=imn; iz<imx; iz++) {
-        iacc = iacc + *(it+ipt);
-        ipt++;
+    if (jmx>=param.jss) jmx = param.jss-1;
+    jt0 = param.it + jmn*param.iss;
+    di = 0;
+    i = ilim>>1;
+    ilo = i-1;
+    ihi = i+1;
+    if (param.iss&1) goto start_hi1;
+    else goto start_lo1;
+    do {
+      if (i==ilo) {
+        ilo--;
+      start_lo1:
+        jt = jt0;
+        for (jz=jmn; jz<=jmx; jz++) {
+          idis = pd[jz-j];
+          imn = i-idis;
+          imx = i+idis;
+          if (imx>=ilim) {
+            if (param.wrap) {
+              imn -= ilim;
+              imx -= ilim;
+            } else imx = ilim;
+          }
+          for (iz=imx; iz>=imn; ) {
+            if (*(jt+iz)) goto cont_hi1;
+            if (iz==0) {
+              if (param.wrap) {
+                iz = ilim;
+                imn += ilim;
+              } else break;
+            } else iz--;
+          }
+          jt += param.iss;
+        }
+        x = (i+0.5)/param.rscl;
+        y = (j+0.5)/param.rscl;
+        ii = j*param.iss + i;
+        cnt = strmhalf(x,y,ii,steps,&param);
+        *(param.it+ii) = 0;
+        cnt += strmhalf(x,y,ii,-steps,&param);
+        if (cnt) {
+          pt = (galint *)param.it;
+          et = pt + lsiz;
+          do *pt++ &= MASK;
+          while (pt!=et);
+          ilo -= pd[0];
+          if (!di) ihi += pd[0];
+        }
       }
-    }
-    if (iacc==0){
-      x = ((gadouble)i2)/rscl;
-      y = ((gadouble)j2)/rscl;
-      xsav = x;
-      ysav = y;
-      gxconv (x+1.0,y+1.0,&xx,&yy,3);
-      gxplot (xx,yy,3);
-      xxsv = xx; yysv = yy;
-      iisav = -999;
-      iacc = 0;
-      acnt = 0;
-      bflg = 0;
-      while (x>=0.0 && x<(gadouble)(is-1) && y>=0.0 && y<(gadouble)(js-1)) {
-        ii = (gaint)x;
-        jj = (gaint)y;
-        xx = x - (gadouble)ii;
-        yy = y - (gadouble)jj;
-        up = u + jj*is+ii; upmask = umask + jj*is+ii;  
-        vp = v + jj*is+ii; vpmask = vmask + jj*is+ii; 
-        if (*upmask==0 || 
-	    *(upmask+1)==0 ||
-            *(upmask+is)==0 || 
-	    *(upmask+is+1)==0) break;
-        if (*vpmask==0 || 
-	    *(vpmask+1)==0 ||
-            *(vpmask+is)==0 || 
-	    *(vpmask+is+1)==0) break;
-        if (flag) {
-          cp = c + jj*is+ii; cpmask = cmask + jj*is+ii;
-          if (*cpmask==0 || 
-	      *(cpmask+1)==0 ||
-              *(cpmask+is)==0 || 
-	      *(cpmask+is+1)==0) icol = 15;
-          else {
-            cv1 = *cp + (*(cp+1)-*cp)*xx;
-            cv2 = *(cp+is) + (*(cp+is+1)-*(cp+is))*xx;
-            cv = cv1 + (cv2-cv1)*yy;
-            icol = gxshdc(shdlvs,shdcls,shdcnt,cv);
+    cont_hi1:
+      i += ++di;
+      if (i==ihi) {
+        ihi++;
+      start_hi1:
+        jt = jt0;
+        for (jz=jmn; jz<=jmx; jz++) {
+          idis = pd[jz-j];
+          imn = i-idis;
+          imx = i+idis;
+          if (imn<=0) {
+            if (param.wrap) {
+              imn += ilim;
+              imx += ilim;
+            } else imn = 0;
           }
-          if (icol!=scol && icol>-1) gxcolr(icol);
-          scol = icol;
+          for (iz=imn; iz<=imx; ) {
+            if (*(jt+iz)) goto cont_lo1;
+            if (iz==ilim) {
+              if (param.wrap) {
+                iz = 0;
+                imx -= ilim;
+              } else break;
+            } else iz++;
+          }
+          jt += param.iss;
         }
-        uv1 = *up + (*(up+1)-*up)*xx;
-        uv2 = *(up+is) + (*(up+is+1)-*(up+is))*xx;
-        uv = uv1 + (uv2-uv1)*yy;
-        vv1 = *vp + (*(vp+1)-*vp)*xx;
-        vv2 = *(vp+is) + (*(vp+is+1)-*(vp+is))*xx;
-        vv = vv1 + (vv2-vv1)*yy;
-        auv = fabs(uv); avv=fabs(vv);
-        if (auv<0.1 && avv<0.1) break;
-        if (auv>avv) {
-          vv = vv*fact/auv;
-          uv = uv*fact/auv;
-        } else {
-          uv = uv*fact/avv;
-          vv = vv*fact/avv;
+        x = (i+0.5)/param.rscl;
+        y = (j+0.5)/param.rscl;
+        ii = j*param.iss + i;
+        cnt = strmhalf(x,y,ii,steps,&param);
+        *(param.it+ii) = 0;
+        cnt += strmhalf(x,y,ii,-steps,&param);
+        if (cnt) {
+          pt = (galint *)param.it;
+          et = pt + lsiz;
+          do *pt++ &= MASK;
+          while (pt!=et);
+          ihi += pd[0];
+          if (!di) ilo -= pd[0];
         }
-        x = x + uv;
-        y = y + vv;
-        ii1 = (gaint)(x*rscl);
-        ij1 = (gaint)(y*rscl);
-        ii1 = ij1*iss + ii1;
-        if (ii1<0 || ii1>=siz) break;
-        if (*(it+ii1)==1) break;
-        if (ii1!=iisav && iisav>-1) *(it+iisav) = 1;
-        if (ii1==iisav) iacc++;
-        else iacc = 0;
-        if (iacc>10) break;
-        iisav = ii1;
-        gxconv (x+1.0,y+1.0,&xx,&yy,3);
-        if (icol>-1) {
-          if (bflg) {gxplot(xold,yold,3); bflg=0;}
-          gxplot (xx,yy,2);
-        } else bflg = 1;
-        xold = xx;
-        yold = yy;
-        acnt++;
-        if (acnt>20) {
-          if (icol>-1) strmar (xxsv,yysv,xx,yy);
-          acnt = 0;
-        }
-        xxsv = xx; yysv = yy;
-      }
-      bflg = 0;
-      x = xsav; y = ysav;
-      gxconv (x+1.0,y+1.0,&xx,&yy,3);
-      gxplot (xx,yy,3);
-      xxsv = xx;
-      yysv = yy;
-      iisav = -999;
-      iacc = 0;
-      acnt = 19;
-      while (x>=0.0 && x<(gadouble)(is-1) && y>=0.0 && y<(gadouble)(js-1)) {
-        ii = (gaint)x;
-        jj = (gaint)y;
-        xx = x - (gadouble)ii;
-        yy = y - (gadouble)jj;
-        up = u + jj*is+ii; upmask = umask + jj*is+ii;  
-        vp = v + jj*is+ii; vpmask = vmask + jj*is+ii; 
-        if (*upmask==0 || 
-	    *(upmask+1)==0 ||
-            *(upmask+is)==0 || 
-	    *(upmask+is+1)==0) break;
-        if (*vpmask==0 || 
-	    *(vpmask+1)==0 ||
-            *(vpmask+is)==0 || 
-	    *(vpmask+is+1)==0) break;
-        if (flag) {
-          cp = c + jj*is+ii; cpmask = cmask + jj*is+ii; 
-          if (*cpmask==0 || 
-	      *(cpmask+1)==0 ||
-              *(cpmask+is)==0 || 
-	      *(cpmask+is+1)==0) icol = 15;
-          else {
-            cv1 = *cp + (*(cp+1)-*cp)*xx;
-            cv2 = *(cp+is) + (*(cp+is+1)-*(cp+is))*xx;
-            cv = cv1 + (cv2-cv1)*yy;
-            icol = gxshdc(shdlvs,shdcls,shdcnt,cv);
+      }
+    cont_lo1:
+      i -= ++di;
+    } while (i>=0);
+    j -= ++dj;
+  start2:
+    pd = pdis[j];
+    jmn = j-jdis;
+    jmx = j+jdis;
+    if (jmn<0) jmn = 0;
+    if (jmx>=param.jss) jmx = param.jss-1;
+    jt0 = param.it + jmx*param.iss;
+    di = 0;
+    i = ilim>>1;
+    ilo = i-1;
+    ihi = i+1;
+    if (param.iss&1) goto start_hi2;
+    else goto start_lo2;
+    do {
+      if (i==ilo) {
+        ilo--;
+      start_lo2:
+        jt = jt0;
+        for (jz=jmx; jz>=jmn; jz--) {
+          idis = pd[jz-j];
+          imn = i-idis;
+          imx = i+idis;
+          if (imx>=ilim) {
+            if (param.wrap) {
+              imn -= ilim;
+              imx -= ilim;
+            } else imx = ilim;
+          }
+          for (iz=imx; iz>=imn; ) {
+            if (*(jt+iz)) goto cont_hi2;
+            if (iz==0) {
+              if (param.wrap) {
+                iz = ilim;
+                imn += ilim;
+              } else break;
+            } else iz--;
           }
-          if (icol!=scol && icol>-1) gxcolr(icol);
-          scol = icol;
+          jt -= param.iss;
         }
-        uv1 = *up + (*(up+1)-*up)*xx;
-        uv2 = *(up+is) + (*(up+is+1)-*(up+is))*xx;
-        uv = uv1 + (uv2-uv1)*yy;
-        vv1 = *vp + (*(vp+1)-*vp)*xx;
-        vv2 = *(vp+is) + (*(vp+is+1)-*(vp+is))*xx;
-        vv = vv1 + (vv2-vv1)*yy;
-        auv = fabs(uv); avv=fabs(vv);
-        if (auv<0.1 && avv<0.1) break;
-        if (auv>avv) {
-          vv = vv*fact/auv;
-          uv = uv*fact/auv;
-        } else {
-          uv = uv*fact/avv;
-          vv = vv*fact/avv;
+        x = (i+0.5)/param.rscl;
+        y = (j+0.5)/param.rscl;
+        ii = j*param.iss + i;
+        cnt = strmhalf(x,y,ii,steps,&param);
+        *(param.it+ii) = 0;
+        cnt += strmhalf(x,y,ii,-steps,&param);
+        if (cnt) {
+          pt = (galint *)param.it;
+          et = pt + lsiz;
+          do *pt++ &= MASK;
+          while (pt!=et);
+          ilo -= pd[0];
+          if (!di) ihi += pd[0];
+        }
+      }
+    cont_hi2:
+      i += ++di;
+      if (i==ihi) {
+        ihi++;
+      start_hi2:
+        jt = jt0;
+        for (jz=jmx; jz>=jmn; jz--) {
+          idis = pd[jz-j];
+          imn = i-idis;
+          imx = i+idis;
+          if (imn<=0) {
+            if (param.wrap) {
+              imn += ilim;
+              imx += ilim;
+            } else imn = 0;
+          }
+          for (iz=imn; iz<=imx; ) {
+            if (*(jt+iz)) goto cont_lo2;
+            if (iz==ilim) {
+              if (param.wrap) {
+                iz = 0;
+                imx -= ilim;
+              } else break;
+            } else iz++;
+          }
+          jt -= param.iss;
+        }
+        x = (i+0.5)/param.rscl;
+        y = (j+0.5)/param.rscl;
+        ii = j*param.iss + i;
+        cnt = strmhalf(x,y,ii,steps,&param);
+        *(param.it+ii) = 0;
+        cnt += strmhalf(x,y,ii,-steps,&param);
+        if (cnt) {
+          pt = (galint *)param.it;
+          et = pt + lsiz;
+          do *pt++ &= MASK;
+          while (pt!=et);
+          ihi += pd[0];
+          if (!di) ilo -= pd[0];
+        }
+      }
+    cont_lo2:
+      i -= ++di;
+    } while (i>=0);
+  } while (j);
+
+  free (param.it);
+}
+
+gaint strmhalf (gadouble x, gadouble y, gaint i, gaint steps,
+                const struct strmparam *param) {
+gadouble xx,yy,xxsv,yysv,xmax,ymax,dx,dy,d,len,dis,dis2,xdis;
+gadouble *up,*vp,*cp,uv,uv2,vv,vv2,cv,cv2,auv,avv,*xfact;
+gaint isav,ii,jj,iz,jz,imn,imx,jmn,jmx;
+gaint is,iss,ilim,bflg,iacc,cnt,icol,scol,wrap,brk;
+char *upmask,*vpmask,*cpmask,*jt;
+
+  is = param->is;
+  iss = param->iss;
+  ilim = iss-1;
+  dis = param->dis;
+  dis2 = dis*dis;
+  xfact = param->xfact;
+  wrap = param->wrap;
+  cnt = 0;
+  scol = -9;
+  icol = 1;
+  isav = i;
+  iacc = 0;
+  len = steps>0 ? ARRMINLEN/2.0 : ARRDIST-ARRMINLEN/2.0;
+  bflg = 1;
+  brk = 0;
+  xmax = (gadouble)(is-1);
+  ymax = (gadouble)(param->js-1);
+  gxconv (x+1.0,y+1.0,&xx,&yy,3);
+  while (!brk) {
+    xxsv = xx; yysv = yy;
+    ii = (gaint)x;
+    jj = (gaint)y;
+    xx = x - (gadouble)ii;
+    yy = y - (gadouble)jj;
+    i = jj*is+ii;
+    up = param->u + i; upmask = param->umask + i;
+    vp = param->v + i; vpmask = param->vmask + i;
+    if (*upmask==0 || *vpmask==0) break;
+    uv = *up;
+    vv = *vp;
+    if (xx!=0.0) {
+      if (*(upmask+1)==0 || *(vpmask+1)==0) break;
+      uv += (*(up+1)-uv)*xx;
+      vv += (*(vp+1)-vv)*xx;
+    }
+    if (yy!=0.0) {
+      if (*(upmask+is)==0 || *(vpmask+is)==0) break;
+      uv2 = *(up+is);
+      vv2 = *(vp+is);
+      if (xx!=0.0) {
+        if (*(upmask+is+1)==0 || *(vpmask+is+1)==0) break;
+        uv2 += (*(up+is+1)-uv2)*xx;
+        vv2 += (*(vp+is+1)-vv2)*xx;
+      }
+      uv += (uv2-uv)*yy;
+      vv += (vv2-vv)*yy;
+    }
+    if (param->cflag) {
+      cp = param->c + i; cpmask = param->cmask + i;
+      if (*cpmask==0) { icol = UNDEFCOL; goto setcol; }
+      cv = *cp;
+      if (xx) {
+        if (*(cpmask+1)==0) { icol = UNDEFCOL; goto setcol; }
+        cv += (*(cp+1)-cv)*xx;
+      }
+      if (yy) {
+        if (*(cpmask+is)==0) { icol = UNDEFCOL; goto setcol; }
+        cv2 = *(cp+is);
+        if (xx) {
+          if (*(cpmask+is+1)==0) { icol = UNDEFCOL; goto setcol; }
+          cv2 += (*(cp+is+1)-cv2)*xx;
+        }
+        cv += (cv2-cv)*yy;
+      }
+      icol = gxshdc(param->shdlvs,param->shdcls,param->shdcnt,cv);
+    setcol:
+      if (icol!=scol && icol>-1) gxcolr(icol);
+      scol = icol;
+    }
+    uv *= param->dyt[jj];
+    vv *= param->dxt[ii];
+    if (param->isgeo)
+      vv *= cos(param->jgrab(param->jvals,param->jmin+y)*DEG2RAD);
+    auv = fabs(uv); avv=fabs(vv);
+    d = auv>avv ? auv : avv;
+    if (d==0.0) break;
+    d *= (gadouble)steps;
+    dx = uv/d;
+    dy = vv/d;
+    x += dx;
+    y += dy;
+
+    /* If we got beyond grid limits, plot to the edge or continue at the
+       other side, if we can wrap around */
+
+    if (x<0.0) {
+      if (wrap) x += xmax;
+      else {
+        y -= dy;
+        dy *= (dx-x)/dx;
+        y += dy;
+        x = 0.0;
+        brk = 1;
+      }
+    } else if (x>=xmax) {
+      if (wrap) x -= xmax;
+      else {
+        y -= dy;
+        dy *= (dx-(x-xmax))/dx;
+        y += dy;
+        x = xmax;
+        xx = (gadouble)iss;
+        ii = ilim;
+        brk = 1;
+        goto xdone;
+      }
+    }
+    xx = x*param->rscl;
+    ii = (gaint)xx;
+  xdone:
+    if (y<0.0) {
+      x -= dx;
+      dx *= (dy-y)/dy;
+      x += dx;
+      y = 0.0;
+      brk = 1;
+    }
+    else if (y>=ymax) {
+      x -= dx;
+      dx *= (dy-(y-ymax))/dy;
+      x += dx;
+      y = ymax;
+      yy = (gadouble)param->jss;
+      jj = param->jss-1;
+      brk = 1;
+      goto ydone;
+    }
+    yy = y*param->rscl;
+    jj = (gaint)yy;
+  ydone:
+    i = jj*iss + ii;
+
+    /* Break if we hit ourselves or if we make no progress */
+
+    if (i!=isav) {
+      if (*(param->it+i)) break;
+      *(param->it+isav) = -1;
+      isav = i;
+      iacc = 0;
+    } else {
+      if (++iacc==2*STEPSPERBOX) break;
+    }
+
+    /* Break if we are too near to another stream line */
+
+    xx += 0.5;
+    yy += 0.5;
+    jmn = (gaint)(yy-dis);
+    jmx = (gaint)(yy+dis);
+    if (jmn<0) jmn = 0;
+    if (jmx>param->jss) jmx = param->jss;
+    jt = param->it + jmn*iss;
+    dy = (gadouble)(jmn+1) - yy;
+    for (jz=jmn; jz<jmx; jz++) {
+      xdis = dis2 - dy*dy;
+      if (xdis>0.0) { /* May be negative due to rounding errors */
+        xdis = sqrt(xdis)*xfact[jz];
+        imn = (gaint)(xx-xdis);
+        imx = (gaint)(xx+xdis);
+        if (imn<=0) {
+          if (wrap) {
+            imn += ilim;
+            imx += ilim;
+          } else imn = 0;
         }
-        x = x - uv;
-        y = y - vv;
-        ii1 = (gaint)(x*rscl);
-        ij1 = (gaint)(y*rscl);
-        ii1 = ij1*iss + ii1;
-        if (ii1<0 || ii1>=siz) break;
-        if (*(it+ii1)==1) break;
-        if (ii1!=iisav && iisav>-1) *(it+iisav) = 1;
-        if (ii1==iisav) iacc++;
-        else iacc = 0;
-        if (iacc>10) break;
-        iisav = ii1;
-        gxconv (x+1.0,y+1.0,&xx,&yy,3);
-        if (icol>-1) {
-          if (bflg) {gxplot(xold,yold,3); bflg=0;}
-          gxplot (xx,yy,2);
-        } else bflg = 1;
-        xold = xx;
-        yold = yy;
-        acnt++;
-        if (acnt>20) {
-          if (icol>-1) strmar(xx,yy,xxsv,yysv);
-          acnt = 0;
+        for (iz=imn; iz<imx; ) {
+          if (*(jt+iz)>0)
+            return cnt;
+          if (iz==ilim) {
+            if (wrap) {
+              iz = 0;
+              imx -= ilim;
+            } else break;
+          } else iz++;
         }
-        xxsv = xx; yysv = yy;
       }
+      jt += iss;
+      dy += 1.0;
+    }
+
+    /* Plot line segment and arrow if appropriate */
+
+    gxconv (x+1.0,y+1.0,&xx,&yy,3);
+    if (icol>-1) {
+      if (bflg) {gxplot(xxsv,yysv,3); bflg=0;}
+      gxplot (xx,yy,2);
+      cnt++;
+    } else bflg = 1;
+    dx = xx-xxsv;
+    dy = yy-yysv;
+    d = sqrt(dx*dx + dy*dy);
+    if (steps>0) {
+      len += d;
+      if (len>=ARRDIST) {
+        if (icol>-1) { strmar(xxsv,yysv,xx,yy,param->asiz); bflg=1; }
+        len = 0.0;
+      }
+    } else {
+      if (len>=ARRDIST) {
+        if (icol>-1) { strmar(xx,yy,xxsv,yysv,param->asiz); bflg=1; }
+        len = 0.0;
+      }
+      len += d;
     }
-    i2++;
-    if (i2==iss) { i2 = 0; j2++; }
   }
-  free (it);
+
+  return cnt;
 }
 
-static gadouble a150 = 150.0*3.1416/180.0;
+static gadouble a150 = 150.0*pi/180.0;
 
-void strmar (gadouble xx1, gadouble yy1, gadouble xx2, gadouble yy2) {
+void strmar (gadouble xx1, gadouble yy1, gadouble xx2, gadouble yy2,
+             gadouble asiz) {
 gadouble dir;
 
+  if (asiz==0.0) return;
   dir = atan2(yy2-yy1,xx2-xx1);
   gxplot (xx2,yy2,3);
-  gxplot (xx2+0.05*cos(dir+a150),yy2+0.05*sin(dir+a150),2);
-  gxplot (xx2,yy2,3);
-  gxplot (xx2+0.05*cos(dir-a150),yy2+0.05*sin(dir-a150),2);
+  gxplot (xx2+asiz*cos(dir+a150),yy2+asiz*sin(dir+a150),2);
   gxplot (xx2,yy2,3);
+  gxplot (xx2+asiz*cos(dir-a150),yy2+asiz*sin(dir-a150),2);
 }
 
 /* Given a shade value, return the relevent color */
 
 gaint gxshdc (gadouble *shdlvs, gaint *shdcls, gaint shdcnt, gadouble val) {