kzmig.c.9.2.99

seismic software,very useful

     	if (!getparfloat("oofo",&ofomin)) ofomin = 0.;
    	if (!getparfloat("dofo",&dofo)) dofo = 999999.;
    	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 = 6;
    	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 = nzo;
	bh.fold = nofo;
	if(dzo>32.) {
		bh.hdt = dzo;
	} else {
		bh.hdt = dzo*1000;
	}
	if(nofo>1) {
		bh.tsort = 2;
	} else {
		bh.tsort = 4;
	}
 	if(traceout==1) fputhdr(outfp,&ch,&bh);


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

	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;

	if (getparstring("backupi",&backupi)) {
		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;
	if(ntrdsk>0) isave=1;

	fprintf(jpfp,"\n");
	fprintf(jpfp," -------- KZMIG 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 the reference traveltimes and lateral slowness	*/
	dr = (dds<ddl)?dds:ddl;
	rmax = sqrt((exst-fst)*(exst-fst)+(eyst-flt)*(eyst-flt));
	if(rmax<sqrt((est-fxst)*(est-fxst)+(elt-fyst)*(elt-fyst)))
		rmax = sqrt((est-fxst)*(est-fxst)+(elt-fyst)*(elt-fyst));
	if(rmax>0.5*offmax+sqrt(apers*apers+aperl*aperl)+dxst+dyst)
		rmax = 0.5*offmax+sqrt(apers*apers+aperl*aperl)+dxst+dyst;

	nr = 1+(int)(rmax/dr+0.99);

/*
	fprintf(stderr," before timeb dds=%g ddl=%g \n",dds,ddl);
	fprintf(stderr," before timeb rmax=%d dr=%g nr=%d nzt=%d \n",
	rmax,dr,nr,nzt);
*/

 	tb0 = alloc1float(nzt*nr);
 	pb0 = alloc1float(nzt*nr);


 	timeb_(&nr,&nzt,&dr,&dzt,&fzt,&dvz,&v0,tb0,pb0);
	free1float(pb0);
 	pb = alloc1float(nz*nr);
 	tb = alloc1float(nz*nr);
 	timeb_(&nr,&nz,&dr,&dzt,&fz,&dvz,&v0,tb,pb);


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

	nlcore = nl;

	if(ntrend==0) {
		nsout = ns;
		nlout = nl;
		dsout = ds;
		dlout = dl;
	} 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);
	}

	/* check grid paramters if migration output line is slant	*/
	if(dlt>999990) {
		if(fabs(fst-sstart)<0.1) fst = sstart;
		if(fabs(flt-lstart)<0.1) flt = lstart;
		tmp = sqrt((send-sstart)*(send-sstart)
			+(lend-lstart)*(lend-lstart));
		tmp = tmp/dst+1.5;
		itmp = tmp;
		if(nlt!=1 || fst!=sstart || flt!=lstart || itmp!=nst)
err("grid parameters in traveltime table are inconsistent with output!\n");
 	}
		
	lwork = nsout * nlcore;
	lwork = lwork*nzo*nofo*sizeof(float);

	ltmp = nz*nsout*nlout*nofo+nt+ntl+nt*mtrace;
	ltmp = ltmp+2*nsout*nlout+mtrace*4;
	ltmp = ltmp  + nfft*(nf+2) + ntl*nf + nsave;
	ltmp = ltmp * sizeof(float);
	ltmp = ltmp + (nsout*nlout*nofo+mtrace+nf*2)*sizeof(int);
	ltmp += (nz*nsout*nlout*ntab+nzt*nst*nlt)*sizeof(float);

	ltmp += (nzt*nr+2*nz*nr)*sizeof(float);
	ltmp += 2*nz*nsout*nlout*sizeof(float);
	lwork = lwork + ltmp;
	tmp = (lwork+1024*1024-1)/(1024*1024); 
	if(lwork>mlimit) err("mlimit too small; need mlimit=%d \n ",(int)tmp);

	tmp = mlimit - lwork;
	tmp = tmp / (nz*nsout*nlout*sizeof(float));
	lwork = lwork - nz*nsout*nlout*ntab*sizeof(float); 
	ntab =  ntab + tmp;
	if(ntab>nxst*nyst)  ntab = nxst*nyst; 
	lwork = lwork + nz*nsout*nlout*ntab*sizeof(float); 
	
	fprintf(jpfp," Total Memory used (in Byte) = %g \n", (float)lwork);
	fprintf(jpfp,"\n");

	/* allocation of memory */
	mig = ealloc2float(nsout*nlcore*nzo,nofo);	
	fold = ealloc2float(nsout*nlcore*nz,nofo);	
	off = (int*) emalloc(nsout*nlout*nofo*sizeof(int));	
	tgain = (float*) emalloc(nt*sizeof(float));
	zgain = (float*) emalloc(nzo*sizeof(float));
 	trace = (float*) emalloc(ntl*sizeof(float));	
	tras = (float*) emalloc(nt*mtrace*sizeof(float));	
	imutel = (int*) emalloc(mtrace*sizeof(int));
	s = (float*) emalloc(nsout*nlout*sizeof(float));
	l = (float*) emalloc(nsout*nlout*sizeof(float));
	st = (float*) emalloc(nst*sizeof(float));
	lt = (float*) emalloc(nst*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));
 	work1 = (float*) emalloc(nfft*nf*sizeof(float));
 	sqrtf = (float*) emalloc(nfft*sizeof(float));
	tracef = (float*) emalloc(ntl*nf*sizeof(float));
	work2 = (float*) emalloc(nfft*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));
 	ttab = ealloc2float(nz*nsout*nlout,ntab);
	tt = (float*) emalloc(nzt*nst*nlt*sizeof(float));	

	datetime=(char*) emalloc(28*sizeof(char));

	itrace = 0;

	bzero(fold[0],nofo*nsout*nlout*nz*sizeof(float));
	bzero(mig[0],nofo*nsout*nlcore*nzo*sizeof(float));
	if(ntrend>0) {
		for(io=0;io<nsout*nofo;io++) disend[io] =  99999999999.;
		dismax = disaper(apers,aperl,sstart,lstart,send,lend);
		xy2abc(sstart,lstart,send,lend,&aa,&bb,&cc);
	}

	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 {
		for(is=0;is<ntrend;is++) {
			s[is] = sstart+is*(send-sstart)/(ntrend-1.);
			l[is] = lstart+is*(lend-lstart)/(ntrend-1.);
		}
		for(is=0;is<nst;is++) {
			st[is] = sstart+is*(send-sstart)/(nst-1.);
			lt[is] = lstart+is*(lend-lstart)/(nst-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];
	}

/*	check migration output range	*/
	stmin = fst; if (est<fst) stmin = est;
        ltmin = flt; if (elt<flt) ltmin = elt;
        stmax = est; if (est<fst) stmax = fst;
        ltmax = elt; if (elt<flt) ltmax = flt;
        if(stmin>sinmin || stmax<sinmax-0.01
        || ltmin>linmin || ltmax<linmax-0.01) {
                warn("stmin=%g > sinmin=%g ? \n",stmin,sinmin);
                warn("stmax=%g < sinmax-0.01=%g ? \n",stmax,sinmax-0.01);
                warn("ltmin=%g > linmin=%g ? \n",ltmin,linmin);
                warn("ltmax=%g < linmax-0.01=%g ? \n",ltmax,linmax-0.01);
                err(" migration output is out of time tables!\n");
        }

	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 = nzo*nofo*nlout*nsout;
		if(fread(mig[0],sizeof(float),ltmp,imgfp)!=ltmp)
			err(" read error on %s \n",diskimg);
		efclose(imgfp);
		fprintf(jpfp," diskimg file opened  \n");

		fprintf(jpfp," open diskfld file... \n");
		fldfp = efopen(diskfld,"r+");
		fseek2g(fldfp,0,1);
		ltmp = nofo*nsout*nlout*nz;
		if(fread(fold[0],sizeof(float),ltmp,fldfp)!=ltmp)
			err(" read error on %s \n",diskfld);
		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++) {
			if(fread(&tr,sizeof(char),240,hdrfp)!=240)
				err(" read error on %s \n",diskhdr);
			off[ix] = tr.offset;
			if(ntrend>0 ) disend[ix] = tr.ungpow;
		}
		ipre = 1;
	} else if(isave==1) {
		imgfp = efopen(diskimg,"w+r");
		fseek2g(imgfp,0,1);
		bzero(mig[0],nofo*nsout*nlout*nzo*sizeof(float));
		fwrite(mig[0],sizeof(float),nofo*nsout*nlout*nzo,imgfp);
 		fprintf(jpfp," initialize diskimg file... \n");
 		fprintf(jpfp," initialization of diskimg file done \n");
		efclose(imgfp);

		fprintf(jpfp," initialize diskfld file... \n");
		fldfp = efopen(diskfld,"w+r");
		fseek2g(fldfp,0,1);
		bzero(fold[0],nofo*nsout*nlout*nz*sizeof(float));
		fwrite(fold[0],sizeof(float),nofo*nsout*nlout*nz,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.ungpow = 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);