sugoupillaud.c

su 的源代码库

		  prod(pr,d1,n+1,2*tmax+n,mix1,2*tmax+n); 
		  /* pr, d1, mix1 are Z-sampled */
		  prod(qr,u1,n+1,2*tmax+n,mix2,2*tmax+n); 
		  /* qr, u1, mix2 are Z-sampled */
		  pcomb(1.,1.,mix1,mix2,2*tmax+n,2*tmax+n,mix1,2*tmax+n);
		  upsample(mix1,2*tmax+n,mix1,2*tmax+n);
		  pshift(mix1,2*tmax+n,-(k-1),mix1,2*tmax+n);
		  for (i=0; i<2*tmax; ++i)	   
		    dk[i]=mix1[i]/ak;        /* dk is z-sampled here */
	  
		  /* total wave field in layer k */
		  pcomb(1.,1.,uk,dk,2*tmax,2*tmax,x,2*tmax);

		} 
	    }
	} 
      else 
	/********************** II. BURIED SOURCE *************************
	 d0 needed:
	     
	        P_(l-1) - r[0]*Q_(l-1) - pV*qr_(l-1) + pV*r[0]*pr_(l-1)
	   d0 = -------------------------------------------------------
	                 z^(l-1) (1+r[1])*...*(1+r[l-1])                  

	*******************************************************************/
	{ 
	  recurs(l-1,Q,P,n,r);
	  pcomb(1.,-r[0],P,Q,n+1,n+1,t1,n+1);

	  rev(Q,n+1,qr);
	  pshift(qr,n+1,l-1-porder(Q,n+1),qr,n+1);
	  rev(P,n+1,pr);
	  pshift(pr,n+1,l-1-porder(P,n+1),pr,n+1);

	  pcomb(-pV*1.,pV*r[0],qr,pr,n+1,n+1,t2,n+1);
	  pcomb(1.,1.,t1,t2,n+1,n+1,d0,n+1);
	  for (i=1, ak=1; i<l; ++i)  ak *= 1+r[i];
	  for (i=0; i<n+1; ++i)       d0[i] *= 1/ak;
	  /* d0 is ALWAYS Z-sampled;
	     shift by (l-1) z postponed to avoid negative powers of z in d0 */
	  
	  /* for any buried source, the combination d0*u1 is needed */

	  prod(u1,d0,2*tmax+n,n+1,mix1,2*tmax+n); /* u1, d0 are Z-sampled   */
	  upsample(mix1,2*tmax+n,mix,2*(tmax+n)); /* now d0*u1 is z-sampled */
	  
	  /* including the postponed (l-1)z shift */
	  pshift(mix,2*(tmax+n),-(l-1),mix1,2*tmax+n);
	  
	  if( k>=l )
	    {
	      /* II.1. Receiver Below the Buried Source
		 (no need to compute u1b; u1 and d0 are sufficient) */
	      
	      /* d1 needed: d1=1-r[0]*u1 ; Z-sampled */
	      pcomb(1.,-r[0],uni,u1,1,2*tmax+n,d1,2*tmax+n);
	    
	      if(k<=n+1) /* receiver in the layers */
		{
		  kk=k;
		  delay=0;
		} 
	      else   /* receiver below the layers:
			The wavefield is the same as that in layer n+1 but 
			delayed by k-(n+1) one-way traveltime units z. */
		{
		  kk=n+1;
		  delay=k-n-1;
		} 
	      
	      /* Compute the field in layer kk and delay it by z=delay: */
	      
	      for (i=1, ak=1; i<kk; ++i)  ak *= 1-r[i];
	      
	      recurs(kk-1,Q,P,n,r);

	      /* computing uk */
	      /* uk=0 for a receiver below the layers: k>n; */
	      if(k<=n) 
		{
		  prod(Q,d1,n+1,2*tmax+n,mix3,2*(tmax+n));
		  prod(P,u1,n+1,2*tmax+n,mix4,2*(tmax+n));
		  pcomb(1/ak,1/ak,mix3,mix4,2*(tmax+n),2*(tmax+n),mix3,2*(tmax+n));
		  prod(d0,mix3,n+1,2*(tmax+n),mix3,2*(tmax+n));
		  upsample(mix3,2*(tmax+n),mix3,2*(tmax+n));
		  /* for the shift in d0 not accounted before: */
		  pshift(mix3,2*(tmax+n),-(l-1+kk-1),uk,2*tmax);
		} 
	      else
		memset( (void *) uk, 0, 2*tmax*DSIZE);

	      /* computing dk */
	      rev(Q,n+1,qr);
	      pshift(qr,n+1,kk-1-porder(Q,n+1),qr,n+1);
	      rev(P,n+1,pr);
	      pshift(pr,n+1,kk-1-porder(P,n+1),pr,n+1);

	      prod(pr,d1,n+1,2*tmax+n,mix3,2*(tmax+n));
	      prod(qr,u1,n+1,2*tmax+n,mix4,2*(tmax+n));
	      pcomb(1/ak,1/ak,mix3,mix4,2*(tmax+n),2*(tmax+n),mix3,2*(tmax+n));
	      upsample(mix3,2*(tmax+n),mix3,2*(tmax+n));
	      /* include the postponed -(l-1)z shift in d0: */
	      pshift(mix3,2*(tmax+n),-(l-1+kk-1),mix3,2*(tmax+n));
	      pshift(mix3,2*(tmax+n),delay,mix3,2*(tmax+n));
	      for (i=0; i<2*tmax; ++i)  dk[i]=mix3[i];

	      /* total field in layer k */
	      pcomb(1.,1.,uk,dk,2*tmax,2*tmax,x,2*tmax);
	    
	      /* If source and receiver are buried in the same layer:
		 Have to average the results from II.1. RECEIVER BELOW and 
		 II.2. RECEIVER ABOVE the buried source to ensure 
		 displacement=0 or pressure=1 at the source at time 0. */
	      /* Store the result from II.1: */
	      if (k==l) for (i=0; i<2*tmax; ++i)  storx[i]=x[i];
	    
	    }
	  
	  
	  if (k<=l) 
	    {
	    /* II.2. Receiver Above the Buried Source:
	       have to compute u1b (u1*d0 is already in mix1):
	       
	                         Q_(l-1) - pV*pr_(l-1)
	       u1b = u1*d0 + -------------------------------
	                     z^(l-1) (1+r[1])*...*(1+r[l-1])        */
	    
	    
	      recurs(l-1,Q,P,n,r);
	      rev(P,n+1,pr);
	      pshift(pr,n+1,l-1-porder(P,n+1),pr,n+1);
	      pcomb(1.,-pV*1.,Q,pr,n+1,n+1,t1,n+1);
	      upsample(t1,n+1,td1,2*n+2);
	      pshift(td1,2*(n+1),-(l-1),td1,2*(n+1));
	      for (i=1, ak=1; i<l; ++i)  ak *= 1+r[i];
	      for (i=0; i<2*(n+1); ++i)  td1[i] *= 1/ak;
	      pcomb(1.,1.,mix1,td1,2*tmax+n,2*(n+1),u1b,2*tmax+n);  
	      /* u1b is z-sampled */
	    
	      if(k==1)
		/* SURFACE SEISMOGRAM for a buried source
		   u1b is enough to compute it: x = (1-r[0])*u1b    */
		for (i=0; i<2*tmax; ++i)  x[i] = u1b[i]*(1-r[0]);
	      
	      else /* 1< k <= l */
		{		
		  /* BURIED RECEIVER ABOVE the buried source:
	      	     have to propagate the field in the first layer 
		     (upgoing u1b; downgoing -r[0]*u1b) to layer k. */
		  
		  recurs(k-1,Q,P,n,r);
	      
		  /* uk related term of the seismogram: */
		  pcomb(1.,-r[0],P,Q,n+1,n+1,t1,n+1);   /* t1 is Z-sampled */
	      
		  /* dk related term of the seismogram: */
		  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);

		  pcomb(1.,-r[0],qr,pr,n+1,n+1,t2,n+1);	/* t2 is Z-sampled */

		  pcomb(1.,1.,t1,t2,n+1,n+1,t1,n+1);
		  upsample(t1,n+1,td1,2*n+2);
		  prod(td1,u1b,2*n+2,2*tmax+n,mix,2*(tmax+n));
		  pshift(mix,2*(tmax+n),-(k-1),mix2,2*tmax+n);
		  for (i=1, ak=1; i<k; ++i)  ak *= 1-r[i];
		  
		  /* total field in layer k: */
		  for (i=0; i<2*tmax; ++i)  x[i]=mix2[i]/ak;

		  /* If source and receiver are buried in the same layer:
		     Have to average the results from II.1. RECEIVER BELOW 
		     (stored in storx) and II.2. RECEIVER ABOVE the buried 
		     source (just obtained) to ensure displacement=0 or 
		     pressure=1 at the source at time 0.                */
		  /* average seismogram: */
		  if(k==l)  pcomb(.5,.5,storx,x,2*tmax,2*tmax,x,2*tmax);
		  
		}
	    }
	}
    }
  
  /* output synthetic seismogram (resampled to Z) */
  i=(k-l)/2;
  if ( 2*i == k-l ){  
    for (i=0; i<tmax; ++i)  tr.data[i]=x[2*i];
    tr.delrt=0;
  }
  else{
    for (i=0; i<tmax; ++i)  tr.data[i]=x[2*i+1];
    tr.delrt=tr.dt/2000;
  }

  tr.ns=tmax;
  tr.trid=1;

  puttr(&tr);
  
  return(CWP_Exit());

  /* Free space */
  free1double(r);
  free1double(x);
  free1double(u1);
  free1double(d1);
  free1double(u1b);
  free1double(uk); 
  free1double(dk); 
  free1double(storx); 
  free1double(Q);
  free1double(P);
  free1double(qr);
  free1double(pr);
  free1double(cn);
  free1double(t1);
  free1double(t2);
  free1double(d0);
  free1double(td1);
  free1double(mix1);
  free1double(mix2);
  free1double(mix3);
  free1double(mix4);
  free1double(mix);
}


/*--------------------------- functions used ---------------------------*/

/*************************************************************************
imin, imax - return the smaller/larger integer of two integers
--------------------------------------------------------------------------
Notes:  Not used in main(), only in the functions below. 
*************************************************************************/
int imin( int c1, int c2 )
{
  if(c1<=c2)
    return(c1);
  else
    return(c2);
}

int imax( int c1, int c2 )
{
  if(c1>=c2)
    return(c1);
  else
    return(c2);
}


/*************************************************************************
porder - return the order of a polynomial 
         (the highest power with a non-zero coefficient)
--------------------------------------------------------------------------
Input:
 p[]     array of polynomial coefficients
 l       polynomial length (length of p[]) 
--------------------------------------------------------------------------
Notes: If the array p[] contains coefficients of negative powers, 
       the result would be biased by the maximum negative power.
       Example: If p[] = a0*z^(-1) + a1*z^0 + a2*z^1 + a3*z^2 and a3=0, 
       the max. power in p[] is 1 but porder(...) will return 2 instead.	
--------------------------------------------------------------------------
Author:  CWP: Albena Mateeva   2000
-------------------------------------------------------------------------*/
int porder( double p[], int l )
{
  int i=0,k=0,j;
  
  do{
    j=k;
    if(p[l-1-i]==0)  k=k+1;
    ++i;
  } while(k>j);
  
  return(l-1-j);
}


/*************************************************************************
pcomb - linear combination of two polynomials
--------------------------------------------------------------------------
Input:
 alfa    scaling factor
 beta    scaling factor
 p1[]    array of polynomial coefficients
 p2[]    array of polynomial coefficients
 l1      length of p1[] 
 l2      length of p2[] 
 l       length of output array c[]

Output:
 c[]     array of polynomial coefficients c[i] = alfa*p1[i] + beta*p2[i]
--------------------------------------------------------------------------
Note: calling statements of the kind  pcomb(a,b, P ,Q,l1,l2, P ,l) allowed
--------------------------------------------------------------------------
Author:  CWP: Albena Mateeva  2000
------------------------------------------------------------------------*/
void pcomb( double alfa, double beta, double p1[], double p2[], 
	    int l1, int l2, double c[] , int l)
{
  int i,ii;
  double *tmp;
  
  ii=imax(l1,l2);

  /* Allocate and zero temporary array */ 
  tmp=ealloc1double(ii);
  memset( (void *) tmp, 0, ii*DSIZE);
  
  if (l1>=l2)
    for(i=l1-1; i>l2-1; i--)  tmp[i] = alfa*p1[i];
  else
    for(i=l2-1; i>l1-1; i--)  tmp[i] = beta*p2[i];

  for(i=0; i<imin(l1,l2); ++i)  tmp[i] = alfa*p1[i] + beta*p2[i];
  
  /* Clear output array */
  memset( (void *) c, 0, l*DSIZE);

  /* Output */
  for(i=0; i<=imin(ii,l-1); ++i)  c[i]=tmp[i];
  
  /* Free space */