kirmigs.c

seismic software,very useful

/*
  RPC-DECL { kirmigs(in char* datain,out char* dataout,int lseekin,int lseekout,
	             char *timefile,char *ampfile,int iamp,
	   	     int nt,int nx,int nz,int lt,
	             float tmin,float x0,float z0,
	             float dt,float dx,float dz,
	             float smint,float xmint,float zmint,
	             float dst,float dxt,float dzt,
	             int nst,int nxt,int nzt,
	             float ofimin,float ofimax,
	             float ofomin,float dofo,int nofo,
		     inout int *mtrace,int i2p5,
		     int intps,int intpx,int intpz,int intype, float dcdp,
	             int amptype,int mlimit,int lagc,float tpow,
            	     float aper, float apanl, float apanr, 
		     float v0, int cdpx0, int dcdpx, char *dipfile,
		     float dxdip, float dxtol); }
    RPC-MAIN-DECL {
extern "C" {
    void initargs(int argc,char *argv[]) ;
}
    }

    RPC-MAIN-CODE {
	initargs(argc, argv);
    }
*/

#include <stdio.h>
#include <osfcn.h>
#include <math.h>
#include <sys/param.h>
#include "grid.h"
#include "su.h"
#include "segy.h"


segytrace tra;

extern "C" {
void *memcpy(void *trr, void *, int );
void f2p5_(...);
void kirmig_(...);
void agc_(...);
void tp_(...);
void bzero(void *, int );
void dipsgrid_cpp(char *dfile, float *dipl, float *dipr,
        double z0, double dz, int nz, int fcdp, int dcdp, int ncdp);
}



int kirmigs(char *datain, char *dataout, int lseekin, int lseekout, 
	    char *timefile, char *ampfile, int iamp,
	    int nt, int nx, int nz, int lt, 
	    float itmin, float iix0, float iz0, 
	    float idt, float idx, float idz,
	    float ismint, float ixmint, float izmint, 
	    float idst, float idxt, float idzt,
	    int nst, int nxt, int nzt,
	    float iofimin, float iofimax,
	    float iofomin, float idofo, int nofo,
	    int *mtrace, int i2p5,
	    int intps, int intpx, int intpz,int intype, float idcdp,
	    int amptype,int mlimit, int lagc, float itpow,
            float iaper, float iapanl, float iapanr, float iv0, 
	    int cdpx0, int dcdpx, char *dipfile, float idxdip, float idxtol) {

/* subroutine kmigs */

	int ix, iz, it, i;
	float tminl;
	float *trace, *trt, *migs, *fold, *tr, *ts, *ar, *as, *sigxt, *sigzt;
	int *inxt, *inzt;
	short *ttbl, *atbl;
	float ascale, tscale;
	long sx, gx;
	int ntrace, iofo;
	int ix0, iof, latbl;
	float xs, xr, tmp, *trz, *trr,dldt,dl,offout,dlm;
	float *twk1,*twk2,*awk1,*awk2,*trwk;
	FILE *tfp, *afp;
	FILE *infp, *outfp;
	int mute;
	int chkcdp;
	long sy, gy, icdp1;
	float offset;
	char *afile, *tfile;
	int lttbl,llimit;
	int lmem ;
	int lwin, mt, notrace, itmp;
	float rmso, zw;
	float *dips;
	int idip;
	int isize;

        /* convert double to float */
	float tmin = itmin , x0 = iix0 , z0 = iz0 ; 
	float dt = idt , dx = idx , dz =idz ;
	float smint = ismint , xmint = ixmint , zmint = izmint ;
	float dst = idst , dxt = idxt , dzt = idzt ;
	float ofimin = iofimin , ofimax = iofimax ;
	float ofomin = iofomin , dofo = idofo ;
	float dcdp = idcdp;
	float tpow = itpow; 
	float v0 = iv0;
	float aper = iaper;
	float apanl = iapanl;
	float apanr = iapanr;
	float dxdip  = idxdip;
	float dxtol = idxtol;

        char chost[MAXHOSTNAMELEN] ;
	int nxdip;

        gethostname(chost,MAXHOSTNAMELEN) ;
	fprintf(stderr,"starting migration at %s for\t%f < offsets <= %f \n",
			chost, ofomin-dofo*.5,ofomin+(nofo-.5)*dofo);
	
	tmp = nx * dx/dxdip;
	nxdip = (int)tmp;
	if(nxdip<1) nxdip = 1;
	if(nxdip>nx) nxdip = nx;
	dxdip = dx * nx/nxdip;
	 
	lwin = lagc/(int)(dt*1000.);

	/* in case the dp starts the migration again, do not remigrate */ 
	isize = 0;
	if((outfp = fopen(dataout,"r"))!=NULL) {
                 fseek(outfp,0L,SEEK_END);
                 isize= (int) ftell(outfp);
                 fclose(outfp);
        }
	if( isize==(nz*sizeof(float)+HDRBYTES)*nx*nofo ) {
	   fprintf(stderr,"output %s already migrated \n",dataout);
           return 0 ;
	}

	/* set mtrace to be -1 in case of error return to kirmig */
	*mtrace = -1;
	
        infp = fopen(datain,"r");
        if( infp == NULL ) {
	   fprintf(stderr,"kirmigs datain=%s not found\n",datain);
           return 1 ;
        }

        outfp = fopen(dataout,"w");
        if( outfp == NULL ) {
	   fprintf(stderr,"kirmigs dataout=%s not found\n",dataout);
           return 1 ;
        }

	/* check to see if there is any input data */
	efseek(infp,0L,2);
	notrace = 0;
        if((eftell(infp)-lseekin)<=0) {
		warn("no trace in input %s \n",datain); 
		notrace = 1;
		*mtrace = 0;
	}

	if (notrace==1) goto skipmig; 

        /* turn off buffering for outfp */
	/* setbuf(outfp,NULL); */

	/*
	fprintf(stderr,"within kirmigs mlimit=%d \n",mlimit);
	fprintf(stderr,"x0=%f z0=%f nx=%d nz=%d \n",x0,z0,nx,nz);
	fprintf(stderr,"dx=%f dz=%f \n",dx,dz);
	fprintf(stderr,"timefile=%s \n",timefile);
	fprintf(stderr,"ampfile=%s \n",ampfile);
	fprintf(stderr,"nt=%d dt=%f tmin=%f \n",nt,dt,tmin);
	fprintf(stderr,"ofomin=%f dofo=%f nofo=%d \n",ofomin,dofo,nofo);
	fprintf(stderr,"xmint=%f zmint=%f smint=%f \n",xmint,zmint,smint);
	fprintf(stderr,"dxt=%f dzt=%f dst=%f \n",dxt,dzt,dst);
	fprintf(stderr,"nxt=%d nzt=%d nst=%d \n",nxt,nzt,nst);
	fprintf(stderr,"lt=%d \n",lt);
	fprintf(stderr,"dcdp=%f dxtol=%f \n",dcdp,dxtol);
	fprintf(stderr,"lseekin=%d lseekout=%d \n",lseekin,lseekout);
	*/

	idip = 1;
	if(dipfile[0]=='N' && dipfile[1]=='U',
		dipfile[2]=='L' && dipfile[3]=='L') idip=0;

	/* memory allocations */
	lmem = (2*nzt*nxt + nofo*nx*nz + nofo + 4*nx*nz + nx + 2*nz 
		+ nzt + 2*lt + nt)*sizeof(float) 
	      +(nx + nz) * sizeof(int) + idip*2*nxdip*nz*sizeof(float);

	llimit = mlimit * 1024 * 1024;	
	if ( lmem > llimit ) {
	   fprintf(stderr,"Need at least memory size mlimit=%d (Bytes)\n",
	 	   lmem);
	   return 1;
	} 

        twk1 = (float*) malloc(nzt*nxt*sizeof(float));
	if( twk1 == 0 ) return 1 ;
        twk2 = (float*) malloc(nzt*nxt*sizeof(float));
	if( twk2 == 0 ) return 1 ;
	migs = (float*) malloc(nofo*nx*nz*sizeof(float));
	if( migs == 0 ) return 1 ;
	fold = (float*) malloc(nofo*sizeof(float));
	if( fold == 0 ) return 1 ;
	ts = (float*) malloc(nx*nz*sizeof(float));
	if( ts == 0 ) return 1 ;
	tr = (float*) malloc(nx*nz*sizeof(float));
	if( tr == 0 ) return 1 ;
	as = (float*) malloc(nx*nz*sizeof(float));
	if( as == 0 ) return 1 ;
	ar = (float*) malloc(nx*nz*sizeof(float));
	if( ar == 0 ) return 1 ;
	sigxt = (float*) malloc(nx*sizeof(float));
	if( sigxt == 0 ) return 1 ;
	sigzt = (float*) malloc(nz*sizeof(float));
	if( sigzt == 0 ) return 1 ;
	inxt = (int*) malloc(nx*sizeof(int));
	if( inxt == 0 ) return 1 ;
	inzt = (int*) malloc(nz*sizeof(int));
	if( inzt == 0 ) return 1 ;

	trt = (float*) malloc(nzt*sizeof(float));
	if( trt == 0 ) return 1 ;
	trz = (float*) malloc(nz*sizeof(float));
	if( trz == 0 ) return 1 ;
	trace = (float*) malloc(lt*sizeof(float));
	if( trace == 0 ) return 1 ;
        trwk = (float*) malloc(lt*sizeof(float));
	if( trwk == 0 ) return 1 ;
	trr = (float*) malloc(nt*sizeof(float));
	if( trr == 0 ) return 1 ;
	dips = (float*) malloc(idip*2*nxdip*nz*sizeof(float));
	if( dips == 0 ) return 1 ;


	lttbl = nst*nxt*nzt*sizeof(short);
	if ( lmem + lttbl > llimit ) {  
	   lttbl = 1 * sizeof(short);
	   tfile = (char*) malloc(strlen(timefile)+1);
	   sprintf(tfile,"%s\0",timefile);
	}
	else {
	   tfile = "null";
	}
	lmem +=lttbl;
	ttbl = (short*) malloc(lttbl);
	if( ttbl == 0 ) return 1 ;
        if ( iamp == 1 ) {
           latbl = nst*nxt*nzt*sizeof(short);
	   if ( lmem + latbl > llimit ) {  
	      latbl = 1 * sizeof(short);
	      afile = (char*) malloc(strlen(ampfile)+1);
	      sprintf(afile,"%s\0",ampfile);
	   }
	   else {
	      afile = "null";
	   }
           awk2 = (float*) malloc(nzt*nxt*sizeof(float));
	   lmem += nzt*nxt*sizeof(float); 
	   if( awk2 == 0 ) return 1 ;
           awk1 = (float*) malloc(nzt*sizeof(float));
	   lmem += nzt*sizeof(float); 
	   if( awk1 == 0 ) return 1 ;
        } else {
           latbl = 1 * sizeof(short);
	   afile = "null";
           awk2 = (float*) malloc(latbl);
	   lmem += latbl; 
	   if( awk2 == 0 ) return 1 ;
           awk1 = (float*) malloc(latbl);
	   lmem += latbl; 
	   if( awk1 == 0 ) return 1 ;
        }
	atbl = (short*) malloc(latbl);
	lmem += latbl; 
	if( atbl == 0 ) return 1 ;
	fprintf(stderr,"total memory used (Byte) =%d \n",lmem); 
	if(lttbl==1) fprintf(stderr,"travel time table disk i/o used \n"); 
	if(iamp==1 && latbl==1) 
	   fprintf(stderr,"amplitude table disk i/o used \n");

	/* read in travel time table and amplitude table */
	tfp = fopen(timefile,"r");
        if( tfp == NULL ) {
	   fprintf(stderr,"kirmigs: timefile=%s not found\n",timefile);
           return 1 ;
        }

	if ( iamp == 1 && latbl > 1 ) { 
	   afp = fopen(ampfile,"r");
           if( afp == NULL ) {
	       fprintf(stderr,"kirmigs: ampfile=%s not found\n",ampfile);
               return 1 ;
           }
	   /* find scale to scale amplitudes */
	   ascale = 0.;
           for(ix=0;ix<nxt*nst;ix++) {
              fread((char *)trt,sizeof(float),nzt,afp);
	      for(iz=0;iz<nzt;iz++) {
	         if(fabs(trt[iz]) > ascale) ascale=fabs(trt[iz]); 
	      }
	   }
	   /* read in amplitudes, scale and store in short *atbl */ 
	   fseek(afp,0,0);
	   if(ascale>0.) ascale = 32000./ascale;
           for(ix=0;ix<nxt*nst;ix++) {
              fread((char *)trt,sizeof(float),nzt,afp);
	      ix0 = ix*nzt;
	      for(iz=0;iz<nzt;iz++) {
	      tmp = trt[iz]*ascale; 
	      atbl[ix0+iz] = (short)tmp;
	      }
	   }
	}
	else {
	   ascale = 1.;
	}
	
	if(lttbl > 1 ) {
	   /* read in times, scale and store in short *ttbl */ 
           /* times in ms */
	   if (dt>=1.) {
              tscale = 32000./(nt*dt);
           }
           else {
              tscale = 1000. * 32000. / (nt*dt*1000.);
           }