-
joeri.brackenhoff authoredjoeri.brackenhoff authored
defineSource.c 11.42 KiB
#define _FILE_OFFSET_BITS 64
#define _LARGEFILE_SOURCE
#define _LARGEFILE64_SOURCE
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <math.h>
#include <string.h>
#include "fdelmodc.h"
#include "segy.h"
/**
* Computes, or read from file, the source signature
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
float r,i;
} complex;
typedef struct _dcomplexStruct { /* complex number */
double r,i;
} dcomplex;
#endif/* complex */
int optncr(int n);
void rc1fft(float *rdata, complex *cdata, int n, int sign);
void cr1fft(complex *cdata, float *rdata, int n, int sign);
int writesufile(char *filename, float *data, int n1, int n2, float f1, float f2, float d1, float d2);
int writesufilesrcnwav(char *filename, float **src_nwav, wavPar wav, int n1, int n2, float f1, float f2, float d1, float d2);
float gaussGen();
float normal(double x,double mu,double sigma);
int comp (const float *a, const float *b);
void spline3(float x1, float x2, float z1, float z2, float dzdx1, float dzdx2, float *a, float *b, float *c, float *d);
int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend, int verbose);
/* random number generators */
double dcmwc4096();
unsigned long CMWC4096(void);
unsigned long xorshift(void);
void seedCMWC4096(void);
/* #define drand48 dcmwc4096 use for different random number generator */
#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 defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwav, int reverse, int verbose)
{
FILE *fp;
size_t nread;
int optn, nfreq, i, j, k, iwmax, tracesToDo;
int iw, n1, namp, optnscale, nfreqscale;
float scl, d1, df, deltom, om, tshift;
float amp1, amp2, amp3;
float *trace, maxampl, scale;
complex *ctrace, tmp;
segy hdr;
scale = 1.0;
n1 = wav.ns;
if (wav.random) { /* initialize random sequence */
srand48(wav.seed+1); seedCMWC4096();
for (i=0; i<8192; i++) {
amp1 = dcmwc4096();
}
}
else {
/* read first header and last byte to get file size */
fp = fopen( wav.file_src, "r" );
assert( fp != NULL);
nread = fread( &hdr, 1, TRCBYTES, fp );
assert(nread == TRCBYTES);
/* read all traces */
tracesToDo = wav.nx;
i = 0;
while (tracesToDo) {
memset(&src_nwav[i][0],0,wav.nt*sizeof(float));
nread = fread(&src_nwav[i][0], sizeof(float), hdr.ns, fp);
assert (nread == hdr.ns);
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread==0) break;
tracesToDo--;
i++;
}
fclose(fp);
}
optn = optncr(n1);
nfreq = optn/2 + 1;
if (wav.nt != wav.ns) {
vmess("Sampling in wavelet is %e while for modeling is set to %e", wav.ds, mod.dt);
vmess("Wavelet sampling will be FFT-interpolated to sampling of modeling");
vmess("file_src Nt=%d sampling after interpolation=%d", wav.ns, wav.nt);
optnscale = wav.nt;
nfreqscale = optnscale/2 + 1;
}
else {
optnscale = optn;
nfreqscale = optnscale/2 + 1;
}
// fprintf(stderr,"define S optn=%d ns=%d %e nt=%d %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
ctrace = (complex *)calloc(nfreqscale,sizeof(complex));
trace = (float *)calloc(optnscale,sizeof(float));
df = 1.0/(optn*wav.ds);
deltom = 2.*M_PI*df;
scl = 1.0/optn;
iwmax = nfreq;
for (i=0; i<wav.nx; i++) {
if (wav.random) {
randomWavelet(wav, src, &src_nwav[i][0], src.tbeg[i], src.tend[i], verbose);
}
else {
memset(&ctrace[0].r,0,nfreqscale*sizeof(complex));
memset(&trace[0],0,optnscale*sizeof(float));
memcpy(&trace[0],&src_nwav[i][0],n1*sizeof(float));
rc1fft(trace,ctrace,optn,-1);
/* Scale source from file with -j/w (=1/(jw)) for volume source injections
no scaling is applied for volume source injection rates */
if (src.injectionrate==0) {
for (iw=1;iw<iwmax;iw++) {
om = 1.0/(deltom*iw);
tmp.r = om*ctrace[iw].i; tmp.i = -om*ctrace[iw].r;
ctrace[iw].r = tmp.r;
ctrace[iw].i = tmp.i;
}
}
if (src.type < 6) { // shift wavelet with +1/2 DeltaT due to staggered in time
tshift=-(0.5*rec.skipdt+1.5)*wav.dt;
for (iw=1;iw<iwmax;iw++) {
om = deltom*iw*tshift;
tmp.r = ctrace[iw].r*cos(-om) - ctrace[iw].i*sin(-om);
tmp.i = ctrace[iw].i*cos(-om) + ctrace[iw].r*sin(-om);
ctrace[iw].r = tmp.r;
ctrace[iw].i = tmp.i;
}
}
/* zero frequency iw=0 set to 0 if the next sample has amplitude==0*/
amp1 = sqrt(ctrace[1].r*ctrace[1].r+ctrace[1].i*ctrace[1].i);
if (amp1 == 0.0) {
ctrace[0].r = ctrace[0].i = 0.0;
}
else { /* stabilization for w=0: extrapolate amplitudes to 0 */
amp2 = sqrt(ctrace[2].r*ctrace[2].r+ctrace[2].i*ctrace[2].i);
amp3 = sqrt(ctrace[3].r*ctrace[3].r+ctrace[3].i*ctrace[3].i);
ctrace[0].r = amp1+(2.0*(amp1-amp2)-(amp2-amp3));
ctrace[0].i = 0.0;
if (ctrace[1].r < 0.0) {
ctrace[0].r *= -1.0;
}
}
for (iw=iwmax;iw<nfreqscale;iw++) {
ctrace[iw].r = 0.0;
ctrace[iw].i = 0.0;
}
memset(&trace[0],0,optnscale*sizeof(float));
cr1fft(ctrace,trace,optnscale,1);
/* avoid a (small) spike in the last sample
this is done to avoid diffraction from last wavelet sample
which will act as a pulse */
maxampl=0.0;
if (reverse) {
for (j=0; j<wav.nt; j++) {
src_nwav[i][j] = scl*(trace[wav.nt-j-1]-trace[0]);
maxampl = MAX(maxampl,fabs(src_nwav[i][j]));
}
}
else {
for (j=0; j<wav.nt; j++) {
src_nwav[i][j] = scl*(trace[j]-trace[wav.nt-1]);
maxampl = MAX(maxampl,fabs(src_nwav[i][j]));
}
}
if (verbose > 3) vmess("Wavelet sampling (FFT-interpolated) done for trace %d", i);
}
}
/* set values smaller than 1e-5 maxampl to zero */
maxampl *= 1e-5;
for (i=0; i<wav.nx; i++) {
for (j=0; j<wav.nt; j++) {
if (fabs(src_nwav[i][j]) < maxampl) src_nwav[i][j] = 0.0;
}
}
free(ctrace);
free(trace);
/* use random amplitude gain factor for each source */
if (src.amplitude > 0.0) {
namp=wav.nx*10; trace = (float *)calloc(2*namp,sizeof(float));
for (i=0; i<wav.nx; i++) {
if (src.distribution) {
scl = gaussGen()*src.amplitude;
k = (int)MAX(MIN(namp*(scl+5*src.amplitude)/(10*src.amplitude),namp-1),0);
d1 = 10.0*src.amplitude/(namp-1);
}
else {
scl = (float)(drand48()-0.5)*src.amplitude;
k = (int)MAX(MIN(namp*(scl+1*src.amplitude)/(2*src.amplitude),namp-1),0);
d1 = 2.0*src.amplitude/(namp-1);
}
trace[k] += 1.0;
/* trace[i] = scl; */
if (wav.random) n1 = wav.nsamp[i];
else n1 = wav.nt;
for (j=0; j<n1; j++) {
src_nwav[i][j] *= scl;
}
}
if (verbose>2) writesufile("src_ampl.su", trace, namp, 1, -5*src.amplitude, 0.0, d1, 1);
/*
qsort(trace,wav.nx,sizeof(float), comp);
for (i=0; i<wav.nx; i++) {
scl = trace[i];
trace[i] = normal(scl, 0.0, src.amplitude);
}
if (verbose>2) writesufile("src_ampl.su", trace, wav.nx, 1, -5*src.amplitude, 0.0, d1, 1);
*/
free(trace);
}
if (verbose>3) writesufilesrcnwav("src_nwav.su", src_nwav, wav, wav.nt, wav.nx, 0.0, 0.0, wav.dt, 1);
/* set maximum amplitude in source file to 1.0 */
/*
assert(maxampl != 0.0);
scl = wav.dt/(maxampl);
scl = 1.0/(maxampl);
for (i=0; i<wav.nx; i++) {
for (j=0; j<n1; j++) {
src_nwav[i*n1+j] *= scl;
}
}
*/
return 0;
}
int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend, int verbose)
{
int optn, nfreq, j, iwmax;
int iw, n1, itbeg, itmax, nsmth;
float df, amp1;
float *rtrace;
float x1, x2, z1, z2, dzdx1, dzdx2, a, b, c, d, t;
complex *ctrace;
n1 = wav.nt; /* this is set to the maximum length (tlength/dt) */
optn = optncr(2*n1);
nfreq = optn/2 + 1;
ctrace = (complex *)calloc(nfreq,sizeof(complex));
rtrace = (float *)calloc(optn,sizeof(float));
df = 1.0/(optn*wav.dt);
iwmax = MIN(NINT(wav.fmax/df),nfreq);
for (iw=1;iw<iwmax;iw++) {
ctrace[iw].r = (float)(drand48()-0.5);
ctrace[iw].i = (float)(drand48()-0.5);
}
for (iw=iwmax;iw<nfreq;iw++) {
ctrace[iw].r = 0.0;
ctrace[iw].i = 0.0;
}
cr1fft(ctrace,rtrace,optn,1);
/* find first zero crossing in wavelet */
amp1 = rtrace[0];
j = 1;
if (amp1 < 0.0) {
while (rtrace[j] < 0.0) j++;
}
else {
while (rtrace[j] > 0.0) j++;
}
itbeg = j;
/* find last zero crossing in wavelet */
// itmax = itbeg+MIN(NINT((tend-tbeg)/wav.dt),n1);
itmax = MIN(NINT(itbeg+(tend-tbeg)/wav.dt),n1);
amp1 = rtrace[itmax-1];
j = itmax;
if (amp1 < 0.0) {
while (rtrace[j] < 0.0 && j>itbeg) j--;
}
else {
while (rtrace[j] > 0.0 && j>itbeg) j--;
}
itmax = j;
/* make smooth transitions to zero aamplitude */
nsmth=MIN(10,itmax);
x1 = 0.0;
z1 = 0.0;
dzdx1 = 0.0;
x2 = nsmth;
z2 = rtrace[itbeg+nsmth];
// dzdx2 = (rtrace[itbeg+(nsmth+1)]-rtrace[itbeg+(nsmth-1)])/(2.0);
dzdx2 = (rtrace[itbeg+nsmth-2]-8.0*rtrace[itbeg+nsmth-1]+
8.0*rtrace[itbeg+nsmth+1]-rtrace[itbeg+nsmth+2])/(12.0);
spline3(x1, x2, z1, z2, dzdx1, dzdx2, &a, &b, &c, &d);
for (j=0; j<nsmth; j++) {
t = j;
rtrace[itbeg+j] = a*t*t*t+b*t*t+c*t+d;
}
x1 = 0.0;
z1 = rtrace[itmax-nsmth];
// dzdx1 = (rtrace[itmax-(nsmth-1)]-rtrace[itmax-(nsmth+1)])/(2.0);
dzdx1 = (rtrace[itmax-nsmth-2]-8.0*rtrace[itmax-nsmth-1]+
8.0*rtrace[itmax-nsmth+1]-rtrace[itmax-nsmth+2])/(12.0);
x2 = nsmth;
z2 = 0.0;
dzdx2 = 0.0;
// fprintf(stderr,"x1=%f z1=%f d=%f x2=%f, z2=%f d=%f\n",x1,z1,dzdx1,x2,z2,dzdx2);
spline3(x1, x2, z1, z2, dzdx1, dzdx2, &a, &b, &c, &d);
for (j=0; j<nsmth; j++) {
t = j;
rtrace[itmax-nsmth+j] = a*t*t*t+b*t*t+c*t+d;
// fprintf(stderr,"a=%f b=%f c=%f d=%f rtrace%d=%f \n",a,b,c,d,j,rtrace[itmax-nsmth+j]);
}
for (j=itbeg; j<itmax; j++) trace[j-itbeg] = rtrace[j];
free(ctrace);
free(rtrace);
return 0;
}
float normal(double x,double mu,double sigma)
{
return (float)(1.0/(2.0*M_PI*sigma*sigma))*exp(-1.0*(((x-mu)*(x-mu))/(2.0*sigma*sigma)) );
}
int comp (const float *a, const float *b)
{
if (*a==*b)
return 0;
else
if (*a < *b)
return -1;
else
return 1;
}