地震波克西霍夫叠后偏移程序源代码

基本信息

源码名称：地震波克西霍夫叠后偏移程序源代码

源码大小：0.01M

文件格式：.cpp

开发语言：C/C++

更新时间：2015-08-17

友情提示：（无需注册或充值，赞助后即可获取资源下载链接）

嘿，亲！知识可是无价之宝呢，但咱这精心整理的资料也耗费了不少心血呀。小小地破费一下，绝对物超所值哦！如有下载和支付问题，请联系我们QQ(微信同号)：813200300

本次赞助数额为： 2 元　

源码介绍

地震波克西霍夫叠后偏移程序源代码

main(int argc, char* argv[])
{
/*int myrank, numprocs;
int namelen;
char processor_name[MPI_MAX_PROCESSOR_NAME];

    MPI_Init(&argc,&argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
    MPI_Get_processor_name(processor_name,&namelen);*/

int b1, b2, b3, nz, ny, nx,n, esize, vel, i,j,k,kk,order;
float sz, sy, sx, o1, o2, o3, dz, dy, dx, vel0, br1, br2, br3,vel1,vel2,vel3,vel4;
float *v;

int sample,nshot,ngx,ngy,ii,jj,num;
float t,dt,ds,dgx,dgy,*danp,*dan,*AMP,*ts,*tg;

nshot=100;   //�ڵ��
ds=80;       //�ڵ��
sample=1200; //ÿ��Ĳ��
dt=1E-3;     //��
ngx=30;      //x��첨��ĸ��
ngy=30;      //y��첨��ĸ�ʽ
dgx=40;      //x��첨��ļ��
dgy=40;      //y��첨��ļ��

esize=4;
dz=2;       //z��
dy=30;      //y��
dx=30;      //x��
b1=1;
b2=1;
b3=1;
order=3;
vel0=1.0/800;
nz=400;
ny=94;
nx=94;

n=nx*ny*nz;
v = new float[n];
ts=new float[n];
tg=new float[n];
danp=new float[n];
dan=new float[n];
AMP=new float[nshot*ngx*ngy*sample];
///////////////////read velocity///////////////////////////

ifstream vv("vel.dat");
for(i=0;i<n;i )
{
vv>>v[i];
danp[i]=0;
}
///////////read segy///////////////

float nx1,nx2,ny1,ny2;
int nxs,nxe,nys,nye;

int nTrace,fileLength,snum;
float vsx,vsy,vgx,vgy,gx,gy;
float *amp,*SX,*GX,*SY,*GY;

amp=new float[sample];
SX=new float[nshot*ngx*ngy];
SY=new float[nshot*ngx*ngy];
GX=new float[nshot*ngx*ngy];
GY=new float[nshot*ngx*ngy];

int ssx,ssy,ggx,ggy;
FILE *in;
in=fopen("model_1200.sgy","rb");
/*fseek(in,3716L,SEEK_SET);
fread(&dt,2,1,in);
dtt=short(dt);
dttt=dtt/(1E 6);
fseek(in,3220L,SEEK_SET);
fread(&sample,4,1,in);
sample=short(sample);
fseek(in,0,SEEK_END); fileLength=ftell(in);
nTrace=(fileLength-3600)/(240 4*sample);*/

fseek(in,3600L,SEEK_SET);
for(i=0;i<nshot*ngx*ngy;i )
{
  fseek(in,72L,SEEK_CUR);
  fread(&ssx,4,1,in);
  HL((unsigned char *) &ssx,4);
     SX[i]=ssx;

  fseek(in,0L,SEEK_CUR);
  fread(&ssy,4,1,in);
  HL((unsigned char *) &ssy,4);
     SY[i]=ssy;

  fseek(in,0L,SEEK_CUR);
  fread(&ggx,4,1,in);
  HL((unsigned char *) &ggx,4);
  GX[i]=ggx;
  //cout<<ggx<<"\t";

        fseek(in,0L,SEEK_CUR);
  fread(&ggy,4,1,in);
  HL((unsigned char *) &ggy,4);
  GY[i]=ggy;

        fseek(in,152L,SEEK_CUR);
  fread(amp,4,sample,in);
  for(j=0;j<sample;j )
  {
   AMP[i*sample j]=amp[j];
  }
  fseek(in,0L,SEEK_CUR);
}
/////////////////travelltime computer and migration////////////////
/*double starttime,endtime;
starttime=MPI_Wtime();
cout<<starttime<<"\n";*/
for(i=0;i<10;i )//nshot
{
   vsy=SY[i*ngx*ngy] 400;
   vsx=SX[i*ngx*ngy] 400;
   fastmarch(order,0,vsy,vsx,b1,b2,b3,nz,ny,nx,dz,dy,dx,vel0,v,ts);//�ڵ��ʱ�̱�
cout<<vsx<<"\n";
   for(j=0;j<ngx;j )//ngy
   {
    for(k=0;k<ngy;k )
    {
     vgx=GX[(i*ngx j)*ngy k] 400;
     vgy=GY[(i*ngx j)*ngy k] 400;
        if((fabs(vsx-vgx)<30)&&(fabs(vsy-vgy)<30))
     {
      fastmarch(order,0,vgy,vgx,b1,b2,b3,nz,ny,nx,dz,dy,dx,vel0,v,tg);//�첨��ʱ��
         /////////��ڵ�ʱ�̱�/////////
      if(vsx<vgx)
      {
       nxs=int(vsx/dx 1.0/2);
       nxe=int(vgx/dx 1.0/2);
      }
      else
      {
       nxs=int(vgx/dx 1.0/2);
       nxe=int(vsx/dx 1.0/2);
      }
      if(vsy<vgy)
      {
       nye=int(vgy/dy 1.0/2);
       nys=int(vsy/dy 1.0/2);
      }
      else
      {
       nys=int(vgy/dy 1.0/2);
       nye=int(vsy/dy 1.0/2);
      }
         for(ii=nxs;ii<=nxe;ii )
      {
       for(jj=nys;jj<=nye;jj )
       {
        for(kk=0;kk<nz;kk )
        {
         t=ts[(ii*ny jj)*nz kk] tg[(ii*ny jj)*nz kk];
                  /////��Ӧ��ȡ////////
                  num=int(t/dt 1/2);
                  if(num<sample)
         {
          float ax=AMP[((i*ngy j)*ngx k)*sample num];
          danp[(ii*ny jj)*nz kk]=danp[(ii*ny jj)*nz kk] ax;
         }
        }
       }
      }
     }
    }
   }
}
/* MPI_Reduce(danp,dan,n,MPI_FLOAT,MPI_SUM,0,MPI_COMM_WORLD);
if(myrank==0)
{
      ofstream tt("mp.dat");
      for(i=0;i<n;i )
      {
          tt<<danp[i]<<"\t";
      }
   }

endtime=MPI_Wtime();
cout<<endtime<<"\n";
printf("It tooks %f seconds\n",endtime-starttime);
MPI_Finalize();*/
ofstream tt("mp.dat");
for(i=0;i<n;i )
{
tt<<danp[i]<<"\t";
}
delete [] AMP;
delete [] danp;
delete [] ts;
delete [] tg;
delete [] v;
return 0;
}

void pqueue_init (int n)
{
x = (float **) malloc ((n 1)*sizeof (float *));
xn = x;
x1 = x 1;
}

void pqueue_close (void)
{
free (x);
}

void pqueue_insert (float* v)
{
float **xi, **xp;
unsigned int p;

xi = xn;
*xi = v;
p = (unsigned int) (xn-x);
for (p >>= 1; p > 0; p >>= 1)
{
    xp = x p;
    if (*v > **xp) break;
    *xi = *xp; xi = xp;
}
*xi = v;
}

float* pqueue_extract (void)
{
unsigned int c;
int n;
float *v, *t;
float **xi, **xc;

v = *x1;
*(xi = x1) = t = *(xn--);
n = (int) (xn-x);
for (c = 2; c <= n; c <<= 1) {
    xc = x c;
    if (c < n && **xc > **(xc 1)) {
      c ; xc ;
    }
    if (*t <= **xc) break;
    *xi = *xc; xi = xc;
}
*xi = t;
return v;
}

static int grid (int i, int n)
{
if (i >=0 && i < n) return i;
if (i < 0) return 0;
return (n-1);
}

static double a, tp, d1, d2, d3, dd[8], dd1,dd2,dd3;
static int nm, n, n1, n2, n3, n12;

void fastmarch(int order, float s1, float s2, float s3,
       int b1, int b2, int b3,
       int nz, int ny, int nx,
       float dz, float dy, float dx,
       float slow0, float *slow, float *ttime)
{
int i1,i2,i3, start1,start2,start3, end1,end2,end3, i, nh;
int stillconstant;
double delta1, delta2, delta3;
float *pt;
unsigned char *mask, *pm;

/* save to static parameters */
n1 = nz; n2 = ny; n3 = nx;
d1 = dz; d2 = dy; d3 = dx;
n12 = n1*n2;
n = n12*n3;

/* should probably check the return value of malloc */
mask = (unsigned char *) malloc (n*sizeof(unsigned char));

for (i=0; i<n; i ) {
ttime[i] = 1.e 20;
mask[i] = FMM_OUT;
}