sufctanismod.c

su 的源代码库

  ll = alloc2float(nz,nx);	
  nn = alloc2float(nz,nx);	
  rho = alloc2float(nz,nx);	
  fx = alloc1float(nt);
  fy = alloc1float(nt);
  fz = alloc1float(nt);
  t = alloc1float(nt);

  /* allocate optional space for store reflection and VSP seismogram */
  if (*reflxfile != '\0') {  /* allocate space for refl_x  */
    refl_x = alloc2float(nt, nx); 
  }
  if (*reflyfile != '\0') {  /* allocate space for refl_y  */
    refl_y = alloc2float(nt, nx); 
  }
  if (*reflzfile != '\0') {  /* allocate space for refl_z  */
    refl_z = alloc2float(nt, nx); 
  }
  if (*vspxfile != '\0') {  /* allocate space for vsp_x  */
    vsp_x = alloc2float(nt, nx); 
  }
  if (*vspyfile != '\0') {  /* allocate space for vsp_y  */
    vsp_x = alloc2float(nt, nx); 
  }
  if (*vspzfile != '\0') {  /* allocate space for vsp_z  */
    vsp_z = alloc2float(nt, nx); 
  }


  /*   initial condition  */
  for (ix=0; ix < nx; ix++)
    for (iz=0; iz < nz; iz++)
      {
	u[ix][iz]=0.0;
	v[ix][iz]=0.0;
	w[ix][iz]=0.0;
	xzsource[ix][iz]=0.0;
      }
		
  for (ix=0; ix<nx; ix++)
    for (iz=0; iz<nz; iz++)
      {
	e11[ix][iz]=0.0;
	e33[ix][iz]=0.0;
	e23[ix][iz]=0.0;
	e13[ix][iz]=0.0;
	e12[ix][iz]=0.0;
      }

  /*   get time response of the source function */
  for (it=0; it<nt; it++) 
    {
      t[it]=it*dt;
      fx[it]=0.0;
      fy[it]=0.0;
      fz[it]=0.0;
    }

  if (indexux) {		/*  x-component of the force  */
    if (wavelet == 1)
      tforce_akb(nt, fx, dt, fpeak);
    if (wavelet == 2)
      tforce_ricker(nt, fx, dt, fpeak);
    if (wavelet == 3) 	
      tforce_spike(nt, fx, dt, fpeak);
    if (wavelet == 4) 	
      tforce_unit(nt, fx, dt, fpeak);
  }
  if (indexuy) {		/*  y-component of the force  */
    if (wavelet == 1)
      tforce_akb(nt, fy, dt, fpeak);
    if (wavelet == 2)
      tforce_ricker(nt, fy, dt, fpeak);
    if (wavelet == 3) 	
      tforce_spike(nt, fy, dt, fpeak);
    if (wavelet == 4) 	
      tforce_unit(nt, fy, dt, fpeak);
  }
  if (indexuz) {		/*  z-component of the force  */
    if (wavelet == 1)
      tforce_akb(nt, fz, dt, fpeak);
    if (wavelet == 2)
      tforce_ricker(nt, fz, dt, fpeak);
    if (wavelet == 3) 	
      tforce_spike(nt, fz, dt, fpeak);
    if (wavelet == 4) 	
      tforce_unit(nt, fz, dt, fpeak);
  }


  /*    obtain density and elastic parameters  */
  vmax = 0.0;
  read_parameter(nx, nz, dx, dz, rho00, drhodx, drhodz, dfile, rho); 
  read_parameter(nx, nz, dx, dz, aa00, daadx, daadz, afile, aa); 
  read_parameter(nx, nz, dx, dz, cc00, dccdx, dccdz, cfile, cc); 
  read_parameter(nx, nz, dx, dz, ff00, dffdx, dffdz, ffile, ff); 
  read_parameter(nx, nz, dx, dz, ll00, dlldx, dlldz, lfile, ll); 
  read_parameter(nx, nz, dx, dz, nn00, dnndx, dnndz, nfile, nn); 

  /*  compute density inverse and maximum velocity  */	
  for (ix=0; ix < nx; ix++)   
    for (iz=0; iz < nz; iz++)
      {
	rho[ix][iz]=1.0/rho[ix][iz];
	if ( cc[ix][iz]*rho[ix][iz] > vmax ) 
	  vmax = cc[ix][iz]*rho[ix][iz];
      }
  vmax = sqrt(vmax);

  /*  give source location  */
  if (*sfile != '\0') {
    read_parameter(nx, nz, dx, dz, aa00, daadx, daadz, sfile, xzsource); 
  } else {
    locate_source(nx, nz, sx, sz, xzsource, source);
  }

  /* Debugging */
  nxcc2=nx;
  nzcc2=nz;
  /* end debugging */

  /*     evolve in time    */
  outfpx = fopen("snapshotx.data", "w");
  outfpy = fopen("snapshoty.data", "w");
  outfpz = fopen("snapshotz.data", "w");
  for (it=0; it < nt; it++) {
      fprintf (stderr,"check it= %d\n", it);

      moving_bc (it, nx, nz, sx, sz, 
		 dx, dz, impulse, movebc, t, vmax, 
		 &mbx1, &mbz1, &mbx2, &mbz2);

      /* FCT correction is localized to the area bounded by */
      /* (fctxend - fctxbeg) by (fctzend-fctzbeg) */
      /* contain the FCT correction boundary by the model boundary */
      /* computed by moving_bc */
      moving_fctbc (mbx1, mbz1, mbx2, mbz2, 
		    nxcc1, nzcc1, nxcc2, nzcc2, 
		    &fctxbeg, &fctzbeg, &fctxend, &fctzend);

      anis_solver2(it, u, v, w, e11,  
		   e33, e23, e12, e13, aa, cc, ff, ll, nn, 
		   rho, xzsource, fx, fy, fz, 
		   dx, dz, dt, nx, nz, dofct, isurf, 
		   eta0, eta, deta0dx, deta0dz, 
		   detadx, detadz, mbx1, mbx2, mbz1, mbz2);

      /*  get reflection seismogram  */
      if (*reflxfile != '\0') {  /* get refl_x  */
	for (ix=0; ix<nx; ix++) {
	  refl_x[ix][it]=u[ix][receiverdepth];
	}
      }
      if (*reflyfile != '\0') {  /* get refl_y  */
	for (ix=0; ix<nx; ix++) {
	  refl_y[ix][it]=v[ix][receiverdepth];
	}
      }
      if (*reflzfile != '\0') {  /* get refl_z  */
	for (ix=0; ix<nx; ix++) {
	  if (w[ix][receiverdepth]>1){
	    warn("problems in %d x=%d \n",it, ix); 
	  }
	  refl_z[ix][it]=w[ix][receiverdepth];
	}
      }

      /*  get VSP seismogram  */
      if ( vspnx >= 0 && vspnx < nx) {/*position not out of range*/
	if (*vspxfile != '\0') {  /* get vsp_x  */
	  for (iz=0; iz<nz; iz++) {
	    vsp_x[iz][it]=u[vspnx][iz];
	  }
	}
	if (*vspyfile != '\0') {  /* get vsp_x  */
	  for (iz=0; iz<nz; iz++) {
	    vsp_y[iz][it]=v[vspnx][iz];
	  }
	}
	if (*vspzfile != '\0') {  /* get vsp_x  */
	  for (iz=0; iz<nz; iz++) {
	    vsp_z[iz][it]=w[vspnx][iz];
	  }
	}
      }
  fwrite(u[0], sizeof(float), nz*nx, outfpx);	
  fwrite(v[0], sizeof(float), nz*nx, outfpy);	
  fwrite(w[0], sizeof(float), nz*nx, outfpz);	

    }	


  fwrite(u[0], sizeof(float), nz*nx, outfp);	

  /*  output seismogram  */
  if (*reflxfile != '\0') {  /* write refl_x */
    if ((outreflx = fopen(reflxfile, "w"))==NULL) {
      fprintf(stderr, "Cannot open file=%s\n", reflxfile); 
      exit(1);
    }
    if (suhead==1){
      su_output(nx,sx,sz,nt,dx,dt,outreflx,refl_x);
    }
    else
      {
	fwrite(refl_x[0], sizeof(float), nx*nt, outreflx);
      }
    fclose(outreflx);
  }
  if (*reflyfile != '\0') {  /* write refl_y */
    if ((outrefly = fopen(reflyfile, "w"))==NULL) {
      fprintf(stderr, "Cannot open file=%s\n", reflyfile); 
      exit(1);
    }
    if (suhead==1){
      su_output(nx,sx,sz,nt,dx,dt,outrefly,refl_y);
    }
    else
      {
	fwrite(refl_y[0], sizeof(float), nx*nt, outrefly);
      }
    fclose(outrefly);
  }
  if (*reflzfile != '\0') {  /* write refl_z */
    if ((outreflz = fopen(reflzfile, "w"))==NULL) {
      fprintf(stderr, "Cannot open file=%s\n", reflzfile); 
      exit(1);
    }
    if (suhead==1){
      su_output(nx,sx,sz,nt,dx,dt,outreflz,refl_z);
    }
    else
      {
	fwrite(refl_z[0], sizeof(float), nx*nt, outreflz);
      }
    fclose(outreflz);
  }
  if (*vspxfile != '\0') {  /* write vsp_x */
    if ((outvspx = fopen(vspxfile, "w"))==NULL) {
      fprintf(stderr, "Cannot open file=%s\n", vspxfile); 
      exit(1);
    }
    if (suhead==1){
      su_output(nz,sx,sz,nt,dx,dt,outvspx,vsp_x);
    }
    else
      {
	fwrite(vsp_x[0], sizeof(float), nz*nt, outvspx);
      }
    fclose(outvspx);
  }
  if (*vspyfile != '\0') {  /* write vsp_x */
    if ((outvspy = fopen(vspyfile, "w"))==NULL) {
      fprintf(stderr, "Cannot open file=%s\n", vspyfile); 
      exit(1);
    }	
    if (suhead==1){
      su_output(nz,sx,sz,nt,dx,dt,outvspy,vsp_y);
    }
    else
      {
	fwrite(vsp_y[0], sizeof(float), nz*nt, outvspy);
      }
  }
  if (*vspzfile != '\0') {  /* write vsp_z */
    if ((outvspz = fopen(vspzfile, "w"))==NULL) {
      fprintf(stderr, "Cannot open file=%s\n", vspzfile); 
      exit(1);
    }
    if (suhead==1){
      su_output(nz,sx,sz,nt,dx,dt,outvspz,vsp_z);
    }
    else
      {
	fwrite(vsp_z[0], sizeof(float), nz*nt, outvspz);
      }
  }

  return(CWP_Exit());
}
/*------------------ end of main program ----------------------------*/

/************************************************************************
* read_parameter --  Obtain elastic parameter, either read from a file 
*			or assume linear variation
*************************************************************************
*************************************************************************
* Input: 
*	int nx		number of grids in x direction 
*	int nz		number of grids in z direction
*	float dx	spatial step in x-direction 	
*	float dz	spatial step in z-direction 	
*	float p00=2.0	elastic parameter at (0, 0) 
*	float dpdx=0.0	parameter gradient in x-direction  d(p)/dx
*	float dpdz=0.0	parameter gradient in z-direction  d(p)/dz 
*	char *file	file name store the elastic parameter
*	
*************************************************************************
* Output: float **pp	elastic parameter array
*
*************************************************************************
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void 	read_parameter(int nx, int nz, float dx, float dz, 
		       float p00, float dpdx, float dpdz, char *file, float **pp)
{
  int 	ix, iz;
  FILE	*infp;

  /*    obtain  elastic parameter  */
  if (*file != '\0') {  /* open requested parameter file */
    if ((infp = fopen(file, "r"))==NULL) {
      fprintf(stderr, "Cannot open file=%s\n", file); 
      exit(1);
    }

    /* read elastic parameter into the array pp[nz][nx] */
    fread(pp[0], sizeof(float), nz*nx, infp);
  } else { /* assume a linear parameter profile */
    for (ix=0; ix < nx; ix++)
      for (iz=0; iz < nz; iz++)
	{
	  pp[ix][iz]=p00+dpdx*(ix-1)*dx
	    +dpdz*(iz-1)*dz;
	}
  } 
}



/************************************************************************
*
* anis_solver2 ---  2nd-order 2-D finite-difference solver for the 
*		first order elastic wave equation
*
*************************************************************************
*************************************************************************
* Input:
*
*	for the elastic wave equation
*	float **u	x-component, solution at previous time level
*	float **v	y-component, solution at previous time level
*	float **w	z-component, solution at previous time level
*	float **aa	elastic parameter
*	float **cc	elastic parameter
*	float **ff	elastic parameter
*	float **ll	elastic parameter
*	float **nn	elastic parameter
*	float **rho	reverse of density field
*	float *fx	time response of the source function, x-component
*	float *fy	time response of the source function, y-component
*	float *fz	time response of the source function, z-component
*	float **xzsource	source location
*
*	int it		current time step
*	float dx	spacing in x-direction
*	float dz	spacing in z-direction
*	int nx		number of grids in x-direction
*	int mbx1	sample start in x-direction
*	int mbx2	sample end in x-direction
*	int nz		number of grids in z-direction
*	int mbz1	sample start in z-direction
*	int mbz2	sample end in z-direction
*	int fctxbeg	left boundary for the FCT correction
*	int fctxend	right boundary for the FCT correction
*	int fctzbeg	top boundary for the FCT correction
*	int fctzend	bottom boundary for the FCT correction 
*	int dofct	FCT correction index
*			=1 do FCT
*			=0 not do the FCT
*	int isurf	index for surface boundary condition
*	float eta0	diffusion coefficient
*	float deta0dx	gradient of eta0 in x-direction
*	float deta0dz	gradient of eta0 in z-direction
*	float eta	anti-diffusion coefficient
*	float detadx	gradient of eta in x-direction
*	float detadz	gradient of eta in z-direction
*
************************************************************************
* Output:
*	float **u	solution of x-component at next time level
*	float **v	solution of y-component at next time level
*	float **w	solution of z-component at next time level
************************************************************************
* Notes: this is a 2nd order finite-difference routine on staggered 
*	with optional FCT correction, and boundary conditions included.
************************************************************************
* Author:	Tong Fei (1993), Colorado School of Mines
************************************************************************/
void	anis_solver2(int it, float **u, float **v, float **w, 
		     float **e11, float **e33, float **e23, float **e12, float **e13,