rmig.c

seismic software,very useful

				ep,strace,ng,datain,xtfile,xtnew);
                        	system(cmd);
				bzero(cmd,1024);
                                sprintf(cmd,
"xzedit xtfile=%s hfile=%s mmap=%d aper=%f sortz=%d xtnew=%s vcid=%d &\n",
				xtfile,hfile,mmap,aper,sortz,xtnew,vcid);
                        	system(cmd);
			}

*/
			/* edit rays */
			/* SeisView rayedit */
/* */

			if(confirm==1 && nsnext<=1 ) {
				if(irayedit==0) {
					irayedit = 1;
					sdfp = efopen(shotdata,"w");
					bh.hns = nt;
					bh.ntrpr = ng;
					auxputhdr(sdfp, &ch, &bh);
					bcopy((char*)&tr,(char*)&tro,240);
					for(ig=0;ig<ng;ig++) {
						tro.sx = sx[ig];
						tro.gx = gx[ig];
						tro.offset = tro.gx - tro.sx;
						for(it=0;it<nt;it++) {
							tro.data[it]=
								trace[ig*nt+it];
						}
						auxputtr(sdfp,&tro);
					}
					fclose(sdfp);
					bzero(cmd,1024);
					sprintf(cmd, 
	"rayedit rayedit.d hfile=%s csmodel=%s data=%s 2> /dev/null",
						hfile,xtnew,shotdata);
					sdfp = epopen(cmd,"w");
				} else {
					auxputhdr(sdfp, &ch, &bh);
					bcopy((char*)&tr,(char*)&tro,240);
					for(ig=0;ig<ng;ig++) {
						tro.sx = sx[ig];
						tro.gx = gx[ig];
						tro.offset = tro.gx - tro.sx;
						for(it=0;it<nt;it++) {
							tro.data[it]=
								trace[ig*nt+it];
						}
						auxputtr(sdfp,&tro);
					}
					fflush(sdfp);
				}
			}

/* */

			if(confirm==1 && nsnext<=1 ) {
				raywait:
           if(!svUIi(uiDone)) goto raywait;
        /*
                		printf(" Have you finished rayedit? (y/n)->");
                                scanf("%1s",&ayes) ;
        */
				/*
				tty = efopen("/dev/tty","r");
                	fprintf(stderr, " Have you finished rayedit? (y/n)->");
                		ayes = fgetc(tty);
                		efclose(tty);
				*/
        /*
				   if(ayes=='n') goto raywait; 
        */
			}

		}
		
		remigrate:

		bzero(nxt,maxne*sizeof(int));
		if(ishotout==1) {
			bzero(migshot,nz*nx*sizeof(float)); 
			bzero(foldshot,nz*nx*sizeof(float)); 
			bzero(ixshot,nx*sizeof(int)); 
		}
		if(ixt==0) {
			xtnewfp = efopen(xtnew,"r");
        		getxrt(xtnewfp, maxnr, &ne, nxt, xs, ts, xhs, zhs, 1);
			efclose(xtnewfp);
		} else {
        		/* read xt file at this shot */
        		getxrt(xtfp, maxnr, &ne, nxt, xs, ts, xhs, zhs, spos);
		}

		/*
                fprintf(stderr, " after getxrt \n");
		*/

		/* event selections at receiver locations */
		maxnes = ne * maxar; 
		xe = (float*) malloc(maxnes*ng*sizeof(float));
		ze = (float*) malloc(maxnes*ng*sizeof(float));
		te = (float*) malloc(maxnes*ng*sizeof(float));
		iee = (int*) malloc(maxnes*ng*sizeof(int));
		ier = (int*) malloc(maxnes*ng*sizeof(int));

		bzero(iray,maxnr*maxne*sizeof(int));
		for(ie=0;ie<ne;ie++) {
			for(ir=0;ir<nxt[ie];ir++) {
				iray[ie*maxnr+ir] = 1;
			}
		}

		/*
                fprintf(stderr, " before evsel \n");
		*/

		evsel(ng, ne, maxnr, maxnes, mmap,
        		sx, sort,
			xs, ts, xhs, zhs, nxt, iray, 
        		xe, ze, te, nes, ier, iee,
			sortz, aper);

		/*
                fprintf(stderr, " after evsel \n");
		*/

                /*fprintf(stderr, "\n");*/

		/* event-mapping depth migration */
		for(ig=0;ig<ng;ig++) {
			jg = idg[ig];

			taper(nt,t0[jg],dt,trace+jg*nt,nes[ig],te+ig*maxnes,
				lpass,ltaper);

			if(ishotout==0 ) {
			  	emig_(&nx,&nz,&x0,&dx,mig,&nt,&t0[jg],
				&dt,&tmute[jg],trace+jg*nt,scale,&nes[ig],
				xe+ig*maxnes,ze+ig*maxnes,te+ig*maxnes,zmig,
				&tm[jg],fold,index,&mwd,&sx[jg],&gx[jg],
				ixlive,&mapmax,&aper);
			} else if(ishotout==1) {
			  	emig_(&nx,&nz,&x0,&dx,migshot,&nt,&t0[jg],
				&dt,&tmute[jg],trace+jg*nt,scale,&nes[ig],
				xe+ig*maxnes,ze+ig*maxnes,te+ig*maxnes,zmig,
				&tm[jg],foldshot,index,&mwd,&sx[jg],&gx[jg],
				ixshot,&mapmax,&aper);
			}
		}

		free(xe);
		free(ze);
		free(te);
		free(iee);
		free(ier);

		if(ishotout==1) {
			itrps = 0;
			for(ix=0;ix<nx;ix++) {
				if(ixshot[ix]==1) {
					itmp = ix*nz;
					for(iz=0;iz<nz;iz++) {
						mig[iz+itmp] +=
							migshot[iz+itmp];
						fold[iz+itmp] +=
							foldshot[iz+itmp];
					}
					ixlive[ix] = 1;

					for(iz=0;iz<nz;iz++) {
						if(foldshot[iz+itmp]>1.) {
							tro.data[iz] = 
							migshot[iz+itmp]/
							foldshot[iz+itmp];
						} else {
							tro.data[iz] = 
							migshot[iz+itmp];
						}
					}
					icount += 1;
					itrps += 1;
					tro.tracl = icount;
					tro.tracr = itrps;
					tro.ep = ep;
					tro.fldr = ep;
					tro.cdp = fcdp + ix*dcdp;
					tro.sx = sxx;
					tro.gx = 2*(x0+ix*dx)-sxx;
					tro.offset = tr.gx - tr.sx;
					auxputtr(shotfp,&tro);
				}
			}
		}

		/* view the shot migrated, apply mute, etc. */
		if(confirm==1 && ( nscount==0 || is==ns-1) ) {

			if(imigedit==0) {
				imigedit = 1;
				/* update grid headers */
				ughupdate(&ughshot,nz,nx,dz,dx,z0,x0,dcdp,fcdp,
					ixlive,&isht0,&nsht);

				/* write to shot grid */
				datafp = efopen(mmshot,"w");
				efwrite(mig+isht0*nz,sizeof(float),
					nz*nsht,datafp);
				efwrite(fold+isht0*nz,sizeof(float),
					nz*nsht,datafp);
				fputusghdr(datafp,&ughshot);
				efclose(datafp);

				/* pick, apply migration mute and stack */

                        	bzero(cmd,1024);
                        	sprintf(cmd,
		     "migedit data=%s hfilein=%s shot=%s 2>/dev/null",
                                	mmfile,hfile,mmshot);
				datafp = epopen(cmd,"w"); 
			} else {
				/* update grid headers */
				ughupdate(&ughshot,nz,nx,dz,dx,z0,x0,dcdp,fcdp,
					ixlive,&isht0,&nsht);
				fputusghdr(datafp,&ughshot);
				efwrite(mig+isht0*nz,sizeof(float),
					nz*nsht,datafp);
				efwrite(fold+isht0*nz,sizeof(float),
					nz*nsht,datafp);
				fflush(datafp);
			}

			bzero(mig,nx*nz*sizeof(float));
			bzero(fold,nx*nz*sizeof(float));
			bzero(ixlive,nx*sizeof(int));
		
			/* if re-migration is needed */
      			if(svUIi(uiReMigrate)) goto remigrate;

      /*
                	printf(
    		" Remigrate of the current shot, after Rayedit? (y/n) ->");
                        scanf("%1s",&ayes) ;
			if(ayes=='y') goto remigrate;
      */

			/* if re-ray tracing is needed */
      			if(svUIi(uiReRayTrace)) goto reraytrace;
      /*
                	printf(
    		" Ray tracing again for the current shot? (y/n) ->");
                        scanf("%1s",&ayes) ;
			if(ayes=='y') goto reraytrace;
      */




			/* update the history file */
			fseek(hisfp,0,0);
			fprintf(hisfp,"Number of Shots Migrated: %d\n",is+1);
			fprintf(hisfp,"Trace Position of Next Shot: %d\n",
				strace+ng);
			fprintf(hisfp,"Source Point of Next Shot: %d\n",
				ep+dep);
			fflush(hisfp);

			/* check with user about the migedit */
         		nsnext=svUIi(uiNsNext);
      /*
                	printf(
    	" Number of Shots to Migrate Next (AFTER SAVING MIGEDIT RESULT!!) ->");
			scanf("%d",&nsnext);
      */
			/*
                	fprintf(stderr,
    	" Number of Shots to Migrate Next (AFTER SAVING MIGEDIT RESULT!!) ->");
			tty = efopen("/dev/tty","r");
                	gets(nshs);
                	efclose(tty);
			*/

			/* clean up xt files */
			rmfile(xtfile);
			rmfile(xtnew);

			if(nsnext==0) goto quit; 
			nscount = nsnext;
		} else if(confirm==0 && ixt==0) {
			/* clean up xt files */
			rmfile(xtfile);
			rmfile(xtnew);
		}

		if (nsnext==1) {
			nscount = 0;
		} else {
			nscount = nscount - 1;
		}


		free(nes);
		free(idg);
		free(sort);

		fprintf(stderr,"migration done for %d traces at shot %d \n",
				ng,is+1); 

	}

	quit:


	if(confirm==1) {
		askagain:
      if(!svUIi(uiAskAgain)) goto askagain;
     /*
                printf(" Have you exited from migedit? (y/n)->");
		scanf("%1s",&ayes);
    */
		/*
		tty = efopen("/dev/tty","r");
                fprintf(stderr, " Have you exited from migedit? (y/n)->");
                ayes = fgetc(tty);
                efclose(tty);
		*/
    /*
		if(ayes=='n') goto askagain;
    */
		/*
		epclose(sdfp);
		epclose(datafp);
		*/
		datafp = efopen(mmfile,"r");
		efread(mig,sizeof(float),nz*nx,datafp);
		efread(fold,sizeof(float),nz*nx,datafp);
		efclose(datafp);

		unlink(mmshot);
		unlink(shotdata);
	}

	
	for(iz=0;iz<nx*nz;iz++) {
		if(fold[iz]>1.) mig[iz] = mig[iz]/fold[iz];
	}

	if(iout==1) {

		outfp = fopen(dataout,"w");
		bh.fold = 1;
		bh.tsort= 4;
		bh.hns = nz;
		bh.ntrpr = nx;
		fputhdr(outfp, &ch, &bh);
		tr.ns = nz;
		tr.dz = dz;
		tr.fz = z0;
		tr.mute = 0;
		tr.muts = 0;
		tr.trid = 1;
		tr.scalco = 1;
		tr.scalel = 1;
		tr.duse = 1;
		tr.delrt = z0;
		tr.cdpt = 1;
		tr.offset = 0;

		for(ix=0;ix<nx;ix++) {
			tr.tracl = ix + 1;
			tr.tracr = ix + 1;
			tr.ep = ix + 1;
			tr.fldr = ix + 1;
			tr.cdp = fcdp + ix*dcdp;
			tr.sx = x0 + ix*dx;
			tr.gx = tr.sx;
		
			for(iz=0;iz<nz;iz++) tr.data[iz] = mig[iz+ix*nz];
			fputtr(outfp,&tr);
		}
		fclose(outfp);
	}
        return 0 ;
}


/* update the unscaled header of a grid, and return pointer position to the
grid array */

void ughupdate(usghed *ugh,int nz,int nx,float dz,float dx,float z0,float x0,
	int dcdp,int fcdp,int *ixlive,int *i2live,int *n2live) {

		int i2, n2, i20;
		float o2;
			bzero((char*)ugh,100);
                        ugh->scale = 1.e-6;
                        ugh->dtype = 4;
                        ugh->n1 = nz;
                        ugh->n2 = nx;
                        ugh->n3 = 2;
                        ugh->n4 = 1;
                        ugh->n5 = 1;
                        ugh->d1 = dz;
                        ugh->d2 = dx;
                        ugh->d3 = 1.;
                        ugh->o1 = z0;
                        ugh->o2 = x0;
                        ugh->o3 = 1.;
                        ugh->dcdp2 = dcdp;
                        ugh->ocdp2 = fcdp;

			i20 = 0;
			
			for(i2=0;i2<nx;i2++) {
				if(ixlive[i2]==1) {
					i20 = i2;
					break;
				}
			}

			n2 = nx - 1;
			for(i2=nx-1;i2>=0;i2--) {

				if(ixlive[i2]==1) {
					n2 = i2;
					break;
				}
			}

			n2 = n2 - i20 + 1;

			o2 = x0 + i20*dx;
                        ugh->o2 = o2;
                        ugh->n2 = n2;
                        ugh->ocdp2 = fcdp + i20*dcdp;
	
			*i2live = i20;
			*n2live = n2;
}

int svUIi(char * dFile){
  char cmd[1024];
  int iRet;
  FILE * fStdOut;

  sprintf (cmd, "testui %s", dFile);
  fStdOut=popen(cmd,"r");
  fscanf(fStdOut,"%d", &iRet);
  pclose(fStdOut);
  return iRet;
} 


void rmfile(char *fname) {
	char cmd[1024];
	bzero(cmd,1024);
  	sprintf (cmd, "/bin/rm -f %s", fname);
	system(cmd);
}