sumigtopo2d.c

su 的源代码库

		  nr,tsum,nzt,fzt,dzt,nxt,fxt,dxt,szif,nangls,nanglg,
		  sigbs,sigbg);

	        ktr++;
	        if((jtr-1)%mtr ==0 ){
			fprintf(jpfp," migrated trace %d\n",jtr);
			fflush(jpfp);
	    	}
	    }
	    jtr++;
	} while (fgettr(infp,&tr) && jtr<ntr);

	fprintf(jpfp," migrated %d traces in total\n",ktr);

	memset((void *) &tro, 0, 240);
	tro.ns = nzo;
	tro.d1 = dzo;
	tro.dt = 1000*(int)(1000*dt+0.5);
	tro.delrt = 0.0;
	tro.f1 = fzo;
	tro.f2 = fxo;
	tro.d2 = dxo;
	tro.trid = 200;

	scal = 4/sqrt(PI)*dxm/v0;
	for(ixo=0; ixo<nxo; ixo++) {
		for(io=0; io<noff; io++)  {
			memcpy((void *) tro.data,
				(const void *) mig[io][ixo], nzo*sizeof(float));
			tro.offset = off0+io*doff;
			tro.tracr = 1+ixo;
			tro.tracl = 1+io+noff*ixo;
			tro.cdp = fxo+ixo*dxo;
			tro.cdpt = 1+io;

			for(izo=0; izo<nzo; ++izo)
				tro.data[izo] *=scal;

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

		}
	}

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

	efclose(jpfp);
	efclose(outfp);

	    
	free2float(tsum);
	free2float(tt);
        free2float(tbs);
        free2float(tbg);
        free2float(pbs);
        free2float(pbg);
        free2float(sigbs);
        free2float(sigbg);
        free2float(angbs);
        free2float(angbg);
	free3float(pb);
	free3float(tb);
	free3float(sigb);
	free3float(angb);
	free3float(ttab);
	free3float(mig);
	return(CWP_Exit());
}

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

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

void decodeSurfaces(int *nrPtr,int **nxzPtr,float ***xPtr,float ***zPtr)
/*************************************************************************
decodeSurfaces - parse surfaces parameter string
**************************************************************************
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[6144],*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,&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 *nxzPtr, float **xPtr, float **zPtr)
/**************************************************************************
decodeSurface - parse one surface specification
***************************************************************************
Input:
string          string representing surface

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);
                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
dxs             x-coordinate increment of shots
nxs             number of shots
srf             surface structure

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

        ns = srf[0].ns;
        ssx = alloc1float(ns);
        ssz = alloc1float(ns);
        sss = alloc1float(ns);

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

        free1(ssx);
        free1(ssz);
        free1(sss);
}

/* residual traveltime calculation based  on reference   time	*/
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,
		float **tb,float **t,float x0)
{
	int ix,iz,jr;
	float xi,ar,sr,sr0;

	for(ix=0; ix<nx; ++ix){
		xi = fx+ix*dx-x0;
		ar = abs(xi)/dr;
		jr = (int)ar;
		sr = ar-jr;
		sr0 = 1.0-sr;
		if(jr>nr-2) jr = nr-2;
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] -= sr0*tb[jr][iz]+sr*tb[jr+1][iz];
	}
} 

/* lateral interpolation	*/

/* sum of two tables	*/
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t)
{
	int ix,iz;

	for(ix=0; ix<nx; ++ix) 
		for(iz=0; iz<nz; ++iz)