// 包括库	--- Include libraries
#include <stdio.h>
#include <math.h>
#include <gtk/gtk.h>            // GUI, Gtk Lib


// 绘制时空浓度输出的像素贴图尺寸	--- Pixmap dimensions to paint the space-time-concentration output
#define X_SIZE 256
#define Y_SIZE 512

// 状态变量的范围			--- Range of the state variable
// 浓度最大值				--- Max value of concentratiom
#define MAX_CONCENTRATION 100
// 浓度最小值                           --- Min value of concentration
#define MIN_CONTENTRATION 0

// 彩色贴图渐变点数			--- Color map gradient number of points
#define N_GRADIENT_POINTS 1000 //



// 定义		--- Definitions
#define N X_SIZE 		// 空间网格数	--- Number of spatial grid points
#define T Y_SIZE 		// 时间步数   	--- Time steps to be simulated (integrated)
#define dt 1 		// 时间步长	--- Temporal length scale (dt)
#define dx 1 		// 空间步长	--- Spatial length scale (dx)
#define diffusion_rate 0.45 	// 扩散系数	--- Diffusion coefficient



// 全球结构	--- Global structures


// 颜色贴图渐变	--- color-map gradient
struct color_gradient{
    int red[(int) N_GRADIENT_POINTS], green[(int) N_GRADIENT_POINTS], blue[(int) N_GRADIENT_POINTS];
} color;


// 模拟数据	--- simulation data
struct simulation_data{
	double u[N];	// u: 为当前时刻的浓度分
			//      --- To store the concentration distribution at the current moment
	int generation;	// generation: 模拟迭代	--- simulation iteration
	double time_buffer[N][T]; // store T denerations (dt) to display as a pixmap
} state;


// 函数声明	--- Function declarations


// Declare PUT PIXEL function to access individual pixel data on a Pixel Buffer. 
//Implemented at the end of documment under main block
void put_pixel(GdkPixbuf *pixbuf, int x, int y, guchar red, guchar green, guchar blue, guchar alpha);
// Function to actually paint the temporal Buffer to a pixmap
static void paint_time_buffer (gpointer data);

// Constructure-like function to make/create/construct a color structure
// n_points can be at MAX the same as N_GRADIENT_POINTS
struct color_gradient makeColorGradient(int n_points){
    struct color_gradient color;
    for(int value = 1; value <= n_points;value ++)
        {
        float ratio;
        float min = (float) 1;
        float max = (float) n_points;
        ratio = 2 * (value-min) / (max - min);
        if (ratio < 1 ){color.blue[value] = (int)(255*(1 - ratio));}else{color.blue[value]=0;};
        if (ratio > 1 ){color.red[value] = (int)(255*(ratio - 1));}else{color.red[value]=0;};
        color.green[value] = 255 - color.blue[value] - color.red[value];
        printf("making a gradient with %d points\n", n_points);
        printf("%d  --->  [ %d %d %d ] ratio = %f\n", value,color.red[value],color.green[value],color.blue[value], ratio);
        }
    return color;
}




// 打印颜色渐变状态	--- Print color-gradient function
void printGradient(void){
// 根据状态变量计算梯度值	--- Calculate gradient value from state variable
float Delta, delta, value;
Delta = (float) (MAX_CONCENTRATION - MIN_CONTENTRATION);
delta = Delta / (float)N_GRADIENT_POINTS;
	// 代	--- generations
	printf("[color value:\t");
	// 最左边3个格点的状态	--- state of the 3 left-most lattice sites
	for (int i = 0; i < 3;i++){
		if (state.u[i] != MIN_CONTENTRATION)
			{value = ceil((state.u[i]-((float)MIN_CONTENTRATION))/delta);}
			else{value = 1;}
		printf("%d ", (int) value);
                 printf("\t");
		}
	printf("...\t");
	// 三个中心格点的状态	--- state of the 3 center lattice sites
        for (int i = floor(N/2)-1 ; i < floor(N/2)+2;i++){
			if (state.u[i] != MIN_CONTENTRATION)
			{value = ceil((state.u[i]-((float)MIN_CONTENTRATION))/delta);}
			else{value = 1;}
                 printf("%d ", (int) value);
                 printf("\t");
                }
	printf("...\t");
	// 3个最右边格点的状态	--- state of the 3 right-most lattice sites
	for (int i = N-3; i < N;i++){
		if (state.u[i] != MIN_CONTENTRATION)
			{value = ceil((state.u[i]-((float)MIN_CONTENTRATION))/delta);}
			else{value = 1;}
                 printf("%d ", (int) value);
                 printf("\t");
		}
	printf("\n");
}



// 打印功能	--- Print function (truncated)
void printGrid(void){
	// 代	--- generations
	printf("[dt x  %d]\t", state.generation);
	// 最左边3个格点的状态	--- state of the 3 left-most lattice sites
	for (int i = 0; i < 3;i++){
		 printf("%.2f ", state.u[i]);
                 printf("\t");
		}
	printf("...\t");
	// 三个中心格点的状态	--- state of the 3 center lattice sites
        for (int i = floor(N/2)-1 ; i < floor(N/2)+2;i++){
                 printf("%.2f ", state.u[i]);
                 printf("\t");
                }
	printf("...\t");
	// 3个最右边格点的状态	--- state of the 3 right-most lattice sites
	for (int i = N-3; i < N;i++){
                 printf("%.2f ", state.u[i]);
                 printf("\t");
		}
	printf("\n");
}



// 初始化功能	--- Initialization function
void initLattice(void){
	for (int  t = 0; t < T; t++) for (int i = 0; i < N; i++) state.time_buffer[i][t]=0.0;
	// 初始化浓度分布	--- Initialize concentration distribution
	for (int i = 0; i < N; i++) {
        	if (i >= 45*N/100 && i <= 55*N/100) {
		state.u[i] = (float) MAX_CONCENTRATION;
        	} else {
            	state.u[i] = 0.0;
        	}
    	}
	for (int i = 0; i < N; i++) state.time_buffer[i][0]= state.u[i];
	// 初始化时间		--- Initialize time
    	state.generation = 0;
}




// 更新功能	--- Update function (to increment one iteration/generation)
void updateLattice(void){
	// u_new： 为下一个时刻的浓度分布
	double u_new[N];//      --- To store the concentration distribution at the next moment
	// 强制边界条件		--- Enforce boundary condition
	u_new[0] = state.u[0];
    u_new[N] = state.u[N];
	// 内部节点计算         --- Internal node calculation
    for (int i = 1; i < N; i++) {
        u_new[i] = state.u[i] + (double)diffusion_rate * dt / (dx*dx) * (state.u[i+1] - 2*state.u[i] + state.u[i-1]);
        }
	// 更新浓度分布         --- Update concentration distribution
    for (int i = 0; i < N; i++) {
		state.u[i] = u_new[i];
        }
    state.generation ++;
}



// 推送时间缓冲		--- Push time buffer
void pushTimeBuffer(int time){
	for (int i=0; i < N; i++) state.time_buffer[i][time]=state.u[i];
}






// 模拟	--- simulation
void integrateSimulation(void){
	// 初始化模拟	--- Initialize simulation
    initLattice();
	printGrid();
	printGradient();
	// 迭代计算		--- Iterative calculation
    for (int t = 1; t < T; t++) {
        updateLattice();
		pushTimeBuffer(t);
		printGrid();
		printGradient();
        }
}



// Activate Gtk function
static void activate (GtkApplication *app, gpointer user_data){
	integrateSimulation();
	// 声明小工具	--- declare Widgets
	GtkWidget *window, *image_space_time_concentration;
    // Gdk Pixel bufer to create an image
    GdkPixbuf *pixbuf;
    // 创建小部件：一个新的窗口	--- Create widgets: a new WINDOW
    window = gtk_application_window_new (app);
    // 设置窗口标题	--- Set window title
	gtk_window_set_title (GTK_WINDOW (window), "Diffusion (FTCS)");
	// 使窗口不可调整大小	--- Make window not resizable
	gtk_window_set_resizable (GTK_WINDOW(window), FALSE);
	// 创建像素缓冲区		--- Create Pixel buffer
	pixbuf = gdk_pixbuf_new(GDK_COLORSPACE_RGB, 0, 8, X_SIZE, Y_SIZE);
	// 从像素缓冲区创建图像	--- Create Image from empty Pixel buffer
    image_space_time_concentration = gtk_image_new_from_pixbuf(pixbuf);
	// 绘制图像				--- Paint Image (use time-buffer data) 
	paint_time_buffer(image_space_time_concentration);
	// 将图像添加到窗口容器	--- Add Image to window container
	gtk_container_add (GTK_CONTAINER (window), image_space_time_concentration);
	// 显示窗口及其小工具	--- Show the window and its Widgets
    gtk_widget_show_all (window);
}




// 主要功能	--- Main function
int main(int    argc, char **argv){
    GtkApplication *app;
    int status;
    app = gtk_application_new ("keymer.lab.yuechuan.diffusion", G_APPLICATION_FLAGS_NONE);
    g_signal_connect (app, "activate", G_CALLBACK (activate), NULL);
    status = g_application_run (G_APPLICATION (app), argc, argv);
    g_object_unref (app);
    return status;
}



static void paint_time_buffer (gpointer data){
	// 创建像素缓冲区		--- Create Pixel buffer
	GdkPixbuf *p;
	p = gdk_pixbuf_new (GDK_COLORSPACE_RGB, 0, 8, X_SIZE, Y_SIZE);
	// 制作颜色渐变			--- Make a color gradient
	color = makeColorGradient((int) N_GRADIENT_POINTS);
	// 根据状态变量计算梯度值	--- Calculate gradient value from state variable
	float Delta, delta, value;
	Delta = (float) (MAX_CONCENTRATION - MIN_CONTENTRATION);
	delta = Delta / (float)N_GRADIENT_POINTS;
  	// 绘制（使用put_pixel)		--- Paint (using put_pixel)	
 	int x, y;
  	for (x = 0; x < X_SIZE; x++)     // 空间 --- Space
        {
        for (y = 0; y < Y_SIZE; y++) // 时间 --- Time
            {
			// 从状态获取梯度值		--- get gradient value from state
			if (state.time_buffer[x][y] != MIN_CONTENTRATION)
					{value = ceil((state.time_buffer[x][y]-((float)MIN_CONTENTRATION))/delta);}
					else{value = 1;}
			// 油漆					--- paint
            put_pixel (p, (int) x, (int) y,
                       (guchar) color.red[(int)value], 
					   (guchar) color.green[(int)value], 
					   (guchar) color.blue[(int)value], 
					   255);
            }
      }
	// 从缓冲区创建图像		--- Create image (data) from buffer
	gtk_image_set_from_pixbuf (GTK_IMAGE (data), GDK_PIXBUF (p));
	// 释放缓冲存储器		--- Release buffer memory
	g_object_unref (p);
};



// put_pixel的实现			--- Implementation of put_pixel. 
// Thanks to code from:
// https://developer.gnome.org/gdk-pixbuf/stable/gdk-pixbuf-The-GdkPixbuf-Structure.html
void put_pixel(GdkPixbuf *pixbuf, int x, int y, guchar red, guchar green, guchar blue, guchar alpha){
    guchar *pixels, *p;
    int rowstride, numchannels;
    numchannels = gdk_pixbuf_get_n_channels(pixbuf);
    rowstride = gdk_pixbuf_get_rowstride(pixbuf);
    pixels = gdk_pixbuf_get_pixels(pixbuf);
    p = pixels + y * rowstride + x * numchannels;
    p[0] = red;     p[1] = green; p[2] = blue; p[3] = alpha;
}