fxymig.c

seismic software,very useful

	int n1,n2,n3,n4,n5,dtype,ierr;
	float o1,o2,o3,o4,o5,d1,d2,d3,d4,d5,ocdp2,oline3,dcdp2,dline3;
	float gmin,gmax,scale;
	int orient, gtype;

	int ix0, iy0, jy0;
	int nyread, nyy, iyread, ixx, iyy; 

	long long lwork, i64, l64, lmem;
	int mlimit;
	long long llimit;
    long long ixywsize,ixytsize,ihdrsize;
	
	char *envs;

	int nybig;

	struct stat sbuf;
	int itime1, itime2, itime3;


   	/* get parameters */
   	initargs(argc,argv);
   	askdoc(1);

	/* open input, output and velocity data sets */
    	if (!getparstring("datain", &datain)) {
		infp = stdin;
	} else {
    		infp = fopen(datain,"r");
	}
    	if (!getparstring("dataout", &dataout)) {
		outfp = stdout;
	} else {
    		outfp = fopen(dataout,"w");
	}
    	if (!getparstring("velfile",&velfile)) err("must specify velfile");
	
    	if (!getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
		lenjpf = 0;
	} else {
		lenjpf = strlen(jpfile);
		jpfp = efopen(jpfile,"w");
		namejp = (char*) emalloc(lenjpf+1);
		bcopy(jpfile,namejp,lenjpf);
		namejp[lenjpf]='\0';
	}
  	if (!getparstring("durfile",&durfile)) {
		lendur = 0;
		iwdur = 0;
		itau0dur = 0;
		itaudur = 0;
	} else {
		lendur = strlen(durfile);
		if( (durfp = fopen(durfile,"r"))!=NULL ) {
			iwdur = 0;
			itau0dur = -999;
			itaudur = 0;
			do {
				fscanf(durfp,"%s %d %s %d %s %d %s %s %s %s %s\n",
				tau0dur, &itau0dur, taudur, &itaudur, wdur, &iwdur,
				dump1,dump2,dump3,dump4,dump5);
			} while (feof(durfp)==0);
		} else {
			durfp = fopen(durfile,"w");
			iwdur = 0;
			itau0dur = -999;
			itaudur = 0;
		}
		fclose(durfp);
		namedur = (char*) emalloc(lendur+1);
		bcopy(durfile,namedur,lendur);
		namedur[lendur]='\0';
	}

	/* perform a fseek64 to force large file options in fread and fwrite */

	fseek64(infp,0,1);
	fseek64(outfp,0,1);

	/* required parameters */
  	if ( getparstring("tracekey",&tracekey) 
		&& getparstring("linekey",&linekey) ) {
		icdps = -1;
		ikey = 1;
		trktype = hdtype(tracekey);
		lnktype = hdtype(linekey);
		indxtrk = getindex(tracekey);
		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 (!getparfloat("maxamp",&maxamp)) maxamp = 1.e+20;
	if (!getparfloat("minvel",&minvel)) minvel = 500.;
	if (!getparfloat("maxvel",&maxvel)) maxvel = 50000;

 
	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 = fopen(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);

	/* 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 ");
		}
	}
	fseek64(velfp,0,0);
	if ( minvel>0. && maxvel>0. ) {
		vxy = (float*) emalloc(nx*ny*sizeof(float));
		for(itau=0;itau<nvs;itau++) {
			efread(vxy,sizeof(float),nx*ny,velfp);
			for(iyy=0;iyy<ny;iyy++) { 
			for(ixx=0;ixx<nx;ixx++) {
				if(vxy[ixx+iyy*nx]<minvel || vxy[ixx+iyy*nx]>maxvel) {
			fprintf(jpfp," velocity=%g error at ix=%d iy=%d itau=%d \n",
					ixx+1,iyy+1,itau+1);
					err("check velfile values ");
				}
			}
			}
		}
		free(vxy);
	}

    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;
	llimit = mlimit;
	llimit = llimit * 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;