/Forward modeling/
#include<stdio.h>
#include<math.h>
#include<graphics.h>
#include<conio.h>
#include<windows.h>
#include<stdlib.h>
void fun();
int main (void)
{

/*Background board*/
initgraph(600,600);       
setbkcolor(WHITE);         
setlinecolor(BLACK);       
cleardevice();             
setorigin(2,598);

fun();
return 0;

}
void fun()
{
/Initializing Variables/
int x,y,t;
double i;
double j;
double dx,dy,dt,dh,v;
double b[179][179];
double constant;
double ricker = 100;
double m,n;
int m1,n1;
double (p)[179][179]=(double()[179][179])malloc(sizeof(double)175179179);
v = 0.5;
dt = dx = dy = dh = 1;
constant = (vdt/dh)(vdt/dh);

/*Gird depiction*/
for(j=2.0;j<=587.4;j+=3.3)
{
	line(2.0,-j,586.1,-j);   
	line(j, -586.1,j,-2.0);       
}


/*Boundary initialization*/
for(t=0;t<175;t++)
{
	for(x=0;x<179;x++)
	{
		for(y=0;y<179;y++)
		{
			p[t][x][y] = 0;//the value of 'p' is 0 on every moment
		}
	}
}

/*Ricker initialization*/
p[0][89][89] = ricker;//t = -1
p[1][89][89] = ricker;//t = 0

/*Propagation of assignment(interation)*/
for(t=1;t<173;t++)
{

	/*(89,89)—(177,177)对角线矩形网格区域*/
	for(x=0;x<89;x++)
	{
			for(y=0;y<89;y++)
			{
				p[t+1][x+89][y+89] = constant*(p[t][(x+89)+1][y+89] + p[t][(x+89)-1][y+89] + p[t][x+89][y+89+1] + p[t][x+89][y+89-1]) + (2-4*constant)*p[t][x+89][y+89] - p[t-1][x+89][y+89];
			}
	}
	/*(89,1)—(177,88)对角线矩形网格区域*/
	for(x=0;x<89;x++)
	{
			for(y=-88;y<0;y++)
			{
				p[t+1][x+89][y+89] = constant*(p[t][(x+89)+1][y+89] + p[t][(x+89)-1][y+89] + p[t][x+89][y+89+1] + p[t][x+89][y+89-1]) + (2-4*constant)*p[t][x+89][y+89] - p[t-1][x+89][y+89];
			}
	}
	/*(89,89)—(1,1)对角线矩形网格区域*/
	for(x=-88;x<1;x++)
	{
			for(y=-88;y<1;y++)
			{
				p[t+1][x+89][y+89] = constant*(p[t][(x+89)+1][y+89] + p[t][(x+89)-1][y+89] + p[t][(x+89)][y+89+1] + p[t][(x+89)][y+89-1]) + (2-4*constant)*p[t][(x+89)][y+89] - p[t-1][x+89][y+89];
			}
	}
	/*(89,89)—(1,177)对角线矩形网格区域*/
	for(x=-88;x<1;x++)
	{
			for(y=0;y<89;y++)
			{
				p[t+1][x+89][y+89] = constant*(p[t][(x+89)+1][y+89] + p[t][(x+89)-1][y+89] + p[t][(x+89)][y+89+1] + p[t][(x+89)][y+89-1]) + (2-4*constant)*p[t][(x+89)][y+89] - p[t-1][x+89][y+89];
			}
	}
}

/*Assign a value to a valid point*/
for(x=0;x<177;x++)
{
	for(y=0;y<177;y++)
	{
		if(p[170][x][y]<=0)
		{
			b[x][y] = 0;
			/*p[t+1][x][y] = 0;*/
		}
		else
		{
			b[x][y] = y;
			/*p[t+1][x][y] = x;*/
		}
	}
}
/*Draw point*/
for(m=2;m<587.4;m+=3.3)
{
	for(n=2;n<587.4;n+=3.3)
	{
		m1 = (int)(m-2)/3.3;
		n1 = (int)(n-2)/3.3;
		putpixel(m1*3.3,-3.3*b[m1][n1],LIGHTGREEN);
	}
}
system("pause");
free(p);
/*Yes ! get it.*/

}