sutifowler.c

su 的源代码库

				fprintf(fpl,"	    stacking velocity %ld infinity\n",iv);
			}
		}		
	}

/******************************************************************************/
	fprintf(fpl,"sutifowler: Opening and zeroing large block matrix disk file\n");
	fprintf(fpl,"	    This can take a while, but all is fruitless if the \n");
	fprintf(fpl,"	    necessary disk space is not there...\n");
/******************************************************************************/
	N[0]=nxfft+2;
	N[1]=ntout*MAX(nv*neta,nvstack);
	fname=VNDtempname(file);
	vnd = VNDop(2,0,2,N,1,sizeof(float),fname,ndir,dir,1);
	VNDfree(fname,"main: freeing fname 1");
	fprintf(fpl,"sutifowler: large file RAM mem buf = %ld bytes\n",
		vnd->NumBytesMemBuf);
	fprintf(fpl,"sutifowler: large file disk area = %ld bytes\n",
		vnd->NumBytesPerBlock*vnd->NumBlocksPerPanel*vnd->NumPanels);


	if(getcvstacks) {
/******************************************************************************/
		fprintf(fpl,"sutifowler: reading input cvstacks\n");
/******************************************************************************/
		for(icmp=0;icmp<ncmps;icmp++) {
			key[0]=icmp;
			key[1]=0;
			for(iv=0;iv<nvstack;iv++) {
				VNDrw('w',0,vnd,1,key,0,
					(char *) tr.data,iv*ntout,1,ntout,
					1,"writing cvstacks to disk");
				if( !gettr(&tr) ) {
				    if(icmp==ncmps-1 && iv==nvstack-1 ) {
					/* all ok, read all the input data */
				    }else{
					err("sutifowler: error reading input cvstacks");
				    }
				}
			}
		}
		goto xffts;
	}
/******************************************************************************/
	fprintf(fpl,
	"sutifowler: beginning constant velocity stacks of the input cmp gathers\n");
/******************************************************************************/
	fname=VNDtempname(file);
	vnda = V2Dop(2,1000000,sizeof(float),fname,nt,mxfold);
	VNDfree(fname,"main: freeing fname 2");
	fprintf(fpl,"sutifowler: cmp gather RAM mem buf = %ld bytes\n",
		vnda->NumBytesMemBuf);

	icmp=0;
	noff=0;
	do {
	   if(tr.cdp!=oldcmp) {
		cvstack(vnda,vnd,icmp,noff,off,mute,lmute,
			nvstack,p2stack,dt,dtout);
		icmp++;
		if(icmp==ncmps) {
			fprintf(fpl,"sutifowler: more input cdps than ncdps parameter\n");
			fprintf(fpl,"	    Will only process ncdps gathers.\n");
			goto done_with_input;
			}
		oldcmp=tr.cdp;
		noff=0;
	   }
	   if(lbtaper>0 || lstaper>0) taper (lstaper,lbtaper,ncmps,icmp,nt,tr.data);
	   factor=scale;
	   for(it=0;it<nt;it++) tr.data[it]*=factor;
	   V2Dw0(vnda,noff,(char *)tr.data,1);
	   off[noff]=tr.offset;
 	   if(ichoose==1 || ichoose==2) {
 		mute[noff]=fmax*off[noff]*off[noff]*dp2;
 	   }else{
 		mute[noff]=0.;
 	   }
	   if(nonhyp) mute[noff]=MAX(mute[noff],2*off[noff]/vmin);
	   noff++;
	   if(noff>mxfold) err("tifowler: input cdp has more traces than mxfold");
	} while ( gettr(&tr) );
	cvstack(vnda,vnd,icmp,noff,off,mute,lmute,
		nvstack,p2stack,dt,dtout);
	icmp++;
done_with_input:
	ncmps=icmp;
	fprintf(fpl,"sutifowler: read and stacked %d cmp gathers\n",ncmps);
	VNDcl(vnda,1);
xffts:
	VNDflush(vnd);

	if(ichoose<5){
/******************************************************************************/
	    fprintf(fpl,"sutifowler: doing forward x -> k spatial fft's\n");
/******************************************************************************/
	    for(it=0;it<(ntout*nvstack);it++) {
		V2Dr0(vnd,it,ccrt,21);
		for(k=ncmps;k<nxfft+2;k++) rt[k]=0.;
		pfarc(1,nxfft,rt,crt);
		V2Dw0(vnd,it,ccrt,22);
	    }
	    VNDr2c(vnd);
	}

	if(ichoose<=3) {
	    fprintf(fpl,"sutifowler: looping over k\n");
	    if(TI && (ichoose==1 || ichoose==2)) { /* build ti vnmo(p) table */
		vndvnmo=ptabledmo(nv,v,etamin,deta,neta,d,vsvp,NP,dp,dp2,file);
		fprintf(fpl,"sutifowler: dmo index(p) RAM mem buf = %ld bytes\n",
			vndvnmo->NumBytesMemBuf);
	    }
	    if(TI && (ichoose==1 || ichoose==3)){ /* build ti vphase(p) table */
		vndvphase=ptablemig(nv,v,etamin,deta,neta,d,vsvp,NP,file);
		fprintf(fpl,"sutifowler: migration scaler(p) RAM mem buf = %ld bytes\n",
			vndvphase->NumBytesMemBuf);
	    }
	    if(ichoose==1 || ichoose==2){
	    	iv=MAX(nv*neta,nvstack);
		fname=VNDtempname(file);
	    	vndb = V2Dop(2,750000,sizeof(complex),
			fname,(long)ntfft,iv);
	    		fprintf(fpl,"sutifowler: (w,v) RAM mem buf = %ld bytes\n",
				vndb->NumBytesMemBuf);
		VNDfree(fname,"main: freeing fname 3");
	    }

/******************************************************************************/
	    for(k=0;k<nxfftny;k++){ 	/* loop over spatial wavenumbers */
/******************************************************************************/
		if(k==(20*(k/20))) {
			fprintf(fpl,"sutifowler: k index = %d out of %d\n",
				k,nxfftny);
		}
		ak=k*dk;
		key[0]=k;
		key[1]=0;
/******************************************************************************/
		if(ichoose==1 || ichoose==2) { /* do Fowler DMO */
/******************************************************************************/
			for(iv=0;iv<nvstack;iv++) {	/* loop over input velocities */
				VNDrw('r',0,vnd,1,key,0,ccrt,iv*ntout,1,ntout,
				31,"Fowler DMO t -> w fft read");
				for(it=ntout;it<ntfft;it++) crt[it]=czero;
				pfacc(-1,ntfft,crt);
				V2Dw0(vndb,iv,ccrt,32);
			}

			for(iw=0;iw<ntfft;iw++) {
				p=0.5*ak/fabs(w[iw]);
				if(TI) {	/* anisotropic TI*/
				    ip=p/dp;
				    if(ip<NP) {
					V2Dr0(vndvnmo,ip,(char *)rindex,40);
				    }else{
					for(iv=0;iv<(nv*neta);iv++) rindex[iv]=-1.;
				    }
				}else{			/* isotropic */
				    p2=p*p;
				    for(iv=0;iv<nv;iv++){
					v2=v[iv]*v[iv];
					rindex[iv]=(1-v2*p2)/(v2*dp2);
				    }
				}	
				V2Dr1(vndb,iw,ccrt,41);
				for(iv=0;iv<nvstack;iv++) ctemp[iv]=crt[iv];
				ints8c(nvstack,1.0,0.0,ctemp,czero,czero,nv*neta,rindex,crt);
				V2Dw1(vndb,iw,ccrt,42);
			}
			for(iv=0;iv<(nv*neta);iv++) {	/* loop over output vel */
				V2Dr0(vndb,iv,ccrt,51);
				pfacc(1,ntfft,crt);
				VNDrw('w',0,vnd,1,key,0,ccrt,iv*ntout,1,ntout,
				52,"Fowler DMO w -> t fft write");		
			}
		}
/******************************************************************************/
		if( ichoose==3 && neta>1 ) {  /* fix up disk order if only doing TI migrations */
/******************************************************************************/
			for(iv=0;iv<nv;iv++) {
				VNDrw('r',0,vnd,1,key,0,ccrt,iv*ntout,1,ntout,
				57,"option 3 fixup for multiple eta read");
				for(ieta=1;ieta<neta;ieta++) {
					VNDrw('w',0,vnd,1,key,0,ccrt,
					iv*ntout+ieta*nv*ntout,1,ntout,
					58,"option 3 fixup for multiple eta write");
				}
			}
		}
/******************************************************************************/
		if( (ichoose==1 || ichoose==3 ) ) { 	/* do Stolt migration */
/******************************************************************************/
			for(iv=0;iv<(nv*neta);iv++) {
				if(TI) {	/* anisotropic TI */
				    V2Dr0(vndvphase,iv,ccrt,50);
				    dpm=rt[0];
				    dw=2.*PI/(ntfft*dtout);
				    iwmin=0.5*ak/( (NP-3)*dpm*dw);
				    for(iw=iwmin+1;iw<ntfftny;iw++) {
					p=0.5*ak/fabs(w[iw]);
					rp=1.0+p/dpm;
					ip=rp;
					wgt=rp-ip;
					factor=wgt*rt[ip+1]+(1.-wgt)*rt[ip];
					rindex[iw]=w[iw]*factor;
					rindex[ntfft-iw]=w[ntfft-iw]*factor;
				    }
				    fw=-2.*PI/dtout;
				    rindex[0]=fw;
				    for(iw=1;iw<iwmin+1;iw++) {
					rindex[iw]=fw;
					rindex[ntfft-iw]=fw;
				    }
				}else{			/* isotropic */
					scale=0.5*v[iv]*ak;
				    	for(iw=0;iw<ntfft;iw++) {
					    if(fabs(w[iw])>scale) {
						factor=scale/w[iw];
						factor=sqrt(1+factor*factor);
						rindex[iw]=w[iw]*factor;
					    }else{
						rindex[iw]=-2.*PI/dtout;
					    }
					}
				}

				VNDrw('r',0,vnd,1,key,0,ccrt,iv*ntout,1,ntout,
					61,"Stolt t -> w fft read");
				for(it=1;it<ntout;it+=2){
					crt[it].r=-crt[it].r;
					crt[it].i=-crt[it].i;
				}
				for(it=ntout;it<ntffts;it++) crt[it]=czero;
				pfacc(1,ntffts,crt);
				dw=2.*PI/(ntffts*dtout);
				fw=-PI/dtout;
				ints8c(ntffts,dw,fw,crt,czero,czero,
					ntfft,rindex,ctemp);
				/* obliquity factor code */
 				for(iw=0;iw<ntfft;iw++){
 					factor=fabs(w[iw]/rindex[iw]);
 					crt[iw].r=factor*ctemp[iw].r;
 					crt[iw].i=factor*ctemp[iw].i;
 				}
				pfacc(-1,ntfft,crt);
				VNDrw('w',0,vnd,1,key,0,ccrt,iv*ntout,1,ntout,
					62,"Stolt w->t fft write");		
			}
		}

	    }
	    fprintf(fpl,"sutifowler: completed loop over wavenumbers\n");
	    if(ichoose==1 || ichoose==2) VNDcl(vndb,1);
	    if(TI && (ichoose==1 || ichoose==2)) VNDcl(vndvnmo,1);
	    if(TI && (ichoose==1 || ichoose==3)) VNDcl(vndvphase,1);
	}

	if(ichoose<5) {
/******************************************************************************/
	    fprintf(fpl,"sutifowler: doing inverse spatial fft's k->x\n");
/******************************************************************************/
	    for(it=0;it<(ntout*nv*neta);it++) {
		V2Dr0(vnd,it,ccrt,71);
		pfacr(-1,nxfft,crt,rt);
		V2Dw0(vnd,it,ccrt,72);
	    }
	    VNDc2r(vnd);
	}

/*****************************************************************/
	fprintf(fpl,"sutifowler: outputting results\n");
/******************************************************************/
	it=0;
	for(icmp=0;icmp<ncmps;icmp++) {
		key[0]=icmp;
		key[1]=0;
		for(ieta=0;ieta<neta;ieta++) {
			for(iv=0;iv<nv;iv++) {
				VNDrw('r',0,vnd,1,key,0,(char *)tr.data,
					iv*ntout+ieta*nv*ntout,1,ntout,82,
					"outputting all velocities for each cmp");
				tr.ns=ntout;
				tr.dt=1000000*dtout;
				tr.cdp=icmp;
				tr.tracf=iv;
				tr.offset=v[iv];
				tr.cdpt=iv;
				tr.sx=icmp*dx;
				tr.gx=icmp*dx;
				it++;
				tr.tracl=it;
				tr.tracr=it;
				tr.fldr=icmp/ngroup;
				tr.ep=10+tr.fldr*ngroup;
				tr.igc=ieta;
				tr.igi=100*(etamin+ieta*deta);
				tr.d1=dtout;
				tr.f1=0.;
				tr.d2=1.;
				tr.f2=0.;
				puttr(&tr);
			}
		}
	}

/* close files and return */
	VNDcl(vnd,1);
	VNDfree(crt,"main: freeing crt");
	VNDfree(ctemp,"main: freeing ctemp");
	VNDfree(v,"main: freeing v");
	VNDfree(p2stack,"main: freeing p2stack");
	VNDfree(mute,"main: freeing mute");
	VNDfree(off,"main: freeing off");
	VNDfree(rindex,"main: freeing rindex");
	VNDfree(w,"main: freeing w");
	if(VNDtotalmem()!=0) {
		fprintf(stderr,"total VND memory at end = %ld\n",
		VNDtotalmem());
	}
	return EXIT_SUCCESS;
}

static void cvstack(VND *vnda, VND *vnd, int icmp, int noff, float *off,
		float *mute, int lmute, int nv, float *p2,
		float dt, float dtout)
/************************************************************************
	compute constant velocity stacks for this cmp and
	load out to disk
*************************************************************************
vnda	VND file with input cmp gather
vnd	VND file for output constant velocity stacks
icmp	cmp number (zero based)
noff	number of offsets for this cmp
off[]	list of offsets for this cmp
mute[]	list of mute times for this cmp
lmute	length of mute taper
nv	number of velocities at which to compute stacks
p2[]	list of 1/(v*v) at which to compute stacks
dt	input sample rate
dtout 	output sample rate
***************************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
***************************************************************************/
{
	int iv;
	int it;
	int itmute;
	long ioff;
	long nt;
	long ntout;
	long key[2];	
	float *rt;
	float *buf;