I am trying to play with the approach in the following Shadertoy
https://www.shadertoy.com/view/lsKGRW
When I tried to do this myself using my own fragment shader, my Runge Kutta just kind of "explodes", there's nothing of the nice Schrodinger wave transition happening. Something non-linear is occuring that isn't occuring in the Shadertoy example.
Now, after a lot of experiments, I just tried to basically copy the code of the Shadertoy to my own glsl fragment shader.
Note that this uses a small bit of preprocessor magic to make sure the code appears in the right buffer, by using GlslCanvas.js. This is explained here https://marketplace.visualstudio.com/items?itemName=circledev.glsl-canvas
as you can see from the code here, the buffers it creates are interpolated "NEAREST"
This is the final code:
#version 300 es
precision highp float;
uniform vec2 u_resolution;
uniform vec2 u_mouse;
uniform float u_time;
#define frame (u_time*60.0)
uniform sampler2D u_buffer0;
uniform sampler2D u_buffer1;
uniform sampler2D u_buffer2;
uniform sampler2D u_buffer3;
#define REFLECTOR_SIZE 5.
#define dt 0.25
float potential(vec2 p) {
bool isTopReflector = p.y < REFLECTOR_SIZE;
bool isBottomReflector = p.y > u_resolution.y - REFLECTOR_SIZE;
bool isSplit = (p.x < u_resolution.x/2. && int(p.y) == int(u_resolution.y/2.));
return float(isTopReflector || isBottomReflector || isSplit);
}
#define cis(theta) vec2(cos(theta),sin(theta))
#define length2(p) dot(p,p)
vec2 psi0(vec2 p) {
p = p / u_resolution.xy - 0.5;
p.x *= u_resolution.x / u_resolution.y;
return exp(-100.*length2(p+vec2(0.675,0.225))) * cis(250.*(p.x+p.y));
}
vec2 divi(vec2 c) {
/* divide by sqrt(-1), ie. rotate 270 deg */
return vec2(c.y, -c.x);
}
#define psi(p) texture(u_buffer0,(p)/u_resolution.xy).xy
#define k1Buf(p) texture(u_buffer1,(p)/u_resolution.xy).xy
#define k2Buf(p) texture(u_buffer2,(p)/u_resolution.xy).xy
#define k3Buf(p) texture(u_buffer3,(p)/u_resolution.xy).xy
#define y1(p) (psi(p))
#define y2(p) (psi(p) + 0.5*dt*k1Buf(p))
#define y3(p) (psi(p) + 0.5*dt*k2Buf(p))
#define y4(p) (psi(p) + dt*k3Buf(p))
vec2 k1P(vec2 p) {
vec2 c = y1(p);
vec2 n = y1(p + vec2(0.,1.));
vec2 e = y1(p + vec2(1,0));
vec2 s = y1(p - vec2(0,1));
vec2 w = y1(p - vec2(1,0));
vec2 laplacian = n + e + s + w - 4.*c;
return divi(-laplacian + potential(p) * c);
}
vec2 k2P(vec2 p) {
vec2 c = y2(p);
vec2 n = y2(p + vec2(0,1));
vec2 e = y2(p + vec2(1,0));
vec2 s = y2(p - vec2(0,1));
vec2 w = y2(p - vec2(1,0));
vec2 laplacian = n + e + s + w - 4.*c;
return divi(-laplacian + potential(p) * c);
}
vec2 k3P(vec2 p) {
vec2 c = y3(p);
vec2 n = y3(p + vec2(0,1));
vec2 e = y3(p + vec2(1,0));
vec2 s = y3(p - vec2(0,1));
vec2 w = y3(p - vec2(1,0));
vec2 laplacian = n + e + s + w - 4.*c;
return divi(-laplacian + potential(p) * c);
}
vec2 k4P(vec2 p) {
vec2 c = y4(p);
vec2 n = y4(p + vec2(0,1));
vec2 e = y4(p + vec2(1,0));
vec2 s = y4(p - vec2(0,1));
vec2 w = y4(p - vec2(1,0));
vec2 laplacian = n + e + s + w - 4.*c;
return divi(-laplacian + potential(p) * c);
}
#if defined(BUFFER_0)
out vec4 o_color;
#define rk4(p) (psi(p) + dt * (1.*k1Buf(p) + 2.*k2Buf(p) + 2.*k3Buf(p) + k4P(p))/6.)
void main() {
o_color = vec4(1.,0.,0.,1.);
vec2 st = gl_FragCoord.xy;
o_color.xy = frame < 10. ? psi0(st) : rk4(st);
}
#elif defined(BUFFER_1)
out vec4 o_color;
void main() {
o_color = vec4(1.,0.,0.,1.);
vec2 st = gl_FragCoord.xy;
o_color.xy = frame < 10. ? psi0(st) : k1P(st);
}
#elif defined(BUFFER_2)
out vec4 o_color;
void main() {
o_color = vec4(1.,0.,0.,1.);
vec2 st = gl_FragCoord.xy;
o_color.xy = frame < 10. ? psi0(st): k2P(st);
}
#elif defined(BUFFER_3)
out vec4 o_color;
void main() {
vec2 st = gl_FragCoord.xy;
o_color.xy = frame < 10. ? psi0(st) : k3P(st);
}
#else
out vec4 o_color;
#define PI 3.14159265358979323846
#define hue2rgb(h) clamp(abs(mod(6.*(h)+vec3(0,4,2),6.)-3.)-1.,0.,1.)
void main() {
vec2 p = gl_FragCoord.xy;
vec2 v = k3Buf(p);
o_color = vec4(1.5 * length(v) * hue2rgb(atan(v.y,v.x)/(2.*PI)) + 0.25 * potential(p), 1);
}
#endif
I have also attached a movie here of how it looks like when I run it. https://gfycat.com/possiblesinfulbantamrooster
Note that the video starts in the middle of a run, and then another run starts a brief second after it - that's just due to the setup of my video recording.
Now, I know a thing or two about GLSL shaders, and about the Runge Kutta technique. I simply cannot for the life of me find why my code is "exploding" and the Shadertoy example isn't. I spent at least four hours debugging this in all sorts of ways, and the only thing I can tell you is that it somehow misbehaves.
There aren't many differences between the Shadertoy and my code anymore. The obvious one is the "Frame" counter, but that looks pretty innocent. It's also not like this whole thing is very fragile, the behaviour is quite different.
The only thing I learned from debugging is, that if I replace the runge kutta 4 method to just return k1P(..) instead of using the four different k values, and I let this k1P return "e" instead of the laplacian interpolation between n/e/w/s, so "return e", then there is a subtle difference in the behaviour: my glsl shader just expands and then "freezes" into discrete values of r/g/b, whereas the Shadertoy grows continuously, until glitching out.
As you can understand, these kind of differential techniques where everything is tied together, are very difficult to debug...
