Diffusion simulation of ink dripping on the surface of the water
I made it using the following diffusion equation.
I'm trying to do path tracing, so I'm outputting an image file in ppm format for practice.
InkFlow.java
package flow;
import java.io.FileWriter;
import java.io.IOException;
public class InkFlow{
//Simulate the diffusion of ink dropped on the surface of the water
static double INK = 255.; //Maximum RGB value
static int mx; //Number of horizontal meshes(Number of pixels in the image? become)
static int my; //Vertical
static int fc=0; //Filename variable
public static void main(String[] args) throws IOException {
/*******
*Setting(Change here)
*******/
//Water surface setting
int x = 10; //Lateral length[cm]
int y = 10; //Vertical length[cm]
//Places where ink is dropped (upper left, upper right, left, center, right, lower left, lower, lower right)
boolean CPlace[][] = { // Cyan
{false,false,false},
{true,false,false},
{false,false,false}
};
boolean MPlace[][] = { // Magenta
{false,true,false},
{false,false,false},
{false,false,false}
};
boolean YPlace[][] = { // Yellow
{false,false,false},
{false,false,false},
{false,true,false}
};
//The size of the ink to hang down
int InkSize = 13;
//Time setting
int t=0; //Initial time 0[msec]
int T=10000; //End time[msec]
int ut = 1000; //File output unit time[msec]
int inkt = 5000; //Time to drip ink[msec]
//Variables for diffusion
double ul = 0.1; //Unit distance[cm]
mx = (int) (x / ul);
my = (int) (y / ul);
double D = 0.25; //Diffusion coefficient(<0.25)
LogSettings(x,y,InkSize,T,ut, inkt ,ul,mx,my,D); //Configuration log
/*******
*Initialization
*******/
//Color initialization (white)
double R[][] = InitialInk();
double G[][] = InitialInk();
double B[][] = InitialInk();
//Drip ink
DipInk(CPlace, MPlace, YPlace, InkSize, R, G, B);
Output(R, G, B);
/*******
*Main loop
*******/
while(t++ <= T) { //Loop until end time T
//diffusion
R = CalcDiffusion(R, D);
G = CalcDiffusion(G, D);
B = CalcDiffusion(B, D);
//Drip ink
if(t < inkt) { //inkt msec hang down
DipInk(CPlace, MPlace, YPlace, InkSize, R, G, B);
}
//output
if(t % ut == 0){
Output(R, G, B);
}
}
System.out.println("done");
}
/*******
*Various methods
*******/
static double[][] InitialInk(){ //Ink initialization
double a[][] = new double[mx][my];
for(int i = 0; i < mx; i++) {
for(int j = 0; j < my; j++) {
a[i][j] = INK;
}
}
return a;
}
static double[][] CalcDiffusion(double[][] a, double D){ //Diffusion calculation
double b[][] = new double[mx][my];
//Inner calculation
for(int i=1;i<(mx-1);i++){
for(int j=1;j<(my-1);j++){ //Diffusion equation
b[i][j] = a[i][j] + D * (a[i+1][j] + a[i-1][j] - 2*a[i][j]) + D * (a[i][j+1] + a[i][j-1] - 2*a[i][j]);
}
}
//Outer calculation
for(int i=0;i<mx;i++){
b[i][my-1] = b[i][my-2];
b[i][0] = b[i][1];
}
for(int j=0;j<my;j++){
b[mx-1][j] = b[mx-2][j];
b[0][j] = b[1][j];
}
return b;
}
static void DipInk(boolean[][] CPlace, boolean[][] MPlace, boolean[][] YPlace, int InkSize, double[][] R, double[][] G, double[][] B) {
for(int i = 0; i < CPlace.length; i++) {
for(int j = 0; j < CPlace[0].length; j++) {
if(CPlace[i][j] || MPlace[i][j] || YPlace[i][j]) {
if(CPlace[i][j]) {
DipWithInkSize(R, InkSize, i, j, 0);
} else {
DipWithInkSize(R, InkSize, i, j, 212);
}
if(MPlace[i][j]) {
DipWithInkSize(G, InkSize, i, j, 0);
} else {
DipWithInkSize(G, InkSize, i, j, 212);
}
if(YPlace[i][j]) {
DipWithInkSize(B, InkSize, i, j, 0);
} else {
DipWithInkSize(B, InkSize, i, j, 212);
}
}
}
}
}
static void DipWithInkSize(double[][] C, int InkSize, int a, int b, double val) {
//Increase the ink mass in a spiral shape according to the ink size
C[(int) (mx*(2*a+1)/6)][(int) (my*(2*b+1)/6)] = val; //First ink
int itr = 1;
int i = 1;
while(true) {
for(int x = itr, y = 0; x > 0; x--,y++) {
C[(int) (mx*(2*a+1)/6)+x][(int) (my*(2*b+1)/6)+y] = val;
i++;
if(i >= InkSize) {
return;
}
}
for(int x = 0, y = itr; y > 0; x--, y--) {
C[(int) (mx*(2*a+1)/6)+x][(int) (my*(2*b+1)/6)+y] = val;
i++;
if(i >= InkSize) {
return;
}
}
for(int x = -itr, y = 0; x < 0; x++,y--) {
C[(int) (mx*(2*a+1)/6)+x][(int) (my*(2*b+1)/6)+y] = val;
i++;
if(i >= InkSize) {
return;
}
}
for(int x = 0, y = -itr; y < 0; x++, y++) {
C[(int) (mx*(2*a+1)/6)+x][(int) (my*(2*b+1)/6)+y] = val;
i++;
if(i >= InkSize) {
return;
}
}
itr++;
}
}
static void Output(double[][] R, double[][] G, double[][] B) throws IOException { //Image file(ppm format)Output of
String fname = "3_" + fc++ + ".ppm";
FileWriter fw = new FileWriter(fname); // Difine file name
// Header
fw.write("P3 \r\n");
fw.write("#The P3 means colors are in ASCII, then columns and rows, then 255 for max color, then RGB triplets \r\n");
fw.write(mx + " " + my + " \r\n");
fw.write("" + (int) INK + " \r\n");
// Body
for(int i = 0; i < mx; i++) {
for(int j = 0; j < my; j++) {
fw.write("" + (int) R[i][j] + " " + (int) G[i][j] + " " + (int) B[i][j]);
fw.write(" \r\n");
}
}
fw.close();
System.out.println("output: " + fname);
}
static void LogSettings(int x, int y, int InkSize, int T, int ut, int inkt, double ul, int mx, int my, double D) {
System.out.println("Water surface(Horizontal x vertical) : " + x + " cm × " + y + " cm");
System.out.println("Ink size: " + ul*ul*InkSize + " cm^2");
System.out.println("Time to hang:" + inkt + " msec");
System.out.println("ending time: " + T + " msec");
System.out.println("Output unit time: " + ut + " msec");
System.out.println("Unit distance: " + ul + " cm");
System.out.println("Number of horizontal meshes: " + mx);
System.out.println("Number of vertical meshes: " + my);
System.out.println("Diffusion coefficient: " + D);
System.out.println("-*-*-*-");
if(mx / 6 < InkSize) {
System.out.println("InkSize is too large");
}
}
}
In console, the output is as follows.
Water surface(Horizontal x vertical) : 10 cm × 10 cm
Ink size: 0.13000000000000003 cm^2
Time to hang:5000 msec
ending time: 10000 msec
Output unit time: 1000 msec
Unit distance: 0.1 cm
Number of horizontal meshes: 100
Number of vertical meshes: 100
Diffusion coefficient: 0.25
-*-*-*-
output: 3_0.ppm
output: 3_1.ppm
output: 3_2.ppm
output: 3_3.ppm
output: 3_4.ppm
output: 3_5.ppm
output: 3_6.ppm
output: 3_7.ppm
output: 3_8.ppm
output: 3_9.ppm
output: 3_10.ppm
done
The output image should look like the one below
Recommended Posts