kzmig.c

seismic software,very useful

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

char *sdoc[] = {
"KZMIG - Kirchhoff depth Migration of 3D poststack/prestack data 	",
" ",
"kzmig [parameters] <input-data >migrated-data				", 
" ",
"Required parameters:							",
"3D master grid parameters:						", 
"x1=                    x coordinate of 1st corner of the 3D master grid",
"y1=                    y coordinate of 1st corner of the 3D master grid",
"s1=                    s coordinate of 1st corner of the 3D master grid",
"l1=                    l coordinate of 1st corner of the 3D master grid",
"cdp1=                  cdp number of 1st corner of the 3D master grid  ",
"x2=                    x coordinate of 2nd corner of the 3D master grid",
"y2=                    y coordinate of 2nd corner of the 3D master grid",
"s2=                    s coordinate of 2nd corner of the 3D master grid",
"l2=                    l coordinate of 2nd corner of the 3D master grid",
"cdp2=                  cdp number of 2nd corner of the 3D master grid  ",
"x3=                    x coordinate of 3rd corner of the 3D master grid",
"y3=                    y coordinate of 3rd corner of the 3D master grid",
"s3=                    s coordinate of 3rd corner of the 3D master grid",
"l3=                    l coordinate of 3rd corner of the 3D master grid",
"cdp3=                  cdp number of 3rd corner of the 3D master grid  ",
"                       (See Note 1 for details)			",
" ",
"dds=                   cdp spacing in inline direction 		",
"ddl=                   Line spacing in crossline 			",
 " ",
"Optional parameters:							",
"ds=dds                 cdp group interval to output			",
"dl=ddl                 Line group  interval to output 			",
"datain=                input data set name (instead of standard input) ",
"dataout=               output data set name (instead of standard output)",
"ttfile=ttfile          file name of input traveltime tables		",
"fzt=from-ttfile        depth of first sample in traveltime tables 	" ,
"nzt=from-ttfile       	Number of depth samples in traveltime tables 	",
"dzt=from-ttfile       	depth interval in traveltime tables		",
"fst=from-ttfile       	first s coordinate in traveltime tables 	", 
"nst=from-ttfile       	Number of s coordinates in traveltime tables	",
"dst=from-ttfile        s interval in traveltime tables			",
"flt=from-ttfile       	first l coordinate in traveltime tables 	", 
"nlt=from-ttfile       	Number of l coordinates in traveltime tables	",
"dlt=from-ttfile        l interval in traveltime tables			",
"fxst=from-ttfile       first s coordinate of source locations in 	",	
"				traveltime tables 	", 
"nxst=from-ttfile       Number of s coordinates of source locations in 	",	
"				traveltime tables 	", 
"dxst=from-ttfile       s interval of source locations in traveltime 	",
"				tables			",
"fyst=from-ttfile       first l coordinate of source locations in 	",	
"				traveltime tables 	", 
"nyst=from-ttfile       Number of l coordinates of source locations in 	",	
"				traveltime tables 	", 
"dyst=from-ttfile       l interval of source locations in traveltime 	",
"				tables			",
"THE ABOVE 15 PARAMETERS ARE READ IN FROM TTFILE ONLY			",
"fzo=fzt                first output depth for migration 		",
"nzo=10*nzt             Number of output migrated samples per trace	", 
"dzo=0.1*dzo            Output migrated depth sample interval  		",
 "                       						",
"flexbin=0              =1 Input data have been flexbinned to normalize	",
"                        amplitudes before migration, the migration will",
"                        not scale the summation by the migration fold	",
"                       =0 Otherwise, use migration fold to normalize 	",
"                        output 					",
"oofo=0.                output minimum offset (>0.)			", 
"dofo=999999.           output offset increment				", 
"nofo=1                 number of output offsets			",
"ofo=0.                 output offsets  (specified as 200,250.,350,600,...) ",
"                       if specified, dofo and oofo will be ignored	", 
"offnear=0              1=any offset (absolute) outside first offset bin ",
"                         is used in migration and summed to min offset", 
"                       0=otherwise will be rejected 			", 
"offfar=0               1=any offset (absolute) outside last offset bin ",
"                         is used in migration and summed to max offset", 
"                       0=otherwise will be rejected 			", 
"sstart=s1              Minimum s (trace position within line) of 	",
"                       migration output				",
"                       (NOTE: this is coordinates, i.e., first trace	",
"                       of a line at s1, 2nd at s1+ds, 3rd at s1+2*ds, ...)",
"ns=(s2-s1)/ds+1        Number of cdp (group) per line of migration output ",
"lstart=l1              Minimum l (line position) of migration output 	",
"                       (NOTE: this is coordinates, i.e., first line	",
"                       l1, 2nd at l1+dl, 3rd at l1+2*dl, ...)		",
"nl=(l3-l1)/dl+1        Number of lines (group) of migration output 	",
"send=                  Ending s position of migration output		",
"lend=                  Ending l position of migration output           ", 
"ntrend=                Number of cdp bins to output from (sstart,lstart) ",
"                       to (send,lend)					", 
"                       (specify send, lend and ntrend will output only ",
"                       one slant line from (sstart,lstart) to 		",
"                       (send,lend) of total ntrend cdps; if not 	",
"                       specified, will output to ns by nl cdps)	", 
"nds=1                  Number of cdp locations per output cdp group	",
"ndl=1                  Number of line locations per output line group  ",
"mlimit=256             Memory limit (in megabytes) to run the program  ",
"gath=0			type of gather: 1--common shot; 2--commom 	",
"				receiver; 0 --other types 		",
"mtrace=1000            maximum number of traces per migration step 	",
"offmax=13000           upper limit of offset (default in ft)		", 
"apers=smax/4.          Migration aperature in inline (s) direction (the",
"                       maximum lateral distance from s position of 	",
"                       the input midpoint to s position of output midpoint),",
"                       where smax is the maximum of lateral distance 	",
"                       in inline direction				",
"aperl=lmax/4.          Migration aperature in crossline (l) direction (the",
"                       maximum lateral distance from l position of 	",
"                       the input midpoint to l position of output midpoint),",
"                       where lmax is the maximum of lateral distance 	",
"                       in crossline direction				",
"angmax=60.             maximum angle (from vertical) in degrees to output ",
"                       migration from an input trace, i.e.,            ",
"                          tmin = 2*dis/(v0*tan(angmax*3.141592654/180.)) ",
"                       where tmin is the minimum time of an output trace, ",
"                       dis is the distance between the output midpoint ",
"                       and input midpoint, v0 is the surface velocity at ",
"                       input midpoint position 			", 
"v0=5000(ft/s)         	reference velocity at the surface.   		",	
"dvz=0.0                reference velocity gradient in depth direction. ",
"tpow=-1.0              power of time gain to be appled to the input    ",
"                       before migration   ",
"zpow=0.0               power of depth gain to be appled to the output  ",
"                       after migration   ",
"isave=0                Mode of saving diskimg, diskhdr and diskfld	",
"                       0=remove the datasets after completion        	",
"                       1=save the datasets after completion          	", 
"                       (will be sent to 1 when backupo is specified)	",
"diskimg=DISKIMAGE      Disk name for storing imaged time sections 	",
"diskhdr=DISKHEADER     Disk name for storing output trace headers  	",
"diskfld=DISKFOLD       Disk name for storing output trace folds  	",
"backupi=               Name of backup tape name of the previous run	",
"                       (if given, e.g., stssxyz.kzmig.bck1, the backup ",
"                        will be read in to overwrite the disk file of	",
"                        diskimg, diskhdr and diskfld; 			",
"                        otherwise, the program will use diskimg, 	",
"                        diskhdr and diskfld, if present)		",
"backupo=               Name of backup tape name of the current run	",
"                       (if given, e.g., stssxyz.kzmig.bck2, the backup ",
"                        tape will be made after the job is done;	",
"                        otherwise, only disk files will be created 	",
"                        if isave=1) 					",
"tapecntl=0             Input tape control (for tape or pipe input)	",
"                       0=tape mounted from begining VSN tpae of a DSN	",
"                         The program will read and not process until 	",
"                         the input trace counter exceeds the number of ",
"                         traces process in the previous run.		",
"                         This option should be used when input is PIPED",
"                         to kzmig					",
"                       1=tape mounted from 1st tape to be processed in	",
"                         the current run. The program will set the 	",
"                         input trace counter to the number of traces 	",
"                         processed in the last run and start to process",
"                         traces right away. Used when input are more   ",
"                         than one disk files (no pipe).		",
"                      -1=input dataset is from disk (not tape and not	",
"                         pipe). The program will them perform a seek	",
"                         to the proper trace position, before reading	",
"                         and processing the data. Use only when a 	",
"                         single disk file is going to be used for 	",
"                         migration.            			",
"ntras=2000000000       Maximum number of traces to process in the 	",
"                       current run, before end of the input		",
"hisfile=KZMIG.HISFILE  histroy file name to indicate how many traces   ",
"                       has been migrated in total upto the last run	",
"                       (warning: remove this file before starting 	",
"                        migration if this is not a retart job, i.e.,   ",
"                        the migration starts at the first input trace) ",
"                       If user does not specify hisfile, KZMIG.HISFILE ",
"                       will be removed after normal completion of job. ",
"f0=10.0                first high cut frequecy (Hz) in anti-aliasing 	",
"                       filters						",
"df=10.0                high cut frequecy (Hz) increment in anti-aliasing",
"                       filters						",
"nf=6                   number of anti-aliasing filters			",
"ftaper=5.0             frequency taper length (Hz) in the filters	",
"ksmax=0.5/dds          maximum wavenumber in line			",
"klmax=0.5/ddl          maximum wavenumber cross line			",  
"rmsmar=50.             maximum allowable ratio between sample amplitude",
"                       and the rms value of amplitudes of the whole 	",
"                       trace. When input amplitude/rms exceeds this 	",
"                       value, the trace will be treated as noise trace,",
"                       and not used in migration. A warning message 	", 
"                       will then be printed to a file (default to 	",
"                       BAD_TRACE_FILE)					",
"badtracefile=          name of file to record bad traces in the input	",
"                       (default to BAD_TRACE_FILE)			",
"jpfile=                Job print file name                             ",
"                          default: the standard error unit             ",
"                            (a) create prefix.exxxx file if using qsub ",
"                                on convex, where prefix is the first   ",
"                                7 characters of the script file name   ",
"                                and xxxx is the job number;            ",
"                            (b) send e-mail when using at command;     ",
"                            (c) send printed message to screen when    ",
"                                running the job without qsub or at     ",
"                                commands.                              ",
"                          specified: the message will be printed       ",
"                            in the named file as the job goes          ",
"ntrdsk=                number of traces migrated per disk file update  ",
"                       (diskimg and diskfld will be updated once for   ",
"                        every migrated ntrdsk traces;                  ",
"                        default to the end of the job                  ",
"ncpu=1                 number of cpu to use in the run 		",
"                       (used in kzmigp program, ignored in kzmig)	",
" ",
"Notes:		",
"I/O managements are mostly copied from ktmig.c by Zhiming Li.	",	
"1. Input data must be SEGY (IEEE) data format				",
"2. Input ttfile consists of nzt*nst*nlt*nxst*nyst floating point samples ",
"   with 100-byte of grid header					", 
"3. s is the lateral position along the inline direction, l is the 	",
"   lateral position along the crossline direction. They may be different",
"   from the (x,y) position of the 3D survey. 				",
"4. source (sx,sy) and receiver (gx,gy) coordinates are used to determine",
"   the midpoint x and y position on the 3D grid:			",
"		x = (gx+sx)/2 * scalco					",
"		y = (gy+sy)/2 * scalco		when scalco>1		", 
"    OR:								",
"		x = (gx+sx)/2 / (-scalco)	when scalco<0		",
"		y = (gy+sy)/2 / (-scalco)				",
"   where gx is the receiver x coordinate from the trace header 	",
"    	  sx is the source x coordinate from the trace header 		",
"         gy is the receiver y coordinate from the trace header 	",
"    	  sy is the source y coordinate from the trace header 		",
"         scalco is the scaler to be applied to gx,gy,sx,sy		",
"5. To migrate the 3-D dataset in stages (so that you get backups of    ",
"   each stage), use ntras, diskimg, diskhdr, diskfld papameters	",
"   if the input dataset is stored as disk file,			",
"   for example,							",
"	stage 1: migrate 50000 traces 					",
"	   kzmig ... ntras=50000 isave=1 diskimg=/huge/id/tmpimg.data 	",
"                    diskhdr=/huge/id/tmphdr.data			", 
"                    diskfld=/huge/id/tmpfld.data			", 
"                    backupo=stssxyz.kzmig.backup1 hisfile=hf ...	",
"	stage 2: migrate next 50000 traces 				",
"	   kzmig ... ntras=50000 isave=1 diskimg=/huge/id/tmpimg.data 	",
"                    diskhdr=/huge/id/tmphdr.data			",
"                    diskfld=/huge/id/tmpfld.data			", 
"                    backupi=stssxyz.kzmig.backup1 hisfile=hf		",
"                    backupo=stssxyz.kzmig.backup2 ...			",
"	stage 3: migrate  the final 10000 traces 			",
"	   kzmig ... ntras=10000 isave=0 diskimg=/huge/id/tmpimg.data 	",
"                    diskhdr=/huge/id/tmphdr.data			",
"                    diskfld=/huge/id/tmpfld.data			", 
"                    backupi=stssxyz.kzmig.backup2 hisfile=hf		",
" ",
"   If a job was completed normally at last stage, backupi should not	",
"   be specified unless diskimg, diskhdr and diskfld are bad, to save	",
"   the time of reading these three datasets into disk.			",
"               							",
"6. To migrate the 3-D dataset in stages (so that you get backups of    ",
"   each stage), use diskimg, diskhdr, diskfld papameters and taperead	",
"   if the input dataset is stored as tapes,				",
"   for example,							",
"	stage 1: migrate first 5 tapes 					",
"          taperead sti=1 ntp=5 ... |					", 
"	   kzmig ... isave=1 diskimg=/huge/id/tmpimg.data 		",
"                    diskhdr=/huge/id/tmphdr.data			", 
"                    diskfld=/huge/id/tmpfld.data			", 
"                    backupo=stssxyz.kzmig.backup1 ...			",
"	stage 2: migrate next 5 tapes 					",
"          taperead sti=6 ntp=5 ... |					", 
"	   kzmig ... isave=1 diskimg=/huge/id/tmpimg.data 		",
"                    diskhdr=/huge/id/tmphdr.data			", 
"                    diskfld=/huge/id/tmpfld.data			", 
"                    backupi=stssxyz.kzmig.backup1 			",
"                    backupo=stssxyz.kzmig.backup2 ...			",
"	stage 3: migrate  the final 10 tapes 				",
"          taperead sti=11 ntp=10 ... |					",
"	   kzmig ... isave=0 diskimg=/huge/id/tmpimg.data 		",
"                    diskhdr=/huge/id/tmphdr.data			", 
"                    diskfld=/huge/id/tmpfld.data			",
"                    backupi=stssxyz.kzmig.backup2 			",
" ",
"   If a job was completed normally at last stage, backupi should not	",
"   be specified unless diskimg, diskhdr and diskfld are bad, to save	",
"   the time of reading these three datasets into disk.			",
"   Do not specify ntras and hisfile parameters in this case.		",
" ",
"7. s and l positions of an input trace are computed using the three 	",
"   master-grid corner positions					",
"									",
"        | y								",
"        |   .l        * (x4,y4)					",
"        |    .     .    .						",
"        |     .  .       .         . s					",
"       (x3,y3) *          .      .					",
"        |        .          .  .					",
"        |         .          * (x2,y2)					",
"        |          .       .						",
"        |            .  .						",
"        |             * (x1,y1)					",
"        |								",
"        |--------------------------------- x				",
"									",
"   (x1,y1) has the smallest s value and the samllest l value		",
"              (s1 is usually =0.0, l1 is usually =0.0)			",
"   (x2,y2) has the largest s value and the smallest l value		",
"              (s2 is usually =(ncdppl-1)*dds, l2 is usually =0.0)	",
"   (x3,y3) has the smallest s value and the largest l value		",
"              (s3 is usually =0.0, l3 is usually =(nlines-1)*ddl 	",
"   where 								",
"         ncdppl is the number of traces per line in the master grid	",
"         dds is the trace spacing (within a line) in the master grid	",
"         nlines is the number of lines in the master grid		",
"         ddl is the line spacing in the master grid 			",
"	  s is the coordinate (in m or ft --- NOT just an integer number)",
"              of trace position within a 3d line 			",
"         l is the coordinate (in m or ft --- NOT just an integer number)",
"              of line position within a 3d survey 			",
" 8. When doing velocity analysis, the output locations are		",
"       the output inline locations (s):				",
" sstart      (ss1),      ss1+dds,      ...,       ss1+(nds-1)*dds,	",
" sstart+ds   (ss2),      ss2+dds,      ...,       ss2+(nds-1)*dds 	",
" sstart+2*ds (ss3),      ss3+dds,      ...,       ss3+(nds-1)*dds 	",
"    ...       								",
" sstart+(ns-1)*ds (ssn), ssn+dds,      ...,       ssn+(nds-1)*dds      ",
"       the output crossline locations (l):				",
" lstart      (ll1),      ll1+ddl,      ...,       ll1+(ndl-1)*ddl,	",
" lstart+dl   (ll2),      ll2+ddl,      ...,       ll2+(ndl-1)*ddl 	",
" lstart+2*dl (ll3),      ll3+ddl,      ...,       ll3+(ndl-1)*ddl 	",
"    ...       								",
" lstart+(nl-1)*ds (lln), lln+ddl,      ...,       lln+(ndl-1)*ddl      ",
" 9. The primary memory requirement is 					",
"            ns*nds*nlcore*nofo*nzo*4 (byte)				", 
"    nlcore is set to nl*ndl.						",
" 10. when nlines is less than 10, a 2.5-D rho filter is applied to 	",
"     compensate for 2D data input					",
" ",
"AUTHOR:		Z. Liu and Z. Li        9/94   			",
NULL};

void migs(int nl, int ns, float **mig, int nz,
	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 *imutel, float *s, float *l,
	int *iofs, float apers, float aperl, float **fold,
	float *work1, float *work2, float *wsave, float *tracef, float *sqrtf,
	int tahd, int *ifcut, int *ltaper, float ksmax, float klmax,
	float f0, float df, int nf, float ftaper, int nsave, int nfft,
	float dz, float fz, float angmax,
	int nxs, int nys, float fxs, float fys, float dxs, float dys,
	int nzo, float fzo, float dzo,
	int nr, float dr, float *tb, float *pb, float **ttab, int igh); 

void pszm3d_(float *trace,int *nt,float *t0,float *dt,float *sx,float *sy,
       float *gx,float *gy,int *imute,float *mig,float *x,float *y,
       int *nxy,int *nzo,float *dzo,float *fzo,float *aperx,float *apery,
       float *fold,float *f0,float *df,int *nf,float *ftaper,
       int *ifcut,int *ltaper,float *tracef,float *wsave,int *nsave,
       float *work1,float *work2,float *sqrtf,int *tahd,int *nfft,
       float *kxmax,float *kymax,float *angmax,float *ts,float *tg,
       float *tit,float *wt,float *ampt,int *nz,float *dz,float *fz,
       float *tb,float *pb,int *nr,float *dr);    

float disp2l(float a, float b, float c,float x, float y);
void xy2abc(float x1, float y1, float x2, float y2,
        float *a, float *b, float *c);
float disaper(float aperx,float apery,float x1,float y1,float x2,float y2);

segytrace tra, tr;
segybhdr bh;
segychdr ch;


main(int argc, char **argv)
{
	int nlines, ncdppl;
   	int nt,nzo,nz,nz0,it,ix,iy;
    	FILE *infp,*outfp,*ttfp,*imgfp,*hdrfp,*hffp,*fldfp;
	float x1,x2,x3,y1,y2,y3,s1,s2,s3,l1,l2,l3;
	char *backupi="", *backupo="";
	char *diskimg="DISKIMAGE", *diskhdr="DISKHEADER";