fxymig.c

seismic software,very useful

	  	efwrite(qbuf,sizeof(float),nx*ny,xytfp);
	}
	free(qbuf);
	fclose(xytfp);

	/* trace fft and transpose data */
	/* initialize fft */
	rffti_(&nfft,wsave);
	jy = -1;
	jx = 0;
	ntrace = 0;

	tmp = (x0m - s1)/dx;
	ix0 = tmp; 
	tmp = (y0m - l1)/dy;
	iy0 = tmp; 

	do {
		if(icdps==-1) {
			gethval(&tra,indxtrk,&trkval);
			ix = vtoi(trktype,trkval);
			ix =  (ix - trstart);
			if(ix%traceinc!=0 || ix<0) continue;
			ix  = ix/traceinc;
			gethval(&tra,indxlnk,&lnkval);
			iy = vtoi(lnktype,lnkval);
			iy =  (iy - lnstart);
			if(iy%lineinc!=0 || iy<0) continue;
			iy = iy/lineinc;
		} else if(icdps==1) { 
			if(cdpnum==0) {
				iy = (tra.cdp - cdp1)/ncdppl/cdpinc;
				ix = (tra.cdp-cdp1-iy*ncdppl*cdpinc)/cdpinc 
					-ix0;
				iy = iy - iy0;
			} else {
				ix = (tra.cdp - cdp1)/nlines/cdpinc;
				iy = (tra.cdp-cdp1-ix*nlines*cdpinc)/cdpinc 
					-iy0;
				ix = ix - ix0;
			}
		} else { 
			/* x and y positions of midpoint */
			if(tra.scalco>1) {
				xm = (tra.sx + tra.gx)*0.5*tra.scalco;
				ym = (tra.sy + tra.gy)*0.5*tra.scalco;
			} else if(tra.scalco<0) {
				xm = (tra.sx + tra.gx)*0.5/(-tra.scalco);
				ym = (tra.sy + tra.gy)*0.5/(-tra.scalco);
			} else {
				xm = (tra.sx + tra.gx)*0.5;
				ym = (tra.sy + tra.gy)*0.5;
			}
			/* migration-grid line and trace positions */
			xy2sl(x1,y1,s1,l1,x2,y2,s2,l2,x3,y3,s3,l3,
				xm,ym,&sm,&lm);
			tmp = (lm - y0m)/dy + 0.5;
			iy = tmp;
			tmp = (sm - x0m)/dx + 0.5;
			ix  = tmp; 
		}


		/*
		fprintf(jpfp,"is=%d il=%d xm=%g ym=%g\n",ix+1,iy+1,xm,ym);
		*/

		/* trace outside migration grid */
		if (iy<0 || iy >= ny || ix<0 || ix>=nx || iy<jy ) continue;
		if (jy == -1) jy = iy;

		/* fft and save to a buffer */
		if (iy==jy) {
			for(it=0;it<nfft;it++) trfft[it] = 0.;
			if(tra.trid==1) bcopy((char*)tra.data,
				(char*)trfft,nt*sizeof(float));
			bcopy((char*)&tra,trhdr+ix*HDRBYTES,HDRBYTES);
			fold[ix] = fold[ix] + 1;
			for(it=0;it<nt;it++) {
				if(fabs(trfft[it])>maxamp) {
			fprintf(jpfp," bad amplitude %g (>%g) found at it=%d ix=%d iy=%d\n",
					trfft[it],maxamp,it+1,ix+1,iy+1);
			err(" bad amplitude %g (>%g) found at it=%d ix=%d iy=%d\n",
					trfft[it],maxamp,it+1,ix+1,iy+1);
				}	
			}
			rfftf_(&nfft,trfft,wsave);
			for(iw=ifmin;iw<=ifmax;iw++) {
				ir = 2 * iw - 1;
				ii = 2 * iw;
				cp[ix+(iw-ifmin)*nx].r = cp[ix+(iw-ifmin)*nx].r
							 + trfft[ir];
				cp[ix+(iw-ifmin)*nx].i = cp[ix+(iw-ifmin)*nx].i
							 + trfft[ii];
			}
			jx = jx + 1;
		/* change line */
		} else {
			fprintf(jpfp,
				" Input %d live traces for line %d\n",jx,jy+1);
			fflush(jpfp);
			ntrace += jx;
			/* write this line to disk */
			wl2d(nx,ny,nf,jy,cp,trhdr,fold,xywfp,hdrfp,jpfp,
			 	trktype,lnktype,
				trkval,lnkval,indxtrk,indxlnk,trstart,lnstart,
				traceinc, lineinc, ikey, 
				cpbig,nybig,0); 
			/* reintialize for the next line */
			bzero(cp,nx*nf*sizeof(complex));
			bzero(trhdr,nx*HDRBYTES*sizeof(char));
			bzero(fold,nx*sizeof(float));
			jx = 0;
			jy = iy;
			/* save first trace of next line */ 
			for(it=0;it<nfft;it++) trfft[it] = 0.;
			if(tra.trid==1) bcopy((char*)tra.data, 
				(char*)trfft, nt*sizeof(float));		
			bcopy((char*)&tra,trhdr+ix*HDRBYTES,HDRBYTES);
			fold[ix] = fold[ix] + 1;
			rfftf_(&nfft,trfft,wsave);
			for(iw=ifmin;iw<=ifmax;iw++) {
				ir = 2 * iw - 1;
				ii = 2 * iw;
				cp[ix+(iw-ifmin)*nx].r = cp[ix+(iw-ifmin)*nx].r
							 + trfft[ir];
				cp[ix+(iw-ifmin)*nx].i = cp[ix+(iw-ifmin)*nx].i
							 + trfft[ii];
			}
			jx = jx + 1;
		} 
	} while(fgettr(infp,&tra));


	/* for the last line */
	if (jy<ny && jx > 0 ) {
		fprintf(jpfp,
			" Input %d Live Traces For Line %d\n",jx,jy+1);
		fflush(jpfp);
		ntrace += jx;
		/* write this line to disk */
		wl2d(nx,ny,nf,jy,cp,trhdr,fold,xywfp,hdrfp,jpfp,
			trktype,lnktype,
			trkval,lnkval,indxtrk,indxlnk,trstart,lnstart,
			traceinc,lineinc,ikey, 
			cpbig,nybig,1); 
	}
	fprintf(jpfp, 
		" There Are %d Total Live Traces for Migration \n",ntrace);
	fprintf(jpfp, " Input fft Of Time Completed And Saved \n");
	fflush(jpfp);

	if(ntrace==0) err(" Check 3D setup parameters: no trace input ntrace=0 ");

	no_fft_needed:
	/* We want to close standard input now, because input is done.
	   However, the SUN f77 compiler does not allow first fortran open
	   after the standard input is closed. So we issue a dummy
	   fortran open call before we close the standard input */ 
/*
	dummyopen_(&dummyunit);
*/
	fclose(infp);

	/* free space used in trace fft */

    free(cp);
    free(cpbig);
	free(trfft); 
	free(wsave);
	free(fold);
	free(trhdr); 

	/* allocate more disk space for diskxyt if needed */
	if(ntau>ntaulast) {
		fprintf(jpfp, " Add more disk space for %s \n",diskxyt);
		fprintf(jpfp, "  ... last ntau = %d \n",ntaulast);
		fprintf(jpfp, "  ... currrent ntau = %d \n",ntau);
		xytfp = fopen(diskxyt,"r+");
		i64 = ntaulast;
		i64 = i64 * nx * ny;
		i64 = i64  * sizeof(float);
		fseek64(xytfp,i64,0);
		qbuf = (float*) emalloc(nx*ny*sizeof(float));
		bzero(qbuf,nx*ny*sizeof(float));

		for(it=ntaulast;it<ntau;it++) {
	  		fwrite(qbuf,sizeof(float),nx*ny,xytfp);
		}
		free(qbuf);
		fclose(xytfp);
	}

	/* allocate space for migration */

	lmem = 0;
    	cp = (complex*)emalloc(nx*ny*ncp*sizeof(complex));
	lmem += nx*ny*ncp*sizeof(complex);
       	bcp = (complex*)emalloc(nkx*nky*ncpu*sizeof(complex));
	lmem += nkx*nky*ncpu*sizeof(complex);
       	cpp = (complex*)emalloc(nx*ny*nwmig*sizeof(complex));
	lmem += nx*ny*nwmig*sizeof(complex);
       	caux = (complex*)emalloc(naux*ncpu*sizeof(complex));
	lmem += naux*ncpu*sizeof(complex);
       	shift = (complex*)emalloc(nf*sizeof(complex));
	lmem += nf*sizeof(complex);
    	aa = (complex*)emalloc(lplane*ncpu*sizeof(complex));
	lmem += lplane*ncpu*sizeof(complex);
    	a = (complex*)emalloc(lplane*ncpu*sizeof(complex));
	lmem += lplane*ncpu*sizeof(complex);
    	r = (complex*)emalloc(lplane*ncpu*sizeof(complex));
	lmem += lplane*ncpu*sizeof(complex);
    	cpt = (float*)emalloc(lplane*2*ncpu*sizeof(float));
	lmem += lplane*ncpu*2*sizeof(float);
    	al = (complex*)emalloc(lvec*ncpu*sizeof(complex));
	lmem += lvec*ncpu*sizeof(complex);
    	ar = (complex*)emalloc(lvec*ncpu*sizeof(complex));
	lmem += lvec*ncpu*sizeof(complex);
    	cpbuf = (complex*)emalloc(lplane*sizeof(complex));
	lmem += lplane*sizeof(complex);
    	cph = (complex*)emalloc(ncph*sizeof(complex));
	lmem += ncph*sizeof(complex);
	if(iisave>0) {
    		aasave = (complex*)emalloc(iisave*lplane*ncpu*sizeof(complex));
		lmem += iisave*lplane*ncpu*sizeof(complex);
    		asave = (complex*)emalloc(iisave*lplane*ncpu*sizeof(complex));
		lmem += iisave*lplane*ncpu*sizeof(complex);
	} else {
    		aasave = (complex*)emalloc(sizeof(complex));
		lmem += sizeof(complex);
    		asave = (complex*)emalloc(sizeof(complex));
		lmem += sizeof(complex);
	}
    	aux1 = (float*)emalloc(naux1*ncpu*sizeof(float));
	lmem += naux1*ncpu*sizeof(float);
    	aux2 = (float*)emalloc(naux2*ncpu*sizeof(float));
	lmem += naux2*ncpu*sizeof(float);
/*
    	aux1 = (float*)emalloc((naux1+naux2)*ncpu*sizeof(float));
	lmem += (naux1+naux2)*ncpu*sizeof(float);
    	aux2 = aux1+naux1*ncpu;
*/
    	om = (float*)emalloc(nf*sizeof(float));
	lmem += nf*sizeof(float);
    	oms = (float*)emalloc(lvec*ncpu*sizeof(float));
	lmem += lvec*ncpu*sizeof(float);
    	vyx = (float*)emalloc(lplane*sizeof(float));
	lmem += lplane*sizeof(float);
    	vxy = (float*)emalloc(nx*ny*sizeof(float));
	lmem += nx*ny*sizeof(float);
    	v = (float*)emalloc(nv*ny*nx*sizeof(float));
	lmem += nv*nx*ny*sizeof(float);
		q = (float*)emalloc(nq*nx*ny*sizeof(float));
	lmem += nq*nx*ny*sizeof(float);
    	va = (float*)emalloc(nv*sizeof(float));
	lmem += nv*sizeof(float);
    	w = (float*)emalloc(ncpu*sizeof(float));
	lmem += ncpu*sizeof(float);
    	qbuf = (float*)emalloc(lplane*sizeof(float));
	lmem += lplane*sizeof(float);
    	dzov = (float*)emalloc(nx*ny*sizeof(float));
	lmem += nx*ny*sizeof(float);

	/*
	err("allocation done lmem=%f \n",(float)lmem);
	*/
	fprintf(jpfp," Total Memory used (in MB) = %g \n", 
		(float)lmem/1024./1024.);
	fprintf(jpfp,"\n");

/* initialize backup disks*/
	if(lendur!=0 && itau0dur>0 && ntau==ntaulast) {
		stat(diskxyw,&sbuf);
		itime1 = (int) sbuf.st_mtime;
		stat(diskxyt,&sbuf);
		itime2 = (int) sbuf.st_mtime;
		stat(durfile,&sbuf);
		itime3 = (int) sbuf.st_mtime;
		if(itime3<itime1 || itime3<itime2) {
			if(ibackd!=1) {
			fprintf(jpfp, 
			" Updates of disks diskxyw or diskxyt failed during last run \n");
			err(" Updates of disks diskxyw or diskxyt failed during last run");
			} 
		}
	} else {
		itime1 = 1;
		itime2 = 2;
		itime3 = 3;
	}
	if(ibackd==1) {

		lenxytb = strlen(diskxytb);
		lenxywb = strlen(diskxywb);
		namexytb = (char*) emalloc(lenxytb+1);
		namexywb = (char*) emalloc(lenxywb+1);
		bcopy(diskxytb,namexytb,lenxytb);
		bcopy(diskxywb,namexywb,lenxywb);
		namexytb[lenxytb]='\0';
		namexywb[lenxywb]='\0';

		hdrbfp = fopen(diskhdrb,"w");
		fseek64(hdrfp,0,0);
		fprintf(jpfp," Copy %s to %s \n",diskhdr,diskhdrb); 
		for(ix=0;ix<nx*ny;ix++) {
	  		fread(&tra,sizeof(char),HDRBYTES,hdrfp);
	  		fwrite(&tra,sizeof(char),HDRBYTES,hdrbfp);
		}
		efclose(hdrbfp);
		
		if(itime3>=itime1 && itime3>=itime2) {

			xywbfp = fopen(diskxywb,"w");
			fseek64(xywbfp,0,0);
			bzero(cp,nx*ny*sizeof(complex));
			i64 = 0;
			fseek64(xywfp,i64,0);
			fprintf(jpfp," Copy %s to %s \n",diskxyw,diskxywb); 
			for(iw=0;iw<nf;iw++) {
	  			efread(cp,sizeof(complex),nx*ny,xywfp);
	  			efwrite(cp,sizeof(complex),nx*ny,xywbfp);
			}
			efclose(xywbfp);

			xytbfp = fopen(diskxytb,"w");
			fseek64(xytbfp,0,0);
			if(idataxyt==1) {
				fprintf(jpfp," Copy %s to %s \n",diskxyt,diskxytb); 
				if(itau0>1) {
                    xytfp = efopen(diskxyt,"r");
                    fseek64(xytfp,0,0);
                    for(it=0;it<itau0-1;it++) {
                          efread(q,sizeof(float),nx*ny,xytfp);
                          efwrite(q,sizeof(float),nx*ny,xytbfp);
                    }
                    efclose(xytfp);
                }
				bzero(q,nx*ny*sizeof(float));
				for(it=itau0-1;it<ntau;it++) {
	  				efwrite(q,sizeof(float),nx*ny,xytbfp);
				}
		  	} else {
				bzero(q,nx*ny*sizeof(float));
				for(it=0;it<ntau;it++) {
	  				efwrite(q,sizeof(float),nx*ny,xytbfp);
				}
			}
			efclose(xytbfp);
		} else {
			xywbfp = fopen(diskxywb,"r");
			fseek64(xywbfp,0,0);
			bzero(cp,nx*ny*sizeof(complex));
			i64 = 0;
			fseek64(xywfp,i64,0);

			fprintf(jpfp,
			" Updates of disks diskxyw or diskxyt failed during last run \n");
			fprintf(jpfp," Copy %s to %s \n",diskxywb,diskxyw); 
			for(iw=0;iw<nf;iw++) {
	  			efread(cp,sizeof(complex),nx*ny,xywbfp);
	  			efwrite(cp,sizeof(complex),nx*ny,xywfp);
			}
			efclose(xywbfp);

			xytbfp = fopen(diskxytb,"r");
			fseek64(xytbfp,0,0);
            xytfp = efopen(diskxyt,"r+");
            fseek64(xytfp,0,0);
			fprintf(jpfp," Copy %s to %s \n",diskxytb,diskxyt); 
            for(it=0;it<ntau;it++) {
            	efread(q,sizeof(float),nx*ny,xytbfp);
                efwrite(q,sizeof(float),nx*ny,xytfp);
            }
            efclose(xytbfp);
			iwdur = 0;
			fprintf(jpfp, " Remigrate all frequency slices: set iwdur=0 \n");
		}
	}

	/* read the ffted wavefield into cp if needed */ 
	if(ncp==nf) {
		i64 = 0;
		fseek64(xywfp,i64,0);
		bzero(cp,nx*ny*nf*sizeof(complex));
	  	fread(cp,sizeof(complex),nx*ny*ncp,xywfp);
		if(isave==0) {
			fclose(xywfp);
			eremove(diskxyw);
		}
	} else {
		fclose(xywfp);
	}