/****************************************************************                         2dcmp_presdm                      ****                         =============                     ****                         Cheng Jiubing                     ****                        2001.01-2002.02                    ****                                                           ****  Copyright (c) TongJi University, ShangHai, China, 1999.  ****                   All rights reserved.                    ****************************************************************/#include "fft.h"#include "mute_direct.h"#include "cjbsegy.h"#include "alloc.h"#include "mpi.h"/***** SEG-Y header *********/       Y_3200   y3200;       bhed     head_400;       cjbsegy  tr, vtr;/***** SEG-Y header *********/void keyin_par(char *iwfile, char *ivfile, char *owfile1, char *owfile2,      char *owfile3, int *seismic, int *operator, int *nxv, int *nzv,     int *ncdp0, float *dmx, float *dhx,  float *depth, float *dz,      float *f1, float *f2, float *f3, float *f4, int *iflag_fft,     float *dxv, float *dzv, float *xv0, int *ncdpstep_cig,     float *gamamin, float *gamamax, int *ngama);void read_par(char *iwfile, char *ivfile, char *owfile1, char *owfile2,       char *owfile3, int *seismic, int *operator, int *nxv, int *nzv,      int *ncdp0, float *dmx, float *dhx,  float *depth, float *dz,      float *f1, float *f2, float *f3, float *f4, int *iflag_fft,      float *dxv, float *dzv, float *xv0, int *ncdpstep_cig,     float *gamamin, float *gamamax, int *ngama);void scan_statistics(FILE *iwfp, int nt, float *xmmin, float *xmmax,      float *offsetmin, float *offsetmax, int *ntrace);void Input_FFT_Transp(FILE *iwfp, FILE *fp, int nmx, int nhx, int ntrace,      int nt, int ntfft, int nw, int nf1, int nf2, int nf3, int nf4,      float xmmin, float offsetmin, float dmx, float dhx, float dt, float *tapert);void velocity(float **vv, float *vz, float ***dss, float ***dsg, int nxv,      int nzv, int nmx, int nhx, int nz, float xmmin, float offsetmin,     float xv0, float dxv, float dzv, float dz, float dmx, float dhx);void filter(complex *cq,int nw,int nf1,int nf2,int nf3,int nf4);void presdm(complex **wave, float ***cimage, float **image, float ***dss,       float ***dsg, float *vz, float *taperx, float *taperh,     int nmx, int nhx, int nz, float dmx, float dhx, float dz, float w,     int operator, float xmmin, float offsetmin, int ncig, int ncdpstep_cig, int myid);/****** split-step fourier operator ******/void ssf(complex **wave, float ***cimage, float **image, float ***dss,      float ***dsg, float *vz, float *taperx, float *taperh, int nmx,     int nhx, int iz, float dmx, float dhx, float dz, float w,      float xmmin, float offsetmin, int ncig, int ncdpstep_cig, int myid);/****** quasi-linear fourier/born operator ******/void LQELBF(complex **wave, float ***cimage, float **image, float ***dss,      float ***dsg, float *vz, float *taperx, float *taperh, int nmx,     int nhx, int iz, float dmx, float dhx, float dz, float w,      float xmmin, float offsetmin, int ncig, int ncdpstep_cig, int myid);void Stack_ImageCIGs_allnodes(char *owfile1, char *owfile2, float ***cimage,      float **image, int seismic, int nmx, int nhx, int nz, int ncdp0,     float offsetmin, float dhx, float xmmin, float dmx, float dz,      int ncig, int ncdpstep_cig, int myid);/* taper boundary processing */void tapfunc(float *taper, int nx, int nxleft, int nxright, int iflag);/* delay depth between MPI nodes to avoid I/O traffic  */void myid_node_delay(int myid);void memory_used(int nxv, int nmx, int nhx, int nz, int ncig, int operator);/* estimate the memory need for this program */void odcig2adcig_gama(char *owfile2, char *owfile3, int seismic,      int nmx, int nhx, int ngama, int nz, int ncdp0,      int ncig, int ncdpstep_cig, float xmmin, float dmx,     float offsetmin, float dhx, float fgama, float dgama, float dz);void fwd_FK_sstack_cjb (float dt, int nt, int nx, float xmin, float dx, int ngama,	float fgama, float dgama, float fmin, float **traces, float **out_traces);void xindex (int nx, float ax[], float x, int *index);void intlin (int nin, float xin[], float yin[], float yinl, float yinr, 	int nout, float xout[], float yout[]);main ( argc, argv)int  argc;char **argv;{        /************************ Keyin parameters *******************************/        char      iwfile[70];                /* input pre-stack gather */        char      ivfile[70];                /* input velocity file */        char      owfile1[70];               /* output image by "t=0, h=0" before outputing CIGs*/        char      owfile2[70];               /* output image-point gather*/        char      owfile3[70];               /* output image by stacking CIGs*/        int       kflag;                           int       seismic;                   /* segy data: seismic=1; su data: seismic=2*/        int       operator;                  /* imaging operator type*/        int       nxv;                       /* velocity lateral gridpoint number*/        int       nzv;                       /* velocity vertical gridpoint number */         int       ncdp0;                     /* wavefield first cdp number */         float     dmx;                       /* midpoint interval */        float     dhx;                       /* half-offset interval */        float     depth;                     /* extrapolate depth */        float     dz;                        /* extrapolate depth interval*/        float     f1;                        /* f1 of frequency range: f1-f2--f3-f4*/        float     f2;                        /* f2 of frequency range: f1-f2--f3-f4*/        float     f3;                        /* f3 of frequency range: f1-f2--f3-f4*/        float     f4;                        /* f4 of frequency range: f1-f2--f3-f4*/        float     dxv;                       /* velocity lateral grid interval */         float     dzv;                       /* velocity vertical grid interval */         float     xv0;                       /* velocity lateral starting coordinate */	int       ncdpstep_cig;              /* CIGs cdp number interval*/	float     gamamin;                   /* minimin incidence angle for CIGs*/	float     gamamax;                   /* minimax incidence angle for CIGs*/        int       ngama=81;                  /* Number of Gama in CIGs */        /************************ Keyin parameters *******************************/        int      ncig, np, myid, num, ierr, i, ip, nnp, iww, mpinode;        int      ntrace, nt, nz, nf1, nf2, nf3, nf4, ntfft, nmx, nxvv,nhx, nkh, nw, nww;        int      ix, ixv, imx, iz, it, iw, ikx, ipp, iflag_fft, unit=1;	float    dt, w, kx, dw, dkx, fw, fkx, zero, dkh;          float    xmmin, xmmax, xsgmin, xsgmax, offsetmin, offsetmax;        float    fgama, dgama;	float    *taperx, *taperh, *tapert, **vv, **v, *vz, **image, ***dss,***dsg;        complex  **wave;        float    ***cimage;	FILE     *fp, *fpv, *ivfp, *iwfp;               char     *tmpfile;               /*temporay frequency domain wavefields*/	MPI_Status status;        MPI_Init(&argc, &argv);        MPI_Comm_size(MPI_COMM_WORLD, &np);        MPI_Comm_rank(MPI_COMM_WORLD, &myid);        if(myid==0){        printf("===================================================\n");        printf("                                                   \n");        printf("             2D Midpoint-Offset Domain             \n");        printf("              PreStack Depth Migration             \n");        printf("                                                   \n");        printf("              -------------------------            \n");        printf("                  SSF,Fourier-Born                 \n");        printf("              -------------------------            \n");        printf("                                                   \n");        printf("            Tongji University,Shanghai,China       \n");        printf("                    Cheng Jiubing                  \n");        printf("===================================================\n");        printf("Notices:                                           \n");        printf("  1. The tool can process single-side shotting,    \n");        printf(" and for the mid-shotting data, we only use one    \n");        printf(" side of the data.                                 \n");        printf("  2. Input wavefields can be any prestack data     \n");        printf(" with or without midpoint-offset gather sorting.   \n");        printf("  3. Input velocity data need SU or SEG-Y header,  \n");        printf(" but the starting point coordinate, grid-interval, \n");        printf(" and grid number must be given.                    \n");         printf("===================================================\n");        printf("Output:                                            \n");        printf("1. Image result by 't=0,h=0' before outputing CIGs \n");        printf("2. Ph-CIGs by 't=0' with 7-point SINC Interpolation\n");        printf("3. New Image results stacked by muted CIGs         \n");        printf("===================================================\n");        }        kflag=2;        if(kflag==1&&np==1){          keyin_par(iwfile, ivfile, owfile1, owfile2, owfile3,                     &seismic, &operator, &nxv, &nzv, &ncdp0,                     &dmx, &dhx, &depth, &dz, &f1, &f2, &f3, &f4,                    &iflag_fft, &dxv, &dzv, &xv0, 		    &ncdpstep_cig, &gamamin, &gamamax, &ngama);        }else{          read_par(iwfile, ivfile, owfile1, owfile2, owfile3,                    &seismic, &operator, &nxv, &nzv, &ncdp0,                    &dmx, &dhx, &depth, &dz, &f1, &f2, &f3, &f4,                   &iflag_fft, &dxv, &dzv, &xv0, 		   &ncdpstep_cig, &gamamin, &gamamax, &ngama);        }        if(myid==0){        printf("iwfile= %s\n",iwfile);        printf("ivfile= %s\n",ivfile);        printf("owfile1= %s\n",owfile1);        printf("owfile2= %s\n",owfile2);        printf("owfile3= %s\n",owfile3);        printf("seismic= %d\n",seismic);        printf("operator= %d\n",operator);        printf("nxv, nzv: %d %d\n",nxv, nzv);        printf("ncdp0= %d\n",ncdp0);        printf("dmx, dhx: %f %f\n",dmx, dhx);        printf("dz= %f\n",dz);        printf("f1, f2, f3, f4: %f %f %f %f\n",f1, f2, f3, f4);        printf("iflag_fft= %d\n",iflag_fft);        printf("dxv, dzv: %f %f\n",dxv, dzv);        printf("xv0= %f\n",xv0);	printf("ncdpstep_cig= %d\n",ncdpstep_cig);        printf("gamamin= %f gamamax=%f\n",gamamin,gamamax);        printf("ngama= %d\n",ngama);        }        nz=(depth+0.0005)/dz;        /*#1: read trace point number and depth sampling rate from seismic header */         if((iwfp = fopen(iwfile,"rb"))==NULL)        { printf("Open iwfile error !\n"); exit(1); }        if(seismic==1){           fread(&y3200, sizeof(Y_3200), 1, iwfp);           fread(&head_400, sizeof(bhed), 1, iwfp);           fread(&tr, sizeof(cjbsegy), 1, iwfp);           fseek(iwfp, (long)(-sizeof(cjbsegy)), 1);        }        else{           fread(&tr, sizeof(cjbsegy), 1, iwfp);           rewind(iwfp);        }        nt = tr.ns;        dt=tr.dt/1000000.0;	if(myid==0){        printf("***************************************************\n");        printf("Seismic Header Information: nt= %d, dt= %f(s)\n",nt,dt);	}        /*#2: determine FFT/frequency parameters */         ntfft = npfa(nt);        nww = ntfft/2+1;        dw = 2.0*PI/(ntfft*dt);        fw=2.0*PI*f1;        nf1=fw/dw+0.5;        fw=2.0*PI*f2;        nf2=fw/dw+0.5;        fw=2.0*PI*f3;        nf3=fw/dw+0.5;        fw=2.0*PI*f4;        nf4=fw/dw+0.5;         nw=nf4-nf1+1;        fw=nf1*dw;        vv=alloc2float(nzv, nxv);        zero2float(vv, nzv, nxv);                if((ivfp=fopen(ivfile,"rb"))==NULL)          { printf("Open ivfile error !\n"); exit(1);}        for(ix=0;ix<nxv;ix++){            fread(vv[ix], FSIZE, nzv, ivfp);	  // for(iz=0;iz<nzv;iz++) printf("%f ",vv[ix][iz]);	}        fclose(ivfp);        if(myid==0){           printf("***************************************************\n");          printf("*             Begin to Scan SEG-Y Header          *\n");          printf("***************************************************\n");          scan_statistics(iwfp, nt, &xmmin, &xmmax, &offsetmin, &offsetmax, &ntrace);          printf("@@@@: xmmin= %f, xmmax= %f\n",xmmin, xmmax);          printf("@@@@: offsetmin= %f,  offsetmax= %f\n",offsetmin, offsetmax);        }        MPI_Bcast(&xmmin, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);        MPI_Bcast(&xmmax, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);        MPI_Bcast(&offsetmin, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);        MPI_Bcast(&offsetmax, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);        MPI_Bcast(&ntrace, 1, MPI_INT, 0, MPI_COMM_WORLD);	offsetmax=fabs(offsetmin)>fabs(offsetmax)?fabs(offsetmin):fabs(offsetmax);	offsetmin=-offsetmax;        nhx=(offsetmax-offsetmin)/(2*dhx)+1;        nmx=(xmmax-xmmin)/dmx+1;         nkh= npfa(nhx);        dkh= 2.0*PI/(nkh*dhx);        /*** Wavefield is ranged in order of fabs(offset), so we use CIGs with Ph>=0.0 ***/        if(myid==0){          printf("@@@@: nw= %d  nmx= %d   nhx= %d\n",nw, nmx, nhx);        }        tmpfile="Temp_FFT_Transp_DepthMig";        if(!iflag_fft&&myid==0){          printf("***************************************************\n");          printf("               Input-FFT-Transp.                   \n");          printf("***************************************************\n");          fp=fopen(tmpfile, "wb");          fseek(iwfp, (long)(-ntrace*(sizeof(cjbsegy)+nt*FSIZE)), 1);	  tapert=alloc1float(nt);	  tapfunc(tapert, nt, 15, 50, +1);           Input_FFT_Transp(iwfp, fp, nmx, nhx, ntrace, nt, ntfft, nw,                            nf1, nf2, nf3, nf4, xmmin, offsetmin,                            dmx, dhx, dt, tapert);	  free1float(tapert);