diff --git a/flat.js b/flat.js
index c96d21c..cd0b4cf 100644
--- a/flat.js
+++ b/flat.js
@@ -27,6 +27,20 @@ function smoothing(radius, distance)
   return value * value * value / volume;
 }
 
+function smoothingSlope(radius, distance)
+{
+  if (distance >= radius)
+    return 0;
+  const value = radius * radius - distance * distance;
+  const scale = -24.0 / (Math.PI * Math.pow(radius, 8));
+  return scale * distance * value * value;
+}
+
+function densityToPressure(density, targetDensity, pressure)
+{
+  return (density - targetDensity) * pressure;
+}
+
 class FlatRenderer {
   positionProperty(position, circleRadius, particleCount)
   {
@@ -45,6 +59,69 @@ class FlatRenderer {
     return property;
   }
 
+  positionGradient(position, circleRadius, particleCount)
+  {
+    var property = new Float32Array([0, 0]);
+    for (var i = 0; i < particleCount; i++) {
+      const particlePosition = [this.particles[i * this.particleStride + 0],
+                                this.particles[i * this.particleStride + 1]];
+      const signedDistance = sub(position, particlePosition);
+      const distance = length(signedDistance);
+
+      const direction = [signedDistance[0] / distance, signedDistance[1] / distance];
+      const particleProperty = this.particles[i * this.particleStride + 4];
+      const particleDensity = this.particles[i * this.particleStride + 5];
+
+      const slope = smoothingSlope(circleRadius, distance);
+
+      const weight = particleProperty * slope / particleDensity;
+      property[0] += direction[0] * weight;
+      property[1] += direction[1] * weight;
+    }
+    return property;
+  }
+
+  positionPressureForce(index, circleRadius, particleCount, targetDensity, pressure, mass)
+  {
+    var pressureForce = new Float32Array([0, 0]);
+
+    const position = [this.particles[index * this.particleStride + 0],
+                      this.particles[index * this.particleStride + 1]];
+
+    for (var i = 0; i < particleCount; i++) {
+      if (i === index) continue;
+
+      const particlePosition = [this.particles[i * this.particleStride + 0],
+                                this.particles[i * this.particleStride + 1]];
+      const signedDistance = sub(position, particlePosition);
+      const distance = length(signedDistance);
+
+      var direction;
+      if (distance === 0) {
+        const x = Math.random();
+        const y = Math.random();
+        const l = length([x, y]);
+        direction = [x / l, y / l];
+      } else {
+        direction = [signedDistance[0] / distance, signedDistance[1] / distance];
+      }
+      const particleDensity = this.particles[i * this.particleStride + 5];
+
+      const slope = smoothingSlope(circleRadius, distance);
+
+      const sharedPressure = (
+        densityToPressure(this.particles[index * this.particleStride * 5], targetDensity, pressure) +
+        densityToPressure(particleDensity, targetDensity, pressure)
+      ) / 2;
+
+      const weight = sharedPressure * slope * mass / particleDensity;
+
+      pressureForce[0] += direction[0] * weight;
+      pressureForce[1] += direction[1] * weight;
+    }
+    return pressureForce;
+  }
+
   particleDensity(aPosition, circleRadius, particleCount)
   {
     var density = 0.0;
@@ -62,7 +139,7 @@ class FlatRenderer {
   {
     this.maxDim = 32;
     this.maxParticles = this.maxDim * this.maxDim;
-    this.particleStride = 4 * 2; // in elements
+    this.particleStride = 4 * 3; // in elements
     this.particles = new Float32Array(this.maxParticles * this.particleStride);
 
     this.maxGradientsDim = 32;
@@ -74,7 +151,7 @@ class FlatRenderer {
     this.gradientBuffers = [];
     this.particleBindGroups = [];
 
-    this.particleConfiguration = new Float32Array(3 * 4);
+    this.particleConfiguration = new Float32Array(4 * 4);
     this.particleConfigurationBuffers = [];
 
     for (let i = 0; i < 2; i++) {
@@ -249,20 +326,13 @@ class FlatRenderer {
         const position = [(x / this.maxGradientsDim) + 0.5 / this.maxGradientsDim,
                           (y / this.maxGradientsDim) + 0.5 / this.maxGradientsDim];
 
-        const stepSize = 0.001;
-        const px = [position[0] + stepSize, position[1]];
-        const py = [position[0], position[1] + stepSize];
-
-        const property = this.positionProperty(position, radius, count);
-
-        const dx = this.positionProperty(px, radius, count) - property;
-        const dy = this.positionProperty(py, radius, count) - property;
+        const gradient = this.positionGradient(position, radius, count);
 
         this.gradients[i * this.gradientStride + 0] = position[0];
         this.gradients[i * this.gradientStride + 1] = position[1];
 
-        this.gradients[i * this.gradientStride + 4] = dx / stepSize;
-        this.gradients[i * this.gradientStride + 5] = dy / stepSize;
+        this.gradients[i * this.gradientStride + 4] = gradient[0];
+        this.gradients[i * this.gradientStride + 5] = gradient[1];
       }
     }
 
@@ -286,11 +356,55 @@ class FlatRenderer {
                              this.particles, 0, this.particles.length);
   }
 
+  updateSimulation(device, frameNumber, configuration)
+  {
+    for (var i = 0; i < configuration.particleCount; i++) {
+      const force = this.positionPressureForce(i,
+                                               configuration.circleRadius,
+                                               configuration.particleCount,
+                                               configuration.targetDensity,
+                                               configuration.pressure,
+                                               configuration.particleMass1);
+      if (force[0] == NaN)
+        throw new Error("force");
+
+      const particleDensity = this.particles[i * this.particleStride + 5];
+
+      const acceleration = [force[0] / particleDensity, force[1] / particleDensity];
+
+      this.particles[i * this.particleStride + 8] += acceleration[0];
+      this.particles[i * this.particleStride + 9] += acceleration[1] + 0.001;
+
+      // collisions
+    }
+
+    for (var i = 0; i < configuration.particleCount; i++) {
+      this.particles[i * this.particleStride + 0] -= this.particles[i * this.particleStride + 8];
+      this.particles[i * this.particleStride + 1] -= this.particles[i * this.particleStride + 9];
+
+      const position = [this.particles[i * this.particleStride + 0],
+                        this.particles[i * this.particleStride + 1]];
+
+      for (var c = 0; c < 2; c++) {
+        if (position[c] < 0) {
+          this.particles[i * this.particleStride + c] = 0.01;
+          this.particles[i * this.particleStride + c + 8] *= -1 * 0.9;
+        }
+
+        if (position[c] > 1) {
+          this.particles[i * this.particleStride + c] = 0.99;
+          this.particles[i * this.particleStride + c + 8] *= -1 * 0.9;
+        }
+      }
+    }
+  }
+
   updateConfiguration(device, frameNumber, configuration)
   {
     const PARTICLE = 0;
     const CIRCLE = 1;
     const MOUSE = 2;
+    const D = 3;
 
     const config = this.particleConfiguration;
     config[PARTICLE * 4 + 0] = configuration.particleCount;
@@ -308,6 +422,11 @@ class FlatRenderer {
     config[MOUSE * 4 + 2] = configuration.lineThickness;
     config[MOUSE * 4 + 3] = configuration.lineLength;
 
+    config[D * 4 + 0] = configuration.targetDensity;
+    config[D * 4 + 1] = configuration.pressure;
+    config[D * 4 + 2] = 0;
+    config[D * 4 + 3] = 0;
+
     device.queue.writeBuffer(this.particleConfigurationBuffers[frameNumber], 0,
                              this.particleConfiguration, 0, this.particleConfiguration.length)
   }
@@ -317,6 +436,7 @@ class FlatRenderer {
     // particle update must be before gradient update
     this.updateParticles(device, frameNumber, configuration);
     this.updateGradients(device, frameNumber, configuration);
+    this.updateSimulation(device, frameNumber, configuration);
     this.updateConfiguration(device, frameNumber, configuration);
 
     renderPass.setPipeline(this.renderPipeline);
diff --git a/flat.wgsl b/flat.wgsl
index 4316cac..b2a4623 100644
--- a/flat.wgsl
+++ b/flat.wgsl
@@ -8,12 +8,14 @@ struct Configuration {
 struct Particle {
   position: vec4f,
   value: vec4f,
+  velocity: vec4f,
 };
 
 struct ParticleConfiguration {
   particle: vec4f,
   circle: vec4f,
   mouse: vec4f,
+  d: vec4f,
 };
 
 struct Gradient {
@@ -88,6 +90,11 @@ fn signedDistanceSegment(p: vec2f, a: vec2f, b: vec2f, r: f32) -> f32
   return length(pa - h * ba) - r;
 }
 
+fn densityToPressure(density: f32, targetDensity: f32, pressure: f32) -> f32
+{
+  return (density - targetDensity) * pressure;
+}
+
 @fragment
 fn fragmentMain(input: VertexOutput) -> @location(0) vec4f
 {
@@ -107,12 +114,18 @@ fn fragmentMain(input: VertexOutput) -> @location(0) vec4f
   let lineThickness = particleConfiguration.mouse.z;
   let lineLength = particleConfiguration.mouse.w;
 
+  let targetDensity = particleConfiguration.d.x;
+  let pressure = particleConfiguration.d.y;
+
+  let particleMass = mix(particleMass1, particleMass2, mass12Mix);
+
   var property: f32 = 0.0;
 
   for (var i: u32 = 0; i < particleCount; i++) {
     let distance = length(input.texture.xy - particle[i].position.xy);
     if (distance < particleSize) {
-      color.x = particle[i].value.x;
+      //color.x = particle[i].value.x;
+      color.x = 1.0;
     }
     //let circleDistance = length(input.texture.xy - circlePosition.xy);
     let influence = smoothing(circleRadius, distance);
@@ -120,16 +133,21 @@ fn fragmentMain(input: VertexOutput) -> @location(0) vec4f
     let particleProperty = particle[i].value.x;
     let particleDensity = particle[i].value.y;
 
-    let p1 = particleProperty * influence;
+    //let p1 = particleProperty * influence;
+    let p1 = influence;
     let p2 = p1 / particleDensity;
 
-    property += mix(p1, p2, mass12Mix);
+    //property += mix(p1, p2, mass12Mix);
+    property += influence;
   }
-  color.z = property * mix(particleMass1, particleMass2, mass12Mix);
+
+  let density = property * particleMass;
+  let fragmentPressure = densityToPressure(density, targetDensity, pressure) * 0.01;
+  color.z = density;
 
   for (var i: u32 = 0; i < 1024; i++) {
     let a = gradient[i].position.xy;
-    let b = a + gradient[i].delta.xy * lineLength;
+    let b = a + gradient[i].delta.xy * particleMass * lineLength;
     let distance = signedDistanceSegment(input.texture.xy, a, b, 0);
     let headDistance = length(input.texture.xy - b);
     if ((distance < lineThickness) || (headDistance < lineThickness * 3.5)) {
@@ -140,6 +158,27 @@ fn fragmentMain(input: VertexOutput) -> @location(0) vec4f
   let testColor = test(input.texture.xy);
   color.z += testColor * testIntensity;
 
-  return vec4(color.xyz, 1.0);
+  let r = vec3f(1, 0, 0);
+  let g = vec3f(0, 1, 0);
+  let b = vec3f(0, 0, 1);
+
+  let n = smoothstep(0, -0.01, fragmentPressure);
+  let p = smoothstep(0, 0.01, fragmentPressure);
+  let z = smoothstep(0.01, 0, abs(fragmentPressure));
+  //return vec4(n, z, p, 1);
+
+  let pressureColor = (r * n + g * z + b * p) * density;
+  return vec4(pressureColor.xyz * (1.0 - color.x) + vec3(1, 1, 1) * color.y, 1.0);
+  //return vec4(color.xyz, 1.0);
   //return vec4(color.x, 0, 0, 1.0);
+
+  /*
+  if (fragmentPressure < -0.1) {
+    return vec4(1, 0, 0, 1);
+  } else if (fragmentPressure > 0.1) {
+    return vec4(0, 1, 0, 1);
+  } else {
+    return vec4(0, 0, 1, 1);
+  }
+  */
 }
diff --git a/index.js b/index.js
index 9d0eebe..7d16781 100644
--- a/index.js
+++ b/index.js
@@ -206,16 +206,7 @@ function remap(low1, high1, low2, high2, value)
   return low2 + (value - low1) * (high2 - low2) / (high1 - low1);
 }
 
-const sliderDefaults = {
-  "circle radius": 385,
-  "circle thickness": 284,
-  "mass12 mix": 1,
-  "particle count": 360,
-  "particle mass1": 126,
-  "particle mass2": 417,
-  "particle size": 51,
-  "test intensity": 0,
-};
+const sliderDefaults = {"test intensity":"0","mass12 mix":"0","circle thickness":"171","line thickness":"107","particle mass2":"13","particle count":"19","target density":"954","particle size":"70","pressure":"19","line length":"6","particle mass1":"1024","circle radius":"379"};
 
 class Slider {
   constructor(label, low, high, max = 1024)
@@ -254,6 +245,7 @@ class Slider {
     } else if (this.label in sliderDefaults) {
       this.input.value = sliderDefaults[this.label];
     }
+
     this.update();
   }
 
@@ -276,7 +268,7 @@ class Slider {
 
 const sliders = {
   particleSize: new Slider("particle size", 0.0, 0.05),
-  particleMass1: new Slider("particle mass1", 0.0, 0.01),
+  particleMass1: new Slider("particle mass1", 0.0, 0.05),
   particleMass2: new Slider("particle mass2", 0.0, 1.0),
   particleCount: new Slider("particle count", 0.0, 1024.0),
   circleRadius: new Slider("circle radius", 0.0, 0.5),
@@ -284,7 +276,9 @@ const sliders = {
   testIntensity: new Slider("test intensity", 0.0, 2.0),
   mass12Mix: new Slider("mass12 mix", 0.0, 1.0, 1),
   lineThickness: new Slider("line thickness", 0.0, 0.005),
-  lineLength: new Slider("line length", 0.0, 0.025),
+  lineLength: new Slider("line length", 0.0, 9.0),
+  targetDensity: new Slider("target density", -10.0, 10.0),
+  pressure: new Slider("pressure", 0.0, 1.0),
 };
 
 const mousePosition = new Float32Array([0, 0, 0, 0]);
@@ -303,6 +297,8 @@ function updateSliders()
     mass12Mix: sliders.mass12Mix.update(),
     lineThickness: sliders.lineThickness.update(),
     lineLength: sliders.lineLength.update(),
+    targetDensity: sliders.targetDensity.update(),
+    pressure: sliders.pressure.update(),
   };
 }