head     1.11;
branch   ;
access   jkc john cartley dhale;
symbols  ;
locks    ; strict;
comment  @ * @;


1.11
date     94.04.06.10.30.55;  author jkc;  state Exp;
branches ;
next     1.10;

1.10
date     94.04.05.15.11.02;  author jkc;  state Exp;
branches ;
next     1.9;

1.9
date     94.04.05.10.50.03;  author jkc;  state Exp;
branches ;
next     1.8;

1.8
date     94.01.10.12.53.18;  author john;  state Exp;
branches ;
next     1.7;

1.7
date     93.11.24.15.40.09;  author john;  state Exp;
branches ;
next     1.6;

1.6
date     92.10.22.16.19.18;  author john;  state Exp;
branches ;
next     1.5;

1.5
date     91.09.05.08.50.50;  author jkcohen;  state Exp;
branches ;
next     1.4;

1.4
date     90.11.15.10.43.31;  author jkc;  state Exp;
branches ;
next     1.3;

1.3
date     90.06.12.10.44.42;  author jkc;  state Exp;
branches ;
next     1.2;

1.2
date     90.05.25.19.25.56;  author jkc;  state Exp;
branches ;
next     1.1;

1.1
date     90.04.18.20.35.04;  author jkc;  state Exp;
branches ;
next     ;


desc
@Create synthetic traces for line scatterer
@


1.11
log
@Fixed setting of trace fields tr.selev, tr.gelev and commented
obscure int truncation.
@
text
@/* SUIMP2D: $Revision: 1.10 $ ; $Date: 94/04/05 15:11:02 $	*/

/*----------------------------------------------------------------------
 * Copyright (c) Colorado School of Mines, 1990.
 * All rights reserved.
 *
 * This code is part of SU.  SU stands for Seismic Unix, a processing line
 * developed at the Colorado School of Mines, partially based on Stanford
 * Exploration Project (SEP) software.  Inquiries should be addressed to:
 *
 *  Jack K. Cohen, Center for Wave Phenomena, Colorado School of Mines,
 *  Golden, CO 80401  (jkc@@dix.mines.colorado)
 *----------------------------------------------------------------------
 */

#include "su.h"
#include "segy.h"

/*********************** self documentation **************/
char *sdoc[] = {
"							",
" SUIMP2D - generate shot records for a line scatterer	",
"           embedded in three dimensions using the Born	",
"	    integral equation				",							
"							",
" suimp2d [optional parameters] >stdout			",
"							",
" Optional parameters					",
"	nshot=1		number of shots			",
"	nrec=1		number of receivers		",
"	c=5000		speed				",
"	dt=.004		sampling rate			",
"	nt=256		number of samples		",
"	x0=1000		point scatterer location	",
"	z0=1000		point scatterer location	",
"	sxmin=0		first shot location		",
"	szmin=0		first shot location		",
"	gxmin=0		first receiver location		",
"	gzmin=0		first receiver location		",
"	dsx=100		x-step in shot location		",
"	dsz=0	 	z-step in shot location		",
"	dgx=100		x-step in receiver location	",
"	dgz=0		z-step in receiver location	",
"							",
NULL};

/* Credits:
 *	CWP: Norm, Jack
 */

/* Theory: Use the 3D Born integral equation (e.g., Geophysics,
 * v51, n8, p1554(7)). Use 2-D delta function for alpha and do
 * remaining y-integral by stationary phase.
 */
/**************** end self doc ***************************/


#define LOOKFAC	2	/* Look ahead factor for npfao	  */
#define PFA_MAX	720720	/* Largest allowed nfft	          */

segy tr;

main(int argc, char **argv)
{
	float c;		/* speed			*/
	float dt;		/* sampling rate		*/
	int nt;			/* number of samples		*/
	int ntsize;		/* ... in bytes			*/
	int nshot;		/* number of shots		*/
	int nrec;		/* number of receivers		*/
	float x0, z0;		/* point scatterer location	*/
	float sxmin, szmin;	/* first shot location		*/
	float gxmin, gzmin;	/* first receiver location	*/
	float dsx;		/* x-step in shot location	*/
	float dsz;		/* z-step in shot location	*/
	float dgx;		/* x-step in receiver location	*/
	float dgz;		/* z-step in receiver location	*/

	float sx, sz;		/* shot location		*/
	float gx, gz;		/* receiver location		*/
	float rs;		/* distance to shot		*/
	float rg;		/* distance to receiver		*/
	float d;		/* rs + rg			*/
	float t;		/* total travel time		*/
	int jt;			/* (int) arrival time on trace	*/
	float k;		/* constant part of response	*/

	register float *rt;	/* real trace			*/
	register complex *ct;	/* complex transformed trace	*/
	int nfft;		/* size of fft 			*/
	int nfftsize;		/* ... in bytes			*/
	int nfby2;		/* nfft/2			*/
	int nfby2p1;		/* nfft/2 + 1			*/
	int nzeros;		/* padded zeroes in bytes	*/
	float spread;		/* geometric spreading factor	*/
	float i32;		/* temp for omega to the 3/2	*/

	register int i;		/* counter			*/
	register int s;		/* shot counter			*/
	register int g;		/* receiver counter		*/



	/* Initialize */
	initargs(argc, argv);
	requestdoc(0);


	/* Get parameters */
	if (!getparint("nshot", &nshot))	nshot = 1;
	if (!getparint("nrec", &nrec))		nrec  = 1;
	if (!getparint("nt", &nt))		nt    = 256;
	if (!getparfloat("c", &c))		c     = 5000.0;
	if (!getparfloat("dt", &dt))		dt    = 0.004;
	if (!getparfloat("x0", &x0))		x0    = 1000.0;
	if (!getparfloat("z0", &z0))		z0    = 1000.0;
	if (!getparfloat("sxmin", &sxmin))	sxmin = 0.0;
	if (!getparfloat("szmin", &szmin))	szmin = 0.0;
	if (!getparfloat("gxmin", &gxmin))	gxmin = 0.0;
	if (!getparfloat("gzmin", &gzmin))	gzmin = 0.0;
	if (!getparfloat("dsx", &dsx))		dsx   = 100.0;
	if (!getparfloat("dsz", &dsz))		dsz   = 0.0;
	if (!getparfloat("dgx", &dgx))		dgx   = 100.0;
	if (!getparfloat("dgz", &dgz))		dgz   = 0.0;


	/* Set the constant header fields */
	tr.ns = nt;
	tr.dt = dt * 1000000.0;
	ntsize = nt * FSIZE;



	/* Set up for fft */
	nfft = npfaro(nt, LOOKFAC * nt);
	if (nfft >= MIN(SU_NFLTS, PFA_MAX))
		err("Padded nt=%d -- too big", nfft);

	nfby2 = nfft / 2;
	nfby2p1 = nfby2 + 1;
	nzeros = (nfft - nt) * FSIZE;


	/* Allocate fft arrays */
	rt   = ealloc1float(nfft);
	ct   = ealloc1complex(nfby2p1);


	/* Set the constant in the response amplitude
	   including scale for inverse fft below      */
	k = 1.0 / (4.0 * sqrt(2.0*c*dt*nfft) * c * dt * dt * nfft * nfft);


	/* Create the traces */
	for (s = 0; s < nshot; ++s) {	/* loop over shots */
		sx = sxmin + s * dsx;
		sz = szmin + s * dsz;
		tr.sx = sx;
		tr.selev = -sz; /* above sea level is positive */
		rs = sqrt((sx - x0)*(sx - x0) + (sz - z0)*(sz - z0));

		for (g = 0; g < nrec; ++g) {	/* loop over receivers */
			memset((void *) tr.data, (int) '\0', ntsize);
			gx = gxmin + g * dgx;
			gz = gzmin + g * dgz;
			tr.gx = gx;
			tr.gelev = -gz; /* above sea level is positive */
			rg = sqrt((gx - x0)*(gx - x0) + (gz - z0)*(gz - z0));
			d = rs + rg;
			t = d/c;
			jt = (int) t/dt; /* truncate to mesh */
			spread = sqrt(rs*rg*d);

			/* Distribute response over two surrounding samples */
			tr.data[jt] = k * ((jt + 1)*dt - t) / spread;
			tr.data[jt + 1] = k * (t - jt*dt) / spread;

			/* Load trace into rt (zero-padded) */
			memcpy((void *) rt, (const void *) tr.data, ntsize);
			memset((void *)(rt + nt), (int) '\0', nzeros);

			/* FFT */
			pfarc(1, nfft, rt, ct);

			/* Formula requires multiplication by
			   abs(omega) to 3/2 power times
			   exp(i pi / 4) sgn omega)    */
			for (i = 0; i < nfby2p1; ++i) {
				i32 = i * sqrt((double) i);
				ct[i] = crmul(ct[i], i32);
			}

			/* Invert and take real part */
			pfacr(-1, nfft, ct, rt);

			/* Load traces back in */
			memcpy( (void *) tr.data, (const void *) rt, ntsize);

			tr.gx = gx;
			tr.gy = gz;
			tr.tracl = s*nrec + g + 1;
			tr.offset = gx - sx;	/* for surface surveys */
			tr.cdp = 0.5*(gx + sx);	/* for surface surveys */

			puttr(&tr);
		} /* end loop on receivers */
	} /* end loop on shots */

	return EXIT_SUCCESS;
}
@


1.10
log
@Apparently the powers of nfft were wrong.  Fixed to the correct(?)
values.
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.9 $ ; $Date: 94/04/05 10:50:03 $	*/
d159 1
a159 1
		tr.sy = sz;
d166 2
d171 1
a171 1
			jt = t/dt;
@


1.9
log
@Added better description of what the code does.
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.8 $ ; $Date: 94/01/10 12:53:18 $	*/
d151 1
a151 1
	k = sqrt((double) nfft) / (4.0 * sqrt(2.0*c*dt) * c * dt * dt);
@


1.8
log
@removed an errant bzero
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.7 $ ; $Date: 93/11/24 15:40:09 $	*/
d23 2
a45 1
/**************** end self doc ***************************/
d51 5
d57 1
d95 1
a95 1
	float spread;		/* 2.5-D spreading factor	*/
@


1.7
log
@added explicit casts to memcpy ; changed bzero to memset
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.6 $ ; $Date: 92/10/22 16:19:18 $	*/
d171 1
a171 1
			bzero(rt + nt, nzeros);
@


1.6
log
@new selfdoc installed
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.5 $ ; $Date: 91/09/05 08:50:50 $	*/
d156 1
a156 1
			bzero(tr.data, ntsize);
d170 1
a170 1
			memcpy(rt, tr.data, ntsize);
d188 1
a188 1
			memcpy(tr.data, rt, ntsize);
@


1.5
log
@string -> String and so forth
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.4 $ ; $Date: 90/11/15 10:43:31 $	*/
d20 24
a43 24
String sdoc = "\
							\n\
SUIMP2D - generate shot records for a line scatterer	\n\
							\n\
suimp2d [optional parameters] >stdout			\n\
							\n\
Optional parameters					\n\
	nshot=1		number of shots			\n\
	nrec=1		number of receivers		\n\
	c=5000		speed				\n\
	dt=.004		sampling rate			\n\
	nt=256		number of samples		\n\
	x0=1000		point scatterer location	\n\
	z0=1000		point scatterer location	\n\
	sxmin=0		first shot location		\n\
	szmin=0		first shot location		\n\
	gxmin=0		first receiver location		\n\
	gzmin=0		first receiver location		\n\
	dsx=100		x-step in shot location		\n\
	dsz=0	 	z-step in shot location		\n\
	dgx=100		x-step in receiver location	\n\
	dgz=0		z-step in receiver location	\n\
							\n\
";
d99 1
a99 1
	askdoc(0);
@


1.4
log
@bcopy -> memcpy
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.3 $ ; $Date: 90/06/12 10:44:42 $	*/
d20 1
a20 1
string sdoc = "\
d103 15
a117 15
	if (!igetpar("nshot", &nshot))		nshot = 1;
	if (!igetpar("nrec", &nrec))		nrec  = 1;
	if (!igetpar("nt", &nt))		nt    = 256;
	if (!fgetpar("c", &c))			c     = 5000.0;
	if (!fgetpar("dt", &dt))		dt    = 0.004;
	if (!fgetpar("x0", &x0))		x0    = 1000.0;
	if (!fgetpar("z0", &z0))		z0    = 1000.0;
	if (!fgetpar("sxmin", &sxmin))		sxmin = 0.0;
	if (!fgetpar("szmin", &szmin))		szmin = 0.0;
	if (!fgetpar("gxmin", &gxmin))		gxmin = 0.0;
	if (!fgetpar("gzmin", &gzmin))		gzmin = 0.0;
	if (!fgetpar("dsx", &dsx))		dsx   = 100.0;
	if (!fgetpar("dsz", &dsz))		dsz   = 0.0;
	if (!fgetpar("dgx", &dgx))		dgx   = 100.0;
	if (!fgetpar("dgz", &dgz))		dgz   = 0.0;
@


1.3
log
@Using pfa
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.2 $ ; $Date: 90/05/25 19:25:56 $	*/
d170 1
a170 1
			bcopy(tr.data, rt, ntsize);
d188 1
a188 1
			bcopy(rt, tr.data, ntsize);
@


1.2
log
@vec --> ealloc1float, etc.
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.1 $ ; $Date: 90/04/18 20:35:04 $	*/
d51 2
a52 1
#define	FACMAX		12		/* For FFTPACK	*/
a53 1

d81 2
a82 6
	float *hcos;		/* hold cosines for fft		*/
	float *hsin;		/* hold sines for fft		*/
	register float *xr;	/* real part of trace		*/
	register float *xi;	/* imaginary part of trace	*/
	float *wr;		/* work area for fft		*/
	float *wi;		/* work area for fft		*/
d86 1
a86 2
	int nfac;		/* number of factors of nfft	*/
	int facs[FACMAX];	/* contains factors of nfft	*/
a87 2
	float xrold;		/* temporary for xr[i]		*/
	float xiold;		/* temporary for xi[i]		*/
d126 1
d128 4
a131 2
	fftfac(&nt, &nfft, &nfac, facs);
	if (nfft > SU_NFLTS)	err("Padded nt=%d -- too big", nfft);
d133 2
a134 2
	nfftsize = nfft * FSIZE;
	nzeros = nfftsize - ntsize;
d138 2
a139 3
	hcos = ealloc1float(nfft);  hsin = ealloc1float(nfft);
	xr   = ealloc1float(nfft);  xi   = ealloc1float(nfft);
	wr   = ealloc1float(nfft);  wi   = ealloc1float(nfft);
a141 4
	/* Set up tables */
	ffttab(&nfft, hcos, hsin);


d169 3
a171 4
			/* Load traces into xr, xi (zero-padded) */
			bcopy(tr.data, xr, ntsize);
			bzero(xr + nt, nzeros);
			bzero(xi, nfftsize);
d173 2
a174 2
			/* Fft */
			fft(xr, xi, wr, wi, &nfft, facs, &nfac, hcos, hsin);
d178 2
a179 6
			   exp(i pi / 4) sgn omega)
			   Take conjugate for inverse fft     */
			xi[0] = xr[0] = 0.0;
			for (i = 1; i <= nfby2; ++i) {
				xrold = xr[i];
				xiold = xi[i];
d181 1
a181 4
				xr[i] =  (xrold - xiold) * i32;
				xi[i] = -(xrold + xiold) * i32;
				xr[nfft - i] =  xr[i];
				xi[nfft - i] = -xi[i];
d185 4
a188 2
			fft(xr, xi, wr, wi, &nfft, facs, &nfac, hcos, hsin);
			bcopy(xr, tr.data, ntsize);
@


1.1
log
@Initial revision
@
text
@d1 1
a1 1
/* SUIMP2D: $Revision: 1.17 $ ; $Date: 89/09/20 19:36:21 $	*/
d142 3
a144 3
	hcos = vec(nfft);  hsin = vec(nfft);
	xr   = vec(nfft);  xi   = vec(nfft);
	wr   = vec(nfft);  wi   = vec(nfft);
@
