sugoupillaud.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* SUGOUPILLAUD: $Revision: 1.3 $ ; $Date: 2006/11/07 22:58:42 $	*/

#include "su.h"
#include "segy.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"									     ",
" SUGOUPILLAUD - calculate 1D impulse response of                            ",
"                non-absorbing Goupillaud medium                             ",
"									     ",
" sugoupillaud < stdin > stdout [optional parameters]			     ",
"									     ",
" Required parameters:                                                       ",
"        none                                                                ",
"                                                                            ",
" Optional parameters:                                                       ",
"        l=1           source layer number; 1 <= l <= tr.ns                  ",
"                      Source is located at the top of layer l.              ",
"        k=1           receiver layer number; 1 <= k                         ",
"                      Receiver is located at the top of layer k.	     ",
"        tmax          number of output time-samples; default:               ",
"                      tmax=NINT((2*tr.ns-(l-1)-(k-1))/2)  if k < tr.ns      ",
"                      tmax=k                              if k >=tr.ns      ",
"        pV=1          flag for vector field seismogram                      ",
"                      (displacement, velocity, acceleration);               ",
"                      =-1 for pressure seismogram.                          ",
"        verbose=0     silent operation, =1 list warnings                    ",
"									     ",
" Input: Reflection coefficient series:                                      ",
"                                                                            ",
"		               impedance[i]-impedance[i+1]                   ",
"		       r[i] = -----------------------------                  ",
"		               impedance[i]+impedance[i+1]                   ",
"                                                                            ",
"        r[0]= surface refl. coef. (as seen from above)                      ",
"        r[n]= refl. coef. of the deepest interface                          ",
"		                                        		     ",
" Input file is to be in SU format, i.e., binary floats with a SU header.    ",
"									     ",
" Remarks:                                                                   ",
" 1. For vector fields, a buried source produces a spike of amplitude 1      ",
" propagating downwards and a spike of amplitude -1 propagating upwards.     ",
" A buried pressure source produces spikes of amplitude 1 both in the up-    ",
" and downward directions.                                                   ",
"    A surface source induces only a downgoing spike of amplitude 1 at the   ",
" top of the first layer (both for vector and pressure fields).              ",
" 2. The sampling interval dt in the header of the input reflectivity file   ",
" is interpreted as a two-way traveltime thicknes of the layers. The sampling",
" interval of the output seismogram is the same as that of the input file.   ",
NULL};

/* 
 * Credits:
 *	CWP: Albena Mateeva, May 2000, a summer project at Western Geophysical
 *
 *
 * ANOTATION used in the code comments [arises from the use of z-transforms]:
 *           Z-sampled: sampling interval equal to the TWO-way 
 *                      traveltime of the layers; 
 *           z-sampled: sampling interval equal to the ONE-way
 *                      traveltime of the layers;
 *
 * REFERENCES:
 *
 *	1. Ganley, D. C., 1981, A method for calculating synthetic seismograms 
 *	which include the effects of absorption and dispersion. 
 *	Geophysics, Vol.46, No. 8, p. 1100-1107.
 * 
 *	The burial of the source is based on the Appendix of that article.
 *
 *	2. Robinson, E. A., Multichannel Time Series Analysis with Digital 
 *	Computer Programs: 1983 Goose Pond Press, 2nd edition.
 *
 *	The recursive polynomials Q, P used in this code are described
 *	in Chapter 3 of the book: Wave Propagation in Layered Media.
 *
 *	My polynomial multiplication and division functions "prod" and
 *	"pratio" are based on Robinson's Fortran subroutines in Chapter 1.
 *
 *	4. Clearbout, J. F., Fundamentals of Geophysical Data Processing with
 *	Applications to Petroleum Prospecting: 1985 Blackwell Scientific 
 *	Publications.
 *
 *	Chapter 8, Section 3: Introduces recursive polynomials F, G in a 
 *	more intuitive way than Robinson.
 *	
 *	The connection between the Robinson's P_k, Q_k and Clearbout's 
 *	F_k, G_k is:
 *				P_k(Z) = F_k(Z)
 *				Q_k(Z) = - Z^(k) G_k(1/Z)
 *
 */
/**************** end self doc ***********************************/

/* Prototypes of functions used internally */
int  porder(double p[], int l);
void pcomb(double alfa, double beta, double p1[], double p2[], 
	   int l1, int l2, double c[], int l);
void prod(double p1[], double p2[], int l1, int l2, double pp[], int ll);
void pratio(double p1[], double p2[], int l1, int l2, int L, 
	    double q[], int lq);
void pshift(double p[], int l, int k, double psh[], int ll);
void recurs(int k, double q[], double p[], int n, double r[]);
void rev(double p[], int l, double pr[] );
void upsample(double p4[], int l, double p2[], int ll);


segy tr;

int
main(int argc, char **argv)
{
  int n;	         /* number of subsurface interfaces 	        */
  int l, k;	         /* source and receiver layers		        */
  int tmax;	         /* number of output samples 		        */
  int pV;	         /* field-type flag			        */
  double *r;	         /* reflection coefficient series 		*/
  double *x;	         /* output seismogram			        */

  int i;		 /* loop counter			        */
  int verbose;   	 /* verbose flag			        */
  int kk, delay;	 /* help accomodate receiver in the lower 
			    homogeneous half-space (k>n+1)	        */
  double *u1, *d1;	 /* up- and down-going waves in	               
			    the 1st layer for surface source            */
  double *u1b;		 /* up-going waves in the first layer	
			    for buried source			        */
  double *uk, *dk;	 /* up- and down-going waves in layer k;
			    output seismogram:  x=uk+dk                 */ 
  double *Q, *P;	 /* recursive polynomials; see Robinson         */
  double *qr, *pr;     	 /* reverse recursive polynomials	        */
  double *cn;		 /* r[0]*Qn - Pn; for convenience	        */
  double uni[1]={1.};	 /* unit "polynomial" - for convenience	        */
  double *t1, *t2;       /* temporary storage for Q,P var-s (Z-sampled) */
  double *td1;           /* temp. storage for z-sampled Q,P             */
  double *mix1, *mix2;   /* temporary storage for Q*U-type terms        */
  double *mix3, *mix4;	 /* temporary storage for Q*U-type terms        */
  double *mix;           /* temporary storage for z-sampled mix1,2      */
  double ak;		 /* product of k-1 transmission coef-s	        */
  double *d0;		 /* shaping for a buried source; see Ganley     */
  double *storx;	 /* temporary storage for x when k=l>1	        */


  /* Initialize */
  initargs(argc, argv);
  requestdoc(1);
	
  /* Get subsurface model */
  if (!gettr(&tr)) err("Can't get reflectivity");
  n = tr.ns - 1 ;

  /* Get parameters */
  if (!getparint("k",&k)) 	        k=1;
  if (!getparint("l",&l))	        l=1;
  if (!getparint("tmax",&tmax)){
    if (k<n+1)                          tmax=(2*n+2-(l-1)-(k-1))/2;
    else                                tmax=k;
  }
  if (!getparint("pV",&pV))	        pV=1;
  if (!getparint("verbose",&verbose))	verbose=0;

  /* Check parameters */
  if (n<=0) err("The number of subsurface interfaces n must be >=1!");
  if (k<1)  err("Receiver layer k must be >=1 (k=1 corresponds to surface seismogram).");
  if (l<1)  err("Source layer l must be >= 1 (l=1 corresponds to a surface source).");
  if (l>(n+1)) err("The current version of the program requires l<=n+1.");
  if (tmax<0)  err("The number of the output time samples tmax cannot be negative.");
  if(!( pV==1 || pV==-1 )) err("The field-type flag pV should be either 1 or -1.");

  /* Verbose */
  if (verbose) {
    warn("Number of layers n=%d", n);
    warn("Source layer l=%d", l);
    warn("Receiver layer k=%d", k);
    warn("Output time samples tmax=%d", tmax);
    warn("Field type flag pV=%d", pV);
  }

  /* Allocate Memory */
  r = ealloc1double(n+1);
  x = ealloc1double(2*tmax);
  u1 = ealloc1double(2*tmax+n);
  d1 = ealloc1double(2*tmax+n);
  u1b = ealloc1double(2*tmax+n);
  uk = ealloc1double(2*tmax);
  dk = ealloc1double(2*tmax);
  storx = ealloc1double(2*tmax);
  Q = ealloc1double(n+1); 
  P = ealloc1double(n+1);	
  qr = ealloc1double(n+1); 
  pr = ealloc1double(n+1); 
  cn = ealloc1double(n+1); 
  t1 = ealloc1double(n+1);
  t2 = ealloc1double(n+1); 
  d0 = ealloc1double(n+1);
  td1 = ealloc1double(2*(n+1));
  mix1 = ealloc1double(2*tmax+n);
  mix2 = ealloc1double(2*tmax+n);
  mix3 = ealloc1double(2*(tmax+n));
  mix4 = ealloc1double(2*(tmax+n));
  mix = ealloc1double(2*(tmax+n));
  
  /* Read and check reflectivity values */
  for (i=0; i<=n; ++i) 
    {
      r[i] = tr.data[i];
      if(r[i]>1. || r[i]<-1.)
	err("Invalid reflection coefficient encountered.");
    }
  
  /* Account for field type (pressure/displacement) */
  for (i=0; i<=n; ++i)  r[i]*=pV;
  
  /********************* Begin computation **************************/
  if(n==0)
    /* Trivial case - implies source in l==1, receiver in k>=1 */
    {
      memset( (void *) x, 0, 2*tmax*DSIZE);
      x[0] = 1.;
      pshift(x,2*tmax,k-1,x,2*tmax);
    } 
  else 
    {
      /* Computing u1: upgoing field in layer 1
         u1 is needed for all source-receiver configurations	*/
      
      recurs(n,Q,P,n,r);		          /* compute Qn, Pn    */
      pcomb(r[0],-1.,Q,P,n+1,n+1,cn,n+1);         /* cn = r[0]*Qn - Pn */  
      pratio(Q,cn,n+1,n+1,2*tmax+n,u1,2*tmax+n);  /* u1 = Qn/cn	       */

      /* sampling interval of u1: Z (two-way time-thickness of layers) */
    
    
      if(l==1)
	/****************** I. SURFACE SOURCE *******************/
	{ 
	  
	  if(k>n)
	    /* I.1. Receiver in the homogen. half-space below the layers */
	    /* Only downgoing waves come to the receiver => x = dk       */
	    /* The wavefield in layer k>=(n+1) is the same as in layer 
	       n+1 but delayed by k-(n+1) one-way traveltime units z.    */
	    {  
	      /* compute wavefield in layer n+1:
                 -z^(n)(1+r[1])...(1+r[n])/cn */
	
	      /* ak = - (1+r[1])*...*(1+r[n]) */
	      for (i=1, ak=1; i<=n; ++i)  ak *= 1+r[i];
	      ak *= -1.;

	      /* compute 1/cn; Z-sampled so far	*/
	      pratio(uni,cn,1,n+1,2*tmax,dk,2*tmax);
	      	
	      /* resample 1/cn to z */
	      upsample(dk,2*tmax,dk,2*tmax);

	      /* z^(n)/cn; z-sampled */
	      pshift(dk,2*tmax,n,dk,2*tmax);
	
	      /* final value of dk at the top 
		 of layer n+1: dk=ak*z^(n)/cn */
	      for (i=0; i<2*tmax; ++i)  dk[i]*=ak;
	      
	      /* account for the one-way delay 
		 between layer n+1 and k:     */
	      pshift(dk,2*tmax,k-n-1,x,2*tmax);

	    } 
	  else 
	    /* I.2. Receiver in the layers k:[1;n]; surface source    */
	    /* To compute the up- and down-going fields in layer k, 
	       first compute the up- and down-going fields in layer 1 
	       and propagate them to layer k.  u1 is already found.   */ 
	    {
	      /* Computing d1 = -Pn/cn (P=Pn computed before) */
	      pratio(P,cn,n+1,n+1,2*tmax+n,d1,2*tmax+n); 
	      for (i=0; i<2*tmax+n; ++i)  
		d1[i] *= -1.;		   /* d1 is Z-sampled */
	      
	      if(k==1)
		{ 
		  /* I.2.1. SURFACE SEISMOGRAM from a surface source
		            (no need to propagate u1 and d1)          */
		  
		  /* total field x = d1 + u1 */
		  pcomb(1.,1.,u1,d1,2*tmax+n,2*tmax+n,x,2*tmax);
		  /* x is Z-sampled here */
		  /* resample x to z */
		  upsample(x,2*tmax,x,2*tmax);
		} 
	      else 
		{
		  /* I.2.2. RECEIVER AT DEPTH k:(1;n]; surface source; 
		            have to propagate u1 and d1 to layer k   */
		  
		  /*	     d1*Q_(k-1) + u1*P_(k-1)
		     uk = -------------------------------
			  z^(k-1) (1-r[1])*...*(1-r[k-1])		       

			  
			     d1*pr_(k-1) + u1*qr_(k-1)
		     dk = -------------------------------
			  z^(k-1) (1-r[1])*...*(1-r[k-1])   	     */
		  
	  
		  /* scaling coeficient ak=(1-r[1])...(1-r[k-1]) */
		  for (i=1, ak=1; i<k; ++i)  ak *= 1-r[i];
		  
		  /*______ upgoing field uk: ______*/ 
		  recurs(k-1,Q,P,n,r);
		  prod(Q,d1,n+1,2*tmax+n,mix1,2*tmax+n);
		  prod(P,u1,n+1,2*tmax+n,mix2,2*tmax+n);
		  pcomb(1.,1.,mix1,mix2,2*tmax+n,2*tmax+n,mix1,2*tmax+n);
		  upsample(mix1,2*tmax+n,mix,2*(tmax+n));
		  pshift(mix,2*(tmax+n),-(k-1),mix1,2*tmax+n);
		  for (i=0; i<2*tmax; ++i)  
		    uk[i]=mix1[i]/ak;        /* uk is z-sampled here */
		  
		  /*____ downgoing field dk: _______*/
		  rev(Q,n+1,qr);
		  pshift(qr,n+1,k-1-porder(Q,n+1),qr,n+1);
		  rev(P,n+1,pr);
		  pshift(pr,n+1,k-1-porder(P,n+1),pr,n+1);