二维弹性介质波动方程数值模拟动态图像

/*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(BLACK);         
	setlinecolor(WHITE);       
	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)*175*179*179);
	v = 0.5;
	dt = dx = dy = dh = 1;	
	constant = (v*dt/dh)*(v*dt/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=2;t<170;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/12)*(16*(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]) - (p[t][(x+89)+2][y+89] + p[t][(x+89)-2][y+89] + p[t][x+89][y+89+2] + p[t][x+89][y+89-2])) + (2-5*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/12)*(16*(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]) - (p[t][(x+89)+2][y+89] + p[t][(x+89)-2][y+89] + p[t][x+89][y+89+2] + p[t][x+89][y+89-2])) + (2-5*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/12)*(16*(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]) - (p[t][(x+89)+2][y+89] + p[t][(x+89)-2][y+89] + p[t][x+89][y+89+2] + p[t][x+89][y+89-2])) + (2-5*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/12)*(16*(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]) - (p[t][(x+89)+2][y+89] + p[t][(x+89)-2][y+89] + p[t][x+89][y+89+2] + p[t][x+89][y+89-2])) + (2-5*constant)*p[t][x+89][y+89] - p[t-1][x+89][y+89];
				}
		}
	}
	for(t=4;t<170;t++)
	{
		/*Assign a value to a valid point*/
		for(x=0;x<177;x++)
		{
			for(y=0;y<177;y++)
			{
				if(p[t][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],WHITE);
			}
		}

		Sleep(60);

		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],BLACK);
			}
		}
	}
	system("pause");
	free(p);
	/*Yes ! get it.*/
}

你也可以认为这是水波，总之是二维各向同性介质中的波传播情况。