kzmig.c

seismic software,very useful

                        fldfp = efopen(diskfld,"w");
                        fseek2g(fldfp,0,1);
                        efwrite(fold[0],sizeof(float),
				nz*nlout*nsout*nofo,fldfp);
                        efclose(fldfp);
                        imgfp = efopen(diskimg,"w");
                        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); 

 
	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) {
		fldfp = efopen(diskfld,"w");
		fseek2g(fldfp,0,1);
		efwrite(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)  {
		imgfp = efopen(diskimg,"w");
		fseek2g(imgfp,0,1);
		fwrite(mig[0],sizeof(float),nofo*nsout*nlout*nzo,imgfp);
		efclose(imgfp);
	}

	if(isave==1) fprintf(jpfp," backup to disk done \n");

	ntrpre = ntrpre + ntrace;
	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) {
	     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);
				efread(&tra,sizeof(char),240,hdrfp);

				tra.offset = ofo[io];
				tra.tracl = ix + 1;
				tra.tracr = iy + 1;
				tra.cdpt = io + 1;
				tra.dz = dzo;
				tra.fz = fzo;
				ofs = ofo[io];
				if(tra.trid==1) {
					mx = (tra.gx+tra.sx)/2;
                        		my = (tra.gy+tra.sy)/2;
					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;
                        		}
				}
 
				tra.f1 = fzo;
				tra.ns = nzo;
				tra.d1 = dzo;
				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;
					}
				}

				/* 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, float fys, float dxs, float dys,
	int nzo, float fzo, float dzo,
	int nr, float dr, float *tb, float *pb, float **ttab, int gath) 
{


	int itr, i, ir0, it, io;
	float tmp;

	int ixs,iys,itab;
	float ax,sx,ay,sy,*ts,*tg,*tit,*wt,*ampt;
	int nsl;
	float ssold,slold,gsold,glold;

	nsl = ns * nl;
	/* allocate work arrays	*/
	ts = alloc1float(nsl*nz);
	tg = alloc1float(nsl*nz);
	tit = alloc1float(nz+1);
	wt = alloc1float(nz+1);
	ampt = alloc1float(nz+1);
/*
	tit = alloc1float((nz+1)*nsl);
	wt = alloc1float((nz+1)*nsl);
	ampt = alloc1float((nz+1)*nsl);
*/
 
 
 	/* linearly interpolate input trace */
	for(itr=0;itr<ktrace;itr++) {
		ir0 = itr*nt;
		if(ntl!=nt) {
               		for(it=0;it<ntl;it++) {
                   		tmp = it*dldt;
                       		i = (int)tmp;
                       		tmp = tmp - i;
                   		if(i>=0 && i<nt-1) {
                       			trace[it] = 
						(1.-tmp)*tras[i+ir0]+
						tmp*tras[i+1+ir0];
                       		}
			}	
		} else {
               		for(it=0;it<ntl;it++)
                     		trace[it] = tras[it+ir0];
		}

	   if(itr==0) {
		ssold = ss[0] - 10000;
		slold = sl[0] - 10000;
		gsold = gs[0] - 10000;
		glold = gl[0] - 10000;
	   }

/*
	fprintf(stderr," ss=%g ssold=%g sl=%g slold=%g itr=%d \n",
		ss[itr],ssold,sl[itr],slold,itr); 
*/

	   /* linearly inpterpolate traveltimes to source and receiver */
	    if(fabs(ss[itr]-ssold)>0.1 || fabs(sl[itr]-slold)>0.1) {
		
		ssold = ss[itr];
		slold = sl[itr];

  		ax = (ss[itr]-fxs)/dxs;
		ixs = (int)ax;
   		if(ixs==nxs-1 && nxs>1) ixs = nxs-2;
 		sx = ax-ixs;
		if(sx<=0.01) sx = 0.0;
		if(sx>=0.99) sx = 1.0;
		ay = (sl[itr]-fys)/dys;
		iys = (int)ay;
 		if(iys==nys-1 && nys>1) iys = nys-2;
 		sy = ay-iys;
		if(sy<=0.01) sy = 0.0;
		if(sy>=0.99) sy = 1.0;
		itab = iys*nxs+ixs;
 		if(nxs==1 && nys==1)
	    	    blinint_(&ns,&nl,&nz,&sx,&sy,ttab[itab],
  		      ttab[itab],ttab[itab],ttab[itab],ts);
		else if(nxs==1 && nys>1)
	    	    blinint_(&ns,&nl,&nz,&sx,&sy,ttab[itab],
  		      ttab[itab],ttab[itab+nxs],ttab[itab+nxs],ts);
 		else if(nxs>1 && nys==1)
	    	    blinint_(&ns,&nl,&nz,&sx,&sy,ttab[itab],
  		      ttab[itab+1],ttab[itab],ttab[itab+1],ts);
		else
	    	    blinint_(&ns,&nl,&nz,&sx,&sy,ttab[itab],
  		      ttab[itab+1],ttab[itab+nxs],ttab[itab+nxs+1],ts);
 
	    } 
/*
	fprintf(stderr," gs=%g gsold=%g gl=%g glold=%g itr=%d \n",
		gs[itr],gsold,gl[itr],glold,itr); 
*/
	    if(fabs(gs[itr]-gsold)>0.1 || fabs(gl[itr]-glold)>0.1) {
		
		gsold = gs[itr];
		glold = gl[itr];
	
  		ax = (gs[itr]-fxs)/dxs;
		ixs = (int)ax;
   		if(ixs==nxs-1 && nxs>1) ixs = nxs-2;
 		sx = ax-ixs;
		if(sx<=0.01) sx = 0.0;
		if(sx>=0.99) sx = 1.0;
		ay = (gl[itr]-fys)/dys;
		iys = (int)ay;
 		if(iys==nys-1 && nys>1) iys = nys-2;
 		sy = ay-iys;
		if(sy<=0.01) sy = 0.0;
		if(sy>=0.99) sy = 1.0;
		itab = iys*nxs+ixs;

 		if(nxs==1 && nys==1)
	    	    blinint_(&ns,&nl,&nz,&sx,&sy,ttab[itab],
  		      ttab[itab],ttab[itab],ttab[itab],tg);
		else if(nxs==1 && nys>1)
	    	    blinint_(&ns,&nl,&nz,&sx,&sy,ttab[itab],
  		      ttab[itab],ttab[itab+nxs],ttab[itab+nxs],tg);
 		else if(nxs>1 && nys==1)
	    	    blinint_(&ns,&nl,&nz,&sx,&sy,ttab[itab],
  		      ttab[itab+1],ttab[itab],ttab[itab+1],tg);
		else
	    	    blinint_(&ns,&nl,&nz,&sx,&sy,ttab[itab],
  		      ttab[itab+1],ttab[itab+nxs],ttab[itab+nxs+1],tg);
 	    }
 
 			/* 3d prestack depth migration */
			io = iofs[itr]-1;

/*
   	fprintf(stderr,"io=%i itr=%d imute=%d \n",io,itr,imutel[itr]);  
*/

		pszm3d_(trace,&ntl,&tminl,&dtl,
			&ss[itr],&sl[itr],&gs[itr],&gl[itr],
 			&imutel[itr],mig[io],s,l,&nsl,
			&nzo,&dzo,&fzo,
 			&apers,&aperl,fold[io],
			&f0,&df,&nf,&ftaper,
			ifcut,ltaper,tracef,
			wsave,&nsave,work1,work2,sqrtf,&tahd,&nfft,
			&ksmax,&klmax,&angmax,
			ts,tg,tit,wt,ampt,&nz,&dz,&fz,
			tb,pb,&nr,&dr);
	}

 	free1float(ts);
	free1float(tg);
	free1float(tit);
	free1float(wt);
	free1float(ampt);
 
}

/* compute distance between a point to a line */
/* (x,y) to ax+by+c=0 */

float disp2l(float a, float b, float c,float x, float y) {
        float tmp;
                tmp = fabs(a*x+b*y+c)/sqrt(a*a+b*b);
                return tmp;
}

/* compute a, b, c for (x1,y1) (x2,y2) */
void xy2abc(float x1, float y1, float x2, float y2,
        float *a, float *b, float *c) {
        *a = - (y2 - y1);
        *b = (x2 - x1);
        *c = x1*y2-y1*x2;
}

float disaper(float aperx,float apery,float x1,float y1,float x2,float y2) {

        float a, b, c;
        float x, y;

        xy2abc(x1,y1,x2,y2,&a,&b,&c);
        x = x1 + aperx;
        if((x2>x1 && y2>y1) || (x2<x1 && y2<y1)) {
                y = y1 - apery;
        } else {
                y = y1 + apery;
        }
        return disp2l(a,b,c,x,y);
}