sudatumk2dr.c

su 的源代码库

        mtmax = 2*dxgo*sin(angmax*PI/180.)/(v0*dt);
        if(mtmax<1) mtmax = 1;

	jtr = 1;
	s = fxso;
	ktr = 0;
	fxin = fxgo;
	for(is=0; is<nxso; ++is)
		ng[is]=0;

	/* Store headers in tmpfile while getting a count */
	hdrfp = etmpfile(); 

	do {

	    float as,res;
	    int is;
	    sx = tr.sx;
	    gx = tr.gx;

	    if(sx!=s) {
		s = sx;
                if(gx>fxgo) fxin = gx;
	    }
	
	    io = (int)((sx-fxso)/dxso);

	    if(MIN(sx,gx)>=fs && MAX(sx,gx)<=es && 
	       MAX(gx-sx,sx-gx)<=ABS(offmax) ){


                /* Number of traces in receiver's gathers  */
                ng[io]++;

                efwrite(&tr, 1, HDRBYTES, hdrfp);

		/*  Down or up ward continuation of the receivers  */
	    	as = (gx-fs)/ds;
	    	is = (int)as; 

		if(is==ns-1) is=ns-2;
		res = as-is;
		if(res<=0.01) res = 0.0;
		if(res>=0.99) res = 1.0;
		sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tsum);
                sum2(nr,nzt,1-res,res,pb[is],pb[is+1],pbg);
                sum2(nr,nzt,1-res,res,tb[is],tb[is+1],tbg);
                sum2(nr,nzt,1-res,res,sigb[is],sigb[is+1],sigbg);
                sum2(nr,nzt,1-res,res,angb[is],angb[is+1],angbg);

                as = (gx-fxi)/dxi;
                is = (int)as;
                if(is==nxi-1) is=nxi-2;
                res = as-is;
                if(res<=0.01) res = 0.0;
                if(res>=0.99) res = 1.0;
                nanglg = (1-res)*nangl[is] + res*nangl[is+1];

		dat2d(tr.data,nt,ft,dt,sx,gx,dats[io],aperx,nxgo,fxin,
		      dxgo,nzo,fzo,dzo,mtmax,gx,fmax,nxi,fxi,dxi,angmax,tbg,
                      pbg,angbg,nr,tsum,nzt,fzt,dzt,nxt,fxt,dxt,
                      antiali,sgn,szif,nanglg,sigbg,verbose);

	        ++ktr;
	        if((jtr-1)%mtr ==0 ){
			fprintf(jpfp," Datumed receiver %d\n",jtr);
			fflush(jpfp);
	    	}

	    }
	    ++jtr;
	} while (fgettr(infp,&tr) && jtr<=ntr);

	fprintf(jpfp," Datumed %d receivers in total\n",ktr);

	/* Seg-Y header					*/

	rewind(hdrfp);

	scal = 1/sqrt(2*PI)*dxgo;
	for(io=0; io<nxso; io++)  {
		for(ixo=0; ixo<ng[io]; ixo++) {

			efread(&tro, 1, HDRBYTES, hdrfp);
			memcpy((void *) tro.data,
			   (const void *) dats[io][ixo], nt*sizeof(float));

			for(it=0; it<nt; it++)
			   tro.data[it] *=scal;

			/* write out */
			fputtr(outfp,&tro); 

		}
	}

	fprintf(jpfp," \n");
	fprintf(jpfp," output done\n");
	fflush(jpfp);

	efclose(jpfp);
	efclose(outfp);
	efclose(hdrfp);

	free1int(ng);	    
	free2float(tsum);
	free2float(tt);
        free2float(tbg);
        free2float(pbg);
        free2float(sigbg);
        free2float(angbg);
        free3float(pb);
        free3float(tb);
        free3float(sigb);
        free3float(angb);
	free3float(ttab);
	free3float(dats);
	return(CWP_Exit());
}

/**********************************************************************
        Subroutines adapted from Dave Hale's modeling library

        decodeReflectors
        decodeReflector
        makeref                 Autor: Dave Hale, CSM, 09/17/91
**********************************************************************/

void decodeSurfaces(int sgn, int *nrPtr, int **nxzPtr, float ***xPtr,
		    float ***zPtr)
/*************************************************************************
decodeSurfaces - parse surfaces parameter string
**************************************************************************
Input :
sgn		Sign of the datuming process

Output:
nrPtr           pointer to nr an int specifying number of surfaces = 2
nxzPtr          pointer to nxz specifying number of (x,z) pairs defining the
                surfaces
xPtr            pointer to array[x][nr] of x values for each surface
zPtr            array[z][nr] of z values for each surface
**************************************************************************/
{
        int nr,*nxz,ir;
        float **x,**z;
        char t[4096],*s;

        /* count surfaces */
        nr = countparname("surf");
        if (nr==0) nr = 1;

        /* allocate space */
        nxz = ealloc1(nr,sizeof(int));
        x = ealloc1(nr,sizeof(float*));
        z = ealloc1(nr,sizeof(float*));

        /* get surfaces */
        for (ir=0; ir<nr; ++ir) {
                if (!getnparstring(ir+1,"surf",&s)) s = "0,0;99999,0";
                strcpy(t,s);
                if (!decodeSurface(t,sgn,&nxz[ir],&x[ir],&z[ir]))
                        err("Surface number %d specified "
                                "incorrectly!\n",ir+1);
        }

        /* set output parameters before returning */
        *nrPtr = nr;
        *nxzPtr = nxz;
        *xPtr = x;
        *zPtr = z;
}

/* parse one surface specification; return 1 if valid, 0 otherwise */

int decodeSurface (char *string, int sgn, int *nxzPtr, float **xPtr,
		   float **zPtr)
/**************************************************************************
decodeSurface - parse one surface specification
***************************************************************************
Input:
string          string representing surface
sgn		Sign of the datuming process

Output:
nxzPtr          pointer to number of x,z pairs
xPtr            array of x values for one surface
zPtr            array of z values for one surface
**************************************************************************/
{
        int nxz,ixz;
        float *x,*z;
        char *s,*t;

        s = string;
        s = strtok(s,",;\0");
        /* count x and z values, while splitting string into tokens */
        for (t=s,nxz=0; t!=NULL; ++nxz)
                t = strtok(NULL,",;\0");

        /* total number of values must be even */
        if (nxz%2) return 0;

        /* number of (x,z) pairs */
        nxz /= 2;

        /* 2 or more (x,z) pairs are required */
        if (nxz<2) return 0;

        /* allocate space */
        x = ealloc1(nxz,sizeof(float));
        z = ealloc1(nxz,sizeof(float));

        /* convert (x,z) values */
        for (ixz=0; ixz<nxz; ++ixz) {
                x[ixz] = atof(s);
                s += strlen(s)+1;
                z[ixz] = atof(s)*sgn;
                s += strlen(s)+1;
        }

        /* set output parameters before returning */
        *nxzPtr = nxz;
        *xPtr = x;
        *zPtr = z;
        return 1;
}

void makesurf(float dsmax,int nr,int *nu,float **xu,float **zu,
        Surface **r)
/*****************************************************************************
Make piecewise cubic surfaces
******************************************************************************
Input:
dsmax           maximum length of surface segment
nr              number of surfaces = 2
nu              array[nr] of numbers of (x,z) pairs; u = 0, 1, ..., nu[ir]
xu              array[nr][nu[ir]] of surface x coordinates x(u)
zu              array[nr][nu[ir]] of surface z coordinates z(u)

Output:
r               array[nr] of surfaces
******************************************************************************
Notes:
Space for the nu, xu, and zu arrays is freed by this function, since
they are only used to construct the surfaces.
*****************************************************************************/
{
        int ir,iu,nuu,iuu,ns,is;
        float x,z,xlast,zlast,dx,dz,duu,uu,ds,fs,rsx,rsz,rsxd,rszd,
                *u,*s,(*xud)[4],(*zud)[4],*us;
        SurfaceSegment *ss;
        Surface *rr;

        /* allocate space for surfaces */
        *r = rr = ealloc1(nr,sizeof(Surface));

        /* loop over surfaces */
        for (ir=0; ir<nr; ++ir) {

                /* compute cubic spline coefficients for uniformly sampled u */
                u = ealloc1float(nu[ir]);
                for (iu=0; iu<nu[ir]; ++iu)
                        u[iu] = iu;
                xud = (float(*)[4])ealloc1float(4*nu[ir]);
                csplin(nu[ir],u,xu[ir],xud);
                zud = (float(*)[4])ealloc1float(4*nu[ir]);
                csplin(nu[ir],u,zu[ir],zud);

                /* finely sample x(u) and z(u) and compute length s(u) */
                nuu = 20*nu[ir];
                duu = (u[nu[ir]-1]-u[0])/(nuu-1);
                s = ealloc1float(nuu);
                s[0] = 0.0;
                xlast = xu[ir][0];
                zlast = zu[ir][0];
                for (iuu=1,uu=duu; iuu<nuu; ++iuu,uu+=duu) {
                        intcub(0,nu[ir],u,xud,1,&uu,&x);
                        intcub(0,nu[ir],u,zud,1,&uu,&z);
                        dx = x-xlast;
                        dz = z-zlast;
                        s[iuu] = s[iuu-1]+sqrt(dx*dx+dz*dz);
                        xlast = x;
                        zlast = z;
                }

                /* compute u(s) from s(u) */
                ns = 1+s[nuu-1]/dsmax;
                ds = s[nuu-1]/ns;
                fs = 0.5*ds;
                us = ealloc1float(ns);
                yxtoxy(nuu,duu,0.0,s,ns,ds,fs,0.0,(float)(nu[ir]-1),us);

                /* compute surface segments uniformly sampled in s */
                ss = ealloc1(ns,sizeof(SurfaceSegment));
                for (is=0; is<ns; ++is) {
                        intcub(0,nu[ir],u,xud,1,&us[is],&rsx);
                        intcub(0,nu[ir],u,zud,1,&us[is],&rsz);
                        intcub(1,nu[ir],u,xud,1,&us[is],&rsxd);
                        intcub(1,nu[ir],u,zud,1,&us[is],&rszd);
                        ss[is].x = rsx;
                        ss[is].z = rsz;
                        ss[is].c = rsxd/sqrt(rsxd*rsxd+rszd*rszd);
                        ss[is].s = -rszd/sqrt(rsxd*rsxd+rszd*rszd);
                }

                /* fill in surface structure */
                rr[ir].ns = ns;
                rr[ir].ds = ds;
                rr[ir].ss = ss;

             /* free workspace */
                free1float(us);
                free1float(s);
                free1float(u);
                free1float((float*)xud);
                free1float((float*)zud);

                /* free space replaced by surface segments */
                free1(xu[ir]);
                free1(zu[ir]);
        }

        /* free space replaced by surface segments */
        free1(nu);
        free1(xu);
        free1(zu);
}

/* Estimation of the Z coordinate of the surfaces  */
  void zcoorSurfaces(float fx,float dx,int nx,float fxs,float dxs,int nxs,
        Surface *srf,float **szif,float *sz,float *nangl)
/*****************************************************************************
From the cubic spline calculation, the Z coor. of the surface are estimated
******************************************************************************
Input:
fx              First x position in the surfaces
dx              x sampling interval
nx              number of output location points
fxs             x-coordinate of the first shot (travel-time tables)
dxs             x-coordinate increment of shots (travel-time tables)
nxs             number of shots (travel-time tables)
srf             surface structure

Output:
szif            z[] coordinates of the surfaces (output locations)
sz              z[] coordinates of the current surface (shot positions)
nangl           array of angles that the normal form with the vertical
*****************************************************************************/
{
        int   nss,ik,jx,nsi,nsf,ns,ixi,is,ix;
        float x,ax,ax0,az,temp;
        float **ssx,**ssz,*sss;
        SurfaceSegment *ss;

        nsi = srf[0].ns;
        nsf = srf[1].ns;
        ns = MAX(nsi,nsf);
        ssx = alloc2float(2,ns);
        ssz = alloc2float(2,ns);
        sss = alloc1float(nsi);

        /* loop over surface  */
        for(is=0; is<2; ++is){
            /* number of segments, segment length */
            nss = srf[is].ns;
            ss = srf[is].ss;
            for (jx=0; jx<nss ; jx++) {
                ssx[jx][is] = ss[jx].x;
                ssz[jx][is] = ss[jx].z;
                if(is==0) sss[jx] = ss[jx].s;
            }
            ixi = 0;
            for (ix=0; ix<nx; ++ix) {
                x = fx + ix*dx;
                for (ik=ixi; ik<nss-1; ++ik) {
                    if (ssx[ik][is]<=x && ssx[ik+1][is]>=x) {
                       az = ssz[ik][is] - ssz[ik+1][is];
                       ax0 = ssx[ik+1][is] - x;
                       ax = ssx[ik+1][is] - ssx[ik][is];
                       szif[ix][is] = ax0*az/ax+ssz[ik+1][is];
                       if (is==0){
                          az = sss[ik] - sss[ik+1];
                          temp = ax0*az/ax + sss[ik+1];
                          nangl[ix] = asin(temp);
                       }
                       ixi = ik;
                    }
                }
            }
            if(is==0){
                ixi = 0;
                for (ix=0; ix<nxs; ++ix) {
                    x = fxs + ix*dxs;
                    for (ik=ixi; ik<nss-1; ++ik) {
                        if (ssx[ik][is]<=x && ssx[ik+1][is]>=x) {
                            az = ssz[ik][is] - ssz[ik+1][is];
                            ax0 = ssx[ik+1][is] - x;
                            ax = ssx[ik+1][is] - ssx[ik][is];
                            sz[ix] = ax0*az/ax+ssz[ik+1][is];
                            ixi = ik;
                        }
                    }
                }
            }
        }
        for(is=0;is<2;is++){
            szif[0][is] = 2*szif[1][is] - szif[2][is];
            szif[nx-1][is] = 2*szif[nx-2][is] - szif[nx-3][is];
        }
        sz[0] = 2*sz[1] - sz[2];
        nangl[0] = 2*nangl[1] - nangl[2];
        sz[nxs-1] = 2*sz[nxs-2] - sz[nxs-3];
        nangl[nx-1] = 2*nangl[nx-2] - nangl[nx-3];