fxymig.c

seismic software,very useful

/*     features:
		up-to 90 degree accuracy
		absorbing boundaries
		dispersion suppression
		stable
		implicit f.d. scheme
		splitting method used 
		f.d. errors compensation
		parallel computation

*/

#include "usu.h"
#include "ghdr.h"
#include "gridhd.h"
#include "grid.h"
#include "comva.h"
#include "su.h"
#include "segy.h"
#include "header.h"
#include <sys/stat.h>

char *sdoc = 
"FXYMIG  - frequency-space domain 3-d migration			 	\n"
"\n"
"fxymig [parameters] <3d-stack >migrated-3d-data			\n" 
"\n"
"Required parameters:							\n"
"tracekey=              segy key word defining trace number within line \n"
"linekey=               segy key work defining line number 		\n"
"Or the following Four parameters defining 3D master grid:		\n"
"cdp1=                  first cdp number of 3D master grid              \n"
"cdpinc=                cdp number increment of 3D master grid          \n"
"ncdppl=                number of cdp (traces) per line of 3D master grid \n"
"nlines=                number of lines of 3D master grid		\n" 
"                       (See Note 6 for details)			\n"
"Or the following 12 3D master grid parameters:				\n" 
"x1=                    x coordinate of 1st corner of the 3D master grid\n"
"y1=                    y coordinate of 1st corner of the 3D master grid\n"
"s1=                    s coordinate of 1st corner of the 3D master grid\n"
"l1=                    l coordinate of 1st corner of the 3D master grid\n"
"x2=                    x coordinate of 2nd corner of the 3D master grid\n"
"y2=                    y coordinate of 2nd corner of the 3D master grid\n"
"s2=                    s coordinate of 2nd corner of the 3D master grid\n"
"l2=                    l coordinate of 2nd corner of the 3D master grid\n"
"x3=                    x coordinate of 3rd corner of the 3D master grid\n"
"y3=                    y coordinate of 3rd corner of the 3D master grid\n"
"s3=                    s coordinate of 3rd corner of the 3D master grid\n"
"l3=                    l coordinate of 3rd corner of the 3D master grid\n"
"                       (See Note 7 for details)			\n"
"                       If cdp1, cdpinc, ncdppl and nliness are given,  \n"
"                       s1=l1=0. 					\n"
"\n"
"ds=                    trace to trace (inline) interval to migrate	\n"
"dl=                    line to line interval to migrate (required when nl>1)\n"
"velfile=               INTERVAL velocity file name (velocity is stored as \n"
"                       v(ns,nl,nvs), i.e., nvs planes of (ns*nl))	\n"
"\n"
"Optional parameters:							\n"
"datain=                input dataset name (default to standard input) 	\n" 
"dataout=               output dataset name (default to standard output)\n" 
"cdpnum=0               cdp numbering type (0=inline then crossline)	\n" 
"                                           1=crossline then inline)	\n" 
"                       used only when cdp1, cdpinc, ncdppl and nlines  \n"
"                       are specified for 3d master grid.   		\n"
"The following two parameters define the migration grid origin:	\n"
"trstart=            	starting trace (inline) number to migrate	\n" 
"traceinc=1             trace number increment to migrate  \n"
"lnstart=            	starting line number to migrate	\n" 
"lineinc=1              line number increment to migrate 		\n"
"Or the following two parameters define the migration grid origin,	\n"
"                       when tracekey and linekey are not specified:	\n"
"sstart=s1              Minimum s (trace position within line) to migrate \n"
"                       (NOTE: this is coordinates, i.e., first trace	\n"
"                       of a line at s1, 2nd at s1+ds, 3rd at s1+2*ds, ...)\n"
"lstart=l1              Minimum l (line position) to migrate data	\n"
"                       (NOTE: this is coordinates, i.e., first line	\n"
"                       l1, 2nd at l1+dl, 3rd at l1+2*dl, ...)		\n"
"ns=ncdppl              Number of cdp (traces) per line to migrate	\n"
"                       (when ncdppl is specified, otherwise required) 	\n"
"nl=nlines              Number of lines in 3-d volume to migrate	\n"
"                       (when nlines is not specified, 1 is assumed)	\n"
"nvs=from-velfile       Number of depth (time) slices in velocity file	\n"
"dvs=from-velfile       Depth (time) slice interval in velocity file	\n"
"                       in METERS or FT if dz is NOT zero (DEPTH MIG)	\n"
"                       in MS if dz is zero (TIME MIG)			\n"
"ntau=nt                Number of output migrated samples per trace	\n" 
"dtau=dt                Output migrated time sample interval (in ms)\n"
"                       (dt is the input time sample interval from 	\n"
"                       trace header)					\n"
"dz=0.                  Output depth interval (in METERS or FT) 	\n"
"                       if not zero, each output trace will contain	\n"
"                       ntau samplles of DEPTH-migrated samples with    \n"
"                       sample interval of dz.	(DEPTH MIGRATION)	\n"
"                       if zero, each output trace will contain         \n"
"                       ntau samples of TIME-migrated samples with	\n"
"                       sample interval of dtau. (TIME MIGRATION)	\n"
"iestep=5               Extrapolation step in output samples			\n"
"iorder=1               Order of migration operator			\n"
"                       (0=45 1=65 2=80 3=87 4=89.99 5=90 degrees)	\n"
"nfft=trace_length*1.5  fft length of trace 				\n"
"dfc=0.03               Dip-filtering constant (0.<=dfc<=1)		\n"
"fmin=(0.5/dt)/20       Minimum frequency to migrate (hz)		\n"
"fmax=(0.5/dt)*1/2      Maximum frequency to migrate (hz) at z=0 or t=0	\n"
"fmaxend=fmax           Maximum frequency to migrate (hz) at tzend \n"
"tzend=99999            Maximum time (in ms) or maximum depth (in m or ft) \n"
"                       where fmaxend is specified \n"
"icstep=25              F.D. error compensation step in output samples (0=no)\n"
"nffs=ns*1.5            Fft length of inline direction when icstep>0	\n"
"nffl=nl*1.5            Fft length of crossline direction when icstep>0	\n"
"lvec=                  Vectorization length (computed internally)	\n"
"                       for 3D migration, must be the maximum of ns and nl \n"
"                       for 2D migration, must be less than or equal to	\n"
"                       nfft/2						\n"
"mlimit=512             Memory limit (in megabytes) to run the program  \n"
"isave=0                Mode of saving disk datasets of wavefield, 	\n"
"                       trace headers, and the imaged 3d volume         \n"
"                       0=remove these datasets after completion        \n" 
"                       1=save these datasets after completion          \n" 
"diskxyw=DISKWAVE       Disk name for storing wavefield			\n" 
"                       (frequency slices) during downward extrapolation \n"
"diskhdr=DISKHEADER     Disk name for storing input trace headers	\n"
"diskxyt=DISKIMAGE      Disk name for storing imaged time/depth slices  \n"
"ncpu=2                 number of CPUs to use (0<= ncpu <=32)		\n"
"durfile=               Disk update record file				\n"
"                       If specified, diskxyw and diskxyt will be updated \n"
"			             once nwmig frequencies have been downward \n"
"                        migrated nsteps. This allows restart of the job \n"
"                        from where the job was stopped by unexpected \n"
"                        system interference at the cost of more elapse time.\n"
"                        If not given (default), the disks will be backup \n"
"                        and reusable for the next step only after the \n"
"                        current job finishes normally. \n"
"nwmig=ncpu             number of frequency slices to migrate nsteps before \n"
"                        disk backup (should be integer multiple of ncpu ) \n"
"nsteps=                number of depth/time slices to hold in memory \n"
"                        recommended value 100, 200, ....    \n"
"                        (default to a number computed internally 	\n"
"                        according the memory limit)			\n" 
"jpfile=                Job print file name				\n"
"                          default: the standard error unit 		\n"
"                            (a) create prefix.exxxx file if using qsub \n"
"                                on convex, where prefix is the first 	\n"
"                                7 characters of the script file name 	\n"
"                                and xxxx is the job number;		\n"
"                            (b) send e-mail when using at command;   	\n" 
"                            (c) send printed message to screen when 	\n"
"                                running the job without qsub or at 	\n"
"                                commands.				\n" 
"                          specified: the message will be printed 	\n"
"                            in the named file as the job goes		\n" 
"backupi=               Name of backup tape name of the previous run	\n"
"                       (if given, e.g., stssxyz.fxymig.bck1, the backup \n"
"                        will be read in to overwrite the disk files	\n"
"                        diskxyw, diskhdr, and diskxyt; 		\n"
"                        otherwise, the program will use diskxyw, 	\n"
"                        diskhdr, diskxyt if they are present)		\n" 
"backupo=               Name of backup tape name of the current run	\n"
"                       (if given, e.g., stssxyz.fxymig.bck2, the backup \n"
"                        tape will be made after the job is done;	\n"
"                        otherwise, only diskxyw, diskhdr, and diskxyt 	\n"
"                        will be created if isave=1) 			\n"
"diskxywb=              Name of disk to backup diskxyw during migration \n" 
"diskhdrb=              Name of disk to backup diskhdr during migration \n"
"diskxytb=              Name of disk to backup diskxyt during migration \n"
"                       (Only when the above three disks are specified, \n"
"                       the program will backup three files every nsteps) \n"
"traceout=1             output migrated seismic trace (1=yes 0=no) \n"
"                       (traceout=0 may be used when running distributed \n"
"                       fxymig on multiple machines). Use fxymerge to   \n"
"                       output migrated traces from diskfiles of the \n"
"                       multiple machines  \n"
"ifmin=                 minimum frequency index to migrate (calulated from \n"
"                       the fxymcal program) for this job \n"
"                       (frequency 0Hz is at index 0, nyquist at nfft/2) \n"
"ifmax=                 maximum frequency index to migrate (calulated from \n"
"                       the fxymcal program) for this job \n"
"                       if both ifmin and ifmax is specified, only the \n"
"                       frequency slices specified between ifmin and ifmax \n"
"                       will be migrated in the job \n"
"maxamp=1.0e+20         maximum absolute amplitude allowed in input \n"
"                       when an amplitude of input trace is larger than \n"
"                       maxamp, the job will be canceled \n"
"minvel=500             minimum velocity allowed in input velocity \n"
"maxvel=50000           maximum velocity allowed in input velocity \n"
"                       when the input velocity is outside (minvel,maxvel) \n"
"                       the job will be canceled \n"
"                       if minvel or maxvel is specified as negative number \n"
"                       no check will be performed \n"
"\n"
"Notes:									\n"
"1. Input stack must be ASCII-IEEE SEGY data format				\n"
"2. Input interval velocity file must be transposed:			\n"
"   i.e., stored as velocity time/depth slices for time/depth migration \n"
"3. s is the lateral position along the inline direction, l is the 	\n"
"   lateral position along the crossline direction. They may be different\n"
"   from the (x,y) position of the 3D survey. 				\n"
"4. When cdp1, cdpinc, ncdppl, and nlines are specified, the cell position\n"
"   of the inpute trace is determined by its cdp number. Otherwise, the	\n"
"   source (x,y) and receiver (x,y) coordinates are used to determine	\n"
"   the midpoint x and y position on the 3D grid:			\n"
"		x = (gx+sx)/2 * scalco					\n" 
"		y = (gy+sy)/2 * scalco		when scalco>1		\n" 
"    OR:								\n"	
"		x = (gx+sx)/2 / (-scalco)	when scalco<0		\n" 
"		y = (gy+sy)/2 / (-scalco)				\n" 
"   where gx is the receiver x coordinate from the trace header 	\n"
"    	  sx is the source x coordinate from the trace header 		\n"
"         gy is the receiver y coordinate from the trace header 	\n"
"    	  sy is the source y coordinate from the trace header 		\n"
"         scalco is the scaler to be applied to gx,gy,sx,sy		\n" 
"5. To migrate the 3-D dataset in stages (so that you get backups of    \n"
"   each stage), use ntau, isave, diskxyw, diskhdr, diskxyt papameters	\n"
"   for example,							\n"
"	stage 1: extrapolate to 500-th time/depth sample		\n"
"	   fxymig ... ntau=500 isave=1 diskxyw=/huge/id/tmpxyw.data 	\n"
"	              diskhdr=/huge/id/tmphdr.data 			\n"
"                     diskxyt=/huge/id/tmpimg.data 			\n"
"                     backupo=stssxyz.fxymig.backup1 ...		\n"
"	stage 2: extrapolate to the 1000-th time/depth sample		\n"
"	   fxymig ... ntau=1000 isave=1 diskxyw=/huge/id/tmpxyw.data 	\n"
"	              diskhdr=/huge/id/tmphdr.data 			\n"
"                     diskxyt=/huge/id/tmpimg.data 			\n"
"                     backupi=stssxyz.fxymig.backup1 ...		\n"
"                     backupo=stssxyz.fxymig.backup2 ...		\n"
"	stage 3: extrapolate to the final (1500-th) time/depth sample	\n"
"	   fxymig ... ntau=1500 isave=0 diskxyw=/huge/id/tmpxyw.data 	\n"
"	              diskhdr=/huge/id/tmphdr.data 			\n"
"                     backupi=stssxyz.fxymig.backup2 ...		\n"
"6. cdp1, cdpinc, ncdppl, nlines are used to define master 3D grid:	\n"
"               							\n"
"  cdp1+(nlines-1)*ncdppl*cdpinc  cdp1+(nlines*ncdppl-1)*cdpinc		\n"
"                    *------------------*    <----------line nlines	\n"
"                    |                  |				\n"
"                    |                  |				\n"
"                    |                  |				\n"
"                    |                  |				\n"
"cdp1+2*ncdppl*cdpinc|                  |				\n"
" cdp1+ncdppl*cdpinc |                  |    <----------line 2		\n"
"                    *------------------*    <----------line 1 		\n"
"                  cdp1            cdp1+(ncdppl-1)*cdpinc 		\n"
"                                                                       \n"
" within the master 3D grid, cdp number is incremented with cdpinc, i.e.,\n"
" first line starts with cdp1, second line starts with cdp1+ncdppl*cdpinc,\n"
" etc., the last line starts with cdp1+(nlines-1)*ncdppl*cdpinc.	\n"
" Input data may have some missing traces at some cdps, but the cdp 	\n"
" numbers of those input traces must match the cdp number definition 	\n"
" shown above. The s1 and l1 values at the starting cell (cdp1) position\n"
" will then be zero.							\n"
"7. s and l positions of an input trace are computed using the three 	\n"
"   master-grid corner positions					\n" 
"									\n"
"        | y								\n"
"        |   .l        * (x4,y4)					\n"
"        |    .     .    .						\n"
"        |     .  .       .         . s					\n"
"       (x3,y3) *          .      .					\n"
"        |        .          .  .					\n"
"        |         .          * (x2,y2)					\n"
"        |          .       .						\n"
"        |            .  .						\n"
"        |             * (x1,y1)					\n"
"        |								\n"
"        |--------------------------------- x				\n"
"									\n"
"   (x1,y1) has the smallest s value and the samllest l value		\n"
"              (s1 is usually =0.0, l1 is usually =0.0)			\n"
"   (x2,y2) has the largest s value and the smallest l value		\n"
"              (s2 is usually =(ncdppl-1)*ds, l2 is usually =0.0)	\n"
"   (x3,y3) has the smallest s value and the largest l value		\n"
"              (s3 is usually =0.0, l3 is usually =(nlines-1)*dl 	\n"
"   where 								\n"
"         ncdppl is the number of traces per line in the master grid	\n"
"         ds is the trace spacing (within a line) in the master grid	\n"
"         nlines is the number of lines in the master grid		\n"
"         dl is the line spacing in the master grid 			\n"
"	  s is the coordinate (in m or ft --- NOT just an integer number)\n"
"              of trace position within a 3d line 			\n"
"         l is the coordinate (in m or ft --- NOT just an integer number)\n"
"              of line position within a 3d survey 			\n"
"8. There are ns by nl traces output, regardless of input traces.	\n"
"   Trace number within line (starting from 1 with increment of 1) will \n"
"   be updated in trace header word 'tracl', while line number (starting\n"
"   from 1 with increment of 1) will be updated in trace header word	\n" 
"   'tracr'. Input missing trace position will be filled with zero traces \n"
"   before migration, and the resulting migrated traces at these locations\n" 
"   will be outputed with header word 'trid' equal 2 (dead trace). Other\n"
"   output header values at these missing trace locations may not be 	\n"
"   correct, since they are simply copied from the nearest live input 	\n"
"   traces.								\n" 
"9. Memory requirement can be estimated by, 				\n"
"  ns*nl*nwmig*8+ns*nl*7*ncpu*8+ns*nl*(nsteps+4+nsteps*dtau/dvs)*4 (BYTES)\n"
"\n"
 
"AUTHOR:		Zhiming Li,        	6/92   			\n"
;

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

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

segytrace tra;
segybhdr bh;
segychdr ch;

main(int argc, char **argv)
{
	int cdp1, nlines, ncdppl, cdpinc, icdps;
   	int nt,nx,ny,ntau,nw,nfft,nfftq,it,iw,ix,iy,iorder,itau0,nsteps;
    	int iestep,icstep,nffx,nffy,naux1,naux2,naux;
    	int nvs,nv,nq,ncp,lvec,lplane,isave,iq;
    	int nkx,nky,ntmp,estep,cstep,ntrace;
	int idataxyt,idataxyw,idatahdr,jx,jy,ir,ii;
	int iwmin, iwmax;
   	float dt,dx,dy,dtau,dvs,dz,dfc,fmin,fmax,x0m,y0m,vref; 
	int ifmin, ifmax, traceout,ifrange,nf, ifminr;
	float fmaxend, tzend;
	int iwend, ntauend;
	float pi,dw,wmin,wmax,xm,ym,tmp;
	float *q,*v,*vxy,*vyx,*om,*oms,*aux1,*aux2,*va;
	float *fold,*wsave,*trfft;
   	complex *cp,*aa,*a,*r,*bcp,*caux,*shift,*al,*ar,*cpp;
	complex *cpbig;
	float *cpt;
   	char *diskxyt,*diskxyw,*diskhdr,*velfile;
   	char *diskxytb,*diskxywb,*diskhdrb;
   	char *datain,*dataout,*trhdr;
   	FILE *infp,*outfp,*velfp,*xywfp,*xytfp,*hdrfp;
   	FILE *xywbfp,*xytbfp,*hdrbfp;
	char *namevel, *namexyt, *namexyw;
	char *namexytb, *namexywb;
	int lenxywb=0, lenxytb=0;
	int lenvel, lenxyw, lenxyt, dummyunit=77;
	float x1,x2,x3,y1,y2,y3,s1,s2,s3,l1,l2,l3,sm,lm;
	complex *asave, *aasave; 
	int iisave=0;
	char *backupi, *backupo;
	int ibacki, ibacko, ibackd;
	int cdpnum;
	char *jpfile;
	char *durfile;
	char *namejp;
	char *namedur;
	int lenjpf;
	int lendur;
	int itau0dur, itaudur, iwdur;
	char tau0dur[10], taudur[9], wdur[7];
	char dump1[3], dump2[3], dump3[2], dump4[8], dump5[4];
	FILE *jpfp;
	FILE *durfp;
	int ncpu, nwmig;
	float *w;
	complex *cpbuf;
	float *qbuf, *dzov;
	complex *cph;
	int ncph=32768;
	float temp;
	int ntaulast;
	float maxamp, minvel, maxvel;
	int itau;

	String tracekey="tracl", linekey="tracr", trktype, lnktype;
	Value trkval, lnkval;
	int indxtrk, indxlnk, trstart, lnstart;
	int traceinc, lineinc, ikey=0, iline, itrace;
	float ftrace, fline;
	
	ghed gh;