added SI appplication for the calculation of virtual sources on receiver postions 3D

corrvir/Makefile

+# Makefile for corrvir
+
+include ../Make_include
+
+ALLINC  = -I.
+LIBS    += -L$L -lgenfft -lm 
+#OPTC += -g
+
+all: corrvir
+
+PRG = corrvir
+
+SRCC	= $(PRG).c \
+		getFileInfo.c \
+		correlate.c \
+		atopkge.c \
+		docpkge.c \
+		wallclock_time.c  \
+		getpars.c
+
+OBJC	= $(SRCC:%.c=%.o)
+
+$(PRG):	$(OBJC) corr.h
+	$(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o $(PRG)  $(OBJC) $(LIBS)
+
+install: $(PRG) 
+	cp $(PRG) $B
+
+clean:
+		rm -f core $(OBJC) $(OBJM) $(PRG)
+
+realclean:
+		rm -f core $(OBJC) $(OBJM) $(PRG) $B/$(PRG) 
+
+
+print:	Makefile $(SRC)
+	$(PRINT) $?
+	@touch print
+
+count:
+	@wc $(SRC)
+
+tar:
+	@tar cf $(PRG).tar Makefile $(SRC) && compress $(PRG).tar
+
+
+

corrvir/atopkge.c

+../utils/atopkge.c
\ No newline at end of file

corrvir/corr.h

+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+
+typedef struct _tracePos { /* Type */
+	int gx;
+	int gy;
+	int gelev;
+} tracePos;
+
+typedef struct _traceCoord { /* Type */
+	int x;
+	int y;
+	int peg;
+	int fpos;
+} traceCoord;
+
+typedef struct _crgPos { /* Type */
+	int gx;
+	int gy;
+	int rpeg;
+	int gelev;
+	int nsrc;
+	traceCoord *src;
+} crgPos;
+

corrvir/correlate.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+
+double wallclock_time(void);
+
+typedef struct { /* complex number */
+	float r,i;
+} complex;
+
+int correlate(complex *cmaster, complex *cslaves, int nfreq, int ncor, double *tcorr, int verbose)
+{
+	int j, istation, icc;
+	double  t0, t1;
+	complex cC;
+
+	t0 = wallclock_time();
+	for (istation=0; istation<ncor; istation++) {
+		icc = istation*nfreq;
+
+#pragma ivdep
+		for (j=0; j<nfreq; j++) {
+			/* A*B */
+			cC.r = cmaster[j+icc].r*cslaves[icc+j].r+cmaster[j+icc].i*cslaves[icc+j].i;
+			cC.i = cmaster[j+icc].r*cslaves[icc+j].i-cmaster[j+icc].i*cslaves[icc+j].r;
+			/* AB* */
+//			cC.r = cmaster[j+icc].r*cslaves[icc+j].r+cmaster[j+icc].i*cslaves[icc+j].i;
+//			cC.i = -cmaster[j+icc].r*cslaves[icc+j].i+cmaster[j+icc].i*cslaves[icc+j].r;
+			/* A*B + AB* */
+//			cC.r = cmaster[j+icc].r*cslaves[icc+j].r+cmaster[j+icc].i*cslaves[icc+j].i;
+//			cC.i = 0.0;
+
+			cslaves[icc+j] = cC;
+		}
+
+	}
+	t1 = wallclock_time();
+	*tcorr += t1-t0;
+
+	return 0;
+}
+
+int coherence(complex *cmaster, complex *cslaves, int nfreq, int ncor, float reps, float epsmax, double *tcorr, int verbose)
+{
+	int j, istation, icc;
+	float maxden, *den, leps, am1, am2, scl;
+	double  t0, t1;
+	complex cC;
+
+	t0 = wallclock_time();
+    den = (float *)malloc(nfreq*ncor*sizeof(float));
+    assert(den != NULL);
+
+	for (istation=0; istation<ncor; istation++) {
+		icc = istation*nfreq;
+		for (j=0; j<nfreq; j++) {
+			am1 = sqrt(cmaster[j+icc].r*cmaster[j+icc].r+cmaster[j+icc].i*cmaster[j+icc].i);
+			am2 = sqrt(cslaves[j+icc].r*cslaves[j+icc].r+cslaves[j+icc].i*cslaves[j+icc].i);
+			den[j+icc] = am1*am2;
+		}
+	}
+    leps = reps*maxden;
+
+	for (istation=0; istation<ncor; istation++) {
+		icc = istation*nfreq;
+
+		maxden=0.0;
+		for (j=0; j<nfreq; j++) {
+			maxden = MAX(maxden, den[j+icc]);
+		}
+    	leps = reps*maxden;
+#pragma ivdep
+		for (j=0; j<nfreq; j++) {
+			/* A*B */
+	        if (den[j+icc]>epsmax*maxden) scl = 1.0/(den[j+icc]); 
+			else if (den[j+icc]<epsmax*maxden && den[j+icc]!=0) scl = 1.0/(den[j+icc]+leps);
+   	        else if (den[j+icc]==0) scl = 1.0;
+
+			cC.r = (cmaster[j+icc].r*cslaves[icc+j].r+cmaster[j+icc].i*cslaves[icc+j].i)*scl;
+			cC.i = (cmaster[j+icc].r*cslaves[icc+j].i-cmaster[j+icc].i*cslaves[icc+j].r)*scl;
+
+			cslaves[icc+j] = cC;
+		}
+
+	}
+	free(den);
+
+	t1 = wallclock_time();
+	*tcorr += t1-t0;
+
+	return 0;
+}

corrvir/demo/compvir_all.sh

+#!/bin/bash
+
+filename_in=/vardim/home/boulleng/Canada_Data/Lalor_shots_clean_4000ms_agc_sortgdelsdel_sdel136out.su
+filename_out=virshots_fmax160_FROM_all.su
+
+
+ulimit -c unlimited
+../corrvir \
+	file_shots=${filename_in} \
+	file_out=${filename_out} \
+	nsources=908 \
+	src_sel=0 \
+	fmax=160 \
+	nbm=1 \
+	nreceivers=2685 \
+	normsrc=0 \
+	normalize=0 \
+	cohr=0 \
+	causal=1 \
+	verbose=4

corrvir/docpkge.c

+../utils/docpkge.c
\ No newline at end of file

corrvir/getFileInfo.c

+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include "par.h"
+#include "segy.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+int getFileInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, int verbose)
+{
+    FILE    *fp;
+    size_t  nread;
+	off_t bytes, ret, trace_sz, ntraces;
+    int i, itrace, one_shot, igath;
+    float *trace, cmin;
+    segy hdr;
+    
+    fp = fopen( file_name, "r" );
+    assert( fp != NULL);
+    nread = fread( &hdr, 1, TRCBYTES, fp );
+    assert(nread == TRCBYTES);
+    ret = fseeko( fp, 0, SEEK_END );
+	if (ret<0) perror("fseeko");
+    bytes = ftello( fp );
+
+	*n1 = hdr.ns;
+	*d1 =1e-6*hdr.dt;
+	*d2 = hdr.d2;
+
+    trace_sz = sizeof(float)*(*n1)+TRCBYTES;
+    ntraces  = (int) (bytes/trace_sz);
+	*n2 = ntraces;
+
+	if (verbose) {
+		fprintf(stderr,"file_base=%s has nt=%d dt=%f and %d traces\n", file_name, *n1, *d1, *n2);
+	}
+
+    return 0;
+}
+

corrvir/getpars.c

+../utils/getpars.c
\ No newline at end of file

corrvir/par.h

+../utils/par.h
\ No newline at end of file

corrvir/segy.h

+/* Copyright (c) Colorado School of Mines, 2001.*/
+/* All rights reserved.                       */
+
+/* segy.h - include file for SEGY traces
+ *
+ * declarations for:
+ *	typedef struct {} segy - the trace identification header
+ *	typedef struct {} bhed - binary header
+ *
+ * Note:
+ *	If header words are added, run the makefile in this directory
+ *	to recreate hdr.h.
+ *
+ * Reference:
+ *	K. M. Barry, D. A. Cavers and C. W. Kneale, "Special Report:
+ *		Recommended Standards for Digital Tape Formats",
+ *		Geophysics, vol. 40, no. 2 (April 1975), P. 344-352.
+ *	
+ * $Author: john $
+ * $Source: /usr/local/cwp/src/su/include/RCS/segy.h,v $
+ * $Revision: 1.23 $ ; $Date: 1998/03/26 23:48:18 $
+ */ 
+
+#define TRCBYTES                240
+
+
+/* TYPEDEFS */
+typedef struct {	/* segy - trace identification header */
+
+	int tracl;	/* trace sequence number within line */
+
+	int tracr;	/* trace sequence number within reel */
+
+	int fldr;	/* field record number */
+
+	int tracf;	/* trace number within field record */
+
+	int ep;	/* energy source point number */
+
+	int cdp;	/* CDP ensemble number */
+
+	int cdpt;	/* trace number within CDP ensemble */
+
+	short trid;	/* trace identification code:
+			1 = seismic data
+			2 = dead
+			3 = dummy
+			4 = time break
+			5 = uphole
+			6 = sweep
+			7 = timing
+			8 = water break
+			9---, N = optional use (N = 32,767)
+
+			Following are CWP id flags:
+
+			 9 = autocorrelation
+
+			10 = Fourier transformed - no packing
+			     xr[0],xi[0], ..., xr[N-1],xi[N-1]
+
+			11 = Fourier transformed - unpacked Nyquist
+			     xr[0],xi[0],...,xr[N/2],xi[N/2]
+
+			12 = Fourier transformed - packed Nyquist
+	 		     even N:
+			     xr[0],xr[N/2],xr[1],xi[1], ...,
+				xr[N/2 -1],xi[N/2 -1]
+				(note the exceptional second entry)
+			     odd N:
+			     xr[0],xr[(N-1)/2],xr[1],xi[1], ...,
+				xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2]
+				(note the exceptional second & last entries)
+
+			13 = Complex signal in the time domain
+			     xr[0],xi[0], ..., xr[N-1],xi[N-1]
+
+			14 = Fourier transformed - amplitude/phase
+			     a[0],p[0], ..., a[N-1],p[N-1]
+
+			15 = Complex time signal - amplitude/phase
+			     a[0],p[0], ..., a[N-1],p[N-1]
+
+			16 = Real part of complex trace from 0 to Nyquist
+
+			17 = Imag part of complex trace from 0 to Nyquist
+
+			18 = Amplitude of complex trace from 0 to Nyquist
+
+			19 = Phase of complex trace from 0 to Nyquist
+
+			21 = Wavenumber time domain (k-t)
+
+			22 = Wavenumber frequency (k-omega)
+
+			23 = Envelope of the complex time trace
+
+			24 = Phase of the complex time trace
+
+			25 = Frequency of the complex time trace
+
+			30 = Depth-Range (z-x) traces
+
+			43 = Seismic Data, Vertical Component 
+
+			44 = Seismic Data, Horizontal Component 1 
+
+			45 = Seismic Data, Horizontal Component 2 
+
+			46 = Seismic Data, Radial Component
+
+			47 = Seismic Data, Transverse Component  
+
+			101 = Seismic data packed to bytes (by supack1)
+			
+			102 = Seismic data packed to 2 bytes (by supack2)
+			*/
+
+	short nvs;	/* number of vertically summed traces (see vscode
+			   in bhed structure) */
+
+	short nhs;	/* number of horizontally summed traces (see vscode
+			   in bhed structure) */
+
+	short duse;	/* data use:
+				1 = production
+				2 = test */
+
+	int offset;	/* distance from source point to receiver
+			   group (negative if opposite to direction
+			   in which the line was shot) */
+
+	int gelev;	/* receiver group elevation from sea level
+			   (above sea level is positive) */
+
+	int selev;	/* source elevation from sea level
+			   (above sea level is positive) */
+
+	int sdepth;	/* source depth (positive) */
+
+	int gdel;	/* datum elevation at receiver group */
+
+	int sdel;	/* datum elevation at source */
+
+	int swdep;	/* water depth at source */
+
+	int gwdep;	/* water depth at receiver group */
+
+	short scalel;	/* scale factor for previous 7 entries
+			   with value plus or minus 10 to the
+			   power 0, 1, 2, 3, or 4 (if positive,
+			   multiply, if negative divide) */
+
+	short scalco;	/* scale factor for next 4 entries
+			   with value plus or minus 10 to the
+			   power 0, 1, 2, 3, or 4 (if positive,
+			   multiply, if negative divide) */
+
+	int  sx;	/* X source coordinate */
+
+	int  sy;	/* Y source coordinate */
+
+	int  gx;	/* X group coordinate */
+
+	int  gy;	/* Y group coordinate */
+
+	short counit;	/* coordinate units code:
+				for previous four entries
+				1 = length (meters or feet)
+				2 = seconds of arc (in this case, the
+				X values are longitude and the Y values
+				are latitude, a positive value designates
+				the number of seconds east of Greenwich
+				or north of the equator */
+
+	short wevel;	/* weathering velocity */
+
+	short swevel;	/* subweathering velocity */
+
+	short sut;	/* uphole time at source */
+
+	short gut;	/* uphole time at receiver group */
+
+	short sstat;	/* source static correction */
+
+	short gstat;	/* group static correction */
+
+	short tstat;	/* total static applied */
+
+	short laga;	/* lag time A, time in ms between end of 240-
+			   byte trace identification header and time
+			   break, positive if time break occurs after
+			   end of header, time break is defined as
+			   the initiation pulse which maybe recorded
+			   on an auxiliary trace or as otherwise
+			   specified by the recording system */
+
+	short lagb;	/* lag time B, time in ms between the time break
+			   and the initiation time of the energy source,
+			   may be positive or negative */
+
+	short delrt;	/* delay recording time, time in ms between
+			   initiation time of energy source and time
+			   when recording of data samples begins
+			   (for deep water work if recording does not
+			   start at zero time) */
+
+	short muts;	/* mute time--start */
+
+	short mute;	/* mute time--end */
+
+	unsigned short ns;	/* number of samples in this trace */
+
+	unsigned short dt;	/* sample interval; in micro-seconds */
+
+	short gain;	/* gain type of field instruments code:
+				1 = fixed
+				2 = binary
+				3 = floating point
+				4 ---- N = optional use */
+
+	short igc;	/* instrument gain constant */
+
+	short igi;	/* instrument early or initial gain */
+
+	short corr;	/* correlated:
+				1 = no
+				2 = yes */
+
+	short sfs;	/* sweep frequency at start */
+
+	short sfe;	/* sweep frequency at end */
+
+	short slen;	/* sweep length in ms */
+
+	short styp;	/* sweep type code:
+				1 = linear
+				2 = cos-squared
+				3 = other */
+
+	short stas;	/* sweep trace length at start in ms */
+
+	short stae;	/* sweep trace length at end in ms */
+
+	short tatyp;	/* taper type: 1=linear, 2=cos^2, 3=other */
+
+	short afilf;	/* alias filter frequency if used */
+
+	short afils;	/* alias filter slope */
+
+	short nofilf;	/* notch filter frequency if used */
+
+	short nofils;	/* notch filter slope */
+
+	short lcf;	/* low cut frequency if used */
+
+	short hcf;	/* high cut frequncy if used */
+
+	short lcs;	/* low cut slope */
+
+	short hcs;	/* high cut slope */
+
+	short year;	/* year data recorded */
+
+	short day;	/* day of year */
+
+	short hour;	/* hour of day (24 hour clock) */
+
+	short minute;	/* minute of hour */
+
+	short sec;	/* second of minute */
+
+	short timbas;	/* time basis code:
+				1 = local
+				2 = GMT
+				3 = other */
+
+	short trwf;	/* trace weighting factor, defined as 1/2^N
+			   volts for the least sigificant bit */
+
+	short grnors;	/* geophone group number of roll switch
+			   position one */
+
+	short grnofr;	/* geophone group number of trace one within
+			   original field record */
+
+	short grnlof;	/* geophone group number of last trace within
+			   original field record */
+
+	short gaps;	/* gap size (total number of groups dropped) */
+
+	short otrav;	/* overtravel taper code:
+				1 = down (or behind)
+				2 = up (or ahead) */
+
+	/* local assignments */
+	float d1;	/* sample spacing for non-seismic data */
+
+	float f1;	/* first sample location for non-seismic data */
+
+	float d2;	/* sample spacing between traces */
+
+	float f2;	/* first trace location */
+
+	float ungpow;	/* negative of power used for dynamic
+			   range compression */
+
+	float unscale;	/* reciprocal of scaling factor to normalize
+			   range */
+
+	int ntr; 	/* number of traces */
+
+	short mark;	/* mark selected traces */
+
+        short shortpad; /* alignment padding */
+
+
+	short unass[14];	/* unassigned--NOTE: last entry causes 
+			   a break in the word alignment, if we REALLY
+			   want to maintain 240 bytes, the following
+			   entry should be an odd number of short/UINT2
+			   OR do the insertion above the "mark" keyword
+			   entry */
+
+} segy;
+

corrvir/wallclock_time.c

+../utils/wallclock_time.c
\ No newline at end of file