kzmig.c

seismic software,very useful

			}
		}
		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.);
		}
	}

	/* bilinear interpolation of output depth from hzgrid */
	if(ihzgrid==1) {
		for(itr=0;itr<nsout*nlout;itr++) {
			sss = s[itr];
			lll = l[itr];
			bilint_(&one,&nshz,&nlhz,&oshz,&olhz,&dshz,&dlhz,&sss,&lll,hz,&hzz);
			hzz = hzz/dzo;
			itmp = hzz + 0.5;
			hzz = itmp * dzo;
			fzout[itr] = hzz;
			/*
			fprintf(stderr,"fzout=%g itr=%d \n",fzout[itr],itr);
			*/
			if(itr==0) {hzmin = hzz; hzmax = hzz;}
			if(hzmin>hzz) hzmin = hzz;
			if(hzmax<hzz) hzmax = hzz;
		}
		if(hzmin<fzt || hzmax >ezt) {
	fprintf(stderr," hzgrid min=%g max=%g outside of tracvel time depths\n",
		hzmin,hzmax);
		}
	}

	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+1.0 || stmax<sinmax-1.0
        || ltmin>linmin+1.0 || ltmax<linmax-1.0) {
                warn(" migration output is out of time table ranges !\n");
    }

    if(stmin>sinmin+1.0) err(" migration output smmin=%g < travel time stmin=%g \n",sinmin,stmin);
    if(stmax<sinmax-1.0) err(" migration output smmax=%g > travel time stmax=%g \n",sinmax,stmax);
    if(ltmin>linmin+1.0) err(" migration output lmmin=%g < travel time ltmin=%g \n",linmin,ltmin);
    if(stmax<sinmax-1.0) err(" migration output lmmax=%g > travel time ltmax=%g \n",linmax,ltmax);
    
    /*
    if(stmin>sinmin+1.0 || stmax<sinmax-1.0
        || ltmin>linmin+1.0 || ltmax<linmax-1.0) {
                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);
					puthval(&tr,indxtrk,&trkval);
					puthval(&tr,indxlnk,&lnkval);
				}
				if(dzo<=32.) { 
					tr.dt = dzo*1000+0.5;
				} else {
					tr.dt = dzo;
				}
				tr.dt = dzo*1000.+0.5;
				tmp = fzout[ix+iy*nsout]+0.5;
				itmp = tmp;
				tr.delrt = (short)itmp;
				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;
					fwrite(&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;
	m1000 = 0;
	oldheadx = newheadx = 0;
	oldheady = newheady = 0;
	ssmin = slmin = 99999999.0;
	ssmax = slmax = -99999999.0;
 
	if(itgain==1) {
		for(it=0;it<nt;it++) {
			tmp = t0 + it*dt;
			if(tmp!=0.0) {
				tgain[it] = pow(tmp,tpow);
			} else {
				tgain[it] = 0.;
			}
		}
 	} else {
		for(it=0;it<nt;it++) tgain[it] = 1.0;
	}
		
	fprintf(jpfp," \n");
	fprintf(jpfp," Migration Parameters \n");
	fprintf(jpfp," -------------------- \n");
	fprintf(jpfp," x1=%g y1=%g s1=%g l1=%g cdp1=%d\n",x1,y1,s1,l1,cdp1);
	fprintf(jpfp," x2=%g y2=%g s2=%g l2=%g cdp2=%d\n",x2,y2,s2,l2,cdp2);
	fprintf(jpfp," x3=%g y3=%g s3=%g l3=%g cdp3=%d\n",x3,y3,s3,l3,cdp3);
	fprintf(jpfp," sstart=%g lstart=%g ns=%d nl=%d ds=%g dl=%g \n",
		sstart,lstart,nsread,nlread,ds,dl);
	if(ntrend>0) fprintf(jpfp," send=%g lend=%g ntrend=%d dend=%g \n",
			send,lend,ntrend,dend);
	if(ntrend>0) fprintf(jpfp," dismax=%g aa=%g bb=%g cc=%g \n",
			dismax,aa,bb,cc);
	fprintf(jpfp," fzo=%g nzo=%d dzo=%g\n",fzo,nzo,dzo);
	if(ihzgrid==1) fprintf(jpfp," hzgrid=%s\n",hzgrid);
	fprintf(jpfp," nds=%d ndl=%d dds=%g ddl=%g \n",
		nds,ndl,dds,ddl);
	fprintf(jpfp," nss=%d nll=%d is0=%d il0=%d \n",
		nss,nll,is0,il0);
	fprintf(jpfp," cdppdds=%d cdppds=%d cdppddl=%d cdppdl=%d cdpnum=%d\n",
		cdppdds, cdppds, cdppddl, cdppdl, cdpnum);
	fprintf(jpfp," ntrpre=%d nt=%d ntl=%d dt=%g t0=%g nfft=%d\n",
		ntrpre,nt,ntl,dt,t0,nfft);
fprintf(jpfp," tpow=%g zpow=%g ofomin=%g dofo=%g obtol=%f nofo=%d isave=%d \n",
		tpow,zpow,ofomin,dofo,obtol,nofo,isave);
	fprintf(jpfp," apers=%g aperl=%g angmax=%g tapecntl=%d\n",
		apers,aperl,angmax,tapecntl);
	fprintf(jpfp," nlcore=%d mtrace=%d ntras=%d mlimit=%d ntrdsk=%d \n",
		nlcore,mtrace,ntras,(int)(llimit/(1024*1024)),ntrdsk);
	fprintf(jpfp," f0=%g df=%g nf=%d ftaper=%g ksmax=%g klmax=%g\n",
		f0,df,nf,ftaper,ksmax,klmax);
	fprintf(jpfp," offmax=%g ntab=%d nr=%d \n",offmax,ntab,nr);
 	fprintf(jpfp," fzt=%g dzt=%g nzt=%d \n",fzt,dzt,nzt);
	fprintf(jpfp," fst=%g dst=%g nst=%d \n",fst,dst,nst);
	fprintf(jpfp," flt=%g dlt=%g nlt=%d \n",flt,dlt,nlt);
	fprintf(jpfp," fxst=%g dxst=%g nxst=%d \n",fxst,dxst,nxst);
	fprintf(jpfp," fyst=%g dyst=%g nyst=%d \n",fyst,dyst,nyst);
	fprintf(jpfp," v0=%g dvz=%g gath=%d \n",v0,dvz,gath);
	fprintf(jpfp," diskimg=%s\n",diskimg);
	fprintf(jpfp," diskfld=%s\n",diskfld);
	fprintf(jpfp," diskhdr=%s\n",diskhdr);
	fprintf(jpfp," ttfile=%s\n",ttfile);
	fprintf(jpfp," backupi=%s\n",backupi);
	fprintf(jpfp," backupo=%s\n",backupo);
	fprintf(jpfp," hisfile=%s\n",hisfile);
	fprintf(jpfp," ncpu=%d traceout=%d\n",ncpu,traceout);
	if(on2trace==1) fprintf(jpfp," total_number_of_input_trace=%d \n",ntotal);
	if(ikey==1) {
		fprintf(jpfp," tr1=%g tr2=%g tr3=%g \n",tr1,tr2,tr3);
		fprintf(jpfp," ln1=%g ln2=%g ln3=%g \n",ln1,ln2,ln3);
		fprintf(jpfp," tracekey=%s linekey=%s \n",tracekey,linekey);
		fprintf(jpfp," strace=%d jthtrace=%d sline=%d jthline=%d \n",
			strace,jthtrace,sline,jthline);
		fprintf(jpfp," mintrace=%d maxtrace=%d \n",mintrace,maxtrace);
		fprintf(jpfp," minline=%d maxline=%d  \n",minline,maxline);
	}
	fprintf(jpfp," \n");


/*
	gettime(datetime);
 	fprintf(jpfp," start migration at %s ... \n",datetime);
*/
 	fprintf(jpfp," start migration at ... \n");
 	fflush(jpfp); 

 	m1000 = ntrpre/1000;
	ntrace = 0;

	if(tapecntl==1) {
		if(jtrace>0) {
        	lpos = nt*sizeof(float)+240;
            lpos = lpos * jtrace + 3600;
            fseek2g(infp,lpos,0);
            if(outputonly==0) fgettr(infp,&tra);
        }
		ltrace = jtrace;
		jtrace = ntrpre;
	} else if(tapecntl==-1) {
		ltrace = 0;
		jtrace = ntrpre;
		lpos = nt*sizeof(float)+240;
		lpos = lpos * ntrpre + 3600;
		fseek2g(infp,lpos,0);
		if(outputonly==0) fgettr(infp,&tra);
	} else if(tapecntl==0) {
		ltrace = jtrace;
		jtrace = ntrpre - jtrace;
	} else if(tapecntl==2) {
		ltrace = jtrace;
		jtrace = ntrpre;
	}

	/* set old header for the first trace */
	if(gath==1){
	 	oldheadx = tra.sx;
 	 	oldheady = tra.sy;
 	} else if(gath==2){
	 	oldheadx = tra.gx;
 	 	oldheady = tra.gy;
 	} 
	range = gottrace = 1;
	ready = done = 0;
	nxsb = nysb = nxs0b = nys0b = 0;

	if (outputonly==1) goto no_mig;
    
	do {
 		
		/* if gather has changed or no more input traces */
		if ( ready || (!gottrace && ktrace>0) ) {
  			fxs = fxst+nxs0*dxst;
			fys = fyst+nys0*dyst;

/*
	fprintf(stderr,"fxs=%g nxs=%d fys=%g nys=%d \n",fxs,nxs,fys,nys);
				fprintf(jpfp,"read ttable time=%f \n",walltime());
*/

	    	for(iys=0,ys=fys; iys<nys; ++iys,ys+=dyst) 
  		    for(ixs=0,xs=fxs; ixs<nxs; ++ixs,xs+=dxst){
 				itab = iys*nxs+ixs;
			
 		    	nseek = nxst*(nys0+iys)+nxs0+ixs;
				nseek = nseek*nlt*nst*nzt;
   		    	fseek2g(ttfp,nseek*sizeof(float),SEEK_SET);

 				fread(tt,sizeof(float),nst*nlt*nzt,ttfp);

				/* compute residual t */
				if(dlt<999990) 
      	   			resit_(tt,&nst,&nlt,&nzt,&izw,&nr,&dst,&dlt,
				    	&dr,&fst,&flt,&xs,&ys,tb0);
				else   
      	    		resit0_(tt,&nst,&nzt,&izw,&nr,&dr,
			    		st,lt,&xs,&ys,tb0);
  
				/* intrepolation t to output lateral positions */ 
   				latint_(&nsout,&nlout,&nz,s,l,ttab[itab],
			 		&nst,&nlt,&nzt,&nz0,&fst,&flt,&dst,&dlt,tt);
	      	}
 
			/* migration of ktrace */
			migs(nlout,nsout,mig,nz,
				ktrace,ntl,dldt,trace, 
				tras,dtl,tminl,nt, 
				ss,sl,gs,gl, 
				imutel,s,l,
				iofs,apers,aperl,fold,
				work1,work2,wsave,tracef,sqrtf,
				tahd,ifcut,ltaper,ksmax,klmax,
				f0,df,nf,ftaper,nsave,nfft, 
				dzt,fz,angmax, 
				nxs,nys,fxs,fys,dxst,dyst,
				nzo,fzout,dzo,
				nr,dr,tb,pb,ttab,gath);

			itrace = itrace + ktrace;
			ndskwt = ndskwt + ktrace;
 			ktrace = 0;