suinvco3d.c

su 的源代码库

					    sz*((1.0-sx)*d32xs[iy][ix][iz+1] 
						  + sx*d32xs[iy][ix+1][iz+1]));

  			d32s+=sy*((1.0-sz)*((1.0-sx)*d32xs[iy+1][ix][iz] +
						sx*d32xs[iy+1][ix+1][iz]) + 
					sz*((1.0-sx)*d32xs[iy+1][ix][iz+1]
					     + sx*d32xs[iy+1][ix+1][iz+1]));
			
			d32g=(1.0-syg)*((1.0-sz)*((1.0-sxg)*d32xg[iyg][ixg][iz]
						 + sxg*d32xg[iyg][ixg+1][iz])
					+ sz*((1.0-sxg)*d32xg[iyg][ixg][iz+1] 
						+ sxg*d32xg[iyg][ixg+1][iz+1]));

  			d32g+=syg*((1.0-sz)*((1.0-sxg)*d32xg[iyg+1][ixg][iz] +
						sxg*d32xg[iyg+1][ixg+1][iz]) + 
					sz*((1.0-sxg)*d32xg[iyg+1][ixg][iz+1]
					     + sxg*d32xg[iyg+1][ixg+1][iz+1]));

			d33s = (1.0-sy)*((1.0-sz)*((1.0-sx)*d33xs[iy][ix][iz] +
						     sx*d33xs[iy][ix+1][iz]) + 
					    sz*((1.0-sx)*d33xs[iy][ix][iz+1] 
						  + sx*d33xs[iy][ix+1][iz+1]));

  			d33s+=sy*((1.0-sz)*((1.0-sx)*d33xs[iy+1][ix][iz] +
						sx*d33xs[iy+1][ix+1][iz]) + 
					sz*((1.0-sx)*d33xs[iy+1][ix][iz+1]
					     + sx*d33xs[iy+1][ix+1][iz+1]));
			
			d33g=(1.0-syg)*((1.0-sz)*((1.0-sxg)*d33xg[iyg][ixg][iz]
						 + sxg*d33xg[iyg][ixg+1][iz])
					+ sz*((1.0-sxg)*d33xg[iyg][ixg][iz+1] 
						+ sxg*d33xg[iyg][ixg+1][iz+1]));

  			d33g+=syg*((1.0-sz)*((1.0-sxg)*d33xg[iyg+1][ixg][iz] +
						sxg*d33xg[iyg+1][ixg+1][iz]) + 
					sz*((1.0-sxg)*d33xg[iyg+1][ixg][iz+1]
					     + sxg*d33xg[iyg+1][ixg+1][iz+1]));

		}

		if (com_type==2 || com_type ==3 ) {

			a1s = (1.0-sy)*((1.0-sz)*((1.0-sx)*a1xs[iy][ix][iz] +
						     sx*a1xs[iy][ix+1][iz]) + 
					    sz*((1.0-sx)*a1xs[iy][ix][iz+1] 
						  + sx*a1xs[iy][ix+1][iz+1]));

  			a1s+=sy*((1.0-sz)*((1.0-sx)*a1xs[iy+1][ix][iz] +
						sx*a1xs[iy+1][ix+1][iz]) + 
					sz*((1.0-sx)*a1xs[iy+1][ix][iz+1]
					     + sx*a1xs[iy+1][ix+1][iz+1]));
			
			a1g=(1.0-syg)*((1.0-sz)*((1.0-sxg)*a1xg[iyg][ixg][iz]
						 + sxg*a1xg[iyg][ixg+1][iz])
					+ sz*((1.0-sxg)*a1xg[iyg][ixg][iz+1] 
						+ sxg*a1xg[iyg][ixg+1][iz+1]));

  			a1g+=syg*((1.0-sz)*((1.0-sxg)*a1xg[iyg+1][ixg][iz] +
						sxg*a1xg[iyg+1][ixg+1][iz]) + 
					sz*((1.0-sxg)*a1xg[iyg+1][ixg][iz+1]
					     + sxg*a1xg[iyg+1][ixg+1][iz+1]));

			b1s = (1.0-sy)*((1.0-sz)*((1.0-sx)*b1xs[iy][ix][iz] +
						     sx*b1xs[iy][ix+1][iz]) + 
					    sz*((1.0-sx)*b1xs[iy][ix][iz+1] 
						  + sx*b1xs[iy][ix+1][iz+1]));

  			b1s+=sy*((1.0-sz)*((1.0-sx)*b1xs[iy+1][ix][iz] +
						sx*b1xs[iy+1][ix+1][iz]) + 
					sz*((1.0-sx)*b1xs[iy+1][ix][iz+1]
					     + sx*b1xs[iy+1][ix+1][iz+1]));
			
			b1g=(1.0-syg)*((1.0-sz)*((1.0-sxg)*b1xg[iyg][ixg][iz]
						 + sxg*b1xg[iyg][ixg+1][iz])
					+ sz*((1.0-sxg)*b1xg[iyg][ixg][iz+1] 
						+ sxg*b1xg[iyg][ixg+1][iz+1]));

  			b1g+=syg*((1.0-sz)*((1.0-sxg)*b1xg[iyg+1][ixg][iz] +
						sxg*b1xg[iyg+1][ixg+1][iz]) + 
					sz*((1.0-sxg)*b1xg[iyg+1][ixg][iz+1]
					     + sxg*b1xg[iyg+1][ixg+1][iz+1]));

			pxs = sin(a1s)*cos(b1s);
			pys = sin(a1s)*sin(b1s);
			pzs = cos(a1s);

			pxg = sin(a1g)*cos(b1g);
			pyg = sin(a1g)*sin(b1g);
			pzg = cos(a1g);

			ctheta = pxs*pxg+pys*pyg+pzs*pzg;
		}
		
		if ( com_type ==3 ) {
			d11 = (pxs+pxg);
			d12 = (pys+pyg);
			d13 = (pzs+pzg);
			d21 = d21s+d21g;
			d22 = d22s+d22g;
			d23 = d23s+d23g;
			d31 = d31s+d31g;
			d32 = d32s+d32g;
			d33 = d33s+d33g;

			det = d11*d22*d33 + d21*d32*d13 + d31*d12*d23
			    - d31*d22*d13 - d21*d12*d33 - d11*d32*d23;

		}

		if (com_type==1 ) weight=1.0; 
		if (com_type==2 ) weight = ctheta;
		if (com_type==3 ) 
			weight = det/amps/ampg/(1.0+ctheta);

			samp = (tsd+tgd)*odt;
			nsamp = samp;
			if (nsamp>nt-2) continue;
			samp = samp-nsamp;
			/* check operator aliasing */
			if (alias==1) {
			  if (pxs+pxg > vout[iyo][ixo][izo]*aliasx) weight =0.0;
			  if (pys+pyg > vout[iyo][ixo][izo]*aliasy) weight =0.0;
			} 
			
			/*output summation*/
			outtrace[iyo][ixo][izo] += weight*(samp*tr.data[nsamp+1]
				+(1.0-samp)*tr.data[nsamp]);

		    }
		  }
		}
		if (geo_type ==1 ) gettr(&tr); 
	  }
	  if( geo_type !=1 ) gettr(&tr); 
 	  if(verbose) fprintf(stderr, "\tfinish line%d\n", ixm+1); 
	}

 
	/* write trace */
	temp = dxm*dym/16/PI/PI/PI;  /* beta_1 */
/*	temp = dxm*dym/8/PI/PI/PI; */ /* beta_2 */
	for (iyo=0; iyo<nyo; ++iyo)
		for(ixo=0; ixo<nxo; ++ixo){
		/* scale output traces */
 		for(izo=0; izo<nzo; ++izo) {	
			if (com_type == 3) 
			   outtrace[iyo][ixo][izo] *= vout[iyo][ixo][izo];
			outtrace[iyo][ixo][izo] *= temp;
		}
		/* make headers for output traces */
		tro.offset = offs;
		tro.tracl = tro.tracr = 1+iyo*nxo+ixo;
		tro.ns = nzo;
		tro.d1 = dzo;
		tro.ns = nzo;
		tro.f1 = fzo;
		tro.f2 = fxo;
		tro.d2 = dxo;
		tro.cdp = ixo;
		tro.trid = 200;
		/* copy migrated data to output trace */
		memcpy((void *) tro.data,(const void *) outtrace[iyo][ixo],nzo*sizeof(float));
		puttr(&tro); 
	}

	free3float(txs);
	free3float(txg);
	free3float(vout);
	if (com_type == 2 || com_type == 3) {
		free3float(a1xs);
		free3float(a1xg);
		free3float(b1xs);
		free3float(b1xg);
	}
	if (com_type == 3) {
		free3float(ampxs);
		free3float(ampxg);
		free3float(d21xs);
		free3float(d21xg);
		free3float(d22xs);
		free3float(d22xg);
		free3float(d23xs);
		free3float(d23xg);
		free3float(d31xs);
		free3float(d31xg);
		free3float(d32xs);
		free3float(d32xg);
		free3float(d33xs);
		free3float(d33xg);
	}

	return EXIT_SUCCESS;
}

void GetFileNameType1(char *infile,float xs,float ys,float xs0,float xs1,float ys0,float ys1,char *xsfile) 
{
	char 	xstmp[80], append[20];
	int	iappxs,iappys;
	float	x,y;

	strcpy(xstmp,infile);

	x = xs;
	if (x < xs0) x = xs0;
	if (x > xs1) x = xs1;
	y = ys;
	if (y < ys0) y = ys0;
	if (y > ys1) y = ys1;

/*	if (x >= 0.0 ) iappxs=(x+0.0001)*1000;
	else iappxs=(x-0.0001)*1000;
	if (y >= 0.0) iappys=(y+0.0001)*1000;
	else iappys=(y-0.0001)*1000;*/
	if (x >= 0.0 ) iappxs=x+0.0001;
	else iappxs=x-0.0001;
	if (y >= 0.0) iappys=y+0.0001;
	else iappys=y-0.0001;

	sprintf(append,"%i",iappxs);
	strcat(xstmp, append);
	sprintf(append,"_%i",iappys);
	strcat(xstmp, append);
	strcpy(xsfile, xstmp);

}

void GetFileNameType2(char *infile,float xs,float xs0,float xs1,char *xsfile)
{
	char 	xstmp[80],append[20];
	int	iappxs;
	float   x;

	strcpy(xstmp,infile);

	x = xs;
	if (x < xs0) x = xs0;
	if (x > xs1) x = xs1;

/*	if (x >= 0.0 ) iappxs=(x+0.0001)*1000;
	else iappxs=(x-0.0001)*1000; */
	if (x >= 0.0 ) iappxs=x+0.0001;
	else iappxs=x-0.0001;

	sprintf(append,"%i",iappxs);
	strcat(xstmp, append);
	strcpy(xsfile, xstmp);
}

short ReadTable(char *infile,int nxt,int nyt,int nzt,float ***outArray) 
{
	FILE	*ifp;
	int	iread, ix,iy,iz;
	float	*t3d;

	t3d = ealloc1float(nzt*nxt*nyt);

	if ( (ifp=fopen(infile,"r")) ==NULL ) { 
		fprintf(stderr,"cannot open file %s\n", infile);
		return -1;
	}

	iread=fread(t3d,sizeof(float),nzt*nxt*nyt,ifp);
	if (iread != nxt*nyt*nzt) {
		fprintf(stderr,"Only read %d values of %s\n",iread,infile);
		free1float(t3d);
		return -2;
	}
	fclose(ifp);

   	for(iy=0; iy<nyt; ++iy)
	   for(ix=0; ix<nxt; ++ix)
		for(iz=0; iz<nzt; ++iz)
			outArray[iy][ix][iz]=t3d[iz+ix*nzt+iy*nzt*nxt]; 

	free1float(t3d);
	return 0;
}

short ReadVelocity(char *vfile,int nxv,int nyv,int nzv,float fxv,float dxv,
float fyv,float dyv,float fzv,float dzv,int nxo,int nyo,int nzo,float fxo,
float dxo, float fyo,float dyo,float fzo,float dzo,float ***vout) 
{
	FILE	*ifp;
	int	iread, ixo,iyo,izo,ix,iy,iz;
	int	i1,i2,i3,i4,i5,i6,i7,i8;
	float	xi,sx,yi,sy,zi,sz;
	float   *v;

	if ( (ifp=fopen(vfile,"r")) ==NULL ) { 
		fprintf(stderr,"cannot open velocity file %s", vfile);
		return -1;
	}

	v = ealloc1float(nzv*nxv*nyv);
	iread=fread(v,sizeof(float),nzv*nxv*nyv,ifp);
	if (iread != nxv*nyv*nzv) {
		fprintf(stderr,"Only read %d values of %s\n",iread,vfile);
		free1float(v);
		return -2;
	}
	fclose(ifp);
/*	fprintf(stderr,"read %d values of %s\n",iread,vfile);*/

	if (nxv == 1 && nyv ==1) {
	   for(iyo=0; iyo<nyo; ++iyo) 
	   for(ixo=0; ixo<nxo; ++ixo) 
		for(izo=0; izo<nzo; ++izo) {
		   zi = (fzo + izo*dzo - fzv)/dzv;
		   iz = zi;
		   sz = zi - iz;
		   if (iz <0 ) { iz =0; sz = 0.0; }
		   if (nzv>1 && iz >nzv-2 ) { iz =nzv-2; sz = 1.0; }

		   vout[iyo][ixo][izo] = (1.0-sz)*v[iz] + sz*v[iz+1];
		}
	}
	else if (nxv > 1 && nyv ==1) {
	   for(iyo=0; iyo<nyo; ++iyo) 
	   for(ixo=0; ixo<nxo; ++ixo) {
		xi = (fxo + ixo*dxo - fxv)/dxv;
		ix = xi;
		sx = xi - ix;
		if (ix <0 ) { ix =0; sx = 0.0; }
		if (ix >nxv-2 ) { ix =nxv-2; sx = 1.0; }

		for(izo=0; izo<nzo; ++izo) {
		   zi = (fzo + izo*dzo - fzv)/dzv;
		   iz = zi;
		   sz = zi - iz;
		   if (iz <0 ) { iz =0; sz = 0.0; }
		   if (nzv>1 && iz >nzv-2 ) { iz =nzv-2; sz = 1.0; }

		   i1 = ix*nzv + iz;
		   i2 = (ix+1)*nzv + iz;
		   i3 = ix*nzv + iz +1;
		   i4 = (ix+1)*nzv + iz +1;
		   vout[iyo][ixo][izo] = 
			   (1.0-sz)*((1.0-sx)*v[i1] + sx*v[i2]) 
			   + sz*((1.0-sx)*v[i3] + sx*v[i4]);
		}
	   }
	}
	else if (nxv > 1 && nyv >1) {
   	for(iyo=0; iyo<nyo; ++iyo) {
	   yi = (fyo + iyo*dyo - fyv)/dyv;
	   iy = yi;
	   sy = yi - iy;
	   if (iy <0 ) { iy =0; sy = 0.0; }
	   if (iy >nyv-2 ) { iy =nyv-2; sy = 1.0; }

	   for(ixo=0; ixo<nxo; ++ixo) {
		xi = (fxo + ixo*dxo - fxv)/dxv;
		ix = xi;
		sx = xi - ix;
		if (ix <0 ) { ix =0; sx = 0.0; }
		if (ix >nxv-2 ) { ix =nxv-2; sx = 1.0; }

		for(izo=0; izo<nzo; ++izo) {
		   zi = (fzo + izo*dzo - fzv)/dzv;
		   iz = zi;
		   sz = zi - iz;
		   if (iz <0 ) { iz =0; sz = 0.0; }
		   if (nzv>1 && iz >nzv-2 ) { iz =nzv-2; sz = 1.0; }

		   i1 = iy*nxv*nzv + ix*nzv + iz;
		   i2 = iy*nxv*nzv + (ix+1)*nzv + iz;
		   i3 = iy*nxv*nzv + ix*nzv + iz +1;
		   i4 = iy*nxv*nzv + (ix+1)*nzv + iz +1;
		   i5 = (iy+1)*nxv*nzv + ix*nzv + iz;
		   i6 = (iy+1)*nxv*nzv + (ix+1)*nzv + iz;
		   i7 = (iy+1)*nxv*nzv + ix*nzv + iz +1;
		   i8 = (iy+1)*nxv*nzv + (ix+1)*nzv + iz +1;
		   vout[iyo][ixo][izo] = 
			   (1.0-sy)*((1.0-sz)*((1.0-sx)*v[i1] + sx*v[i2]) 
			   + sz*((1.0-sx)*v[i3] + sx*v[i4]))
  			   + sy*((1.0-sz)*((1.0-sx)*v[i5] + sx*v[i6]) 
			   + sz*((1.0-sx)*v[i7] + sx*v[i8]));     

		}
	   }
	}
	}

	free1float(v);
	return 0;   
}