import java.awt.*;
import java.awt.event.*;
import java.awt.image.*;
import java.applet.*;

public class StableFluids extends Applet
implements Runnable,MouseListener,MouseMotionListener {
	
	final int size = 128;
	final int mask = size-1;
	
	double[][] c = new double[size][size];
	double[][] u = new double[size][size];
	double[][] v = new double[size][size];
	double[][] tmp_c = new double[size][size];
	double[][] tmp_u = new double[size][size];
	double[][] tmp_v = new double[size][size];
	double[][] p = new double[size][size];
	double[][] rec_lap = new double[size][size];
	double[][] im = new double[size][size];
	
	boolean dragging = false;
	double mouseX,mouseY;
	double cursorX,cursorY;
	double cursorU,cursorV;
	
	FFT fft;
	
	int scale = 4;
	int[] pixels = new int[size*size];
	MemoryImageSource mis;
	Image image;
	Thread thread;
	
    public void update(Graphics g) {
		paint(g);
	}
	
    public void paint(Graphics g) {
		for (int i=0,y=0;y<size;y++) 
			for (int x=0;x<size;x++,i++) {
				int c = Math.min((int)this.c[y][x],256);
				pixels[i] = 0xff<<24|c<<16|c<<8|c;
			}
		mis.newPixels();
		g.drawImage(image,0,0,size*scale,size*scale,this);
    }
	
	public void mousePressed(MouseEvent e) {
		dragging = true;
		mouseX = cursorX = (double)e.getX()/scale;
		mouseY = cursorY = (double)e.getY()/scale;
		cursorU = 0;
		cursorV = 0;
	}
	
	public void mouseReleased(MouseEvent e) {
		dragging = false;
	}
	
	public void mouseEntered(MouseEvent e) {
	}
	
	public void mouseExited(MouseEvent e) {
	}
	
	public void mouseClicked(MouseEvent e) {
	}
	
	public void mouseDragged(MouseEvent e) {
		mouseX = (double)e.getX()/scale;
		mouseY = (double)e.getY()/scale;
	}
	
	public void mouseMoved(MouseEvent e) {
	}
	
	private void drag() {
		if (dragging) {
			cursorU += 2*((mouseX-cursorX)/10-cursorU)/10;
			cursorV += 2*((mouseY-cursorY)/10-cursorV)/10;
			cursorX += cursorU;
			cursorY += cursorV;
			for (int y=0;y<size;y++)
				for (int x=0;x<size;x++) {
					int dX = x-(int)cursorX;
					int dY = y-(int)cursorY;
					double w = Math.exp(-(dX*dX+dY*dY)/32);
					u[y][x] = (1-w)*u[y][x]+w*cursorU;
					v[y][x] = (1-w)*v[y][x]+w*cursorV;
					/*
					if (dX*dX+dY*dY<100) {
						u[y][x] = cursorU;
						v[y][x] = cursorV;
					}
					*/
				}
		}
	}
	
	double interpolate(double q[][],double x,double y) {
		int trX = (int)x;
		int trY = (int)y;
		double frX = x-trX;
		double frY = y-trY;
		int x0 = trX&mask;
		int y0 = trY&mask;
		int x1 = (x0+1)&mask;
		int y1 = (y0+1)&mask;
		return
			(1-frX)*(1-frY)*q[y0][x0]+
			(1-frX)*frY*q[y1][x0]+
			frX*(1-frY)*q[y0][x1]+
			frX*frY*q[y1][x1];
	}
	
	private void project() {
		for (int y=0;y<size;y++)
			for (int x=0;x<size;x++) {
				p[y][x] = u[y][x]+v[y][x]-u[y][(x+1)&mask]-v[(y+1)&mask][x];
				im[y][x] = 0;
			}
		fft.forward(p,im);
		for (int y=0;y<size;y++)
			for (int x=0;x<size;x++) {
				p[y][x] *= rec_lap[y][x];
				im[y][x] *= rec_lap[y][x];
			}
		fft.inverse(p,im);
		for (int y=0;y<size;y++)
			for (int x=0;x<size;x++) {
				u[y][x] += p[y][(x-1)&mask]-p[y][x];
				v[y][x] += p[(y-1)&mask][x]-p[y][x];
			}
	}
	
	private void move() {
		int x,y;
		double reX,reY;
		for (y=0,reY=size;y<size;y++,reY++)
			for (x=0,reX=size;x<size;x++,reX++) {
				double dx = (u[y][(x+1)&mask]+u[y][x])/2;
				double dy = (v[(y+1)&mask][x]+v[y][x])/2;
				tmp_c[y][x] = interpolate(c,reX-dx,reY-dy);
			}
		for (y=0,reY=size;y<size;y++,reY++)
			for (x=0,reX=size;x<size;x++,reX++) {
				double dx = u[y][x];
				double dy = (
					v[(y+1)&mask][x]+v[y][x]+
					v[(y+1)&mask][(x-1)&mask]+v[y][(x-1)&mask])/4;
				tmp_u[y][x] = interpolate(u,reX-dx,reY-dy);
			}
		for (y=0,reY=size;y<size;y++,reY++)
			for (x=0,reX=size;x<size;x++,reX++) {
				double dx = (
					u[(x+1)&mask][x]+u[y][x]+
					u[(y-1)&mask][(x+1)&mask]+u[(y-1)&mask][x])/4;
				double dy = v[y][x];
				tmp_v[y][x] = interpolate(v,reX-dx,reY-dy);
			}
		double[][] tmp_tmp_c = tmp_c;
		double[][] tmp_tmp_u = tmp_u;
		double[][] tmp_tmp_v = tmp_v;
		tmp_c = c;
		tmp_u = u;
		tmp_v = v;
		c = tmp_tmp_c;
		u = tmp_tmp_u;
		v = tmp_tmp_v;
	}
	
	private void saturate() {
		double maxc = 0;
		double minc = 1;
		for (int y=0;y<size;y++)
			for (int x=0;x<size;x++) {
				maxc = Math.max(maxc,c[y][x]);
				minc = Math.min(minc,c[y][x]);
			}
		double amp = 255/(maxc-minc);
		double off = amp*minc;
		for (int y=0;y<size;y++)
			for (int x=0;x<size;x++) {
				c[y][x] = amp*c[y][x]-off;
			}
	}

	public void simulate() {
		drag();
		project();
		move();
		saturate();
	}
	
	public void run() {
		try {
			while (!Thread.interrupted()) {
				simulate();
				repaint();
				Thread.sleep(10);
			}
		} catch (InterruptedException e) {
		}
	}
	
	public void start() {
		thread = new Thread(this);
		thread.setDaemon(true);
		thread.start();
	}
	
	public void stop() {
		thread.interrupt();
	}
	
    public void init() {
		
		for (int y=0;y<size;y++)
			for (int x=0;x<size;x++) {
				c[y][x] = ((x^y)&64)!=0?255:0;
				u[y][x] = 0;
				v[y][x] = 0;
			}
		
		fft = new FFT(size);
		
		for (int y=0;y<size;y++)
			for (int x=0;x<size;x++) {
				rec_lap[y][x] = 0;
				im[y][x] = 0;
			}
		
		rec_lap[0][0] = 4;
		rec_lap[0][1] = -1;
		rec_lap[1][0] = -1;
		rec_lap[0][mask] = -1;
		rec_lap[mask][0] = -1;
		
		fft.forward(rec_lap,im);
		
		for (int y=0;y<size;y++)
			for (int x=0;x<size;x++) {
				//rec_lap[y][x] = (x|y)!=0?1/(size*size*rec_lap[y][x]):0;
				rec_lap[y][x] = (x|y)!=0?1/(size*size*rec_lap[y][x]):0;
			}
		
		mis = new MemoryImageSource(size,size,pixels,0,size);
		mis.setAnimated(true);
		image = getToolkit().createImage(mis);
		addMouseListener(this);
		addMouseMotionListener(this);
    }
}

class FFT {
	private int size;
	private double rec_size;
	private int[] perm;
	private double[] coef_re;
	private double[] coef_im;
	
	public FFT(int size) {
		this.size = size;
		rec_size = 1.0/size;
		perm = new int[size];
		coef_re = new double[size+1];
		coef_im = new double[size+1];
		for (int i=0;i<=size;i++) {
			double theta = 2*Math.PI*i/size;
			coef_re[i] = Math.cos(theta);
			coef_im[i] = Math.sin(theta);
		}
		for (int i=0;i<size;i++) {
			int a = size;
			int b = 1;
			int c = 0;
			while (a!=0) {
				a >>= 1;
				if ((a&i)!=0) c |= b;
				b <<= 1;
			}
			perm[i] = c;
		}
	}
	
	private void permute(double[] re,double[] im) {
		for (int i=0;i<size;i++) 
			if (perm[i]>i) {
				double tmp_re = re[i];
				double tmp_im = im[i];
				re[i] = re[perm[i]];
				im[i] = im[perm[i]];
				re[perm[i]] = tmp_re;
				im[perm[i]] = tmp_im;
			}
	}
	
	public void forward(double[] re,double[] im) {
		permute(re,im);
		for (int i=0;i<size;i++) {
			re[i] *= rec_size;
			im[i] *= rec_size;
		}
		for (int n=2,N_n=size/n;n<=size;n<<=1,N_n>>=1)
			for (int i=0;i<size;i+=n) 
				for (int j=0,a=i,b=i+n/2,c=0;j<n/2;j++,a++,b++,c+=N_n) {
					double tmp_re = re[b]*coef_re[c]-im[b]*coef_im[c];
					double tmp_im = im[b]*coef_re[c]+re[b]*coef_im[c];
					re[b] = re[a]-tmp_re;
					im[b] = im[a]-tmp_im;
					re[a] += tmp_re;
					im[a] += tmp_im;
				}
	}
	
	public void inverse(double[] re,double[] im) {
		permute(re,im);
		for (int n=2,N_n=size/n;n<=size;n<<=1,N_n>>=1)
			for (int i=0;i<size;i+=n) 
				for (int j=0,a=i,b=i+n/2,c=size;j<n/2;j++,a++,b++,c-=N_n) {
					double tmp_re = re[b]*coef_re[c]-im[b]*coef_im[c];
					double tmp_im = im[b]*coef_re[c]+re[b]*coef_im[c];
					re[b] = re[a]-tmp_re;
					im[b] = im[a]-tmp_im;
					re[a] += tmp_re;
					im[a] += tmp_im;
				}
	}
	
	private void transpose(double[][] re,double[][] im) {
		for (int y=0;y<size;y++)
			for (int x=0;x<y;x++) {
				double tmp_re = re[y][x];
				double tmp_im = im[y][x];
				re[y][x] = re[x][y];
				im[y][x] = im[x][y];
				re[x][y] = tmp_re;
				im[x][y] = tmp_im;
			}
	}
	
	public void forward(double[][] re,double[][] im) {
		for (int i=0;i<size;i++) forward(re[i],im[i]);
		transpose(re,im);
		for (int i=0;i<size;i++) forward(re[i],im[i]);
	}
	
	public void inverse(double[][] re,double[][] im) {
		for (int i=0;i<size;i++) inverse(re[i],im[i]);
		transpose(re,im);
		for (int i=0;i<size;i++) inverse(re[i],im[i]);
	}
}