sureflpsvsh.c

su 的源代码库

	cl = alloc1float(nlayers+nrand_layers);
	ct = alloc1float(nlayers+nrand_layers);
	ql = alloc1float(nlayers+nrand_layers);
	qt = alloc1float(nlayers+nrand_layers);
	rho = alloc1float(nlayers+nrand_layers);
	t = alloc1float(nlayers+nrand_layers);
	lobs = alloc1int(nor+1);
	lobs[nor]=0;

	/*********************************************************************/
	/* read  input parameters from files or command line */
	if (*clfile !='\0') {			/* read from a file */	
		if ((infp=efopen(clfile,"r"))==NULL)
			err("cannot open file of pwave velocities=%s\n",clfile);
		fread(cl,sizeof(float),nlayers,infp);
		efclose(infp);
	} else getparfloat("cl",cl);		/* get from command line */
	if (*qlfile !='\0') {
		if ((infp=efopen(qlfile,"r"))==NULL)
			err("cannot open file of compressional Q=%s\n",qlfile);
		fread(ql,sizeof(float),nlayers,infp);
		efclose(infp);
	} else getparfloat("ql",ql);
	if (*ctfile !='\0') {
		if ((infp=efopen(ctfile,"r"))==NULL)
			err("cannot open file of swave velocities=%s\n",ctfile);
		fread(ct,sizeof(float),nlayers,infp);
		efclose(infp);
	} else getparfloat("ct",ct);
	if (*qtfile !='\0') {
		if ((infp=efopen(qtfile,"r"))==NULL)
			err("cannot open file of shear Q=%s\n",qtfile);
		fread(qt,sizeof(float),nlayers,infp);
		efclose(infp);
	} else getparfloat("qt",qt);
	if (*rhofile !='\0') {
		if ((infp=efopen(rhofile,"r"))==NULL)
			err("cannot open file of densities=%s\n",rhofile);
		fread(rho,sizeof(float),nlayers,infp);
		efclose(infp);
	} else getparfloat("rho",rho);
	if (*tfile !='\0') {
		if ((infp=efopen(tfile,"r"))==NULL)
			err("cannot open file of thicknesses=%s\n",tfile);
		fread(t,sizeof(float),nlayers,infp);
		efclose(infp);
	} else getparfloat("t",t);
	if (*lobsfile !='\0') {
		if ((infp=efopen(lobsfile,"r"))==NULL)
			err("can't open file of receiver layers=%s\n",lobsfile);
		fread(lobs,sizeof(int),nor,infp);
		efclose(infp);
	} else getparint("lobs",lobs);

	/*********************************************************************/
	/* if requested, do interpolation and/or parameter adjustment */
	if (nlint!=0)
		parameter_interpolation (nlayers, intlayers, nintlayers, 
				intlayth, cl, ql, ct, qt, rho, t);	

	/* if requested, compute random velocity layers */
	if (rand==1) {
		random_velocity_layers (&nlayers, &lsource, nrand_layers, sdcl,
			sdct, layer, zlayer, cl, ql, ct, qt, rho, t);
	}

	/* if requested, apply earth flattening approximation */
	if (flt==1) {
		apply_earth_flattening (nlayers, z0, cl, ct, rho, t);
	}


	/*********************************************************************/
	/* get filter parameters */
	if (*filtypefile !='\0') {
		if ((infp=efopen(filtypefile,"r"))==NULL)
			err("cannot open file=%s\n",filtypefile);
		getparint("nfilters",&nfilters);
		filters_type=alloc1int(nfilters);
		fread (filters_type,sizeof(int),nfilters,infp);
		efclose(infp);
	} else {
		nfilters=countparval("filters_type");
		filters_type=alloc1int(nfilters);
		getparint("filters_type",filters_type);
	}
	if (*fphfile !='\0') {
		if ((infp=efopen(fphfile,"r"))==NULL)
			err("cannot open file=%s\n",fphfile);
		filters_phase=alloc1int(nfilters);
		fread (filters_phase,sizeof(float),nfilters,infp);
		efclose(infp);
	} else if (nfilters == countparval("filters_phase")) {
		filters_phase=alloc1int(nfilters);
		getparint("filters_phase",filters_phase);
	} else err("number of elements infilterstype and phase must be equal");
	if (*dbpofile !='\0') {
		if ((infp=efopen(dbpofile,"r"))==NULL)
			err("cannot open file=%s\n",dbpofile);
		dbpo=alloc1float(nfilters);
		fread (dbpo,sizeof(float),nfilters,infp);
		efclose(infp);
	} else if (nfilters == countparval("dbpo")) {
		dbpo=alloc1float(nfilters);
		getparfloat("dbpo",dbpo);
	} else err("number of elements in filters_type and dbpo must be equal");
	if (*f1file !='\0') {
		if ((infp=efopen(f1file,"r"))==NULL)
			err("cannot open file=%s\n",f1file);
		f1=alloc1float(nfilters);
		fread (f1,sizeof(float),nfilters,infp);
		efclose(infp);
	} else if (nfilters == countparval("f1")) {
		f1=alloc1float(nfilters);
		getparfloat("f1",f1);
	} else err("number of elements in filters_type and f1 must be equal");
	if (*f2file !='\0') {
		if ((infp=efopen(f2file,"r"))==NULL)
			err("cannot open file=%s\n",f2file);
		f2=alloc1float(nfilters);
		fread (f2,sizeof(float),nfilters,infp);
		efclose(infp);
	} else if (nfilters == countparval("f2")) {
		f2=alloc1float(nfilters);
		getparfloat("f2",f2);
	} else err("number of elements in filters_type and f2 must be equal");
		

	/*********************************************************************/
	/* allocate space for wavefield computations */
	wavefield1=alloc2float(nt,nx);
	if (wtype==1) {
		wavefield2=alloc2float(nt,nx);
		wavefield3=alloc2float(nt,nx);
	}
	/* get name of output file for processing information */
	if (verbose==2||verbose==3) {
		if (!getparstring("outf",&outf))	outf="info";
		if ((outfp=efopen(outf,"w"))==NULL) {
			warn("cannot open processing file =%s, no processing\n"
			"information file will be generated\n",outf);
			verbose=1;
		}
	}

	/* initialize wavefields */
	if (wtype==1) {
		for (ix=0;ix<nx;ix++) {
			for (it=0;it<nt;it++) {
				wavefield1[ix][it]=0.0;
				wavefield2[ix][it]=0.0;
				wavefield3[ix][it]=0.0;
			}
		}
	} else if (wtype==2) {
		for (ix=0;ix<nx;ix++) {
			for (it=0;it<nt;it++) {
				wavefield1[ix][it]=0.0;
			}
		}
	}

	/* number of time samples in computed traces */
	ntc=tsec/dt;
	if (int_type==2) bp=0.0;

	/*********************************************************************/
	/* Now, compute the actual reflectivities */
	compute_reflectivities (int_type, verbose, wtype, wfield, vsp, flt,
		win, nx, nt, ntc, nor, nf, nlayers, lsource, layern, nfilters,
		filters_phase, nw, np, bp, tlag, red_vel, w1, w2, fx, dx, bx,
		fs, decay, p2w, tsec, fref, wrefp, wrefs, epsp, epss, sigp,
		sigs, pw1, pw2, pw3, pw4, h1, h2, m1, m2, m3, fref, lobs,
		filters_type, dbpo, f1, f2, cl, ct, ql, qt, rho, t, wavefield1,
		wavefield2, wavefield3, outfp);
	/*********************************************************************/

	/* if open, close processing information file */
	if (verbose==2||verbose==3) efclose(outfp);

	/* convolve with a wavelet and write the results out */
	if (wtype==1) {			/* PSV */
		
		/* convolve with a wavelet to produce the seismograms */
		convolve_wavelet (wavelet_type, nx, nt, dt, fpeak, wavefield1); 
		convolve_wavelet (wavelet_type, nx, nt, dt, fpeak, wavefield2); 
		convolve_wavelet (wavelet_type, nx, nt, dt, fpeak, wavefield3); 

		/* output results in SU format */
		if(*wfp!='\0'){
			if ((wfp_file=efopen(wfp,"w"))==NULL)
				err("cannot open pressure file=%s\n",wfp);
			{	register int ix;
				for (ix=0; ix<nx; ix++) {
					for (it=0; it<nt; it++)
						tr1.data[it]=wavefield1[ix][it];

					/* headers*/
					tr1.ns=nt;
					tr1.dt=1000*(int)(1000*dt);
					tr1.offset=(bx+ix*dx)*1000;
	
					/* output trace */
					fputtr(wfp_file, &tr1);
				}
				efclose (wfp_file);
			}
		}
		if (*wfr !='\0') {
			if ((wfr_file=efopen(wfr,"w"))==NULL)
					err("cannot open radial wfield file=%s\n",wfr);
			{	register int ix;
				for (ix=0; ix<nx; ix++) {
					for (it=0; it<nt; it++)
						tr2.data[it]=wavefield2[ix][it];
					tr2.ns=nt;
					tr2.dt=1000*(int)(1000*dt);
					tr2.offset=(bx+ix*dx)*1000;
					fputtr(wfr_file, &tr2);
				}
				efclose (wfr_file);
			}
		}
		if (*wfz !='\0') {
			if ((wfz_file=efopen(wfz,"w"))==NULL)
				err("canno open vertical field file=%s\n",wfz);
			{	register int ix;
				for (ix=0; ix<nx; ix++) {
					for (it=0; it<nt; it++)
							tr3.data[it]=wavefield3[ix][it];
					tr3.ns=nt;
					tr3.dt=1000*(int)(1000*dt);
					tr3.offset=(bx+ix*dx)*1000;
					fputtr(wfz_file, &tr3);
				}
				efclose (wfz_file);
			}
		}
		
		/* free allocated space */
		free2float(wavefield1);
		free2float(wavefield2);
		free2float(wavefield3);

	} else if (wtype==2) {			/* SH */

		/* convolve with a wavelet to produce the seismogram */
		convolve_wavelet (wavelet_type, nx, nt, dt, fpeak, wavefield1); 

		/* output the result in SU format */
		if (*wft !='\0') {
			if ((wft_file=efopen(wft,"w"))==NULL)
				err("cannot open tangential file=%s\n",wft);
			{	register int ix;
				for (ix=0; ix<nx; ix++) {
					for (it=0; it<nt; it++)
							tr1.data[it]=wavefield1[ix][it];
					tr1.ns=nt;
					tr1.dt=1000*(int)(1000*dt);
					tr1.offset=(bx+ix*dx)*1000;
					fputtr(wft_file, &tr1);
				}
				efclose (wft_file);
			}
		}

		/* free allocated space */
		free2float(wavefield1);
	}

	/* free workspace */
	free1float(cl);
	free1float(ct);
	free1float(ql);
	free1float(qt);
	free1float(rho);
	free1float(t);
	free1int(lobs);
	free1int(filters_type);
	free1int(filters_phase);
	free1float(dbpo);
	free1float(f1);
	free1float(f2);
	return EXIT_SUCCESS;
}