fxymig.c.6.15.99

seismic software,very useful

                nq = tmp;
		if(nq<=0) nq = 1;
                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>mlimit && 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>mlimit && 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 used (in Byte) = %g \n", (float)lwork);
	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(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 nq=%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, mlimit/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);


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

	nybig = mlimit/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++) {
	  	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;
			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,
				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,
			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) {
		xytfp = fopen(diskxyt,"r+");
		i64 = ntaulast * 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++) {
	  		efwrite(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);
	*/

/* initialize backup disks*/
	if(lendur!=0 && itau0dur>0 ) {
		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) {