fxymig.c.6.15.99

seismic software,very useful

		indxlnk = getindex(linekey);
		s1 = 0.; l1 = 0.; 
   		if (!getparint("ns",&nx)) err(" ns missing \n");
   		if (!getparint("nl",&ny)) err(" nl missing \n");
   		if (!getparint("traceinc",&traceinc)) traceinc=1;
   		if (!getparint("lineinc",&lineinc)) lineinc=1;
   	} else if ( getparint("cdp1",&cdp1) && getparint("ncdppl",&ncdppl) &&
   	     getparint("nlines",&nlines) && getparint("cdpinc",&cdpinc) ) {
		s1 = 0.; l1 = 0.; x1 = 0.; y1 = 0.;
		s2 = 0.; l2 = 0.; x2 = 0.; y2 = 0.;
		s3 = 0.; l3 = 0.; x3 = 0.; y3 = 0.;
   		if (!getparint("ns",&nx)) nx = ncdppl;
   		if (!getparint("nl",&ny)) ny = nlines; 
		icdps = 1;
		if (cdpinc==0) cdpinc = 1;
   		if (!getparint("cdpnum",&cdpnum)) cdpnum = 0;
	} else {
		icdps = 0;
   		if (!getparfloat("x1",&x1)) err("must specify x1");
   		if (!getparfloat("y1",&y1)) err("must specify y1");
   		if (!getparfloat("s1",&s1)) err("must specify s1");
   		if (!getparfloat("l1",&l1)) err("must specify l1");
   		if (!getparfloat("x2",&x2)) err("must specify x2");
   		if (!getparfloat("y2",&y2)) err("must specify y2");
   		if (!getparfloat("s2",&s2)) err("must specify s2");
   		if (!getparfloat("l2",&l2)) err("must specify l2");
   		if (!getparfloat("x3",&x3)) err("must specify x3");
   		if (!getparfloat("y3",&y3)) err("must specify y3");
   		if (!getparfloat("s3",&s3)) err("must specify s3");
   		if (!getparfloat("l3",&l3)) err("must specify l3");
   		if (!getparint("ns",&nx)) {
			if(getparint("ncdppl",&ncdppl)) {
				nx = ncdppl;
			} else {
				err("must specify ns");
			}
		}
   		if (!getparint("nl",&ny)) {
			if(getparint("nlines",&nlines)) {
				ny = nlines;
			} else {
				ny = 1;  
			}
		}
	}
 
	if(ny==0) ny=1;  
    if (!getparfloat("ds",&dx)) err("Must specify ds!\n");
	dy = 1.;
    if (!getparfloat("dl",&dy) && ny>1) 
	err("Must specify dl for 3D migration");
	if(dy==0.) dy = 1.0;

	/* read id headers */
    fgethdr(infp,&ch,&bh);

	/* read in first trace for nt and dt */
    if (!fgettr(infp,&tra))  err("can't get first trace");
    nt = tra.ns;
    dt = (float)tra.dt/1000000.;

	if(icdps==0 || icdps==-1 ) {
		cdp1 = tra.cdp;
		cdpinc = 1;
		ncdppl = nx;
		nlines = ny;
	}	
	/* get velociry grid header info */
	velfp = efopen(velfile,"r");
	fseek64(velfp,0,1);
	
	ierr = fgetghdr(velfp,&gh);
	if (ierr==0) fromghdr(&gh,&scale,&dtype,&n1,&n2,&n3,
        		&n4,&n5,&d1,&d2,&d3,&d4,&d5,&o1,&o2,&o3,&o4,&o5,
        		&dcdp2,&dline3,&ocdp2,&oline3,
        		&gmin,&gmax,&orient,&gtype);
	fclose(velfp);

	/* optional parameters */
    if (!getparint("nvs",&nvs)) {
		if(ierr==0) {
			nvs = n3;
		} else {
			err("Must specify nvs ");
		}
	}
    if (!getparfloat("dvs",&dvs)) {
		if(ierr==0) {
			dvs = d3;
		} else {
			err("Must specify dvs ");
		}
	}

    if (!getparint("ntau",&ntau)) ntau = nt;
    if (!getparfloat("dtau",&dtau)) {
		dtau = dt;
	} else {
		dtau = dtau*0.001;
	}
    if (!getparfloat("dz",&dz)) dz = 0.;
	vref = 2.*dz/dtau;
	
	if(icdps>=0) {
    		if (!getparfloat("sstart",&x0m)) x0m = s1;
    		if (!getparfloat("lstart",&y0m)) y0m = l1;
	} else {
    		if (!getparint("trstart",&trstart)) err(" trstart missing");
    		if (!getparint("lnstart",&lnstart)) err(" lnstart missing");
		    s1 = o1;
		    l1 = o2;
    		x0m = o1;
    		y0m = o2;
	}

	/* update id headers and write to output */
	bh.hns = ntau;
	bh.hdt = dtau * 1000000; 
	if(dz>0.) {
		if(dz<=32) {
			bh.hdt = dz * 1000.;
		} else {
			warn(" Binary file hdt and trace header dt in meters or feet \n");
			bh.hdt = dz;
		}
	}

    if (!getparint("traceout",&traceout)) traceout = 1;
	if(traceout==1) fputhdr(outfp,&ch,&bh);

	if(ierr==0) {
		/*
		if(o1!=x0m) err("o1=%g x0m=%g \n",o1,x0m);
		if(o2!=y0m) err("o2=%g y0m=%g \n",o2,y0m);
		if(o3!=0.) err("o3=%g \n",o3);
		if(d1!=dx) err("d1=%g dx=%g \n",d1,dx);
		if(d2!=dy) err("d2=%g dy=%g \n",d2,dy);
		if(d3!=dvs) err("d3=%g dvs=%g \n",d3,dvs);
		if(n1!=nx) err("n1=%d nx=%d \n",n1,nx);
		if(n2!=ny) err("n2=%d ny=%d \n",n2,ny);
		if(n3!=nvs) err("n3=%d nvs=%d \n",n3,nvs);
		*/

		if(fabs(o1-x0m)>0.1 || fabs(o2-y0m)>0.1 || 
			fabs(o3)>0.1 ||
		   	fabs(d1-dx)>0.1 || fabs(d2-dy)>0.1 || 
			fabs(d3-dvs)>0.1 ||
		   	n1!=nx || n2!=ny || n3!=nvs) {

			fprintf(jpfp,
				"o1=%g sstart=%g o2=%g lstart=%g o3=%g \n",
				o1,x0m,o2,y0m,o3);
			fprintf(jpfp,
				"d1=%g ds=%g d2=%g dl=%g d3=%g dvs=%g\n",
				d1,dx,d2,dy,d3,dvs);
			fprintf(jpfp,
				"n1=%d ns=%d n2=%d nl=%d n3=%d nvs=%d\n",
				n1,nx,n2,ny,n3,nvs);
			fprintf(jpfp," check velfile header \n");
			fflush(jpfp);

			err("check velfile header ");
		}
	}

    if (!getparint("nfft",&nfft)) nfft = nt * 3 / 2;
    ntmp = (nfft+1)/2*2;
    radix_(&ntmp,&nfft);
	nfftq = nfft/2+1;
    if (!getparint("iestep",&iestep)) {
    	if (!getparint("estep",&estep)) {
			iestep = 5;
		} else {
			iestep = (estep / (dtau*1000.) + 0.5);
		}
	}
    if (!getparint("iorder",&iorder)) iorder = 1;
    if (!getparfloat("dfc",&dfc)) dfc = 0.03;
    if (!getparfloat("fmin",&fmin)) fmin = 0.05 * .5 / dt;
    if (!getparfloat("fmax",&fmax)) fmax = 1. / 2.0 * .5 / dt;
    if (!getparfloat("fmaxend",&fmaxend)) fmaxend = fmax;
	if(fmaxend>fmax) fmaxend=fmax;
    if (!getparfloat("tzend",&tzend)) tzend = 99999.;
    if (!getparint("icstep",&icstep)) {
    	if (!getparint("cstep",&cstep)) { 
			icstep=25;
		} else {
			icstep = cstep / (dtau*1000.) + 0.5; 
		}
	}
    if (!getparint("nffs",&nffx)) nffx = nx * 3 / 2;
    nffx = (nffx+1)/2*2;
    radix_(&nffx,&nkx);
    if (!getparint("nffl",&nffy)) nffy = ny * 3 / 2;
    nffy = (nffy+1)/2*2;
    radix_(&nffy,&nky);
	if(nx==1) nkx=1;
	if(ny==1) nky=1; 
    if (!getparint("lvec",&lvec)) lvec = 0;
    if (!getparint("mlimit",&mlimit)) mlimit = 512;
	mlimit = mlimit * 1024 * 1024;

	idataxyw = 1;
	if (!getparstring("diskxyw",&diskxyw)) {
		diskxyw = (char*) emalloc(80*sizeof(char));
		sprintf(diskxyw,"DISKWAVE\0");
		idataxyw = 0;
	}
	idatahdr = 1;
	if (!getparstring("diskhdr",&diskhdr)) {
		diskhdr = (char*) emalloc(80*sizeof(char));
		idatahdr = 0;
		sprintf(diskhdr,"DISKHEADER\0");
	}
	idataxyt = 1;
	if (!getparstring("diskxyt",&diskxyt)) {
		diskxyt = (char*) emalloc(80*sizeof(char));
		sprintf(diskxyt,"DISKIMAGE\0");
		idataxyt = 0;
	}

	if ( getparstring("diskxytb",&diskxytb) &&
	     getparstring("diskxywb",&diskxywb) &&
	     getparstring("diskhdrb",&diskhdrb) ) {
		ibackd = 1;
	} else {
		ibackd = 0;
	}

	if (!getparint("isave",&isave)) isave = 0;

	ibacki = 0;
	if (getparstring("backupi",&backupi)) {
		ibacki = 1;
		fprintf(jpfp," Backup from the previous run ");
		fflush(jpfp);
		tar3fr(backupi,diskhdr,diskxyt,diskxyw);
	}


	ibacko = 0;
	if (getparstring("backupo",&backupo)) ibacko = 1;
	if(ibacko==1) isave=1;

	if (!getparint("ncpu",&ncpu)) ncpu = 2;
	if (!getparint("nwmig",&nwmig)) nwmig = ncpu;
	if(nwmig<ncpu) {
	   warn(" nwmig:%d < ncpu:%d;  nwmig reset to %d \n",nwmig, ncpu, ncpu);
	   nwmig = ncpu;
	} else {
	/* make sure nwmig is an integer multiple of ncpu */
	   nwmig = ((nwmig+ncpu-1)/ncpu)*ncpu;
	}
	if (ncpu<1 || ncpu>32) err(" check ncpu value in setup");
	
	envs = (char*) emalloc(80*sizeof(char));
    if(!getenv("PARALLEL")) {
    	sprintf(envs,"%s=%d","PARALLEL",ncpu);
        putenv(envs);
    }

    if (!getparint("ifmin",&ifmin)) ifmin = -1;
    if (!getparint("ifmax",&ifmax)) ifmax = -1;
	ifrange=1;
	if(ifmin==-1 || ifmax==-1) ifrange = 0;

	/* detect starting time of migration */
	ixytsize = 0;
	ntaulast = ntau;
	if(idataxyt==1) {
		if((xytfp = fopen(diskxyt,"r"))!=NULL) {
			i64 = 0;
           	fseek64(xytfp,i64,SEEK_END);
           	ixytsize= ftell64(xytfp);
           	fclose(xytfp);
			ntaulast = ixytsize/(nx*ny*sizeof(float));
    	}
	}
	itau0 = ixytsize/(nx*ny*sizeof(float)) + 1;

	if(lendur!=0) {
		if(itau0dur > itau0) err("check durfile parameter itau0dur"); 
		if(itau0dur > 0) itau0 = itau0dur;
		if(itau0dur==-999) itau0=1;
	}

	if (!getparint("nsteps",&nsteps)) nsteps = ntau-itau0+1;
	if(ny==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 + 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>mlimit && 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>mlimit && 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>mlimit && 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<mlimit) 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,"mlimit=%d lwork=%f \n",mlimit,(float)lwork);
*/
                tmp = (mlimit-lwork)/((1.+temp)*nx*ny*sizeof(float));