ktmig.c~

seismic software,very useful

#include "usu.h"
#include "subc.h"
#include "ghdr.h"
#include "gridhd.h"
#include "grid.h"
#include "su.h"
#include "segy.h"
#include "header.h"

char *sdoc = 
"KTMIG - Kirchhoff Time Migration of 2D/3D poststack/prestack data 	\n"
"\n"
"ktmig [parameters] <input-data >migrated-data				\n" 
"\n"
"Required parameters:							\n"
"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"
"cdp1=                  cdp number 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"
"cdp2=                  cdp number 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"
"cdp3=                  cdp number of 3rd corner of the 3D master grid\n"
"                       (See Note 1 for details)			\n"
"\n"
"dds=                   cdp spacing in inline direction 		\n"
"ddl=                   Line spacing in crossline 			\n"
"velfile=               RMS velocity file name (velocity is stored as 	\n"
"                       v(nvs,nsv,nlv)					\n"
"\n"
"Optional parameters:							\n"
"ds=dds                 cdp group interval to output			\n"
"dl=ddl                 Line group  interval to output 			\n"
"datain=                input data set name (instead of standard input) \n"
"dataout=               output data set name (instead of standard output) \n"
"t0v=from-velfile       time (ms) of first sample in velfile  		\n" 
"ntv=from-velfile       Number of time samples per velocity trace	\n"
"dtv=from-velfile       time interval (ms) in velocity file		\n"
"s0v=from-velfile       smallest s coordinate in velfile  		\n" 
"nsv=from-velfile       Number of s coordinates in velocity file	\n"
"dsv=from-velfile       s interval of velocity file			\n"
"l0v=from-velfile       smallest l coordinate in velfile  		\n" 
"nlv=from-velfile       Number of l coordinates in velocity file	\n"
"dlv=from-velfile       l interval of velocity file			\n" 
"tau0=0.                minimum output time in ms			\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"
"oofo=0.                output minimum offset (>0.)			\n" 
"dofo=999999.           output offset increment				\n" 
"obtol=.5*dofo          offset bin tolerance (distance between the \n"
"                       accepted input offset and the output offset \n"
"                       (default to half of the output offset increment). \n"
"                       useful for migrating every other offset:         \n"
"                       e.g.: When dofo=400 and input offset increment \n"
"                       is 200, specify obtol=100 (instead of 200) \n"
"                       kzmig will migrated only those offsets  \n"
"                       whose distance from the output offset bin center \n"
"                       is equal or less than 100. \n"
"nofo=1                 number of output offsets			\n"
"ofo=0.                 output offsets  (specified as 200,250.,350,600,...) \n"
"                       if specified, dofo and oofo will be ignored	\n" 
"sstart=s1              Minimum s (trace position within line) of 	\n"
"                       migration output				\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"
"ns=(s2-s1)/ds+1        Number of cdp (group) per line of migration output \n"
"lstart=l1              Minimum l (line position) of migration output 	\n"
"                       (NOTE: this is coordinates, i.e., first line	\n"
"                       l1, 2nd at l1+dl, 3rd at l1+2*dl, ...)		\n"
"nl=(l3-l1)/dl+1        Number of lines (group) of migration output 	\n"
"send=                  Ending s position of migration output		\n"
"lend=                  Ending l position of migration output           \n" 
"ntrend=                Number of cdp bins to output from (sstart,lstart) \n"
"                       to (send,lend)					\n" 
"                       (specify send, lend and ntrend will output only \n"
"                       one slant line from (sstart,lstart) to 		\n"
"                       (send,lend) of total ntrend cdps; if not 	\n"
"                       specified, will output to ns by nl cdps)	\n" 
"nds=1                  Number of cdp locations per output cdp group	\n"
"ndl=1                  Number of line locations per output line group  \n"
"mlimit=256             Memory limit (in megabytes) to run the program  \n"
"mtrace=1000            Number of traces read per migration step	\n"
"apers=smax/4.          Migration aperature in inline (s) direction (the\n"
"                       maximum lateral distance from s position of 	\n"
"                       the input midpoint to s position of output midpoint),\n"
"                       where smax is the maximum of lateral distance 	\n"
"                       in inline direction				\n"
"aperl=lmax/4.          Migration aperature in crossline (l) direction (the\n"
"                       maximum lateral distance from l position of 	\n"
"                       the input midpoint to l position of output midpoint),\n"
"                       where lmax is the maximum of lateral distance 	\n"
"                       in crossline direction				\n"
"angmax=60.             maximum angle (from vertical) in degrees to output \n"
"                       migration from an input trace, i.e.,            \n"
"                          tmin = 2*dis/(v0*tan(angmax*3.141592654/180.)) \n"
"                       where tmin is the minimum time of an output trace, \n"
"                       dis is the distance between the output midpoint \n"
"                       and input midpoint, v0 is the surface velocity at \n"
"                       input midpoint position 			\n" 
"vmin=                  minimum velocity (default: read from velocity   \n"
"                       header file); must be specified if value 0 is \n"
"                       obtained from the velocity header. 		\n"
"tpow=-1.0              power of time gain to be applied to the input    \n"
"                       before migration 				\n"
"tpowaf=0.0             power of time gain to be applied to the output \n"
"                       after migration 				\n"
"isave=0                Mode of saving diskimg, diskhdr and diskfld	\n"
"                       0=remove the datasets after completion        	\n" 
"                       1=save the datasets after completion          	\n" 
"                       (will be sent to 1 when backupo is specified)	\n" 
"diskimg=DISKIMAGE      Disk name for storing imaged time sections 	\n"
"diskhdr=DISKHEADER     Disk name for storing output trace headers  	\n"
"diskfld=DISKFOLD       Disk name for storing output trace folds  	\n"
"backupi=               Name of backup tape name of the previous run	\n"
"                       (if given, e.g., stssxyz.ktmig.bck1, the backup \n"
"                        will be read in to overwrite the disk file of	\n"
"                        diskimg, diskhdr and diskfld; 			\n"
"                        otherwise, the program will use diskimg, 	\n"
"                        diskhdr and diskfld, if present)		\n" 
"backupo=               Name of backup tape name of the current run	\n"
"                       (if given, e.g., stssxyz.ktmig.bck2, the backup \n"
"                        tape will be made after the job is done;	\n"
"                        otherwise, only disk files will be created 	\n"
"                        if isave=1) 					\n"
"tapecntl=0             Input control (for tape or pipe  or disk input)	\n"
"                       0=tape or pipe input \n"
"                         The program will read and not process until 	\n"
"                         the input trace counter exceeds the number of \n"
"                         traces processed in the previous run of same input\n"
"                       1=more than one disk input \n"
"                         The program will skip the number of traces \n"
"                         processed in the previous run of the same input \n"
"                         before processing \n"
"                         This option should be used when multiple runs \n"
"                         (each with a single input) are needed	\n"
"                      -1=only one disk file input 	\n"
"                         The program will skip the total number of traces \n"
"                         processed in the previous run before processing \n"
"                       2=input from pipe, but processing will start from \n"
"                         first trace (like sudecomp input) \n"
"ntras=2000000000       Maximum number of traces to process in the 	\n"
"                       current run, before end of the input		\n"
"hisfile=KTMIG.HISFILE  histroy file name to indicate how many traces   \n"
"                       have been migrated in total (all runs) and in \n"
"                       current run	\n"
"                       (warning: remove this file before starting 	\n"
"                        migration if this is not a retart job, i.e.,   \n"
"                        the migration starts at the first input trace) \n" 
"                       If user does not specify hisfile, KTMIG.HISFILE \n"
"                       will be removed after normal completion of job. \n"
"f0=10.0                first high cut frequecy (Hz) in anti-aliasing 	\n"
"                       filters						\n" 
"df=10.0                high cut frequecy (Hz) increment in anti-aliasing\n"
"                       filters						\n" 
"nf=7                   number of anti-aliasing filters			\n"
"ftaper=5.0             frequency taper length (Hz) in the filters	\n"
"ksmax=0.5/dds          maximum wavenumber in line			\n"
"klmax=0.5/ddl          maximum wavenumber cross line			\n"  
"rmsmar=50.             maximum allowable ratio between sample amplitude\n"
"                       and the rms value of amplitudes of the whole 	\n"
"                       trace. When input amplitude/rms exceeds this 	\n"
"                       value, the trace will be treated as noise trace,\n"
"                       and not used in migration. A warning message 	\n" 
"                       will then be printed to a file (default to 	\n"
"                       BAD_TRACE_FILE)					\n"
"badtracefile=          name of file to record bad traces in the input	\n"
"                       (default to BAD_TRACE_FILE)			\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"
"ntrdsk=                number of traces migrated per disk file update	\n"
"                       (diskimg and diskfld will be updated once for 	\n"
"                        every migrated ntrdsk traces;			\n" 
"                        default to the end of the job 			\n"
"ncpu=1                 number of CPUs to be used		\n"
"traceout=1             output migrated seismic trace (1=yes 0=no) \n"
"                       (traceout=0 may be used to skip output, if \n"
"                        more input dataset is going to be migrated \n"
"                        in the next run) \n"
"on2trace=0             total number of traces in the input is \n"
"                        specified at the first 4 bytes of the trace \n"
"                        header of the first trace (0=no 1=yes) \n"
"tr1=                   trace number at 1st corner of the 3D master grid\n"
"tr2=                   trace number at 2nd corner of the 3D master grid\n"
"tr3=                   trace number at 3rd corner of the 3D master grid\n"
"ln1=                   line number at 1st corner of the 3D master grid\n"
"ln2=                   line number at 2nd corner of the 3D master grid\n"
"ln3=                   line number at 3rd corner of the 3D master grid\n"
"tracekey=              segy trace header key word for trace number \n"
"linekey=               segy trace header key word for line number \n"
"strace=0               when (trace#-strace)%jthtrace==0, pass to migration \n"
"                       otherwise skip \n"
"jthtrace=1             pass every jth-trace to migration \n"
"mintrace=-2000000000   minimum input trace number to migrate \n"
"maxtrace=2000000000    maximum input trace number to migrate \n"
"sline=0                when (line#-sline)%jthline==0, pass to migration \n"
"                       otherwise skip \n"
"jthline=1              pass every jth-line to migration \n"
"minline=-2000000000    minimum input line number to migrate \n"
"maxline=2000000000     maximum input line number to migrate \n"
"\n"
"Notes:									\n"
"0. 1st corner of the 3D master grid is at minimum line and minimum trace \n"
"   2nd corner of the 3D master grid is at minimum line and maximum trace \n"
"   3rd corner of the 3D master grid is at maximum line and minimum trace \n"
"   when tracekey, linekey, tr1, tr2, tr3, ln1, ln2 and ln3 are given, \n"
"   the output trace headers will be updated at missing input trace location \n"
"   (trid=2) with line and trace numbers \n"
"1. Input data must be SEGY (IEEE) data format				\n"
"2. Input RMS velocity file consists of ntv*nsv*nlv floating point samples \n"
"   with/without grid header						\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. source (sx,sy) and receiver (gx,gy) 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 ntras, diskimg, diskhdr, diskfld papameters	\n"
"   if the input dataset is stored as disk file,			\n"
"   for example,							\n"
"	stage 1: migrate 50000 traces 					\n"
"	   ktmig ... ntras=50000 isave=1 diskimg=/huge/id/tmpimg.data 	\n"
"                    diskhdr=/huge/id/tmphdr.data			\n" 
"                    diskfld=/huge/id/tmpfld.data			\n" 
"                    backupo=stssxyz.ktmig.backup1 hisfile=hf ...	\n"
"	stage 2: migrate next 50000 traces 				\n"
"	   ktmig ... ntras=50000 isave=1 diskimg=/huge/id/tmpimg.data 	\n"
"                    diskhdr=/huge/id/tmphdr.data			\n" 
"                    diskfld=/huge/id/tmpfld.data			\n" 
"                    backupi=stssxyz.ktmig.backup1 hisfile=hf		\n"
"                    backupo=stssxyz.ktmig.backup2 ...			\n"
"	stage 3: migrate  the final 10000 traces 			\n"
"	   ktmig ... ntras=10000 isave=0 diskimg=/huge/id/tmpimg.data 	\n"
"                    diskhdr=/huge/id/tmphdr.data			\n" 
"                    diskfld=/huge/id/tmpfld.data			\n" 
"                    backupi=stssxyz.ktmig.backup2 hisfile=hf		\n"
"\n"
"   If a job was completed normally at last stage, backupi should not	\n"
"   be specified unless diskimg, diskhdr and diskfld are bad, to save	\n"
"   the time of reading these three datasets into disk.			\n"
"               							\n"
"6. To migrate the 3-D dataset in stages (so that you get backups of    \n"
"   each stage), use diskimg, diskhdr, diskfld papameters and taperead	\n"
"   if the input dataset is stored as tapes,				\n"
"   for example,							\n"
"	stage 1: migrate first 5 tapes 					\n"
"          taperead sti=1 ntp=5 ... |					\n" 
"	   ktmig ... isave=1 diskimg=/huge/id/tmpimg.data 		\n"
"                    diskhdr=/huge/id/tmphdr.data			\n" 
"                    diskfld=/huge/id/tmpfld.data			\n" 
"                    backupo=stssxyz.ktmig.backup1 ...			\n"
"	stage 2: migrate next 5 tapes 					\n"
"          taperead sti=6 ntp=5 ... |					\n" 
"	   ktmig ... isave=1 diskimg=/huge/id/tmpimg.data 		\n"
"                    diskhdr=/huge/id/tmphdr.data			\n" 
"                    diskfld=/huge/id/tmpfld.data			\n" 
"                    backupi=stssxyz.ktmig.backup1 			\n"
"                    backupo=stssxyz.ktmig.backup2 ...			\n"
"	stage 3: migrate  the final 10 tapes 				\n"
"          taperead sti=11 ntp=10 ... |					\n" 
"	   ktmig ... isave=0 diskimg=/huge/id/tmpimg.data 		\n"
"                    diskhdr=/huge/id/tmphdr.data			\n" 
"                    diskfld=/huge/id/tmpfld.data			\n" 
"                    backupi=stssxyz.ktmig.backup2 			\n"
"\n"
"   If a job was completed normally at last stage, backupi should not	\n"
"   be specified unless diskimg, diskhdr and diskfld are bad, to save	\n"
"   the time of reading these three datasets into disk.			\n"
"   Do not specify ntras and hisfile parameters in this case.		\n"
"\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)*dds, 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)*ddl 	\n"
"   where 								\n"
"         ncdppl is the number of traces per line in the master grid	\n"
"         dds is the trace spacing (within a line) in the master grid	\n"
"         nlines is the number of lines in the master grid		\n"
"         ddl 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. When doing velocity analysis, the output locations are		\n"
"       the output inline locations (s):				\n"
" sstart      (ss1),      ss1+dds,      ...,       ss1+(nds-1)*dds,	\n"
" sstart+ds   (ss2),      ss2+dds,      ...,       ss2+(nds-1)*dds 	\n"
" sstart+2*ds (ss3),      ss3+dds,      ...,       ss3+(nds-1)*dds 	\n"
"    ...       								\n"
" sstart+(ns-1)*ds (ssn), ssn+dds,      ...,       ssn+(nds-1)*dds      \n"
"       the output crossline locations (l):				\n"
" lstart      (ll1),      ll1+ddl,      ...,       ll1+(ndl-1)*ddl,	\n"
" lstart+dl   (ll2),      ll2+ddl,      ...,       ll2+(ndl-1)*ddl 	\n"
" lstart+2*dl (ll3),      ll3+ddl,      ...,       ll3+(ndl-1)*ddl 	\n"
"    ...       								\n"
" lstart+(nl-1)*ds (lln), lln+ddl,      ...,       lln+(ndl-1)*ddl      \n"
" 9. The primary memory requirement is 					\n"
"            ns*nds*nlcore*nofo*ntau*4 (byte) + velocity grid size (byte)	\n" 
"            + ns*nds*nlcore*nfft*4 (byte) \n"
"    nlcore is set to nl*ndl.						\n"
" 10. when nlines is less than 10, a 2.5-D rho filter is applied to 	\n"
"     compensate for 2D data input					\n"
"\n"
"AUTHOR:		Zhiming Li,        	3/93   			\n"
;

void migs(int nl, int ns, int nofo, float *mig, int ntau,
	FILE *imgfp, int ktrace, int ntl, float dldt, float *trace, 
	float *tras, float dtl, float tminl, int nt, 
	float *ss, float *sl, float *gs, float *gl, 
	int *imute, float *fovt2, float *s, float *l,
	float *tm, int *iofs, float apers, float aperl, float *fold, 
	float *w1, float *work1, float *work2, float *wsave, float *tracef,
	int *ifcut, int *ltaper, float ksmax, float klmax, 
	float f0, float df, int nf, float ftaper, int nsave, int nfft,
	float dtau, float tau0, float angmax,
	int *indxw,int nindxw, int incdxw,
	float *resamp, int ires, int ntres,
	float *ss2, float *scs, float *scg, int ncpu, 
	float  *vr2, float *vi2, float *vq4, float *tau);

void fhigint_(float *f0,float *df,int *nf,float *ftaper,float *dt,
	int *ifcut,int *ltaper,float * wsave,int *nsave,int *nfft);