test5-c press.txt

地球物理大师Claebout的几个有用的程序

#include "math.h"
#include "stdio.h"
#include<stdlib.h>
#include <time.h>

void main()
{
	int it,nt,ih,nh,is,ns,iz,nz,it0,iy; 
	double p[512],b[512],refl[25][16],z[25],geol[25],random;

	FILE *fp1=fopen("plot.txt","w+");
	FILE *fp2=fopen("plot.dat","w+");
	nt=512;
	nh=48;
	ns=10;
	nz=25;

	srand((unsigned)time(NULL));

	for ( iz=0;iz<nz;iz++ )		// reflector depth
	{	
		random=1.0*(double)rand()/(double)RAND_MAX;
		z[iz]=nt*random;		// random() is on the interval (0,32767)
//		printf("%d  %e  %e\n",iz+1,random,z[iz]);
	}

	for ( iz=0;iz<nz;iz++ )		// reflector strength with depth
	{
		random=1.0*(double)rand()/(double)RAND_MAX;
		geol[iz]=2.0*random-1.0;
//		printf("%d,%e,%e\n",iz+1,random,geol[iz]);
	}
	
	for ( is=0;is<ns;is++ )		// give texture to the geology
		for ( iz=0;iz<nz;iz++ ) //ns==10(!=16)
		{	
			random=1.0*(double)rand()/(double)RAND_MAX;
			refl[iz][is]=(1.0+random)*geol[iz];//saved as colmus
		}
/*	for ( ;is<16;is++ )		
		for ( iz=0;iz<nz;iz++ )
			refl[iz][is]=0.0;
*/	
	for ( it=0;it<nt;it++ )		// prepare a wavelet
		b[it]=exp(-(it+1)*0.08)*sin(0.5*(it+1)-0.5);

	for ( is=ns;is>0;is-- )		// shots run backwards
		for ( ih=0;ih<nh;ih++ )	// down cable h=(g-s)/2
		{	
			iy=(is-1)+ih;		// y=midpoint
			iy=1+(iy-ns*(iy/ns));	// periodic with midpoint
			for ( it=0;it<nt;it++ )
				p[it]=0.0;
			for ( iz=0;iz<nz;iz++ )// add in a hyperbola for each layer
			{	it0=(int)sqrt(pow(z[iz],2.0)+100.0*pow(ih,2.0));
				for ( it=0;it<nt-it0;it++ )	// add in the wavelet
					p[it+it0]=p[it+it0]+refl[iz][iy]*b[it];
			}
			for ( it=0;it<nt;it++ )
				fprintf(fp1,"%d  %e\n",it+1,p[it]);
			fwrite(p,sizeof(double),nt,fp2);
		}
	fclose(fp1);
	fclose(fp2);
}