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

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; //ÿ��IJ�������
  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;
  }