//////////////////頭文件////////////////////////////////////

#include?

#include?

#include?

////////////////開辟內存數組//////////////////////

void?**alloc2?（int?n1?int?n2?int?size）
{
???????
	int?i2;
???????
	void?**p;

????????
	if?（（p=（void**）malloc（n2*sizeof（void*）））==NULL）
??????????????
		return?NULL;
????????
	if?（（p[0]=（void*）malloc（n2*n1*size））==NULL）?{
???????????????
		free（p）;
???????????????
		return?NULL;
???????
	}
????????
	for?（i2=0;?i2?????????????
		p[i2]?=?（char*）p[0]+size*n1*i2;
????????
	return?p;

}

///////////////////子波程序////////////////////////////

void?tforce_ricker（int?nt?float?*tforce?float?dt?float?fpeak）
{

	int	it;

	float?PI=3.1415926;
	
	float	t1?t0;?
	
	t0=1.0/fpeak;
	
	for?（it=0;?it		
		t1=it*dt;
		
		tforce[it]?=?exp（-PI*PI*fpeak*fpeak*（t1-t0）*（t1-t0））*（1.0-2.*PI*PI*fpeak*fpeak*（t1-t0）*（t1-t0））;
	
	}

}



main（）
{

	int?NxNz;////////////////模型大小//////////////////

	float?**Vp;//////////////////速度程序////////////////////

	int?ixiz;////////////////循環變量////////////////////////////

	int?Ntit;//////////////時間采樣點與時間變量/////////////////

	int?f;//////////////////////頻率變量///////////////////

	float?dt;//////////////////時間長采樣率///////////////////////////////////////////

	float?*fsouce;////////////////////震源子波函數//////////////////////

	float?o0o1o2o3o4o5;///////////////////////有限差分系數////////////////////

	float?dxdz;//////////////////空間采樣率////////////////////////

	float?**U**U_left**U_right;/////////////////u（ijk）（n+1）=2*u（ijk）（n）-u（ijk）（n-1）;U為u（ijk）（n）U_left為u（ijk）（n-1）U_right為u（ijk）（n+1）

	int?source_xsource_z;/////////////////震源位置/////////////////////////////

	FILE?*fp;

	o1=3.33;

	o2=-0.47;

	o3=0.07;

	o4=-0.01;

	o5=0;

	o0=-2*（o1+o2+o3+o4+o5）;


	Nx=100;

	Nz=100;

	Nt=1000;

	f=30;

	dt=0.001;

	dx=10;

	dz=10;

	source_x=Nx/2;

	source_z=Nz/2;

	////////////開辟內存////////////////////



	Vp?=?（float?**）alloc2（NxNzsizeof（float））;

	///////////////初始化/////////////////////////////////



	for（ix=0;ix
		for（iz=0;iz
			Vp[iz][ix]=0;

		}

	}

		////////////開辟內存////////////////////

	

	fsouce=（float?*）malloc（Nt*sizeof（float））;

	///////////////初始化/////////////////////////////////

	for（it=0;it
		fsouce[it]=0.0;

	}

		////////////開辟內存////////////////////

	U=（float?**）alloc2（Nx+10Nz+10sizeof（flo