kirmigs.c

seismic software,very useful

           for(ix=0;ix<nxt*nst;ix++) {
              fread((char *)trt,sizeof(float),nzt,tfp);
	      ix0 = ix*nzt;
	      for(iz=0;iz<nzt;iz++) {
	         tmp = trt[iz]*tscale;
	         if ( tmp < 32000 ) {
	            ttbl[ix0+iz] = (short)tmp;
	         }
	         else {
	            ttbl[ix0+iz] = 32500;
	         }
	      }
	   }   

           if (dt>=1.) {
              tmin = tmin * 1000.;
              dl = dt * nt / lt;
	      dlm = dl;
           } else {
              tscale = tscale/1000.;
              dt = dt * 1000.;
              tmin = tmin * 1000.;
              dl = dt * nt / lt;
	      dlm = dl;
           }
	} else {
	   dl = dt*nt/lt;
	   dlm = dl * 1000.;
	   tscale = 1.;
	}
	
	

	/* Main loop over traces */
	ix = 0;
	dldt = dl/dt;
	tminl = tmin / dt; 
	ntrace =0;

/* skip lseekin bytes in infp */
	fseek(infp,lseekin,0);

/* see if cdp and offset are used in migration, instead of sx and gx */
/* read first trace */
	fgettr(infp,&tra);
	sy = tra.sy;
	gy = tra.gy;
	if ( sy == gy && dcdp==0. ) {
	   chkcdp = 0;
	} else {
	   chkcdp = 1;
	   if( dcdp == 0. ) {
		fprintf(stderr,"dcdp must be specified ! \n");
		return 1;
	   }
	   icdp1 = cdpx0;
	}
	
	/* read in dips grid */
	
	if(idip==1) {

		tmp = dcdpx*nx/nxdip;
		itmp = (int) tmp;	
		dipsgrid_cpp(dipfile,dips,dips+nz*nxdip,z0,dz,nz,
			cdpx0,itmp,nxdip);
		for(iz=0;iz<nxdip*nz*2;iz++) {
			tmp=dips[iz];
			if(tmp>=90.) {
		 		tmp = 89.9;
			} else if (tmp<=-90.) {
				tmp = -89.9;
			} 
			dips[iz]=tan(tmp*3.141592654/180.);
		}
	}

	/* initialize mtrace */
	*mtrace = 0;

	do {

                int nonzerotrace ;
		/* Load trace into trace (zero-padded) */
		memcpy(trr, tra.data, nt*sizeof(float));

	        /* check if this is a zero trace */
		nonzerotrace = 0 ;	
		for(it=0;it<nt;it++) {
		  if(trr[it] != 0.0) { nonzerotrace = 1; break ; } 
		} 
		if( nonzerotrace == 0 ) tra.trid = 0;
		if( tra.trid != 0 ) {

/* apply 2.5-D filter */
	            if(i2p5==1) f2p5_(trr,&nt);
/* apply agc */
		    if(lwin>0) {
                        mute = (tra.mute-tra.delrt)/tra.dt;
                        if(mute<0) mute=0;
                        if(mute>nt) mute=nt;
                        mt = nt - mute;
                        rmso = 2000.;
                        agc_(trr+mute,&mt,&lwin,trwk,&rmso);
                    }
/* apply tpow */
                    if(fabs(tpow)>0.0001) tp_(trr,&nt,&tpow,&tmin,&dt);


/* linearly interpolate input trace */
	            for(it=0;it<lt;it++) {
		       tmp = tminl+it*dldt; 
	               i = (int)tmp;
	               tmp = tmp - i;
	               if(i>=0 && i<nt-1) {
		          trace[it] = (1.-tmp)*trr[i]+tmp*trr[i+1];
		       }
		       else {
		          trace[it] = 0.;
		       }	
		    }
	

	            /* obtain source and receiver x's from trace hader */
		    if ( chkcdp == 0 ) {
		       sx = tra.sx;
		       gx = tra.gx;	
		       xs = sx;
		       xr = gx;
		       if(tra.scalco>1) {
				xs = xs*tra.scalco;
				xr = xr*tra.scalco;
		       } else if (tra.scalco<0) {
				xs = xs/(-tra.scalco);
				xr = xr/(-tra.scalco);
		       }		
		    }
		    else {
		       xs = (tra.cdp-icdp1)*dcdp + x0 - tra.offset/2.;
		       xr = (tra.cdp-icdp1)*dcdp + x0 + tra.offset/2.;

		    }
	            tmp = (tra.mute-tra.delrt)/dlm + 1.;
	            mute = (int)tmp ; 
		    tmp = tra.mute * 0.5 * v0 * 0.001;
                    tmp = tmp*tmp - 0.25 * fabs(xs-xr) * fabs(xs-xr);
                    if ( tmp > 0. ) {
                        zw = sqrt(tmp);
                    } else {
			zw = 0.;
		    }
		    ix = ix+ 1;

		    /* migration */
		    tmp = fabs(xs-xr);
		    offset = (tmp-ofomin)/dofo+1.5 ;   
		    iofo = (int)offset;
		    if ( (iofo<1 || iofo>nofo) && intype==1 ) break; 

/*
 fprintf(stderr,"xr=%f xs=%f at cdp=%d dcdp=%f icdp1=%d x0=%f offset=%d\n",
	xr,xs,tra.cdp,dcdp,icdp1,x0,tra.offset);	
*/

		}
		ntrace = ntrace + 1;
		if ( ntrace%1000 == 0 )
   	fprintf(stderr,"input %d traces processed at %s \n",ntrace,chost);

	        if(tmp>=ofimin && tmp<=ofimax && tra.trid==1 && mute<lt) {

		   	kirmig_(trace,&xs,&xr,&tmin,&lt,&dl,
                           	migs,&x0,&z0,&dx,&dz,&nx,&nz,
                           	&nofo,&ofomin,&dofo,fold,
                           	ttbl,atbl,&xmint,&zmint,&dxt,&dzt,&nxt,&nzt,
                           	&nst,&smint,&dst,ts,tr,as,ar,&dxtol,
                           	sigxt,sigzt,inxt,inzt,&iamp,mtrace,
                           	twk1,twk2,awk1,awk2,&tscale,trwk,
                           	&intps,&intpx,&intpz,&mute,
			   	&ascale,&amptype,
			   	tfile,afile,&aper,&apanl,&apanr,&zw,
				dips,&idip,&nxdip,&dxdip);
		}

	} while (fgettr(infp,&tra));


   fprintf(stderr,"input %d traces processed at %s \n",ntrace,chost);
	
        /* free spaces */
	free((char *)trt);
	free((char *)trr);
	free((char *)ttbl);
	free((char *)atbl);
	free((char *)twk1);
	free((char *)twk2);
	free((char *)awk1);
	free((char *)awk2);
	free((char *)trwk);
	free((char *)tr);
	free((char *)ts);
	free((char *)ar);
	free((char *)as);
	free((char *)sigxt);
	free((char *)sigzt);
	free((char *)inxt);
	free((char *)inzt);
	free((char *)trace);
	free((char*)dips);

	if(fabs(tpow)>0.0001) {
                trace = (float *) malloc(nz*sizeof(float));
                for(i=0;i<nz;i++) {
                        tmp = (z0+i*dz)/(z0+nz*0.5*dz);
                        tmp = fabs(tmp);
                        if(tmp==0.) {
                                trace[i] = 0.;
                        } else {
                                trace[i] = pow(tmp,-tpow);
                        }
                }
        }

	ascale=1./ascale;
	if ( amptype == 0 ) {
	   ascale = ascale/2.;
	}
	else {
	   ascale=ascale*ascale;
	}
	if ( iamp == 0 ) ascale = 1.;
	
	skipmig:

/* output traces */

	/* skip lseekout bytes in outfp */
	fseek(outfp,lseekout,0);
	
	/* if no input trace, zero output */
	if(notrace==1) {
		bzero((char*)&tra,(240+nz*sizeof(float)));
	}

	/* update trace headers */
	tra.ns = nz;
	tra.trid = 1;

	if ( dz < 30. ) {
	   tmp = dz * 1000.;
           tra.dt = (unsigned short) tmp;
	   tmp = z0 * 1000.;
	   tra.delrt = (unsigned short) tmp;
	} else if (dz <300.) {
	   tmp = dz * 100.;
           tra.dt = (unsigned short) tmp;
	   tmp = z0 * 100.;
	   tra.delrt = (unsigned short) tmp;
	} else {
           tra.dt = (unsigned short) dz;
	   tra.delrt = (unsigned short) z0;
	}
	tra.sy = 0;
	tra.gy = 0;
	/* scale x coordinate if needed */  
	itmp = (int) dx;
	tmp = dx - itmp;
	if(fabs(tmp)>0.01) {
		tra.scalco = -100;
	} else {
		tra.scalco = 1;
	} 

	for(iof=0;iof<nofo;iof++) {
	   offout = ofomin + iof*dofo;
	   if(notrace==0) {
	   	if(fold[iof]>1.) {
	      		tmp = ascale / fold[iof];
	   	}
	   	else {
	      		tmp = ascale;
	   	}
	   }
	   for(ix=0;ix<nx;ix++) {
/* update trace headers */
	      tra.offset = offout;
              tra.cdp = ix+1 ;
              tra.tracf = ix+1 ;
              xs = x0 + ix*dx - offout/2;  
              xr = x0 + ix*dx + offout/2;  
	      if (tra.scalco==-100) {
              	xs = xs * 100.;
              	xr = xr * 100.;
	      }
	      tra.sx = xs;
	      tra.gx = xr;
	      tra.ep = tra.sx;	
	      tra.fldr = tra.sx;	

	      if(notrace==0) {
	      	ix0 = ix*nz+iof*nx*nz;
	      	for (i=0;i<nz;i++) {
	         	trz[i] = migs[ix0+i]*tmp;
	      	}
	      	if(fabs(tpow)>0.0001)
                        for(i=0;i<nz;i++) trz[i] = trz[i]*trace[i];
	      	memcpy(tra.data, trz, nz*sizeof(float));
	      }
	      fputtr(outfp,&tra);
	      /* fflush(outfp); */
	   }
	} 

	if(notrace==0) {
		free((char *)migs);
		free((char *)fold);
		free((char *)trz);
		if(fabs(tpow)>0.0001) free((char *)trace);
	}

	fclose(infp);
	fclose(outfp);

   fprintf(stderr,
	"kirmig done at %s for\t%f < offsets <= %f  for %d live traces\n",
		  chost,ofomin-0.5*dofo,ofomin+(nofo-0.5)*dofo,*mtrace);

	return 0;
}