/*if the horizontal spacing interval is in feet, convert it to meter*/        if(!flag)        dx*=0.3048;	/* loops over depth */	for(iz=0;iz<nz;++iz){        /*the imaging condition*/        for(ix=0;ix<nx;ix++){        for(iw=0,w=fw;iw<nw;w+=dw,iw++){                   complex tmp;                float ratio=10.0;                                if(fabs(ix+ix2-ixshot)*dx<ratio*iz*dz)                tmp=cmul(cp[iw][ix],cp1[iw][ix]);                else tmp=cmplx(0.0,0.0);                  cresult[ix+ix2][iz]+=tmp.r/ntfft;        }        }/* anothe imaging condition, slightly different from the above one, but quiteslow*/        /*        for(iw=0,w=fw;iw<nw;w+=dw,iw++){                float kk=0.0;                complex tmp;                float ratio=1.5;                 if(dip<80)ratio=1.5;                else ratio=1.5;                        for(ix=0;ix<nx;ix++){                     kk+=(pow(cp1[iw][ix].i,2.0)+pow(cp1[iw][ix].r,2.0))/nx;                }                         for(ix=0;ix<nx;ix++){                tmp=cmul(cp[iw][ix],cp1[iw][ix]);                if(fabs(ix-ixshot)*dx<ratio*iz*dz||ixshot-ix<0 )                tmp=crmul(tmp,1.0/(kk+1.0e-10));                  else tmp=cmplx(0.0,0.0);                                cresult[ix+ix2][iz]+=tmp.r/ntfft;                        }                }*/                                      /*get the average velocity*/                 v1=0.0;                for(ix=0;ix<nx;++ix)                {v1+=v[iz][ix]/nx;}                /*compute time-invariant wavefield*/                for(ix=0;ix<nx;++ix)                for(iw=0,w=fw;iw<nw;w+=dw,++iw) {                        kz2=-(1.0/v1)*w*dz;                        cshift2=cmplx(cos(kz2),sin(kz2));                        cp[iw][ix]=cmul(cp[iw][ix],cshift2);                        cp1[iw][ix]=cmul(cp1[iw][ix],cshift2);                }                /*wave-propagation using finite-difference method*/                fdmig( cp, nx, nw,v[iz],fw,dw,dz,dx,dt,dip);                fdmig( cp1,nx, nw,v[iz],fw,dw,dz,dx,dt,dip);                /*apply thin lens term here*/                for(ix=0;ix<nx;++ix)                for(iw=0,w=fw;iw<nw;w+=dw,++iw){                        kz2=-(1.0/v[iz][ix]-1.0/v1)*w*dz;                        cshift2=cmplx(cos(kz2),sin(kz2));                        cp[iw][ix]=cmul(cp[iw][ix],cshift2);                        cp1[iw][ix]=cmul(cp1[iw][ix],cshift2);                }}free2complex(cp);free2complex(cp1);free2float(v);nxshot--;/* time(&t2);warn("\n %d nxshot has been finished in %f seconds",nxshot,difftime(t2,t1));*/if(nxshot)goto loop;	/* restore header fields and write output */	for(ix=0; ix<nxo; ix++){                tr.ns = nz ;                tr.dt = dz*1000000.0 ;                tr.d2 = dx;                tr.offset = 0;                 tr.cdp = tr.tracl = ix;		memcpy( (void *) tr.data, (const void *) cresult[ix],nz*FSIZE);		puttr(&tr);	}		return(CWP_Exit());	}float * ricker(float Freq,float dt,int *Npoint){int i;/* they are the dummy counter*/float Bpar,t,u,*Amp;int Np1,N;        if(Freq==0.0)Freq=30.0;if(dt==0.0)dt=0.004;Bpar=sqrt(6.0)/(PI*Freq);N=ceil(1.35*Bpar/dt);Np1=N;*Npoint=2*N+1;         Amp=alloc1float(*Npoint);        Amp[Np1]=1.0;  for(i=1;i<=N;i++){t=dt*(float)i;u=2.0*sqrt(6.0)*t/Bpar;Amp[Np1+i]=Amp[Np1-i]=0.5*(2.0-u*u)*exp(-u*u/4.0);}return Amp;}void fdmig( complex **cp, int nx, int nw, float *v,float fw,float        dw,float dz,float dx,float dt,int dip){	int iw,ix,step=1;	float *s1,*s2,w,coefa[5],coefb[5],v1,vn,trick=0.1;	complex cp2,cp3,cpnm1,cpnm2;        complex a1,a2,b1,b2;        complex endl,endr;        complex *data,*d,*a,*b,*c;        s1=alloc1float(nx);        s2=alloc1float(nx);        data=alloc1complex(nx);        d=alloc1complex(nx);        a=alloc1complex(nx);        b=alloc1complex(nx);        c=alloc1complex(nx);        if(dip==45){        coefa[0]=0.5;coefb[0]=0.25;         step=1;        }                if(dip==65){        coefa[0]=0.478242060;coefb[0]=0.376369527;        step=1;        }                if(dip==79){        coefa[0]=coefb[0]=0.4575;        step=1;        }        if(dip==80){        coefa[1]=0.040315157;coefb[1]=0.873981642;        coefa[0]=0.457289566;coefb[0]=0.222691983;        step=2;        }                if(dip==87){        coefa[2]=0.00421042;coefb[2]=0.972926132;        coefa[1]=0.081312882;coefb[1]=0.744418059;        coefa[0]=0.414236605;coefb[0]=0.150843924;        step=3;        }                 if(dip==89){        coefa[3]=0.000523275;coefb[3]=0.994065088;        coefa[2]=0.014853510;coefb[2]=0.919432661;        coefa[1]=0.117592008;coefb[1]=0.614520676;        coefa[0]=0.367013245;coefb[0]=0.105756624;        step=4;        }        if(dip==90){        coefa[4]=0.000153427;coefb[4]=0.997370236;        coefa[3]=0.004172967;coefb[3]=0.964827992;        coefa[2]=0.033860918;coefb[2]=0.824918565;        coefa[1]=0.143798076;coefb[1]=0.483340757;        coefa[0]=0.318013812;coefb[0]=0.073588213;        step=5;	}        v1=v[0];vn=v[nx-1];loop:         step--;        for(iw=0,w=fw;iw<nw;iw++,w+=dw){                if(fabs(w)<=1.0e-10)w=1.0e-10/dt;                 for(ix=0;ix<nx;ix++){                        s1[ix]=(v[ix]*v[ix])*coefb[step]/(dx*dx*w*w)+trick;                        s2[ix]=-v[ix]*dz*coefa[step]/(w*dx*dx)*0.5;		}                for(ix=0;ix<nx;ix++){                        data[ix]=cp[iw][ix];                }                cp2=data[1];                cp3=data[2];                cpnm1=data[nx-2];                cpnm2=data[nx-3];                a1=crmul(cmul(cp2,conjg(cp3)),2.0);                b1=cadd(cmul(cp2,conjg(cp2)),cmul(cp3,conjg(cp3)));                if(b1.r==0.0 && b1.i==0.0)                        a1=cexp(cmplx(0.0,-w*dx*0.5/v1));                else                        a1=cdiv(a1,b1);                if(a1.i>0.0)a1=cexp(cmplx(0.0,-w*dx*0.5/v1));                a2=crmul(cmul(cpnm1,conjg(cpnm2)),2.0);                b2=cadd(cmul(cpnm1,conjg(cpnm1)),cmul(cpnm2,conjg(cpnm2)));                if(b2.r==0.0 && b2.i==0.0)                        a2=cexp(cmplx(0.0,-w*dx*0.5/vn));                else                        a2=cdiv(a2,b2);                if(a2.i>0.0)a2=cexp(cmplx(0.0,-w*dx*0.5/vn));                for(ix=0;ix<nx;ix++){                        a[ix]=cmplx(s1[ix],s2[ix]);                        b[ix]=cmplx(1.0-2.0*s1[ix],-2.0*s2[ix]);                }                for(ix=1;ix<nx-1;ix++){		d[ix]=cadd(cadd(cmul(data[ix+1],a[ix+1]),cmul(data[ix-1],a[ix-1])),		cmul(data[ix],b[ix]));                }                d[0]=cadd(cmul(cadd(b[0],cmul(a[0],a1)),data[0]),cmul(data[1],a[1]));		d[nx-1]=cadd(cmul(cadd(b[nx-1],cmul(a[nx-1],a2)),data[nx-1]),		cmul(data[nx-2],a[nx-2]));                for(ix=0;ix<nx;ix++){                        data[ix]=cmplx(s1[ix],-s2[ix]);                        b[ix]=cmplx(1.0-2.0*s1[ix],2.0*s2[ix]);                }                endl=cadd(b[0],cmul(data[0],a1));                endr=cadd(b[nx-1],cmul(data[nx-1],a2));                                for(ix=1;ix<nx-1;ix++){                        a[ix]=data[ix+1];                        c[ix]=data[ix-1];                }                a[0]=data[1];                c[nx-1]=data[nx-2];                                        retris(data,a,c,b,endl,endr,nx,d);                for(ix=0;ix<nx;ix++){                        cp[iw][ix]=data[ix];                }        }if(step) goto loop;        free1complex(data);        free1complex(d);        free1complex(b);        free1complex(c);        free1complex(a);        free1float(s1);        free1float(s2);                        return;}                 void retris(complex *data,complex *a,complex *c, complex *b,                complex endl,complex endr, int nx, complex *d){                         int ix;        complex *e,den;        complex *f;        e=alloc1complex(nx);        f=alloc1complex(nx);        e[0]=cdiv(cneg(a[0]),endl);        f[0]=cdiv(d[0],endl);        for(ix=1;ix<nx-1;++ix){                den=cadd(b[ix],cmul(c[ix],e[ix-1]));                e[ix]=cdiv(cneg(a[ix]),den);                f[ix]=cdiv(csub(d[ix],cmul(f[ix-1],c[ix])),den);        }                 	data[nx-1]=cdiv(csub(d[nx-1],cmul(f[nx-2],c[nx-2])),cadd(endr,cmul(c[nx-2],e[nx-2])));                        for(ix=nx-2;ix>-1;--ix)        data[ix]=cadd(cmul(data[ix+1],e[ix]),f[ix]);        free1complex(e);        free1complex(f);        return;  }