fxymig.c.6.15.99

seismic software,very useful

			fprintf(jpfp, 
			" Updates of disks diskxyw or diskxyt failed during last run \n");
			err(" Updates of disks diskxyw or diskxyt failed during last run");
			} 
		}
	} else {
		itime1 = 1;
		itime2 = 2;
		itime3 = 3;
	}
	if(ibackd==1) {

		lenxytb = strlen(diskxytb);
		lenxywb = strlen(diskxywb);
		namexytb = (char*) emalloc(lenxytb+1);
		namexywb = (char*) emalloc(lenxywb+1);
		bcopy(diskxytb,namexytb,lenxytb);
		bcopy(diskxywb,namexywb,lenxywb);
		namexytb[lenxytb]='\0';
		namexywb[lenxywb]='\0';

		hdrbfp = fopen(diskhdrb,"w");
		fseek64(hdrfp,0,0);
		fprintf(jpfp," Copy %s to %s \n",diskhdr,diskhdrb); 
		for(ix=0;ix<nx*ny;ix++) {
	  		fread(&tra,sizeof(char),HDRBYTES,hdrfp);
	  		fwrite(&tra,sizeof(char),HDRBYTES,hdrbfp);
		}
		efclose(hdrbfp);
		
		if(itime3>=itime1 && itime3>=itime2) {

			xywbfp = fopen(diskxywb,"w");
			fseek64(xywbfp,0,0);
			bzero(cp,nx*ny*sizeof(complex));
			i64 = 0;
			fseek64(xywfp,i64,0);
			fprintf(jpfp," Copy %s to %s \n",diskxyw,diskxywb); 
			for(iw=0;iw<nf;iw++) {
	  			efread(cp,sizeof(complex),nx*ny,xywfp);
	  			efwrite(cp,sizeof(complex),nx*ny,xywbfp);
			}
			efclose(xywbfp);

			xytbfp = fopen(diskxytb,"w");
			fseek64(xytbfp,0,0);
			if(idataxyt==1) {
				fprintf(jpfp," Copy %s to %s \n",diskxyt,diskxytb); 
				if(itau0>1) {
                    xytfp = efopen(diskxyt,"r");
                    fseek64(xytfp,0,0);
                    for(it=0;it<itau0-1;it++) {
                          efread(q,sizeof(float),nx*ny,xytfp);
                          efwrite(q,sizeof(float),nx*ny,xytbfp);
                    }
                    efclose(xytfp);
                }
				bzero(q,nx*ny*sizeof(float));
				for(it=itau0-1;it<ntau;it++) {
	  				efwrite(q,sizeof(float),nx*ny,xytbfp);
				}
		  	} else {
				bzero(q,nx*ny*sizeof(float));
				for(it=0;it<ntau;it++) {
	  				efwrite(q,sizeof(float),nx*ny,xytbfp);
				}
			}
			efclose(xytbfp);
		} else {
			xywbfp = fopen(diskxywb,"r");
			fseek64(xywbfp,0,0);
			bzero(cp,nx*ny*sizeof(complex));
			i64 = 0;
			fseek64(xywfp,i64,0);

			fprintf(jpfp,
			" Updates of disks diskxyw or diskxyt failed during last run \n");
			fprintf(jpfp," Copy %s to %s \n",diskxywb,diskxyw); 
			for(iw=0;iw<nf;iw++) {
	  			efread(cp,sizeof(complex),nx*ny,xywbfp);
	  			efwrite(cp,sizeof(complex),nx*ny,xywfp);
			}
			efclose(xywbfp);

			xytbfp = fopen(diskxytb,"r");
			fseek64(xytbfp,0,0);
            xytfp = efopen(diskxyt,"r+");
            fseek64(xytfp,0,0);
			fprintf(jpfp," Copy %s to %s \n",diskxytb,diskxyt); 
            for(it=0;it<ntau;it++) {
            	efread(q,sizeof(float),nx*ny,xytbfp);
                efwrite(q,sizeof(float),nx*ny,xytfp);
            }
            efclose(xytbfp);
			iwdur = 0;
			fprintf(jpfp, " Remigrate all frequency slices: set iwdur=0 \n");
		}
	}

	/* read the ffted wavefield into cp if needed */ 
	if(ncp==nf) {
		i64 = 0;
		fseek64(xywfp,i64,0);
		bzero(cp,nx*ny*nf*sizeof(complex));
	  	fread(cp,sizeof(complex),nx*ny*ncp,xywfp);
		if(isave==0) {
			fclose(xywfp);
			eremove(diskxyw);
		}
	} else {
		fclose(xywfp);
	}


	/* 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);	
	}

	/* 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);

	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) {
       			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(traceout==1) 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(mlimit>128*1024*1024) {
			nyy = 128*1024*1024 / (nx*ntau*sizeof(float)) ;
		} else {
			nyy = mlimit / (nx*ntau*sizeof(float)) ;
		}
		if(nyy<1) nyy=1;
		q = (float*)emalloc(nyy*nx*ntau*sizeof(float));
		if( q==0 ) err("mlimit 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(traceout==1) 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,
		complex *cpbig, int nybig, int iw2disk) {

	int ix,kx,ntmp,iw; 
	int itmp, mybig, jy;
	long long i64;

	/* 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);
	fwrite(trhdr,sizeof(char),nx*HDRBYTES,hdrfp);
}