diff --git a/.gitignore b/.gitignore
index e73850b0..d04c5012 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,6 @@
 /.cache
 /.classpath
-/.project
+/**/.project
 /.pydevproject
 /.idea
 /.gradle
@@ -15,3 +15,5 @@ build/
 
 core/out/
 demo/out/
+core/bin/
+core/build/
diff --git a/core/src/main/java/edu/mines/jtk/dsp/LocalSmoothingFilter.java b/core/src/main/java/edu/mines/jtk/dsp/LocalSmoothingFilter.java
index 25159822..98029cea 100644
--- a/core/src/main/java/edu/mines/jtk/dsp/LocalSmoothingFilter.java
+++ b/core/src/main/java/edu/mines/jtk/dsp/LocalSmoothingFilter.java
@@ -77,7 +77,7 @@ public LocalSmoothingFilter() {
    * @param niter stop when number of iterations exceeds this limit.
    */
   public LocalSmoothingFilter(double small, int niter) {
-    _small = (float)small;
+    _small = small;
     _niter = niter;
     _ldk = new LocalDiffusionKernel(LocalDiffusionKernel.Stencil.D22);
   }
@@ -92,7 +92,7 @@ public LocalSmoothingFilter(double small, int niter) {
   public LocalSmoothingFilter(
     double small, int niter, LocalDiffusionKernel ldk)
   {
-    _small = (float)small;
+    _small = small;
     _niter = niter;
     _ldk = ldk;
   }
@@ -361,7 +361,7 @@ public void applySmoothL(double kmax, float[][][] x, float[][][] y) {
   private static Logger log = 
     Logger.getLogger(LocalSmoothingFilter.class.getName());
 
-  private float _small; // stop iterations when residuals are small
+  private double _small; // stop iterations when residuals are small
   private int _niter; // number of iterations
   private boolean _pc; // true, for preconditioned CG iterations
   private LocalDiffusionKernel _ldk; // computes y += (I+G'DG)x
@@ -618,23 +618,23 @@ private void solve(Operator2 a, float[][] b, float[][] x) {
     a.apply(x,q);
     saxpy(-1.0f,q,r); // r = b-Ax
     scopy(r,d); // d = r
-    float delta = sdot(r,r); // delta = r'r
-    float bnorm = sqrt(sdot(b,b));
-    float rnorm = sqrt(delta);
-    float rnormBegin = rnorm;
-    float rnormSmall = bnorm*_small;
+    double delta = sdot(r,r); // delta = r'r
+    double bnorm = sqrt(sdot(b,b));
+    double rnorm = sqrt(delta);
+    double rnormBegin = rnorm;
+    double rnormSmall = bnorm*_small;
     int iter;
     log.fine("solve: bnorm="+bnorm+" rnorm="+rnorm);
     for (iter=0; iter<_niter && rnorm>rnormSmall; ++iter) {
       log.finer("  iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin);
       a.apply(d,q); // q = Ad
-      float dq = sdot(d,q); // d'q = d'Ad
-      float alpha = delta/dq; // alpha = r'r/d'Ad
+      double dq = sdot(d,q); // d'q = d'Ad
+      double alpha = delta/dq; // alpha = r'r/d'Ad
       saxpy( alpha,d,x); // x = x+alpha*d
       saxpy(-alpha,q,r); // r = r-alpha*q
-      float deltaOld = delta;
+      double deltaOld = delta;
       delta = sdot(r,r); // delta = r'r
-      float beta = delta/deltaOld;
+      double beta = delta/deltaOld;
       sxpay(beta,r,d); // d = r+beta*d
       rnorm = sqrt(delta);
     }
@@ -649,27 +649,27 @@ private void solve(Operator3 a, float[][][] b, float[][][] x) {
     float[][][] r = new float[n3][n2][n1];
     scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); // r = b-Ax
     scopy(r,d);
-    float delta = sdot(r,r);
-    float bnorm = sqrt(sdot(b,b));
-    float rnorm = sqrt(delta);
-    float rnormBegin = rnorm;
-    float rnormSmall = bnorm*_small;
+    double delta = sdot(r,r);
+    double bnorm = sqrt(sdot(b,b));
+    double rnorm = sqrt(delta);
+    double rnormBegin = rnorm;
+    double rnormSmall = bnorm*_small;
     int iter;
     log.fine("solve: bnorm="+bnorm+" rnorm="+rnorm);
     for (iter=0; iter<_niter && rnorm>rnormSmall; ++iter) {
       log.finer("  iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin);
       a.apply(d,q);
-      float dq = sdot(d,q);
-      float alpha = delta/dq;
+      double dq = sdot(d,q);
+      double alpha = delta/dq;
       saxpy( alpha,d,x);
       if (iter%100<99) {
         saxpy(-alpha,q,r);
       } else {
         scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r);
       }
-      float deltaOld = delta;
+      double deltaOld = delta;
       delta = sdot(r,r);
-      float beta = delta/deltaOld;
+      double beta = delta/deltaOld;
       sxpay(beta,r,d);
       rnorm = sqrt(delta);
     }
@@ -688,25 +688,25 @@ private void solve(Operator2 a, Operator2 m, float[][] b, float[][] x) {
     scopy(b,r);
     a.apply(x,q);
     saxpy(-1.0f,q,r); // r = b-Ax
-    float bnorm = sqrt(sdot(b,b));
-    float rnorm = sqrt(sdot(r,r));
-    float rnormBegin = rnorm;
-    float rnormSmall = bnorm*_small;
+    double bnorm = sqrt(sdot(b,b));
+    double rnorm = sqrt(sdot(r,r));
+    double rnormBegin = rnorm;
+    double rnormSmall = bnorm*_small;
     m.apply(r,s); // s = Mr
     scopy(s,d); // d = s
-    float delta = sdot(r,s); // r's = r'Mr
+    double delta = sdot(r,s); // r's = r'Mr
     int iter;
     log.fine("msolve: bnorm="+bnorm+" rnorm="+rnorm);
     for (iter=0; iter<_niter && rnorm>rnormSmall; ++iter) {
       log.finer("  iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin);
       a.apply(d,q); // q = Ad
-      float alpha = delta/sdot(d,q); // alpha = r'Mr/d'Ad
+      double alpha = delta/sdot(d,q); // alpha = r'Mr/d'Ad
       saxpy( alpha,d,x); // x = x+alpha*d
       saxpy(-alpha,q,r); // r = r-alpha*q
       m.apply(r,s); // s = Mr
-      float deltaOld = delta;
+      double deltaOld = delta;
       delta = sdot(r,s); // delta = r's = r'Mr
-      float beta = delta/deltaOld;
+      double beta = delta/deltaOld;
       sxpay(beta,s,d); // d = s+beta*d
       rnorm  = sqrt(sdot(r,r));
     }
@@ -721,19 +721,19 @@ private void solve(Operator3 a, Operator3 m, float[][][] b, float[][][] x) {
     float[][][] r = new float[n3][n2][n1];
     float[][][] s = new float[n3][n2][n1];
     scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); // r = b-Ax
-    float bnorm = sqrt(sdot(b,b));
-    float rnorm = sqrt(sdot(r,r));
-    float rnormBegin = rnorm;
-    float rnormSmall = bnorm*_small;
+    double bnorm = sqrt(sdot(b,b));
+    double rnorm = sqrt(sdot(r,r));
+    double rnormBegin = rnorm;
+    double rnormSmall = bnorm*_small;
     m.apply(r,s); // s = Mr
     scopy(s,d); // d = s
-    float delta = sdot(r,s); // r's = r'Mr
+    double delta = sdot(r,s); // r's = r'Mr
     int iter;
     log.fine("msolve: bnorm="+bnorm+" rnorm="+rnorm);
     for (iter=0; iter<_niter && rnorm>rnormSmall; ++iter) {
       log.finer("  iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin);
       a.apply(d,q); // q = Ad
-      float alpha = delta/sdot(d,q); // alpha = r'Mr/d'Ad
+      double alpha = delta/sdot(d,q); // alpha = r'Mr/d'Ad
       saxpy( alpha,d,x); // x = x+alpha*d
       if (iter%100<99) {
         saxpy(-alpha,q,r); // r = r-alpha*q
@@ -741,9 +741,9 @@ private void solve(Operator3 a, Operator3 m, float[][][] b, float[][][] x) {
         scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); // r = b-Ax
       }
       m.apply(r,s); // s = Mr
-      float deltaOld = delta;
+      double deltaOld = delta;
       delta = sdot(r,s); // delta = r's = r'Mr
-      float beta = delta/deltaOld;
+      double beta = delta/deltaOld;
       sxpay(beta,s,d); // d = s+beta*d
       rnorm  = sqrt(sdot(r,r));
     }
@@ -783,34 +783,34 @@ public void compute(int i3) {
   }
 
   // Returns the dot product x'y.
-  private static float sdot(float[][] x, float[][] y) {
+  private static double sdot(float[][] x, float[][] y) {
     int n1 = x[0].length;
     int n2 = x.length;
-    float d = 0.0f;
+    double d = 0.0;
     for (int i2=0; i2<n2; ++i2) {
       float[] x2 = x[i2], y2 = y[i2];
       for (int i1=0; i1<n1; ++i1) {
-        d += x2[i1]*y2[i1];
+        d += (double)x2[i1]*(double)y2[i1];
       }
     }
     return d;
   }
-  private static float sdot(final float[][][] x, final float[][][] y) {
+  private static double sdot(final float[][][] x, final float[][][] y) {
     final int n3 = x.length;
-    final float[] d3 = new float[n3];
+    final double[] d3 = new double[n3];
     Parallel.loop(n3,new Parallel.LoopInt() {
       public void compute(int i3) {
         d3[i3] = sdot(x[i3],y[i3]);
       }
     });
-    float d = 0.0f;
+    double d = 0.0;
     for (int i3=0; i3<n3; ++i3)
       d += d3[i3];
     return d;
   }
 
   // Computes y = y + a*x.
-  private static void saxpy(float a, float[][] x, float[][] y) {
+  private static void saxpy(final double a, float[][] x, float[][] y) {
     int n1 = x[0].length;
     int n2 = x.length;
     for (int i2=0; i2<n2; ++i2) {
@@ -821,7 +821,7 @@ private static void saxpy(float a, float[][] x, float[][] y) {
     }
   }
   private static void saxpy(
-    final float a, final float[][][] x, final float[][][] y)
+    final double a, final float[][][] x, final float[][][] y)
   {
     final int n3 = x.length;
     Parallel.loop(n3,new Parallel.LoopInt() {
@@ -832,18 +832,18 @@ public void compute(int i3) {
   }
 
   // Computes y = x + a*y.
-  private static void sxpay(float a, float[][] x, float[][] y) {
+  private static void sxpay(final double a, float[][] x, float[][] y) {
     int n1 = x[0].length;
     int n2 = x.length;
     for (int i2=0; i2<n2; ++i2) {
       float[] x2 = x[i2], y2 = y[i2];
       for (int i1=0; i1<n1; ++i1) {
-        y2[i1] = a*y2[i1]+x2[i1];
+        y2[i1] = (float)(a*(double)y2[i1])+x2[i1];
       }
     }
   }
   private static void sxpay(
-    final float a, final float[][][] x, final float[][][] y)
+    final double a, final float[][][] x, final float[][][] y)
   {
     final int n3 = x.length;
     Parallel.loop(n3,new Parallel.LoopInt() {