kzmig.c.9.2.99

seismic software,very useful

			*/
			continue;
		}

 		/* update minimum and maximum (s,l) coordinates	*/
		if(ss[ktrace]<ssmin) ssmin = ss[ktrace];
		if(gs[ktrace]<ssmin) ssmin = gs[ktrace];
		if(sl[ktrace]<slmin) slmin = sl[ktrace];
		if(gl[ktrace]<slmin) slmin = gl[ktrace];
		if(ss[ktrace]>ssmax) ssmax = ss[ktrace];
		if(gs[ktrace]>ssmax) ssmax = gs[ktrace];
		if(sl[ktrace]>slmax) slmax = sl[ktrace];
		if(gl[ktrace]>slmax) slmax = gl[ktrace];

		/* find traveltime tables for this gather */
		nxs0 = (ssmin-fxst)/dxst;
		nys0 = (slmin-fyst)/dyst;
		nxs = 1+(int)((ssmax-fxst)/dxst+0.99)-nxs0;
		nys = 1+(int)((slmax-fyst)/dyst+0.99)-nys0;
 
 		if(nys*nxs>ntab){
			range = 0;
			nxs = nxsb;
			nys = nysb;
			nxs0 = nxs0b;
			nys0 = nys0b;
			continue;
		}else{
  			nxsb = nxs;
			nysb = nys;
			nxs0b = nxs0;
			nys0b = nys0;
		}


		iofs[ktrace] = io + 1;

 		/* apply tpow */
                if(itgain==1) for(it=0;it<nt;it++) tra.data[it] *= tgain[it];
		for(it=0;it<nt;it++) tras[it+ktrace*nt] = tra.data[it];

			
		itmp = (tra.mute-tra.delrt);
 		itmp = itmp*1000/(int)tra.dt*ifact + 1; 
		if(itmp<1) itmp = 1;
		if(itmp>ntl) itmp = ntl; 
		imutel[ktrace] = itmp;
 
 
/*
fprintf(stderr,"il=%d is=%d ml=%g ms=%g cdp=%d offset=%d io=%d tracl=%d\n",
		il,is,ml,ms,tra.cdp,tra.offset,io,tra.tracl); 
*/

		if(il>=0 && il<nlout && is>=0 && is<nsout) {
            if(ikey==1) {
				ms = s[il*nsout+is];
				ml = l[il*nsout+is];
				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(&tra,indxtrk,&trkval);
				puthval(&tra,indxlnk,&lnkval);
			}
			if(ntrend==0) {
			   if(abs(abs(tra.offset)-ofo[io])
                                <abs(off[io+is*nofo+il*nofo*ns]-ofo[io])) {
				off[io+is*nofo+il*nofo*ns] = abs(tra.offset);
				lpos=io+is*nofo+il*nofo*ns;
				lpos=lpos*240;
				fseek2g(hdrfp,lpos,0);
 				fwrite(&tra,sizeof(char),240,hdrfp);
 			   } 
			} else {
				if(lend==lstart) {
					ps = ms;
					pl = lstart;
				} else {
					tmp = (send-sstart)/(lend-lstart);
					pl = ml + tmp*(lstart*tmp+ms-sstart); 
					pl = pl/(1.+tmp*tmp);
					ps = tmp*(pl-lstart) + sstart;
				}
				if(dsout!=0.) {
					tmp = (ps - sstart)/dsout + 0.5;
					isend = tmp; 
				} else { 
					tmp = (pl - lstart)/dlout + 0.5;
					isend = tmp;
				}
				if (isend>=0 && isend<ntrend) {
				     tmp = sqrt( (ms-s[isend])*(ms-s[isend])
				         + (ml-l[isend])*(ml-l[isend]) );
				     if (tmp<disend[io+isend*nofo]) { 
					 	disend[io+isend*nofo] = tmp;
					 	tra.ungpow  = tmp;
				     	lpos=io+isend*nofo;
				     	lpos=lpos*240;
					 	fseek2g(hdrfp,lpos,0);
 				     	if(fwrite(&tra,sizeof(char),240,hdrfp)!=240)
							err(" error write %s \n",diskhdr);
				      }
			        }  
			}
		}

 		ktrace += 1;

		if(jtrace/1000>m1000) { 
		/*
		gettime(datetime);
		fprintf(jpfp," processing done for %d input traces at %s\n",
			jtrace,datetime);
		*/
		fprintf(jpfp," processing done for %d input traces \n",
			jtrace);

		m1000 = jtrace/1000;
		fflush(jpfp);
		}
	    }while(done);

/*
fprintf(stderr,"jtrace=%i,range=%i ntrace=%d \n",jtrace,range,ntrace);
*/
	
 		if(range){
			if(ntrace==ntras) gottrace = 0;
			else if(!fgettr(infp,&tra))   
				gottrace = 0;
		} else {
			jtrace -= 1;
			ntrace -= 1;
		}

		if(gottrace) {
			if(gath==1){
	 			newheadx = tra.sx;
	 			newheady = tra.sy;
  			} else if(gath==2){
	 			newheadx = tra.gx;
	 			newheady = tra.gy;
 			}
			if(ktrace==mtrace || !range) ready = 1;
			if(newheadx!=oldheadx || newheady!=oldheady){
				oldheadx = newheadx;
				oldheady = newheady;
				ready = 1;
 			}
			if(ready){  
				range = 1;
   				ssmin = slmin = 99999999.0;
				ssmax = slmax = -99999999.0;
			}
			if(!ktrace) ready = 0;
		}
	
		if(ndskwt>=ntrdsk && ntrdsk>0) {
                /* update disk file */
                fldfp = efopen(diskfld,"r+");
                fseek2g(fldfp,0,1);
                fwrite(fold[0],sizeof(float),
				nz*nlout*nsout*nofo,fldfp);
                efclose(fldfp);
                imgfp = efopen(diskimg,"r+");
                fseek2g(imgfp,0,1);
                fwrite(mig[0],sizeof(float),nofo*nsout*nlout*nzo,imgfp);
                efclose(imgfp);
                fprintf(hffp,"Number of Traces Processed: %d %d \n",
                                jtrace,jtrace-ntrpre+ltrace);
                        fflush(hffp);
                        ndskwt=0;
                }

	} while(!done); 

	no_mig:
 
 /*
	gettime(datetime);
	fprintf(jpfp," %d input traces read in this run at %s\n",ntrace,datetime);
	fprintf(jpfp," %d input traces migrated in this run at %s\n",
		itrace,datetime);
*/
	fprintf(jpfp," %d input traces read in this run \n",ntrace);
	fprintf(jpfp," %d input traces migrated in this run \n",
		itrace);

	fflush(jpfp);

	efclose(hdrfp);

	if(isave==1 && ntrace>0 && outputonly==0 ) {
		fldfp = efopen(diskfld,"r+");
		fseek2g(fldfp,0,1);
		fwrite(fold[0],sizeof(float),nlout*nsout*nz*nofo,fldfp);
		efclose(fldfp);
	} 

	if(isave==0) unlink(diskfld);

	if(ftell(btfp)==0) { 
		efclose(btfp);
		unlink("BAD_TRACE_FILE");
	} else {
		efclose(btfp);
	}

	if(isave==1 && ntrace>0 && outputonly==0 )  {
		imgfp = efopen(diskimg,"r+");
		fseek2g(imgfp,0,1);
		fwrite(mig[0],sizeof(float),nofo*nsout*nlout*nzo,imgfp);
		efclose(imgfp);
	}

	if(isave==1 && ntrace>0 && outputonly==0 ) 
		fprintf(jpfp," backup to disk done \n");

	ntrpre = ntrpre + ntrace;
	if(outputonly==0) 
		fprintf(hffp,"Number of Traces Processed: %d 0 \n",ntrpre);
	efclose(hffp);
	if(ihis==0) unlink(hisfile);


	/* output */
	hdrfp = efopen(diskhdr,"r");
	fseek2g(hdrfp,0,1);

	if((ntrace == ntras || feof(infp)!=0 || outputonly==1 ) && traceout==1)  {
	     fprintf(jpfp," start output ... \n");

	     if(izgain==1) {
	     	for(iz=0; iz<nzo; iz++){
			tmp = fzo+iz*dzo;
			zgain[iz] = pow(tmp,zpow);
	     	} 
	     } else {
	     	for(iz=0; iz<nzo; iz++){
			tmp = fzo+iz*dzo;
			zgain[iz] = 1.0;
	     	} 
             }			

  	    for(iy=0;iy<nlout;iy++) {
		for(ix=0;ix<nsout;ix++) {
			for(io=0;io<nofo;io++) {
				itmp = iy*nsout+ix;
				itmp2 = itmp*nzo;
 				bcopy(mig[io]+itmp2,tra.data,nzo*sizeof(float));
				itmp2 = itmp*nz;
				if(flexbin==0) {
  					for(iz=0;iz<nz;iz++) 
				    		fold[io][itmp2+iz] =
 						sqrt(fold0+fold[io][itmp2+iz]);
  					for(iz=0;iz<nzo;iz++){
				    		zi = (fzo+iz*dzo-fz)/dzt;
				    		izi = zi;
				    		if(izi>nz-2) izi = nz-2;
 				    		scale = (1.0-zi+izi)*
							fold[io][itmp2+izi]+
							(zi-izi)*
							fold[io][itmp2+izi+1]; 
				    		tra.data[iz] = tra.data[iz]/
							scale*zgain[iz];
					}
				} else {
  					for(iz=0;iz<nzo;iz++){
				    		tra.data[iz] = tra.data[iz]*
								zgain[iz];
					}
				}

				/* update headers */
				bzero(&tra,240);
				if(fread(&tra,sizeof(char),240,hdrfp)!=240)
					err(" error reading %s \n",diskhdr);

				tra.offset = ofo[io];
				tra.tracl = ix + 1;
				tra.tracr = iy + 1;
				tra.cdpt = io + 1;
				if(dzo<=32.) { 
					tra.dt = dzo*1000+0.5;
				} else {
					tra.dt = dzo;
				}
				tra.delrt = fzo;
				tra.dz = dzo;
				tra.fz = fzo;
				ofs = ofo[io];
				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(tra.trid==1) {
					tmp = tra.gx-tra.sx;
					dd = tra.gy-tra.sy;
					tmp = tmp*tmp + dd*dd;
              		dd = sqrt(tmp);
                   	if(dd>0.) {
                   		tra.sx = mx+ofs*(tra.sx-mx)/dd;
                       	tra.gx = mx+ofs*(tra.gx-mx)/dd;
                       	tra.sy = my+ofs*(tra.sy-my)/dd;
                       	tra.gy = my+ofs*(tra.gy-my)/dd;
                    } else {
                    	tra.sx = mx-ofs/2;
                       	tra.gx = mx+ofs/2;
                       	tra.sy = my;
                       	tra.gy = my;
                    }
				} else {
                    tra.sx = mx-ofs/2;
                    tra.gx = mx+ofs/2;
                    tra.sy = my;
                    tra.gy = my;
				}
 
				tra.ns = nzo;
				tra.mute = 0;
 
 				if(ntrend>0) {
					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) {
						tra.cdp = iiy*ncdppl+iix+cdp1;
					} else {
						tra.cdp = iix*nlines+iiy+cdp1;
					}
				}

                if(ikey==1) {
					ms = s[iy*nsout+ix];
					ml = l[iy*nsout+ix];
					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(&tra,indxtrk,&trkval);
					puthval(&tra,indxlnk,&lnkval);
				}

				/* output */
				fputtr(outfp,&tra);
			}
		}
	    }
	    fprintf(jpfp," output done \n");
	    fflush(jpfp);
	} 
	
	efclose(outfp);
	

	efclose(hdrfp);
	if(isave==0) unlink(diskhdr);
	if(isave==0) unlink(diskimg);


	free(mig);
	free2float(ttab);

	/* backup to tape */
	if(ibacko==1) {
                fprintf(jpfp," Backup %s , %s and %s to %s \n",
                        diskimg,diskhdr,diskfld,backupo);
		tar3to(backupo,diskimg,diskhdr,diskfld);
        }

	return 0;
}



void migs(int nl, int ns, float **mig, int nz,
	int ktrace, int ntl, float dldt, float *trace, 
	float *tras, float dtl, float tminl, int nt, 
	float *ss, float *sl, float *gs, float *gl, 
	int *imutel, float *s, float *l,
	int *iofs, float apers, float aperl, float **fold,
	float *work1, float *work2, float *wsave, float *tracef, float *sqrtf,
	int tahd, int *ifcut, int *ltaper, float ksmax, float klmax,
	float f0, float df, int nf, float ftaper, int nsave, int nfft,
	float dz, float fz, float angmax,
	int nxs, int nys, float fxs,