地震有限差分正演程序

//////////////////头文件////////////////////////////////////

#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