fxymig.c

seismic software,very useful

		if(itau0dur==-999) itau0=1;
	}

	if (!getparint("nsteps",&nsteps)) nsteps = ntau-itau0+1;
	if(ny==1) nsteps = ntau-itau0+1;
	if(nsteps>ntau-itau0+1) nsteps = ntau - itau0 + 1;

    if ( icstep > 0 ) icstep = ((icstep+iestep-1)/iestep)*iestep;
    if ( icstep ==0 ) {
		nkx = 1;
		nky = 1;
	}
	/* compute number of frequencies to migrated */
	pi = 3.141592654;
	dw = 2.*pi/(nfft*dt);
	wmin = fmin*2*pi/dw;
	iwmin = wmin;
	if (iwmin<1) iwmin = 1;
	if (iwmin>=nfftq) iwmin = nfftq-1;
	wmax = fmax*2*pi/dw;
	iwmax = wmax;
	if (iwmax<1) iwmax = 1;
	if (iwmax>=nfftq) iwmax = nfftq-1;
	nw = iwmax - iwmin + 1;
	wmin = iwmin * dw;
	tmp = fmaxend*2.*pi/dw;
	iwend= tmp;
	if(iwend>iwmax) iwend=iwmax;
	iwend = iwend - iwmin + 1;
	if(dz>0.) {
		tmp = tzend/dz + 1.5;
		ntauend = tmp;
	} else {
		tmp = tzend/(dtau*1000.) + 1.5;
		ntauend = tmp;
	}

	if(ifrange==0) {
		ifmin = iwmin;
		ifmax = iwmax;
		nf = ifmax - ifmin + 1;
	} else {
		if(ifmin<iwmin) ifmin = iwmin;
		if(ifmax>iwmax) ifmax = iwmax;
		wmin = ifmin * dw;
		nf = ifmax - ifmin + 1;
	}


	if(ny>1) {
		lvec = nx; 
		if (lvec<ny) lvec=ny;
	} else {
		if (lvec==0) lvec = nf;
	}

	fprintf(jpfp,"\n");
	fprintf(jpfp," -------- FXYMIG PRINTOUT -------- \n");
	fprintf(jpfp,"\n");
	fflush(jpfp);

	/* compute memory requirement for migration*/
	ncp = nf;
	nq = nsteps;

	if(dz>0.) {
		tmp = nq * dz/dvs + 1.5;
	} else {
		tmp = nq*dtau/(dvs*0.001) + 1.5;
	}
	nv = (int) tmp;
	if(nv<1) nv=1;
	if(ny==1) nv = nvs;
	naux1 = (nkx+15) * 4;
	naux2 = (nky+15) * 4;
	naux = nkx;
	if(nky>nkx) naux = nky;
	lplane = nx*ny;
	if(ny==1) lplane=nx*lvec;

	if(ny==1) { ncpu = 1; nwmig = 1;}
 
	if(ny>1 && iorder<=1 && dx==dy ) iisave = 1;

	lmem = nx * ny;
	lmem = lmem * ncp;
	lwork = (lmem + (4*lplane+naux+nkx*nky+2*lvec)*ncpu+nx*ny*nwmig+ncph+nf)
			*sizeof(complex);
	lmem = nx*ny;
	lmem = lmem * (2+nv+nq);
	lwork = lwork + ( (naux1+naux2+lvec)*ncpu+
		+ 2*lplane+nf+lmem+ncpu)
		*sizeof(float) 
		+ iisave*2*lplane*ncpu*sizeof(complex); 

/*
	fprintf(jpfp,"lwork=%f nx=%d ny=%d ncp=%d \n",(float)lwork,nx,ny,ncp);
	fprintf(jpfp,"mlimit=%d lwork=%lld\n",mlimit,lwork);
	fprintf(jpfp,"lplane=%d naux=%d nw=%d nkx=%d nky=%d\n",
			lplane,naux,nw,nkx,nky);
	fprintf(jpfp,"lvec=%d naux1=%d naux2=%d nv=%d nq=%d iisave=%d\n",
			lvec,naux1,naux2,nv,nq,iisave);

	lwork = (nx*ny*ncp + 6*lplane + naux + nw + nkx*nky + 2*lvec)
		*sizeof(complex);
	fprintf(jpfp,"mlimit=%d lwork=%lld\n",mlimit,lwork);
*/


/*
	if(lwork>llimit && ny>1) {
*/
	if(ny>1) {
		fprintf(jpfp,
		"Wavefield Disk I/O Used To Reduce Memory Requirement \n");
		ncp = 1;
		lwork = (nx*ny*ncp + (4*lplane+naux+nkx*nky+2*lvec)*ncpu+nx*ny*nwmig
			+ ncph + nf)*sizeof(complex);
		lmem = nx*ny;
		lmem = lmem * (2+nv+nq);
		lwork = lwork 
			+((naux1+naux2+lvec)*ncpu+2*lplane
			+nf+lmem+ncpu)
			*sizeof(float) 
			+ iisave*2*lplane*ncpu*sizeof(complex); 
	}
	if(lwork>llimit && ny==1) {
		fprintf(jpfp,
		"Velocity Disk I/O Used To Reduce Memory Requirement \n");
		nv = 1;
		lwork = (nx*ny*ncp + (4*lplane+naux+nkx*nky+nx*ny+2*lvec)
			+ ncph + nf)*sizeof(complex);
		lmem = nx*ny;
		lmem = lmem * (2+nv+nq);
		lwork = lwork 
			+((naux1+naux2+lvec)*ncpu+lplane*2
			+nf+lmem+ncpu)
			*sizeof(float) 
			+ iisave*2*lplane*ncpu*sizeof(complex); 
	}
	if(lwork>llimit && ny>1) {
		fprintf(jpfp,
	"Image/Velocity Disk I/O Used To Reduce Memory Requirement \n");
		/*
		for(iq=nq;iq>=1;iq=iq/2) {
			if(dz>0.) {
				tmp = iq * dz/dvs + 1.5;
			} else {
				tmp = iq*dtau/(dvs*0.001) + 1.5;
			}
			nv = (int) tmp;
			if(nv<1) nv=1;
			lwork = (nx*ny*ncp + 
				(4*lplane+naux+nkx*nky+2*lvec)*ncpu+nx*ny*nwmig+ncph+nf)
				*sizeof(complex);
			lwork = lwork 
				+((naux1+naux2+lvec)*ncpu+lplane*2
				+nf+nx*ny*(2+nv+iq)+ncpu)
				*sizeof(float) 
				+ iisave*2*lplane*ncpu*sizeof(complex); 
			if(lwork<llimit) break;
		}
		nq = iq;
		if(nq == 0) nq = 1;
		if(nq>nsteps) nq = nsteps;
		*/
                if(dz>0.) {
                        temp = dz/dvs;
                } else {
                        temp = dtau/(dvs*0.001);
                }
                lwork = (nx*ny*ncp+(4*lplane+naux+nkx*nky+2*lvec)*ncpu
                        +nx*ny*nwmig + ncph + nf)*sizeof(complex);
                lwork = lwork + ( (naux1+naux2+lvec)*ncpu+
                        + lplane*2+nf+nx*ny*2+ncpu)
                        *sizeof(float)
                        + iisave*2*lplane*ncpu*sizeof(complex);

/*
                fprintf(stderr,"llimit=%f lwork=%f \n",(float)llimit,(float)lwork);
*/

                tmp = (llimit-lwork)/((1.+temp)*nx*ny*sizeof(float));
                nq = tmp;
				if(nq>nsteps || nq<0) nq = nsteps;
                nv = nq * temp + 1.5;
/*
                fprintf(stderr,"nq=%d nv=%d \n",nq,nv);
*/
                lwork = lwork + (nq+nv)*nx*ny*sizeof(float);
	} else if(nq<ntau-itau0+1) {
		fprintf(jpfp,
	"Image/Velocity Disk I/O Used To Reduce Memory Requirement \n");
	}

	if(lwork>llimit && iisave>0 && nq==1) {
	      	fprintf(jpfp,
		"FD Coefficients Recomputed To Reduce Memory Requirement\n");
		iisave = 0;
		lwork = (nx*ny*ncp + (4*lplane+naux+nkx*nky+2*lvec)*ncpu 
			+ nx*ny*nwmig + ncph + nf)*sizeof(complex);
		lwork = lwork 
			+ (naux1+naux2+2*lplane+lvec*ncpu+nf
			+ nx*ny*(2+nv+nq)+ncpu)*sizeof(float) 
			+ iisave*2*lplane*ncpu*sizeof(complex); 
	}
	if(lwork>llimit && nq==1 ) {
		fprintf(jpfp,
			"Memory Limit Too Small. Need At Least %g Bytes \n",
			(float)(lwork)); 
		fprintf(jpfp, "Job Cancelled \n");

		err("Memory Limit Too Small. Need At Least %g Bytes \n",
			(float)(lwork)); 
	}

	fprintf(jpfp," Total Memory needed (in MB) = %g \n", 
		(float)lwork/1024./1024.);
	fprintf(jpfp,"\n");


	if(lendur!=0 && itau0dur>0 && itaudur>1) { 
		if(nq<(itaudur-itau0dur+1)) warn(" nq=%d too small \n",nq); 
		if(nq<(itaudur-itau0dur+1)) warn(" itaudur=%d itau0dur=%d \n",
			itaudur,itau0dur);
		nq = itaudur - itau0dur + 1;
		if(nq>nsteps || ntau>ntaulast) nq = nsteps;
		if(dz>0.) {
			tmp = nq * dz/dvs + 1.5;
		} else {
			tmp = nq*dtau/(dvs*0.001) + 1.5;
		}
		nv = (int) tmp;
		if(nv<1) nv=1;
		if(ny==1) nv = nvs;
	}


/*print setup parameters */

	fprintf(jpfp,"\n");
	fprintf(jpfp," Migration Parameters: \n");
	fprintf(jpfp," ===================== \n");
	if(icdps>=0) {
	fprintf(jpfp," cdp1=%d ncdppl=%d nlines=%d \n", cdp1,ncdppl,nlines);
	fprintf(jpfp," x1=%f y1=%f s1=%f l1=%f \n", x1,y1,s1,l1);
	fprintf(jpfp," x2=%f y2=%f s2=%f l2=%f \n", x2,y2,s2,l2);
	fprintf(jpfp," x3=%f y3=%f s3=%f l3=%f \n", x3,y3,s3,l3);
	fprintf(jpfp," sstart=%f lstart=%f \n", x0m,y0m);
	} else {
	fprintf(jpfp," tracekey=%s linekey=%s \n", tracekey,linekey);
	fprintf(jpfp," traceinc=%d lineinc=%d \n", traceinc,lineinc);
	fprintf(jpfp," trstart=%d lnstart=%d \n", trstart,lnstart);
	}
	fprintf(jpfp," ntau=%d nt=%d nw=%d ns=%d nl=%d \n", ntau,nt,nw,nx,ny);
	fprintf(jpfp," dtau=%f dt=%f dz=%f ds=%f dl=%f \n",
		dtau*1000.,dt*1000.,dz,dx,dy);
	fprintf(jpfp," itau0=%d iestep=%d icstep=%d nvs=%d naux=%d \n",
		itau0,iestep,icstep,nvs,naux);
	fprintf(jpfp," nfft=%d nffs=%d nffl=%d nsteps=%d nv=%d ncp=%d \n",
		nfft,nkx,nky,nq,nv,ncp);
	fprintf(jpfp," lvec=%d lplane=%d iorder=%d naux1=%d naux2=%d \n",
		lvec,lplane,iorder,naux1,naux2);
	fprintf(jpfp," dfc=%f vref=%f dvs=%f isave=%d iisave=%d \n", 
		dfc,vref,dvs,isave,iisave);
	fprintf(jpfp, " velfile=%s \n", velfile);
	fprintf(jpfp, " diskxyw=%s \n", diskxyw);
	fprintf(jpfp, " diskhdr=%s \n", diskhdr);
	fprintf(jpfp, " diskxyt=%s \n", diskxyt);
	if(ibackd==1) {
		fprintf(jpfp, " diskxywb=%s \n", diskxywb);
		fprintf(jpfp, " diskhdrb=%s \n", diskhdrb);
		fprintf(jpfp, " diskxytb=%s \n", diskxytb);
	}
	if(lendur>0) {
		fprintf(jpfp, " durfile=%s \n", durfile);
		fprintf(jpfp, " itau0dur=%d itaudur=%d iwdur=%d\n", 
			itau0dur,itaudur,iwdur);
	}

	fprintf(jpfp," fmin=%g fmax=%g mlimit=%d ncpu=%d nwmig=%d \n", 
		fmin, fmax, (int)(llimit/1024./1024.),ncpu,nwmig);
	fprintf(jpfp," fmaxend=%g tzend=%g \n", fmaxend, tzend);
	fprintf(jpfp," traceout=%d \n", traceout);
	if(ifrange==1) {
		fprintf(jpfp," ifmin=%d ifmax=%d nf=%d iwmin=%d iwmax=%d iwend=%d \n", 
			ifmin,ifmax,nf,iwmin,iwmax,iwend+iwmin-1);
	}

	if(ibacki==1) fprintf(jpfp, " backupi=%s \n", backupi);
	if(ibacko==1) fprintf(jpfp, " backupo=%s \n", backupo);
	fprintf(jpfp," ===================== \n");
	fprintf(jpfp,"\n");
	fflush(jpfp);

/* compute relative lowest frequency index */
	ifminr = ifmin - iwmin + 1;


	/* allocate memory for trace fft */
	lwork = nx*nf*sizeof(complex)+(nfft*2+15+nx)*sizeof(float)
			+nx*HDRBYTES*sizeof(char);
	if(lwork>llimit) 
		err("Memory Limit Too Small. Need At Least %g Mega Bytes \n",
			(lwork+1024*1024-1)/(1024*1024)); 

	nybig = llimit/2/(nx*nf*sizeof(complex));
	if(nybig<1) nybig = 1;

   	cp = (complex*)emalloc(nx*nf*sizeof(complex));
   	cpbig = (complex*)emalloc(nybig*nx*nf*sizeof(complex));
	trfft  = (float*)emalloc(nfft*sizeof(float)); 
	wsave  = (float*)emalloc((2*nfft+15)*sizeof(float));
	fold = (float*)emalloc(nx*sizeof(float));
	trhdr = (char*)emalloc(nx*HDRBYTES*sizeof(char));


	/* alloc disk spaces for wavefield data p(x,y,w) and trace headers */
	ixywsize = 0;
	if(idataxyw==1 ) {
		if((xywfp = fopen(diskxyw,"r"))!=NULL) {
			i64 = 0;
       		fseek64(xywfp,i64,SEEK_END);
       		ixywsize= ftell64(xywfp);
       		fclose(xywfp);
		}
   	}

	i64 = nx*ny;
	i64 = i64*nf;
	i64 = i64*sizeof(complex);
   	if( ixywsize!= i64 ) {
		xywfp = fopen(diskxyw,"w+");
	} else {
		xywfp = fopen(diskxyw,"r+");
	} 
    fseek64(xywfp,0,0);

	ihdrsize = 0;
	if(idatahdr==1 ) {
		if((hdrfp = fopen(diskhdr,"r"))!=NULL) {
			i64 = 0;
        	fseek64(hdrfp,i64,SEEK_END);
           	ihdrsize= ftell64(hdrfp);
       		fclose(hdrfp);
        }
	}
    i64 = nx*ny;
	i64 = i64*HDRBYTES;

    if( ihdrsize!=i64 ) {
		hdrfp = fopen(diskhdr,"w+");
	} else {
		hdrfp = fopen(diskhdr,"r+");
	} 
    fseek64(hdrfp,0,0);


	l64 = nx*ny;
	l64 = l64*nf*sizeof(complex);
	if(ihdrsize==i64 && ixywsize==l64 && itau0dur!=-999) goto no_fft_needed;

	bzero(cp,nx*nf*sizeof(complex));
	bzero(trhdr,nx*HDRBYTES*sizeof(char));
	bzero(cpbig,nybig*nx*nf*sizeof(complex));


	/* see if we can allocate the disk space */
	fclose(xywfp);
	xywfp = fopen(diskxyw,"w+"); 
	fclose(hdrfp);
	hdrfp = fopen(diskhdr,"w+"); 

	for(iy=0;iy<ny;iy++) { 
		i64 = iy*nf*nx;
		i64 = i64*sizeof(complex);
		fseek64(xywfp,i64,0);
		fwrite(cp,sizeof(complex),nx*nf,xywfp);

		i64 = iy*nx*HDRBYTES;
		i64 = i64*sizeof(char);
		fseek64(hdrfp,i64,0);
		fwrite(trhdr,sizeof(char),nx*HDRBYTES,hdrfp);
	}
	i64 = 0;
	fseek64(xywfp,i64,0);
	fseek64(hdrfp,i64,0);

	xytfp = fopen(diskxyt,"w+");
	qbuf = (float*) emalloc(nx*ny*sizeof(float));
	bzero(qbuf,nx*ny*sizeof(float));
	fseek64(xytfp,i64,0);
	for(it=0;it<ntau;it++) {