Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
  • plane-wave
  • 2.1.2
  • 2.1.1
  • 2.1.0
  • 2.0.0
  • 1.4
  • 1.4.0
  • 1.3.3
  • 1.3.2
  • 1.3.1
  • 1.3.0
  • 1.2.3
  • 1.2.2
  • 1.2.1
  • 1.2.0
  • 1.1.4
  • 1.1.3
  • 1.1.2
  • 1.1.1
  • 1.1.0
  • 1.0.0
22 results
Blame Permalink
HomG.c 13.78 KiB 
#include "par.h"
#include "segy.h"
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <genfft.h>

#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#endif
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))

#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
    float r,i;
} complex;
#endif/* complex */

int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
double wallclock_time(void);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
int topdet(float *data, int nt);

void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout);
void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout);
void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift);

char *sdoc[] = {
" ",
" HomG - Calculate a Homogeneous Green's function ",
" ",
" authors  : Joeri Brackenhoff 	(J.A.Brackenhoff@tudelft.nl)",
"		   : Jan Thorbecke		(janth@xs4all.nl)",
" ",
" Required parameters: ",
"",
"   file_in= ................. First file of the array of receivers",
"   file_shot= ............... File containing the shot location",
" ",
" Optional parameters: ",
" ",
"   file_out= ................ Filename of the output",
"   numb= .................... integer number of first snapshot file",
"   dnumb= ................... integer number of increment in snapshot files",
"   zmax= .................... Integer number of maximum depth level",
"   inx= ..................... Number of sources per depth level",
"   zrcv= .................... z-coordinate of first receiver location",
"   xrcv= .................... x-coordinate of first receiver location",
"   dzrcv= ................... z-spacing of receivers",
"   dxrcv= ................... x-spacing of receivers",
NULL};

int main (int argc, char **argv)
{
	FILE *fp_in, *fp_shot, *fp_out;
	char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100];
	float *indata, *Ghom, *shotdata, *rtrace, *costaper, scl, rho;
	float dt, dx, t0, x0, xmin, xmax1, sclsxgx, f1, f2, dxrcv, dzrcv, dxpos, offset, dw, *taper;
	int nshots, nt, nw, nx, ntraces, ret, ix, it, is, ir, pos, ifile, file_det, nxs, nzs, sxmin, sxmax;
	int pos1, xcount, zcount, npos, zmax, file_cl, ht, inx, numb, dnumb, indrcv, shift;
	int rmt, smooth, *tol, tolside, tolset, mode, ntap, maxoffset, offset_tmp, count;
	complex *chom, *cshot, *ctrace;
	segy *hdr_in, *hdr_out, *hdr_shot;

	initargs(argc, argv);
	requestdoc(1);

	if (!getparstring("fin", &fin)) fin = NULL;
	if (!getparstring("fshot", &fshot)) fshot = NULL;
    if (!getparstring("fout", &fout)) fout = "out.su";
	if (!getparint("zmax", &zmax)) zmax = 0;
	if (!getparint("inx", &inx)) inx = 0;
	if (!getparfloat("zrcv", &f1)) f1 = 0;
    if (!getparfloat("xrcv", &f2)) f2 = 0;
	if (!getparfloat("dzrcv", &dzrcv)) dzrcv = -1;
    if (!getparfloat("dxrcv", &dxrcv)) dxrcv = -1;
	if (!getparfloat("rho", &rho)) rho=1000.0;
	if (!getparint("numb", &numb)) numb=0;
    if (!getparint("dnumb", &dnumb)) dnumb=1;
	if (!getparint("tolset", &tolset)) tolset=10;
	if (!getparint("mode", &mode)) mode=0;
	if (!getparint("ntap", &ntap)) ntap=0;
	if (!getparint("maxoffset", &maxoffset)) maxoffset=0;
	if (fin == NULL) verr("Incorrect f2 input");
	if (fshot == NULL) verr("Incorrect Green input");

	if (dnumb == 0) dnumb = 1;

	sprintf(fins,"z%d",numb);

	ptr  = strstr(fin,fins);
	pos1 = ptr - fin + 1;

   	sprintf(fbegin,"%*.*s", pos1-1, pos1-1, fin);
   	sprintf(fend,"%s", fin+pos1+1);

	file_det = 1;
	zcount=0;
	nzs=0;

	while (file_det) {
        sprintf(fins,"z%d",nzs*dnumb+numb);
        sprintf(fin,"%s%s%s",fbegin,fins,fend);
        fp_in = fopen(fin, "r");
        if (fp_in == NULL) {
            if (nzs == 0) {
                verr("error on opening basefile=%s", fin);
            }
            else if (nzs == 1) {
                vmess("1 file detected");
				file_det = 0;
         		break;
            }
            else {
                vmess("%d files detected",nzs);
                file_det = 0;
                break;
            }
        }
        fclose(fp_in);
        nzs++;
    }

	if (inx < 1) { 
		inx = 1;
	}

	if (zmax < 1) zmax=1;
	if (zmax < nzs) nzs=zmax;

	nxs = inx;
	count=0;
	npos = nxs*nzs;

	vmess("nxs: %d, nzs: %d",nxs,nzs);

	nshots = 0;
    getFileInfo(fshot, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax1, &sclsxgx, &ntraces);

	if (dxrcv < 0) dxrcv=dx;
	if (dzrcv < 0) dzrcv=dx;

	// ngath zijn het aantal schoten
	shotdata	= (float *)malloc(nt*nx*nshots*sizeof(float));
	hdr_shot	= (segy *)calloc(nx*nshots,sizeof(segy));

	fp_shot = fopen(fshot,"r");
	if (fp_shot == NULL) {
		verr("Could not open file");
	}
	vmess("nt: %d nx: %d nshots: %d",nt,nx,nshots);
	//nx = readData(fp_shot, shotdata, hdr_shot, nt);
	fclose(fp_shot);
	readSnapData(fshot, &shotdata[0], &hdr_shot[0], 1, nx, nt, 0, nx, 0, nt);

	hdr_out     = (segy *)calloc(nxs,sizeof(segy));	
	Ghom		= (float *)malloc(nt*npos*sizeof(float));
	ht			= (int)ceil(nt/2);
	nw 			= ht+1;
	dw			= 2.0*(M_PI)/(dt*nt);
	cshot		= (complex *)malloc(nw*nx*sizeof(complex));
	tol			= (int *)malloc(nxs*sizeof(float));
	taper       = (float *)malloc(nx*sizeof(float));

    /*for (ix=0; ix<nx; ix++) {
        taper[ix] = 1.0;
    }
    if (ntap > 0) {//Create taper
		vmess("Tapering of %d points applied at edges",ntap);
        for (ix=0; ix<ntap; ix++) {
            taper[ix] = (cos((M_PI)*(ix-ntap)/ntap)+1)/2.0;
            taper[nx-1-ix] = (cos((M_PI)*(ix-ntap)/ntap)+1)/2.0;
        }
    }*/

	for (ix = 0; ix < nx; ix++) {
		/*for (it=0; it<nt; it++) {
			shotdata[ix*nt+it] *= taper[ix];
		}*/
		if (hdr_shot[ix].scalco < 0) {
			offset_tmp = (hdr_shot[ix].sx-hdr_shot[ix].gx)/-hdr_shot[ix].scalco;
		}
		else if (hdr_shot[ix].scalco == 0) {
			offset_tmp = hdr_shot[ix].sx-hdr_shot[ix].gx;
		}
		else {
			offset_tmp = (hdr_shot[ix].sx-hdr_shot[ix].gx)*hdr_shot[ix].scalco;
		}
		if (maxoffset > 0 ) {
			if (abs(offset_tmp) > maxoffset) {
				for (it=0;it<nt;it++) {
					shotdata[ix*nt+it] = 0.0;
				}
				vmess("Removed offset:%d",offset_tmp);
			}
		}
		rc1fft(&shotdata[ix*nt],&cshot[ix*nw],nt,-1);
	}

#pragma omp parallel default(shared) \
  private(offset,ctrace,rtrace,chom,indrcv,rmt,ix,it,is) \
  private(indata, hdr_in,fins,fin2,fp_in,offset_tmp)
{
	chom		= (complex *)calloc(nw,sizeof(complex));
	ctrace		= (complex *)malloc(nw*sizeof(complex));
    rtrace		= (float *)malloc(nt*sizeof(float));
	indata		= (float *)malloc(nt*nx*nxs*sizeof(float));
    hdr_in 		= (segy *)calloc(nx*nxs,sizeof(segy));
#pragma omp for 
	for (ir = 0; ir < nzs; ir++) {
        sprintf(fins,"z%d",ir*dnumb+numb);
		sprintf(fin2,"%s%s%s",fbegin,fins,fend);
        fp_in = fopen(fin2, "r");
		if (fp_in == NULL) {
			verr("Danger Will Robinson");
		}
		fclose(fp_in);
		readSnapData(fin2, &indata[0], &hdr_in[0], nxs, nx, nt, 0, nx, 0, nt);
		for (is = 0; is < nxs; is++) {
			for (ix = 0; ix < nx; ix++) {
				/*for (it=0; it<nt; it++) {
            		indata[is*nt*nx+ix*nt+it] *= taper[ix];
        		}*/
				if (hdr_in[is*nx+ix].scalco < 0) {
            		offset_tmp = (hdr_in[is*nx+ix].sx-hdr_in[is*nx+ix].gx)/-hdr_in[is*nx+ix].scalco;
        		}
        		else if (hdr_in[is*nx+ix].scalco == 0) {
            		offset_tmp = hdr_in[is*nx+ix].sx-hdr_in[is*nx+ix].gx;
        		}
        		else {
            		offset_tmp = (hdr_in[is*nx+ix].sx-hdr_in[is*nx+ix].gx)*hdr_in[is*nx+ix].scalco;
        		}
        		if (maxoffset > 0 ) {
            		if (abs(offset_tmp) > maxoffset) {
                		for (it=0;it<nt;it++) {
                    		indata[is*nt*nx+ix*nt+it] = 0.0;
                		}
                	//vmess("Removed offset:%d",offset_tmp);
            		}
				}
				rc1fft(&indata[is*nt*nx+ix*nt],ctrace,nt,-1);
				if (mode==0) { //Single source
					for (it = 1; it < nw; it++) {
						//chom[it].r +=  (4/(rho*dw*it))*(ctrace[it].r*cshot[ix*nw+it].i + ctrace[it].i*cshot[ix*nw+it].r);
						//chom[it].r +=  (4/(rho*dw*it))*2*(ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i);
						chom[it].r +=  2*(ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i);
						//chom[it].r +=  ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i;
						//chom[it].r +=  ctrace[it].r*cshot[ix*nw+it].i + ctrace[it].i*cshot[ix*nw+it].r;
					}
				}
				else { //Multiple sources
                	for (it = 0; it < nw; it++) {
                        /*chom[it].r -=  (2/(rho*dw*it))*(ctrace[it].r*cshot[ix*nw+it].i + ctrace[it].i*cshot[ix*nw+it].r);
						chom[it].i +=  (2/(rho*dw*it))*(ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i);*/
						chom[it].r -=  (ctrace[it].r*cshot[ix*nw+it].i + ctrace[it].i*cshot[ix*nw+it].r);
                        chom[it].i +=  (ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i);
                    }
				}
			}
			cr1fft(&chom[0],rtrace,nt,1);
			indrcv = 0;
            rmt = MIN(nt-indrcv,indrcv)-shift;
			for (it = 0; it < ht; it++) {
				if (it > ht-rmt) {
					Ghom[it*nxs*nzs+is*nzs+ir] = 0.0;
				}
				else {
					Ghom[it*nxs*nzs+is*nzs+ir] = rtrace[ht+it];
				}
			}
			for (it = ht; it < nt; it++) {
				if (it < ht+rmt) {
					Ghom[it*nxs*nzs+is*nzs+ir] = 0.0;
				}
				else {
					Ghom[it*nxs*nzs+is*nzs+ir] = rtrace[it-ht];
				}
            }
		memset(&chom[0].r, 0, nw*2*sizeof(float));
		}
		//vmess("Creating Homogeneous Green's function at depth %d from %d depths",ir+1,nzs);
		count+=1;
		vmess("Creating Homogeneous Green's function at depth %d from %d depths",count,nzs);
	}
    free(chom);free(ctrace);free(rtrace);
	free(indata);free(hdr_in);
}
	free(shotdata);

	vmess("nxs: %d nxz: %d f1: %.7f",nxs,nzs,f1);

	fp_out = fopen(fout, "w+");
	
	for (ir	= 0; ir < nt; ir++) {
		for (ix = 0; ix < nxs; ix++) {
            	hdr_out[ix].fldr	= ir+1;
            	hdr_out[ix].tracl	= ir*nxs+ix+1;
            	hdr_out[ix].tracf	= ix+1;
				hdr_out[ix].scalco  = hdr_shot[0].scalco;
    			hdr_out[ix].scalel	= hdr_shot[0].scalel;
				hdr_out[ix].sdepth	= hdr_shot[0].sdepth;
				hdr_out[ix].trid	= 1;
				hdr_out[ix].ns		= nzs;
				hdr_out[ix].trwf	= nxs;
				hdr_out[ix].ntr		= hdr_out[ix].fldr*hdr_out[ix].trwf;
				hdr_out[ix].f1		= f1;
				hdr_out[ix].f2		= f2/1000;
				hdr_out[ix].dt      = dt*(1E6);
				hdr_out[ix].d1      = dzrcv;
            	hdr_out[ix].d2      = dxrcv;
				hdr_out[ix].sx      = hdr_shot[0].sx;
				hdr_out[ix].gx      = (int)roundf(f2 + (ix*hdr_out[ix].d2)*1000.0);
            	hdr_out[ix].offset	= (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
		}
		ret = writeData(fp_out, &Ghom[ir*nxs*nzs], hdr_out, nzs, nxs);
		if (ret < 0 ) verr("error on writing output file.");
	}
	
	fclose(fp_out);
	return 0;
}

void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift)
{
    int     i, j, n, optn, nfreq, sign;
    float   df, dw, om, tau, scl;
    float   *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
    complex *cdata1, *cdata2, *ccon, tmp;

    optn = optncr(nsam);
    nfreq = optn/2+1;


    cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
    if (cdata1 == NULL) verr("memory allocation error for cdata1");
    cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
    if (cdata2 == NULL) verr("memory allocation error for cdata2");
    ccon = (complex *)malloc(nfreq*nrec*sizeof(complex));
    if (ccon == NULL) verr("memory allocation error for ccov");

    rdata1 = (float *)malloc(optn*nrec*sizeof(float));
    if (rdata1 == NULL) verr("memory allocation error for rdata1");
    rdata2 = (float *)malloc(optn*nrec*sizeof(float));
    if (rdata2 == NULL) verr("memory allocation error for rdata2");

    /* pad zeroes until Fourier length is reached */
    pad_data(data1, nsam, nrec, optn, rdata1);
    pad_data(data2, nsam, nrec, optn, rdata2);

    /* forward time-frequency FFT */
    sign = -1;
    rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
    rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);

    /* apply convolution */
    p1r = (float *) &cdata1[0];
    p2r = (float *) &cdata2[0];
    qr = (float *) &ccon[0].r;
    p1i = p1r + 1;
    p2i = p2r + 1;
    qi = qr + 1;
    n = nrec*nfreq;
    for (j = 0; j < n; j++) {
        *qr = (*p2r**p1r-*p2i**p1i);
        *qi = (*p2r**p1i+*p2i**p1r);
        qr += 2;
        qi += 2;
        p1r += 2;
        p1i += 2;
        p2r += 2;
        p2i += 2;
    }
    free(cdata1);
    free(cdata2);

    if (shift) {
        df = 1.0/(dt*optn);
        dw = 2*PI*df;
        tau = dt*(nsam/2);
        for (j = 0; j < nrec; j++) {
            om = 0.0;
            for (i = 0; i < nfreq; i++) {
                tmp.r = ccon[j*nfreq+i].r*cos(om*tau) + ccon[j*nfreq+i].i*sin(om*tau);
                tmp.i = ccon[j*nfreq+i].i*cos(om*tau) - ccon[j*nfreq+i].r*sin(om*tau);
                ccon[j*nfreq+i] = tmp;
                om += dw;
            }
        }
    }

        /* inverse frequency-time FFT and scale result */
    sign = 1;
    scl = 1.0/((float)(optn));
    crmfft(&ccon[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
    scl_data(rdata1,optn,nrec,scl,con,nsam);

    free(ccon);
    free(rdata1);
    free(rdata2);
    return;
}

void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout)
{
    int it,ix;
    for (ix=0;ix<nrec;ix++) {
       for (it=0;it<nsam;it++)
        datout[ix*nsamout+it]=data[ix*nsam+it];
       for (it=nsam;it<nsamout;it++)
        datout[ix*nsamout+it]=0.0;
    }
}

void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout)
{
    int it,ix;
    for (ix = 0; ix < nrec; ix++) {
        for (it = 0 ; it < nsamout ; it++)
            datout[ix*nsamout+it] = scl*data[ix*nsam+it];
    }
}