fxymig.c

seismic software,very useful

	/* read in velocity if needed */
	if(nv==nvs) {
		if((velfp=fopen(velfile,"r"))==NULL)
			err("Input Velocity File velfile Error");
		fseek64(velfp,0,0);
	  	fread(v,sizeof(float),nx*ny*nvs,velfp);
		fclose(velfp);	
	} else {
		fclose(velfp);
	}

	/* compute frequencies to migrate */ 
	for(iw=0;iw<nf;iw++) om[iw] = wmin +iw*dw;
 
/*
	dump2xplot(v,nx,nvs,0,"v plot");
	dump2xplot(cp,nx*2,nf,1,"cp plot");
*/

	fprintf(jpfp,"\n");
	fprintf(jpfp," Start Migration \n");
	fprintf(jpfp," =============== \n");
	fflush(jpfp);

	fflush(jpfp);

    	/* call migration routine */

	lenvel = strlen(velfile);
	lenxyt = strlen(diskxyt);
	lenxyw = strlen(diskxyw);
	
	namevel = (char*) emalloc(lenvel+1);
	namexyt = (char*) emalloc(lenxyt+1);
	namexyw = (char*) emalloc(lenxyw+1);

	bcopy(velfile,namevel,lenvel);
	bcopy(diskxyt,namexyt,lenxyt);
	bcopy(diskxyw,namexyw,lenxyw);

	namevel[lenvel]='\0';
	namexyt[lenxyt]='\0';
	namexyw[lenxyw]='\0';


	if(dz==0.) dvs = dvs * 0.001;

        fclose(hdrfp);

	if(lenjpf>1) fclose(jpfp);

    wxymig_(&nx,&ny,&ntau,&nw,&itau0,&iestep,&icstep,
		&nvs,&nkx,&nky,&nq,&nv,&ncp,&lvec,&lplane,&iorder,
		&naux1,&naux2,&naux,
		&dt,&dx,&dy,&dtau,&dfc,&vref,&dvs,
		q,v,vxy,
		om,oms,vyx,w,
		aux1,aux2,
		a,aa,r,cpt,
		caux,shift,al,ar,
		cp,bcp,cpp,
		namexyw,namexyt,namevel,
		&lenxyw,&lenxyt,&lenvel,
		asave,aasave,&iisave,va,
		namejp,&lenjpf,&ncpu,
		namexywb,namexytb,
		&lenxywb,&lenxytb,
		qbuf,cpbuf,cph,dzov,&ncph,
		namedur,&lendur,&itaudur,&iwdur,
		&ntauend,&iwend,&nwmig,&nf,&ifminr);

	if(lenjpf>1) jpfp = efopen(jpfile,"a");

	fprintf(jpfp," =============== \n");
	fprintf(jpfp," Migration Done\n");
	fprintf(jpfp,"\n");

	/* free space */
	if(isave==0) free(cp);
	free(bcp);
	free(aa);
	free(a);
	free(r);
	free(cpt);
	free(caux);
	free(shift);
	free(al);
	free(ar);
	free(aux1);
	free(aux2);
	free(om);
	free(oms);
	free(vyx);
	free(vxy);
	free(v);
	free(va);
	free(asave);
	free(aasave);
	if(ncp==1 && idataxyw==0) eremove(diskxyw);

   	/* output migration result q(x,y,t) */
    hdrfp = fopen(diskhdr,"r");
	i64 = 0;
	fseek64(hdrfp,i64,0);

	fprintf(jpfp," Start Output \n");
	fprintf(jpfp," ============ \n");
	if ( nq==(ntau-itau0+1) && ntau>0 ) {
		if(itau0==1) {
			if(traceout==1) {
       			for (iy=0;iy<ny;iy++) {
       				for (ix=0;ix<nx;ix++) {
       					for (it=0;it<ntau-1;it++) 
						tra.data[it+1] = 
						q[(it*ny+iy)*nx+ix];
						tra.data[0] = 0;
	  					fread(&tra,sizeof(char),HDRBYTES,hdrfp);
						tra.ns = ntau;
						tra.dt = dtau * 1000000;
						tra.d1 = dtau * 1000000;
						if(tra.trid==0) {
							tra.trid = 2;
							if(cdpnum==0) {
                               tra.cdp = cdp1+(iy+iy0)*ncdppl*cdpinc
                                         +(ix+ix0)*cdpinc;
                        	} else {
                               tra.cdp = cdp1+(iy+iy0)*cdpinc
                                         +(ix+ix0)*nlines*cdpinc;
                        	}
						}
						if(dz>0.) {
							if(dz<=32.) {
								tra.dt = dz * 1000.;
							} else {
								tra.dt = dz;
							}
							tra.dz = dz;
							tra.fz = 0.;
							tra.mute = 0;
							tra.mutb = 0;
						}
						tra.tracl = ix + 1;
						tra.tracr = iy + 1; 
						tra.delrt = 0;
						if(ikey==1) {
							ftrace = ix*traceinc + trstart + 0.5;
							fline = iy*lineinc + lnstart + 0.5;
							itrace = ftrace;
							iline = fline;
							itov(trktype,&trkval,itrace);
							itov(lnktype,&lnkval,iline);
							puthval(&tra,indxtrk,&trkval);
							puthval(&tra,indxlnk,&lnkval);
						}
						fputtr(outfp,&tra);
					}
				}
			}			
			if(isave==1) {
				xytfp = fopen(diskxyt,"r+");
				fseek64(xytfp,0,0);
				fwrite(q,sizeof(float),nx*ny*ntau,xytfp);
				fclose(xytfp);
			}
		} else {
			xytfp = fopen(diskxyt,"r+");
			i64 = (itau0-1);
			i64 = i64*nx*ny;
			i64 = i64*sizeof(float);
			fseek64(xytfp,i64,0);
			fwrite(q,sizeof(float),nx*ny*nq,xytfp);
			fclose(xytfp);
		}
	} 
	if ( (itau0 > 1 || nq!=ntau) && ntau>0) {
		free(q);
		if(traceout==1) {
			if(llimit>128*1024*1024) {
				nyy = 128*1024*1024 / (nx*ntau*sizeof(float)) ;
			} else {
				nyy = llimit / (nx*ntau*sizeof(float)) ;
			}
			if(nyy<1) nyy=1;
			q = (float*)emalloc(nyy*nx*ntau*sizeof(float));
			if( q==0 ) err("llimit too small");
			xytfp = fopen(diskxyt,"r");
       		for (iy=0;iy<ny;iy=iy+nyy) {
				nyread = nyy;
				if(iy+nyread>ny) nyread = ny - iy;
				bzero((char*)q,nyy*nx*ntau*sizeof(float));

				for(it=0;it<ntau;it++) {
					i64 = it;
					i64 = i64*nx*ny+iy*nx;
					i64 = i64*sizeof(float);
					fseek64(xytfp,i64,0);
	   				fread(q+it*nx*nyy,sizeof(float),
						nx*nyread,xytfp);
				}
				for (iyread=0;iyread<nyread;iyread++) {
					iyy = iy + iyread;
					jy0 = iyread * nx;
       				for (ix=0;ix<nx;ix++) {
       					for (it=0;it<ntau-1;it++) 
						tra.data[it+1] = 
							q[it*nx*nyy+jy0+ix];
						tra.data[0] = 0.;
	  					fread(&tra,sizeof(char),HDRBYTES,hdrfp);
						tra.ns = ntau;
						tra.dt = dtau * 1000000;
						tra.d1 = dtau * 1000000;
						if(tra.trid==0) {
							tra.trid=2;
							if(cdpnum==0) {
                             tra.cdp = cdp1+(iy+iy0)*ncdppl*cdpinc
                                           +(ix+ix0)*cdpinc;
                        	} else {
                             tra.cdp = cdp1+(iy+iy0)*cdpinc
                                           +(ix+ix0)*nlines*cdpinc;
                        	}
						}
						if(dz>0.) {
							tra.dz = dz;
							tra.fz = 0;
							tra.mute = 0;
							tra.mutb = 0;
							if(dz<=32.) {
								tra.dt = dz * 1000.;
							} else {
								tra.dt = dz;
							}
						}
						tra.tracl = ix + 1;
						tra.tracr = iyy + 1; 
						tra.delrt = 0;
						if(ikey==1) {
							ftrace = ix*traceinc + trstart + 0.5;
							fline = iyy*lineinc + lnstart + 0.5;
							itrace = ftrace;
							iline = fline;
							itov(trktype,&trkval,itrace);
							itov(lnktype,&lnkval,iline);
							puthval(&tra,indxtrk,&trkval);
							puthval(&tra,indxlnk,&lnkval);
						}
						fputtr(outfp,&tra);
					}
				}
			}
			fclose(xytfp);
		}
		if(isave==0 ) {
			eremove(diskxyt);
		}
    }

	fprintf(jpfp," Job Done \n");

	free(q);
	fclose(outfp);
	fclose(hdrfp);
	if(isave==0) eremove(diskhdr);
	if(isave==1) {
		if(ncp==nf) {
			xywfp = fopen(diskxyw,"w");
			fseek64(xywfp,0,0);
			fwrite(cp,sizeof(complex),nx*ny*nf,xywfp);
			fclose(xywfp);
		}
	} else {
		if(ncp!=nf) eremove(diskxyw);
	} 

	if(ibacko==1) {
		fprintf(jpfp," Backup %s , %s and %s to %s \n", 
			diskhdr,diskxyt,diskxyw,backupo); 
		tar3to(backupo,diskhdr,diskxyt,diskxyw);
	}

	return 0;
}

/* write a input fft-transformed line to xyw-disk and write trace headers of 
a line to hdr-disk */

void wl2d(int nx,int ny,int nw,int iy,complex *cp,char *trhdr,float *fold,
		FILE *xywfp, FILE *hdrfp, FILE *jpfp,
		String trktype, String lnktype,
		Value trkval, Value lnkval,
		int indxtrk, int indxlnk, int trstart, int lnstart,
		int traceinc, int lineinc, int ikey, 
		complex *cpbig, int nybig, int iw2disk) {

	int ix,kx,ntmp,iw; 
	int itmp, mybig, jy;
	long long i64;
	int iline, itrace;
	float ftrace, fline;

	/* zero fill in missing traces */
	for(ix=0;ix<nx;ix++) {
		if(fold[ix]==0.) {
			fprintf(jpfp,
			"   Missing Trace At Shotpoint is=%d Line il=%d \n",
				ix+1,iy+1);
			fflush(jpfp);
			ntmp = ix - 1;
			if(ntmp<0) ntmp = 0;
			for(kx=ntmp;kx<nx;kx++) {
				if(fold[kx]!=0.) {
					for(iw=0;iw<nw;iw++) {
				         	cp[ix+iw*nx].r = 0.;
				         	cp[ix+iw*nx].i = 0.;
					}
		           		fold[ix]=1.0;
		      			bcopy(trhdr+kx*HDRBYTES,
                        			trhdr+ix*HDRBYTES,
						HDRBYTES);
					bzero(trhdr+28+ix*HDRBYTES,2);
		      			break;
	   			}
			}
		}
	}
	/* normalize */
	for(ix=0;ix<nx;ix++) {
		if(fold[ix]>1.0) {
			for(iw=0;iw<nw;iw++) {
		         	cp[ix+iw*nx].r = cp[ix+iw*nx].r/fold[ix];
		         	cp[ix+iw*nx].i = cp[ix+iw*nx].i/fold[ix];
			}
		}
	}

	jy = iy%nybig;
	mybig = jy + 1;

/*
	fprintf(stderr," iy=%d nybig=%d jy=%d \n",iy,nybig,jy);
*/

	/* save to cpbig */
	for(ix=0;ix<nx;ix++) {
		for(iw=0;iw<nw;iw++) {
			cpbig[iw*nybig*nx+jy*nx+ix].r = cp[ix+iw*nx].r;
			cpbig[iw*nybig*nx+jy*nx+ix].i = cp[ix+iw*nx].i;
		}
	}

	/* write to disk */
	if( iw2disk==1 || mybig==nybig ) {

		itmp = iy / nybig;
		itmp = itmp * nybig;

/*
		fprintf(stderr," iy=%d nybig=%d itmp=%d mybig=%d \n",
			iy,nybig,itmp,mybig);
*/

		for(iw=0;iw<nw;iw++) {
			i64 = iw;
			i64 = i64*nx*ny+itmp*nx;
			i64 = i64*sizeof(complex);
			fseek64(xywfp,i64,0);
	   		fwrite(cpbig+iw*nx*nybig,sizeof(complex),nx*mybig,xywfp);
		}
		bzero(cpbig,nybig*nw*nx*sizeof(complex));
	}

	/* save headers to disk */
	i64 = iy*nx*HDRBYTES*sizeof(char);
	fseek64(hdrfp,i64,0);
	if(ikey==1) {
		for(ix=0;ix<nx;ix++) {
			ftrace = ix*traceinc + trstart + 0.5;
			fline = iy*lineinc + lnstart + 0.5;
			itrace = ftrace;
			iline = fline;
			itov(trktype,&trkval,itrace);
			itov(lnktype,&lnkval,iline);
			bcopy(trhdr+ix*HDRBYTES,(char*)&tra,HDRBYTES);
			puthval(&tra,indxtrk,&trkval);
			puthval(&tra,indxlnk,&lnkval);
			bcopy((char*)&tra,trhdr+ix*HDRBYTES,HDRBYTES);
		}
	}
	fwrite(trhdr,sizeof(char),nx*HDRBYTES,hdrfp);
}