suinvco3d.c

su 的源代码库

	/* initialize varibles */
	amps = ampg =1.0;
	weight = 1.0;
	pxs = pys = pzs =0.0;
	pxg = pyg = pzg =0.0;
	d21s = d22s = d23s = d31s = d32s = d33s =0.0;
	d21g = d22g = d23g = d31g = d32g = d33g =0.0;
	ctheta = 0.0;
	det = 1.0;

	/* read in tables for special case. */
	if ( geo_type==3 ) {
	  /* traveltime table */
	  ReadTable(tfile,nxt,nyt,nzt,txs); 
	  ReadTable(tfile,nxt,nyt,nzt,txg); 

	  if ( com_type == 3 ) {

		/* amplitude table */
		ReadTable(ampfile,nxt,nyt,nzt,ampxs); 
		ReadTable(ampfile,nxt,nyt,nzt,ampxg); 

		/* d21 table */
		ReadTable(d21file,nxt,nyt,nzt,d21xs); 
		ReadTable(d21file,nxt,nyt,nzt,d21xg); 

		/* d22 table */
		ReadTable(d22file,nxt,nyt,nzt,d22xs); 
		ReadTable(d22file,nxt,nyt,nzt,d22xg); 

		/* d23 table */
		ReadTable(d23file,nxt,nyt,nzt,d23xs); 
		ReadTable(d23file,nxt,nyt,nzt,d23xg); 

		/* d31 table */
		ReadTable(d31file,nxt,nyt,nzt,d31xs); 
		ReadTable(d31file,nxt,nyt,nzt,d31xg); 

		/* d32 table */
		ReadTable(d32file,nxt,nyt,nzt,d32xs); 
		ReadTable(d32file,nxt,nyt,nzt,d32xg); 

		/* d33 table */
		ReadTable(d33file,nxt,nyt,nzt,d33xs); 
		ReadTable(d33file,nxt,nyt,nzt,d33xg); 

	   }

	   if (com_type==2 || com_type==3 ) {
		/* a1 table */
		ReadTable(a1file,nxt,nyt,nzt,a1xs); 
		ReadTable(a1file,nxt,nyt,nzt,a1xg); 

		/* b1 table */
		ReadTable(b1file,nxt,nyt,nzt,b1xs); 
		ReadTable(b1file,nxt,nyt,nzt,b1xg); 
	  }
	}

	/* loop over midpoints */
	for (ixm=0; ixm<nxm; ixm+=1){
	  xm = fxm+ixm*dxm;
   	  xs = xm-0.5*offs;
	  xg = xs+offs;
	  xs_beg = xs - dxt*(nxt-1)/2;
	  xg_beg = xg - dxt*(nxt-1)/2;

	  if (geo_type == 2 ) {
	    /* traveltime table */
	    GetFileNameType2(tfile, xs, xt0, xt1, xsfile);
	    GetFileNameType2(tfile, xg, xt0, xt1, xgfile);
	    ReadTable(xsfile,nxt,nyt,nzt,txs); 
	    ReadTable(xgfile,nxt,nyt,nzt,txg); 

	    if ( com_type == 3 ) {

		/* amplitude table */
		GetFileNameType2(ampfile, xs, xt0, xt1, xsfile);
		GetFileNameType2(ampfile, xg, xt0, xt1, xgfile); 
		ReadTable(xsfile,nxt,nyt,nzt,ampxs); 
		ReadTable(xgfile,nxt,nyt,nzt,ampxg); 

		/* d21 table */
		GetFileNameType2(d21file, xs, xt0, xt1, xsfile);
		GetFileNameType2(d21file, xg, xt0, xt1, xgfile); 
		ReadTable(xsfile,nxt,nyt,nzt,d21xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d21xg); 

		/* d22 table */
		GetFileNameType2(d22file, xs, xt0, xt1, xsfile);
		GetFileNameType2(d22file, xg, xt0, xt1, xgfile); 
		ReadTable(xsfile,nxt,nyt,nzt,d22xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d22xg); 

		/* d23 table */
		GetFileNameType2(d23file, xs, xt0, xt1, xsfile);
		GetFileNameType2(d23file, xg, xt0, xt1, xgfile); 
		ReadTable(xsfile,nxt,nyt,nzt,d23xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d23xg); 

		/* d31 table */
		GetFileNameType2(d31file, xs, xt0, xt1, xsfile);
		GetFileNameType2(d31file, xg, xt0, xt1, xgfile); 
		ReadTable(xsfile,nxt,nyt,nzt,d31xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d31xg); 

		/* d32 table */
		GetFileNameType2(d32file, xs, xt0, xt1, xsfile);
		GetFileNameType2(d32file, xg, xt0, xt1, xgfile); 
		ReadTable(xsfile,nxt,nyt,nzt,d32xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d32xg); 

		/* d33 table */
		GetFileNameType2(d33file, xs, xt0, xt1, xsfile);
		GetFileNameType2(d33file, xg, xt0, xt1, xgfile); 
		ReadTable(xsfile,nxt,nyt,nzt,d33xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d33xg); 

	    }

	    if (com_type==2 || com_type==3 ) {
		/* a1 table */
		GetFileNameType2(a1file, xs, xt0, xt1, xsfile);
		GetFileNameType2(a1file, xg, xt0, xt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,a1xs); 
		ReadTable(xgfile,nxt,nyt,nzt,a1xg); 

		/* b1 table */
		GetFileNameType2(b1file, xs, xt0, xt1, xsfile);
		GetFileNameType2(b1file, xg, xt0, xt1, xgfile); 
		ReadTable(xsfile,nxt,nyt,nzt,b1xs); 
		ReadTable(xgfile,nxt,nyt,nzt,b1xg); 

	    }
	  }

	  for (iym=0; iym<nym; ++iym){
	    ys=fym+ iym*dym;
	    ys_beg=ys - dyt*(nyt-1)/2;

	    if (geo_type == 1 ) {
	    /* traveltime table */
	    GetFileNameType1(tfile, xs, ys,xt0, xt1,yt0,yt1, xsfile);
	    GetFileNameType1(tfile, xg, ys,xt0, xt1,yt0,yt1, xgfile);
	    ReadTable(xsfile,nxt,nyt,nzt,txs); 
	    ReadTable(xgfile,nxt,nyt,nzt,txg); 

	    if ( com_type == 3 ) {

		/* amplitude table */
		GetFileNameType1(ampfile, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(ampfile, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,ampxs); 
		ReadTable(xgfile,nxt,nyt,nzt,ampxg); 

		/* d21 table */
		GetFileNameType1(d21file, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(d21file, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,d21xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d21xg); 

		/* d22 table */
		GetFileNameType1(d22file, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(d22file, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,d22xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d22xg); 

		/* d23 table */
		GetFileNameType1(d23file, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(d23file, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,d23xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d23xg); 

		/* d31 table */
		GetFileNameType1(d31file, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(d31file, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,d31xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d31xg); 

		/* d32 table */
		GetFileNameType1(d32file, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(d32file, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,d32xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d32xg); 

		/* d33 table */
		GetFileNameType1(d33file, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(d33file, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,d33xs); 
		ReadTable(xgfile,nxt,nyt,nzt,d33xg); 

	    }

	    if (com_type==2 || com_type==3 ) {
		/* a1 table */
		GetFileNameType1(a1file, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(a1file, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,a1xs); 
		ReadTable(xgfile,nxt,nyt,nzt,a1xg); 

		/* b1 table */
		GetFileNameType1(b1file, xs,ys, xt0, xt1,yt0,yt1, xsfile);
		GetFileNameType1(b1file, xg,ys, xt0, xt1,yt0,yt1, xgfile);
		ReadTable(xsfile,nxt,nyt,nzt,b1xs); 
		ReadTable(xgfile,nxt,nyt,nzt,b1xg); 

	    }
	    }

		/* loop over output points */
		for (iyo=0,yo=fyo; iyo<nyo; ++iyo, yo+=dyo){
		  yi = (yo-ys_beg)/dyt;
		  iy = yi;
		  sy = yi-iy;

		  yig = (yo-ys_beg)/dyt;
		  iyg = yig;
		  syg = yig-iyg;

		  if (iy<0 ||iy > nyt-2 ) continue;
		  if (iy==0 ||sy < 0 ) continue;
		  if (iyg<0 ||iyg > nyt-2 ) continue;
		  if (iyg==0 ||syg < 0 ) continue;

		  for (ixo=0,xo=fxo; ixo<nxo; ++ixo,xo+=dxo){

		    i=sqrt((xo-xm)*(xo-xm)+(yo-ys)*(yo-ys))/dxm;
		    if(i>nxb) continue; 

		    xi = (xo-xs_beg)/dxt;
		    ix = xi;
		    sx = xi-ix;

		    xig = (xo-xg_beg)/dxt;
		    ixg = xig;
		    sxg = xig-ixg;

		    if (ix<0 ||ix > nxt-2 ) continue;
		    if (ix==0 ||sx < 0 ) continue;
		    if (ixg<0 ||ixg > nxt-2 ) continue;
		    if (ixg==0 ||sxg < 0 ) continue;

		    for (izo=1,zo=fzo+dzo; izo<nzo; ++izo,zo+=dzo){ 
			
			/* check range of output */

			if(i>aperture[izo]) continue; 

			/* determine sample indices */
			zi = zo/dzt;
			iz = zi;
			if(iz<1||iz >nzt-2) continue;

			/* bilinear interpolation */
			sz = zi-iz;

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

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

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

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

		if (com_type ==3 ) {
		
			amps = (1.0-sy)*((1.0-sz)*((1.0-sx)*ampxs[iy][ix][iz] +
						     sx*ampxs[iy][ix+1][iz]) + 
					    sz*((1.0-sx)*ampxs[iy][ix][iz+1] 
						  + sx*ampxs[iy][ix+1][iz+1]));

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

			if (amps <TINY ) continue;

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

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

			if (ampg <TINY ) continue;
		}

		if (com_type ==3 ) {
			d21s = (1.0-sy)*((1.0-sz)*((1.0-sx)*d21xs[iy][ix][iz] +
						     sx*d21xs[iy][ix+1][iz]) + 
					    sz*((1.0-sx)*d21xs[iy][ix][iz+1] 
						  + sx*d21xs[iy][ix+1][iz+1]));

  			d21s+=sy*((1.0-sz)*((1.0-sx)*d21xs[iy+1][ix][iz] +
						sx*d21xs[iy+1][ix+1][iz]) + 
					sz*((1.0-sx)*d21xs[iy+1][ix][iz+1]
					     + sx*d21xs[iy+1][ix+1][iz+1]));
			
			d21g=(1.0-syg)*((1.0-sz)*((1.0-sxg)*d21xg[iyg][ixg][iz]
						 + sxg*d21xg[iyg][ixg+1][iz])
					+ sz*((1.0-sxg)*d21xg[iyg][ixg][iz+1] 
						+ sxg*d21xg[iyg][ixg+1][iz+1]));

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

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

  			d22s+=sy*((1.0-sz)*((1.0-sx)*d22xs[iy+1][ix][iz] +
						sx*d22xs[iy+1][ix+1][iz]) + 
					sz*((1.0-sx)*d22xs[iy+1][ix][iz+1]
					     + sx*d22xs[iy+1][ix+1][iz+1]));
			
			d22g=(1.0-syg)*((1.0-sz)*((1.0-sxg)*d22xg[iyg][ixg][iz]
						 + sxg*d22xg[iyg][ixg+1][iz])
					+ sz*((1.0-sxg)*d22xg[iyg][ixg][iz+1] 
						+ sxg*d22xg[iyg][ixg+1][iz+1]));

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

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

  			d23s+=sy*((1.0-sz)*((1.0-sx)*d23xs[iy+1][ix][iz] +
						sx*d23xs[iy+1][ix+1][iz]) + 
					sz*((1.0-sx)*d23xs[iy+1][ix][iz+1]
					     + sx*d23xs[iy+1][ix+1][iz+1]));
			
			d23g=(1.0-syg)*((1.0-sz)*((1.0-sxg)*d23xg[iyg][ixg][iz]
						 + sxg*d23xg[iyg][ixg+1][iz])
					+ sz*((1.0-sxg)*d23xg[iyg][ixg][iz+1] 
						+ sxg*d23xg[iyg][ixg+1][iz+1]));

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

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

  			d31s+=sy*((1.0-sz)*((1.0-sx)*d31xs[iy+1][ix][iz] +
						sx*d31xs[iy+1][ix+1][iz]) + 
					sz*((1.0-sx)*d31xs[iy+1][ix][iz+1]
					     + sx*d31xs[iy+1][ix+1][iz+1]));
			
			d31g=(1.0-syg)*((1.0-sz)*((1.0-sxg)*d31xg[iyg][ixg][iz]
						 + sxg*d31xg[iyg][ixg+1][iz])
					+ sz*((1.0-sxg)*d31xg[iyg][ixg][iz+1] 
						+ sxg*d31xg[iyg][ixg+1][iz+1]));

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

			d32s = (1.0-sy)*((1.0-sz)*((1.0-sx)*d32xs[iy][ix][iz] +
						     sx*d32xs[iy][ix+1][iz]) +