ktmig.c~

seismic software,very useful

		dtau = dtau * 0.001;
	}
   	if (!getparfloat("tau0",&tau0)) {
		tau0 = 0.;
	} else {
		tau0 = tau0 * 0.001;
	}
		
   	if (!getparfloat("oofo",&ofomin)) ofomin = 0.;
   	if (!getparfloat("dofo",&dofo)) dofo = 999999.;
   	if (!getparfloat("obtol",&obtol)) obtol = 0.5*dofo;
   	if (!getparint("nofo",&nofo)) nofo = 1;
	itmp = countparname("ofo");
	if(itmp>0 && itmp!=nofo) err("number of ofo not match with nofo");
	ofo = (float*) emalloc(nofo*sizeof(float));
	if(itmp>0) {
		getparfloat("ofo",ofo);
		if(nofo>1) {
			ofol = ofo[0] - 0.5*(ofo[1]-ofo[0]);
			ofor = ofo[nofo-1] + 0.5*(ofo[nofo-1]-ofo[nofo-2]);
		} else {
			ofol = ofo[0] - 0.5*dofo;
			ofor = ofo[nofo-1] + 0.5*dofo;
			ofomin = ofo[0];
		}
	} else {
		for(io=0;io<nofo;io++) ofo[io] = ofomin + io*dofo;
	}
	iofoin = itmp;
	
    	if (!getparfloat("f0",&f0)) f0 = 10.;
    	if (!getparfloat("df",&df)) df = 10.;
    	if (!getparint("nf",&nf)) nf = 7;
    	if (!getparfloat("ftaper",&ftaper)) ftaper = 5.;

	/* compute maximum wavenumbers to be used in migration */
    	if (!getparfloat("ksmax",&ksmax)) {
		if(ncdppl>1) {
			ksmax = 0.5/dds;
		} else { 
			ksmax = 0.5/0.001;
		}
	}
    	if (!getparfloat("klmax",&klmax)) {
		if(nlines>1) {
			klmax = 0.5/ddl;
		} else {
			klmax = 0.5/0.001;
		}
	}


   	if (!getparint("tapecntl",&tapecntl)) tapecntl = 0;
	
	nsave = ((ntl * 3 / 2)/2)*2;
	radix_(&nsave,&nfft);
	nsave = nfft * 2 + 30;

	/* update id headers and write to output */
	bh.hns = ntau;
	bh.hdt = dtau * 1000000; 
	if(traceout==1) fputhdr(outfp,&ch,&bh);


   	if (!getparint("mlimit",&tmp_mlimit)) tmp_mlimit = 256;
   	mlimit = tmp_mlimit;
        mlimit *= 1024 * 1024;

	getparstring("diskhdr",&diskhdr);
	getparstring("diskimg",&diskimg);
	getparstring("diskfld",&diskfld);
	if(getparstring("hisfile",&hisfile)) ihis=1;

	if(!getparstring("jpfile",&jpfile)) {
                jpfp = stderr;
        } else {
                jpfp = efopen(jpfile,"w");
        }


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

	ibacki = 0;
	if (getparstring("backupi",&backupi)) {
		ibacki = 1;
		fprintf(jpfp," Backup from the previous run ");
                fprintf(jpfp," Backup %s , %s and %s from %s \n",
                        diskimg,diskhdr,diskfld,backupi);
		tar3fr(backupi,diskimg,diskhdr,diskfld);
	}

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

	fprintf(jpfp,"\n");
	fprintf(jpfp," -------- KTMIG PRINTOUT -------- \n");
	fprintf(jpfp,"\n");

	/* detect starting trace of migration */
	jtrace = 0;
	if( (hffp = fopen(hisfile,"r"))!=NULL ) {
		fclose(hffp);
                hffp = efopen(hisfile,"r+");
		buf = (char *) malloc(81*sizeof(char));
		do {
                                bzero(buf,81);
                                if(fgets(buf,81,hffp)==NULL) break;
                                if(strncmp(buf,
                                "Number of Traces Processed: ",28)==0) {
                                        sscanf(buf+28,"%d %d",&ntrpre,&jtrace);
                                }
                } while(feof(hffp)==0);
		free(buf);
		fprintf(jpfp," From history file: \n");
		fprintf(jpfp,"Number of Traces Processed: %d\n",ntrpre);
		if(jtrace>0) {
		fprintf(jpfp,"Number of Traces Processed for this input: %d\n",
			jtrace);
		}
	} else {
		hffp = efopen(hisfile,"w");
		ntrpre = 0;
	}


	/* compute memory requirement for migration*/
 
	ntmp = nl*ndl;
	nl = ntmp;
	ntmp = ns*nds;
	ns = ntmp;

	nlcore = nl;

	incdxw = 10;
	nindxw = ntau / incdxw + 1;
	ires = 4;
	ntres = nt*ires + 1;
	

	if(ntrend==0) {
		nsout = ns;
		nlout = nl;
	} else {
		nsout = ntrend;
		nlout = 1;
		nlcore = 1;
		dsout = (send-sstart)/(nsout-1);
		dlout = (lend-lstart)/(nsout-1);
		dend = sqrt ( (send-sstart)*(send-sstart)
		      + (lend-lstart)*(lend-lstart) ) / (nsout-1);
	}

	lwork = nsout * nlcore;
	lwork = lwork*ntau*nofo*sizeof(float);
	ltmp = nsout*nlout*nofo+nt + ntau + ntl + nt*mtrace;
	ltmp = ltmp + 2*ntau + 2*ntau*mtrace + 2*nsout*nlout + ntv + mtrace*4;
	ltmp = ltmp + ntau + nfft*(nf+nsout*nlout) + ntl*nf + nsave;
	ltmp = ltmp + nsout*nofo + ntres*nf + ntau;
	ltmp = ltmp + ntau*nsout*nlout*2 + ntau*mtrace + ntau*3;
	ltmp = ltmp * sizeof(float);
	ltmp = ltmp+(nsout*nlout*nofo+mtrace*2+nf*2+nindxw*nsout*nlout)*sizeof(int);
	lwork = lwork + ltmp;
	ltmp = ntv*nsv*nlv*sizeof(float);
	lwork = lwork + ltmp;

	ivdisk = 0;
	if(lwork>mlimit) {
		warn("Velocity Disk I/O Used to Reduce Memory \n");
		lwork = lwork - ltmp;
		ivdisk = 1;
	}

	if(lwork>mlimit) err("mlimit too small ");
	
	fprintf(jpfp," Total Memory used (in Byte) = %g \n", (float)lwork);
	fprintf(jpfp,"\n");

	/* allocation of memory */
	mig = (float*) emalloc(nsout*nlcore*nofo*ntau*sizeof(float));	
	fold = (float*) emalloc(nsout*nlout*nofo*sizeof(float));	
	off = (int*) emalloc(nsout*nlout*nofo*sizeof(int));	
	tgain = (float*) emalloc(nt*sizeof(float));
	tugain = (float*) emalloc(ntau*sizeof(float));
	trace = (float*) emalloc(ntl*sizeof(float));	
	tras = (float*) emalloc(nt*mtrace*sizeof(float));	
	imutel = (int*) emalloc(mtrace*sizeof(int));
	tm = (float*) emalloc(ntau*sizeof(float));
	tmig = (float*) emalloc(ntau*sizeof(float));
	fovt2 = (float*) emalloc(ntau*mtrace*sizeof(float));
	s = (float*) emalloc(nsout*nlout*sizeof(float));
	l = (float*) emalloc(nsout*nlout*sizeof(float));
	vel = (float*) emalloc(ntv*sizeof(float));
	ss = (float*) emalloc(mtrace*sizeof(float));
	sl = (float*) emalloc(mtrace*sizeof(float));
	gs = (float*) emalloc(mtrace*sizeof(float));
	gl = (float*) emalloc(mtrace*sizeof(float));
	iofs = (int*) emalloc(mtrace*sizeof(int));
	w1 = (float*) emalloc(ntau*sizeof(float));
	work1 = (float*) emalloc(nfft*nf*sizeof(float));
	tracef = (float*) emalloc(ntl*nf*sizeof(float));
	work2 = (float*) emalloc(nfft*nsout*nlout*sizeof(float));
	wsave = (float*) emalloc(nsave*sizeof(float));
	ifcut = (int*) emalloc(nf*sizeof(int));
	ltaper = (int*) emalloc(nf*sizeof(int));
	disend = (float*) emalloc(nsout*nofo*sizeof(float));
	resamp = (float*) emalloc(ntres*nf*sizeof(float));
	indxw = (int*) emalloc(nindxw*nsout*nlout*sizeof(int));

	ss2 = (float*) emalloc(ntau*sizeof(float));
	scg = (float*) emalloc(ntau*nsout*nlout*sizeof(float));
	scs = (float*) emalloc(ntau*nsout*nlout*sizeof(float));
	vr2 = (float*) emalloc(ntau*mtrace*sizeof(float));
	vi2 = (float*) emalloc(ntau*sizeof(float));
	tau = (float*) emalloc(ntau*sizeof(float));
	vq4 = (float*) emalloc(ntau*sizeof(float));

	if(ivdisk==0) velos = (float*) emalloc(ntv*nsv*nlv*sizeof(float));

	itrace = 0;

	bzero(fold,nofo*nsout*nlout*sizeof(float));
	bzero(mig,nofo*nsout*nlcore*ntau*sizeof(float));
	if(ntrend>0) {
		for(io=0;io<nsout*nofo;io++) disend[io] =  99999999999.;
	}

	if(ntrend==0) { 
		for(il=0;il<nlout;il++) {
			for(is=0;is<nsout;is++) {
				s[is+il*nsout] = sstart+is/nds*ds+is%nds*dds;
			}
		}
		for(il=0;il<nlout;il++) {
			for(is=0;is<nsout;is++) {
				l[is+il*nsout] = lstart+il/ndl*dl+il%ndl*ddl;
			}
		}
	} else {
		fprintf(jpfp,
		" Output s and l Locations of Target Line\n");
		for(is=0;is<ntrend;is++) {
			s[is] = sstart+is*(send-sstart)/(ntrend-1.);
			l[is] = lstart+is*(lend-lstart)/(ntrend-1.);
			fprintf(jpfp,
				"   s=%g l=%g trace=%d \n",s[is],l[is],is+1); 
		}
	}

	sinmin = s[0];
	sinmax = s[0];
	linmin = l[0];
	linmax = l[0];

	for(is=0;is<nsout*nlout;is++) {
		if(sinmin>s[is]) sinmin = s[is];
		if(sinmax<s[is]) sinmax = s[is];
		if(linmin>l[is]) linmin = l[is];
		if(linmax<l[is]) linmax = l[is];
	}
		
	sinmin = sinmin - apers;
	sinmax = sinmax + apers;
	linmin = linmin - aperl;
	linmax = linmax + aperl;

	/* allocate disk space for image */
	if(ntrpre>0) {

		fprintf(jpfp," open diskimg file... \n");
		imgfp = efopen(diskimg,"r+");
		fseek2g(imgfp,0,1);
		ltmp = ntau*nofo*nlout*nsout;
		efread(mig,sizeof(float),ltmp,imgfp);
		fprintf(jpfp," diskimg file opened  \n");

		fprintf(jpfp," open diskfld file... \n");
		fldfp = efopen(diskfld,"r+");
		fseek2g(fldfp,0,1);
		efread(fold,sizeof(float),nofo*nsout*nlout,fldfp);
		efclose(fldfp);
		fprintf(jpfp," diskfld file opened  \n");

		fprintf(jpfp," open diskhdr file... \n");
		hdrfp = efopen(diskhdr,"r+");
		fseek2g(hdrfp,0,1);
		fprintf(jpfp," diskhdr file opened  \n");
		
		for(ix=0;ix<nsout*nlout*nofo;ix++) {
			efread(&tr,sizeof(char),240,hdrfp);
			off[ix] = tr.offset;
			if(ntrend>0 ) disend[ix] = tr.dz;
		}
		ipre = 1;
	} else if(isave==1) {
		imgfp = efopen(diskimg,"w+r");
		fseek2g(imgfp,0,1);
		bzero(tm,ntau*sizeof(float));
		fprintf(jpfp," initialize diskimg file... \n");
		for(ix=0;ix<nsout*nlout*nofo;ix++) 
			efwrite(tm,sizeof(float),ntau,imgfp);
		fprintf(jpfp," initialization of diskimg file done \n");
		if(nlcore==nlout) efclose(imgfp);

		fprintf(jpfp," initialize diskfld file... \n");
		fldfp = efopen(diskfld,"w+r");
		fseek2g(fldfp,0,1);
		bzero(fold,nofo*nsout*nlout*sizeof(float));
		efwrite(fold,sizeof(float),nofo*nsout*nlout,fldfp);
		efclose(fldfp);
		fprintf(jpfp," initialization of diskfld file done  \n");
	}

	if(ntrpre==0) {
		fprintf(jpfp," initializing diskhdr file... \n");
		hdrfp = efopen(diskhdr,"w");
		bzero(&tr,240);
		tr.offset = 99999999;
		tr.ns = tra.ns;
		tr.dt = tra.dt;
		tr.trid = 2;
		if(ntrend>0) tr.dz = 99999999999.;
		for(iy=0;iy<nlout;iy++) {
			tr.tracr = iy + 1;
			for(ix=0;ix<nsout;ix++) {
				tmp = (l[iy*nsout+ix] - l1)/ddl + 0.5;
				iiy = tmp;
				tmp = (s[iy*nsout+ix] - s1)/dds + 0.5;
				iix = tmp;
				if(cdpnum==0) {
					tr.cdp = iiy*ncdppl + iix + cdp1;
				} else {
					tr.cdp = iix*nlines + iiy + cdp1;
				}
				tr.tracl = ix + 1;
				ms = s[iy*nsout+ix];
				ml = l[iy*nsout+ix];
				sl2xy(s1,l1,x1,y1,s2,l2,x2,y2,s3,l3,x3,y3,ms,ml,&mx,&my);
				if(ikey==1) {
				fline = (ml-l1)*(ln3-ln1)/(l3-l1) + ln1 + 0.5;
				ftrace = (ms-s1)*(tr2-tr1)/(s2-s1) + tr1 + 0.5;
					intline = fline;
					inttrace = ftrace;
					itov(trktype,&trkval,inttrace);
					itov(lnktype,&lnkval,intline);
					puthval(&tr,indxtrk,&trkval);
					puthval(&tr,indxlnk,&lnkval);
				}
				for(io=0;io<nofo;io++) {
					tr.cdpt  = io+1;
					ofs = ofo[io];
					tr.sx = mx-ofs/2;
					tr.gx = mx+ofs/2;
					tr.sy = my;
					tr.gy = my;
					efwrite(&tr,sizeof(char),240,hdrfp);
				}
			}
		}
		fprintf(jpfp," initializing diskhdr file done \n");
		for(ix=0;ix<nofo*nsout*nlout;ix++) off[ix] = 99999999;
		fflush(hdrfp);
		
	}


	ktrace = 0;
	cdppre = tra.cdp - 1;
	ispre = -9999;
	ilpre = -9999;
	m1000 = 0;