sudoc:su_prog_rules
Differences
This shows you the differences between two versions of the page.
Next revision | Previous revision | ||
sudoc:su_prog_rules [2019/07/30 20:52] – created seisunix | sudoc:su_prog_rules [2019/07/30 21:55] (current) – seisunix | ||
---|---|---|---|
Line 1: | Line 1: | ||
'' | '' | ||
- | Creating a new Seismic | + | ====== Development rules for Seismic |
+ | Seismic Unix has remained a stable, easy to install, and transportable package by following a few simple rules. | ||
+ | ===== Don't change the header structure ==== | ||
+ | The first of these is that **changing the SU header definitions is forbidden.** The header structure looks deceptively simple, but contains a hidden trap, that it is possible to easily introduce a gap in the header that would render the SU data made under a new system mutually incompatible with that made by other versions of SU. | ||
+ | |||
+ | Suppose, however, you need another header field? The way to approach this is to use one of the less commonly used fields, define a keyword field so that the user may change this. | ||
+ | |||
+ | ===== Avoid mix-language programming ===== | ||
+ | |||
+ | One of the stumbling blocks of older versions of SU was the use of subroutines written in different languages, such as Fortran, made porting the code quite difficult. Expert programmers likely will scoff at this statement, but remember, SU is designed so that non-expert users may easily install, use, and expand this package. | ||
+ | |||
+ | ===== Avoid making SU dependent on Third Party packages ===== | ||
+ | |||
+ | Many users have wondered why SU does not depend on such widely available packages such as FFTW (Fastest Fourier Transform in the West) or GSL (the GNU scientific library). The basic reason is that SU predates many of the available open source packages that users have available. In this way, SU has always been developed from the point of view of being completely self-contained. | ||
+ | |||
+ | The reality is that the more items that a user must fetch to install your code, the less likely that they will want to use it. You will also be at the mercy of whims of third party developers who may decide arbitrarily to change their code. | ||
+ | |||
+ | ===== Do not include commercial software ===== | ||
+ | |||
+ | There are a number of " | ||
+ | |||
+ | ===== Be conscious of the open source licenses ===== | ||
+ | |||
+ | The SU materials are released under a Berkeley-style open source license. This choice avoids the viral aspect of the GNU Public License. Commercial developers may think twice about making use of GPL Licensed items, because of the requirement that improvements must be returned to the originators of code, and that derivative code automatically inherits the GPL License. | ||
+ | |||
+ | The limitation is that while code released under a Berkeley License may be absorbed into a GPL software collection, the reverse is not true. You may not include GPL licensed code in the SU distribution. | ||
+ | |||
+ | ===== Avoid using global variables ===== | ||
+ | If you need to pass a variable to a subroutine then pass it as an argument to the subroutine, do not define it as a global variable. | ||
+ | |||
+ | ===== Document your code thoroughly ===== | ||
+ | Make sure there are Notes: in the Self Documentation. | ||
+ | Document each declared variable. Put technical information describing the algorithm being used. Remember that " | ||
+ | |||
+ | Document subroutines as though they will be clipped out and put in a library---because if a subroutine is useful, it will be separated from its parent program. | ||
+ | |||
+ | Pretend that the code will not run without the documentation. | ||
+ | |||
+ | ===== Develop from existing code ==== | ||
+ | |||
+ | An example program is **sufiler.c** which is located in **$CWPROOT/ | ||
+ | |||
+ | Proceeding block by block in the code, | ||
+ | |||
+ | **The copyright statement ** | ||
+ | < | ||
+ | /* Copyright (c) Colorado School of Mines, 2011.*/ | ||
+ | /* All rights reserved. | ||
+ | |||
+ | </ | ||
+ | |||
+ | **The RCS (revision control system) version label ** | ||
+ | < | ||
+ | /* SUFILTER: $Revision: 1.23 $ ; $Date: 2011/11/12 00:09:00 $ */ | ||
+ | </ | ||
+ | |||
+ | **The include files that define the SU data type** | ||
+ | < | ||
+ | #include " | ||
+ | #include " | ||
+ | </ | ||
+ | |||
+ | ** The begin selfdoc flag which allows updated to capture the self documentation for tools such as suname ** | ||
+ | < | ||
+ | / | ||
+ | </ | ||
+ | |||
+ | **The self documentation ** | ||
+ | < | ||
+ | char *sdoc[] = { | ||
+ | " | ||
+ | " SUFILTER - applies a zero-phase, sine-squared tapered filter | ||
+ | " | ||
+ | " sufilter <stdin >stdout [optional parameters] | ||
+ | " | ||
+ | " Required parameters: | ||
+ | " | ||
+ | " | ||
+ | " Optional parameters: | ||
+ | " | ||
+ | " | ||
+ | " | ||
+ | " | ||
+ | " | ||
+ | " Defaults: | ||
+ | " | ||
+ | " | ||
+ | " | ||
+ | " Examples of filters: | ||
+ | " Bandpass: | ||
+ | " Bandreject: sufilter <data f=10, | ||
+ | " Lowpass: | ||
+ | " Highpass: | ||
+ | " Notch: | ||
+ | NULL}; | ||
+ | </ | ||
+ | |||
+ | **The attribution block, authors, technical references, and header fields and other technical details go here** | ||
+ | < | ||
+ | /* Credits: | ||
+ | | ||
+ | | ||
+ | * | ||
+ | * Possible optimization: | ||
+ | * filter is 0.0. | ||
+ | * | ||
+ | * Trace header fields accessed: ns, dt, d1 | ||
+ | */ | ||
+ | </ | ||
+ | |||
+ | ** The end self documentation flag, this tells updatedoc that this is the end of the documentation** | ||
+ | < | ||
+ | / | ||
+ | </ | ||
+ | |||
+ | **Function prototypes** | ||
+ | < | ||
+ | /* Prototype of function used internally */ | ||
+ | void polygonalFilter(float *f, float *amps, | ||
+ | int npoly, int nfft, float dt, float *filter); | ||
+ | </ | ||
+ | |||
+ | **internally defined constants** | ||
+ | < | ||
+ | #define PIBY2 | ||
+ | #define FRAC0 | ||
+ | #define FRAC1 | ||
+ | #define FRAC2 | ||
+ | #define FRAC3 | ||
+ | #define LOOKFAC 2 /* Look ahead factor for npfao */ | ||
+ | #define PFA_MAX 720720 | ||
+ | </ | ||
+ | |||
+ | **Declaration of SEGY type objects in global space.** | ||
+ | < | ||
+ | segy tr; | ||
+ | </ | ||
+ | |||
+ | ** The main ** | ||
+ | < | ||
+ | int | ||
+ | main(int argc, char **argv) | ||
+ | { | ||
+ | </ | ||
+ | |||
+ | ** Declarations ** | ||
+ | < | ||
+ | register float *rt; /* real trace */ | ||
+ | register complex *ct; /* complex transformed trace */ | ||
+ | float *filter; | ||
+ | float *f; /* array of filter frequencies | ||
+ | int npoly; | ||
+ | float *amps; | ||
+ | int namps; | ||
+ | int icount, | ||
+ | float dt; /* sample spacing | ||
+ | float nyq; /* nyquist frequency | ||
+ | int nt; /* number of points on input trace */ | ||
+ | int nfft; /* number of points for fft trace */ | ||
+ | int nf; /* number of frequencies (incl Nyq) */ | ||
+ | int verbose; | ||
+ | cwp_Bool seismic; | ||
+ | </ | ||
+ | |||
+ | ** Initialize the command line, or output selfdoc ** | ||
+ | < | ||
+ | /* Initialize */ | ||
+ | initargs(argc, | ||
+ | requestdoc(1); | ||
+ | </ | ||
+ | |||
+ | ** Get parameters from the first trace and from the command line** | ||
+ | < | ||
+ | /* Get info from first trace */ | ||
+ | if (!getparint(" | ||
+ | if (!gettr(& | ||
+ | seismic = ISSEISMIC(tr.trid); | ||
+ | if (seismic) { | ||
+ | if (verbose) | ||
+ | dt = ((double) tr.dt)/ | ||
+ | } | ||
+ | else { | ||
+ | if (verbose) warn(" | ||
+ | |||
+ | dt = tr.d1; | ||
+ | |||
+ | } | ||
+ | |||
+ | /* error trapping so that the user can have a default value of dt */ | ||
+ | if (!(dt || getparfloat(" | ||
+ | dt = .004; | ||
+ | if (verbose) { | ||
+ | warn(" | ||
+ | warn(" | ||
+ | } | ||
+ | } | ||
+ | |||
+ | nt = tr.ns; | ||
+ | nyq = 0.5/dt; | ||
+ | </ | ||
+ | |||
+ | ** read trace by trace and perform the filtering ** | ||
+ | |||
+ | < | ||
+ | /* Set up FFT parameters */ | ||
+ | nfft = npfaro(nt, LOOKFAC * nt); | ||
+ | if (nfft >= SU_NFLTS || nfft >= PFA_MAX) | ||
+ | err(" | ||
+ | |||
+ | nf = nfft/2 + 1; | ||
+ | |||
+ | /* Get frequencies that define the filter */ | ||
+ | if ((npoly = countparval(" | ||
+ | f = ealloc1float(npoly); | ||
+ | getparfloat(" | ||
+ | } else { | ||
+ | npoly = 4; | ||
+ | f = ealloc1float(npoly); | ||
+ | |||
+ | f[0] = FRAC0 * nyq; | ||
+ | f[1] = FRAC1 * nyq; | ||
+ | f[2] = FRAC2 * nyq; | ||
+ | f[3] = FRAC3 * nyq; | ||
+ | } | ||
+ | |||
+ | /* Check f values */ | ||
+ | if(npoly < 2) warn(" | ||
+ | for(ifs=0; ifs < npoly-1; ++ifs) | ||
+ | if(f[ifs] < 0.0 || f[ifs] > f[ifs+1]) | ||
+ | err(" | ||
+ | | ||
+ | /* Get filter amplitude values*/ | ||
+ | if ((namps = countparval(" | ||
+ | amps = ealloc1float(namps); | ||
+ | getparfloat(" | ||
+ | } else { | ||
+ | namps = npoly; | ||
+ | amps = ealloc1float(namps); | ||
+ | |||
+ | /* default is a trapezoidal bandpass filter */ | ||
+ | for(iamps=0; | ||
+ | amps[iamps]=1.; | ||
+ | | ||
+ | amps[0]=0.; amps[namps-1]=0.; | ||
+ | } | ||
+ | if (!(namps==npoly)) | ||
+ | err(" | ||
+ | | ||
+ | /* Check amps values */ | ||
+ | for(iamps = 0, icount=0; iamps < namps ; ++iamps) { | ||
+ | if( amps[iamps] > 0. ) ++icount; | ||
+ | if( amps[iamps] < 0.) err(" | ||
+ | } | ||
+ | if (icount==0) err(" | ||
+ | for(iamps = 0, icount=0; iamps < namps-1 ; ++iamps) { | ||
+ | if(!(amps[iamps]==amps[iamps+1])) ++icount; | ||
+ | } | ||
+ | if (icount==0) warn(" | ||
+ | </ | ||
+ | |||
+ | ** Dynamic allocation which hides " | ||
+ | < | ||
+ | |||
+ | /* Allocate fft arrays */ | ||
+ | rt = ealloc1float(nfft); | ||
+ | ct = ealloc1complex(nf); | ||
+ | filter = ealloc1float(nf); | ||
+ | |||
+ | /* Build the polygonal filter filter[]*/ | ||
+ | polygonalFilter(f, | ||
+ | |||
+ | </ | ||
+ | |||
+ | ** note that most of the time we loop over input traces, these are in "do -- while" loops ** | ||
+ | < | ||
+ | /* Main loop over traces */ | ||
+ | do { | ||
+ | register int i; | ||
+ | |||
+ | /* Load trace into rt (zero-padded) */ | ||
+ | memcpy((void *) rt, (const void *) tr.data, nt*FSIZE); | ||
+ | memset((void *) (rt + nt), 0 , (nfft-nt)*FSIZE); | ||
+ | |||
+ | /* FFT, filter, inverse FFT */ | ||
+ | pfarc(1, nfft, rt, ct); | ||
+ | for (i = 0; i < nf; ++i) ct[i] = crmul(ct[i], | ||
+ | pfacr(-1, nfft, ct, rt); | ||
+ | |||
+ | /* Load traces back in, recall filter had nfft factor */ | ||
+ | for (i = 0; i < nt; ++i) tr.data[i] = rt[i]; | ||
+ | |||
+ | puttr(& | ||
+ | } while (gettr(& | ||
+ | |||
+ | return(CWP_Exit()); | ||
+ | } | ||
+ | |||
+ | </ | ||
+ | |||
+ | ** Definition of the subroutine that builds the filter ** | ||
+ | < | ||
+ | void polygonalFilter(float *f, float *amps, int npoly, | ||
+ | int nfft, float dt, float *filter) | ||
+ | / | ||
+ | polygonalFilter -- polygonal filter with sin^2 tapering | ||
+ | ************************************************************************** | ||
+ | Input: | ||
+ | f | ||
+ | amps array[npoly] of amplitude values | ||
+ | npoly size of input f and amps arrays | ||
+ | dt time sampling interval | ||
+ | nfft number of points in the fft | ||
+ | |||
+ | Output: | ||
+ | filter | ||
+ | ************************************************************************** | ||
+ | Notes: Filter is to be applied in the frequency domain | ||
+ | ************************************************************************** | ||
+ | Author: | ||
+ | *************************************************************************/ | ||
+ | #define PIBY2 | ||
+ | { | ||
+ | int *intfr; | ||
+ | int icount, | ||
+ | int taper=0; | ||
+ | int nf; /* number of frequencies (incl Nyq) */ | ||
+ | int nfm1; /* nf-1 */ | ||
+ | float onfft; | ||
+ | float df; /* frequency spacing (from dt) */ | ||
+ | |||
+ | | ||
+ | intfr=alloc1int(npoly); | ||
+ | |||
+ | nf = nfft/2 + 1; | ||
+ | nfm1 = nf - 1; | ||
+ | onfft = 1.0 / nfft; | ||
+ | |||
+ | /* Compute array of integerized frequencies that define the filter*/ | ||
+ | df = onfft / dt; | ||
+ | for(ifs=0; ifs < npoly ; ++ifs) { | ||
+ | intfr[ifs] = NINT(f[ifs]/ | ||
+ | if (intfr[ifs] > nfm1) intfr[ifs] = nfm1; | ||
+ | } | ||
+ | |||
+ | /* Build filter, with scale, and taper specified by amps[] values*/ | ||
+ | /* Do low frequency end first*/ | ||
+ | for(icount=0; | ||
+ | filter[icount] = amps[0] * onfft; | ||
+ | |||
+ | /* now do the middle frequencies */ | ||
+ | for(ifs=0 ; ifs< | ||
+ | | ||
+ | ++taper; | ||
+ | for(icount=intfr[ifs]; | ||
+ | float c = PIBY2 / (intfr[ifs+1] - intfr[ifs] + 2); | ||
+ | float s = sin(c*(icount - intfr[ifs] + 1)); | ||
+ | float adiff = amps[ifs+1] - amps[ifs]; | ||
+ | filter[icount] = (amps[ifs] + adiff*s*s) * onfft; | ||
+ | } | ||
+ | } else if (amps[ifs] > amps[ifs+1]) { | ||
+ | ++taper; | ||
+ | for(icount=intfr[ifs]; | ||
+ | float c = PIBY2 / (intfr[ifs+1] - intfr[ifs] + 2); | ||
+ | float s = sin(c*(intfr[ifs+1] - icount + 1)); | ||
+ | float adiff = amps[ifs] - amps[ifs+1]; | ||
+ | | ||
+ | } | ||
+ | } else | ||
+ | if(!(taper)){ | ||
+ | for(icount=intfr[ifs]; | ||
+ | | ||
+ | } else { | ||
+ | for(icount=intfr[ifs]+1; | ||
+ | | ||
+ | } | ||
+ | } | ||
+ | |||
+ | /* finally do the high frequency end */ | ||
+ | for(icount=intfr[npoly-1]+1; | ||
+ | filter[icount] = amps[npoly-1] * onfft; | ||
+ | } | ||
+ | |||
+ | } | ||
+ | |||
+ | |||
+ | </ |
sudoc/su_prog_rules.txt · Last modified: 2019/07/30 21:55 by seisunix