sugoupillaud.c

su 的源代码库

  free1double(tmp);
}


/*************************************************************************
rev - returns reverse-polinomial coefficients 
--------------------------------------------------------------------------
Input:
 p[]     array of polynomial coefficients
 l       length of p[]

Output:
 pr[]    array of polynomial coefficients 
         pr[0]=p[k], pr[1]=p[k-1], ... , pr[k]=p[0]
         where k is the order of the input polynomial p[]
--------------------------------------------------------------------------
Note:    calling statements of the kind  rev( P ,l, P )  allowed
--------------------------------------------------------------------------
Author:  CWP: Albena Mateeva  2000
------------------------------------------------------------------------*/
void rev(double p[],int l,double pr[])
{
  int i;
  double *tmp;
  int order;
  
  order = porder(p,l);

  /* Allocate and zero temporary array */  
  tmp = ealloc1double(order+1);
  memset( (void *) tmp, 0, (order+1)*DSIZE);
  
  for(i=0; i<=order; ++i)  tmp[i] = p[order-i];
  
  /* Clear output array */ 
  memset( (void *) pr, 0, l*DSIZE);
  
  /* Output */
  for(i=0; i<=order; ++i)  pr[i] = tmp[i];
  
  /* Free space */
  free1double(tmp);
}


/*************************************************************************
pshift - shift a z-polynomial of order n by k units 
--------------------------------------------------------------------------
Input:
 p[]     array of polynomial coefficients
 l       length of p[] 
 k       integer (has the meaning of time-shift in this code)
 ll      length of output array psh[]

Output:
 psh[]   the coefficients of a polynomial psh[] = z^(k)*p[]
         where p[] is the input polynomial in z
--------------------------------------------------------------------------
Notes: k should be such that the resulting polynomial has only (+) powers;
       calling statements of the kind	pshift( P ,l,k, P ,ll)  allowed.
--------------------------------------------------------------------------
Author:  CWP: Albena Mateeva  2000
------------------------------------------------------------------------*/
void pshift( double p[], int l, int k, double psh[], int ll )
{
  int ii,i;
  double *tmp;
  
  ii=imax(l+k,l);
  
  /* Allocate and zero temporary array */
  tmp = ealloc1double(ii);
  memset( (void *) tmp, 0, ii*DSIZE);
  
  if(k>=0)
    for (i=0; i<l; ++i)   tmp[i+k] = p[i];
  else
    for (i=-k; i<l; ++i)  tmp[i+k] = p[i];
  
  /* Clear output array */
  memset( (void *) psh, 0, ll*DSIZE);
  
  /* Output */
  for(i=0; i<imin(l,ll); ++i)  psh[i] = tmp[i];
  
  /* Free space */
  free1double(tmp);
}


/*************************************************************************
prod - polynomial multiplication
--------------------------------------------------------------------------
Input:
 p1[]    array of polynomial coefficients
 p2[]    array of polynomial coefficients
 l1      length of p1[]
 l2      length of p2[]
 ll      length of output array pp[]

Output:
 pp[]    array of polynomial coefficients
         pp is essentially the convolution of p1 and p2
--------------------------------------------------------------------------
Note:    calling statements of the kind  prod( P ,Q,l1,l2, P ,ll)  allowed
--------------------------------------------------------------------------
Author:  CWP: Albena Mateeva  2000
------------------------------------------------------------------------*/
void prod( double p1[], double p2[], int l1, int l2, double pp[], int ll )
{
  int i,j,k;	
  double *tmp;

  /* Allocate and zero temporary array */ 
  tmp = ealloc1double(l1+l2-1);
  memset( (void *) tmp, 0, (l1+l2-1)*DSIZE);
  
  for (i=0; i<l1; ++i){
    for (j=0; j<l2; ++j){
      k=i+j;
      tmp[k]=tmp[k]+p1[i]*p2[j];
    }
  }
  
  /* Clear output array */ 
  memset( (void *) pp, 0, ll*DSIZE);

  /* Output */
  for(i=0; i<imin(l1+l2-1,ll); ++i)   pp[i]=tmp[i];
  
  /* Free space */
  free1double(tmp);
}


/*************************************************************************
pratio - polynomial division
--------------------------------------------------------------------------
Input:
 p1[]    array of polynomial coefficients
 p2[]    array of polynomial coefficients
 l1      length of p1[]
 l2      length of p2[]
 L       desired number of output coefficients
 lq      length of output array q[]

Output:
 q[]     the first L coefficients of the polynomial q=p1/p2
--------------------------------------------------------------------------
Note:    see E. Robinson, Chap.1 for the algorithm (a Fortran subroutine)
--------------------------------------------------------------------------
Author:  CWP: Albena Mateeva  2000
------------------------------------------------------------------------*/
void pratio( double p1[], double p2[], int l1, int l2, int L, 
	     double q[], int lq )
{
  int i=0,k=0,j,ii,jj;
  int lo,o1,o2;
  double *tmp1, *tmp2;
  
  o1=porder(p1,l1);
  o2=porder(p2,l2);
  
  /* Allocate temporary arrays */
  tmp1 = ealloc1double(L+o2+1);
  tmp2 = ealloc1double(o2+1);
  
  /* Zero temporary arrays */
  memset( (void *) tmp1, 0, (L+o2+1)*DSIZE);
  memset( (void *) tmp2, 0, (o2+1)*DSIZE);
  
  do{
    jj=k;
    if(p2[i]==0) k=k+1;
    ++i;
  } while(k>jj);
  
  pshift(p2,l2,-jj,tmp2,o2+1);
  
  lo=imin(o1,L+o2);
  
  for (i=0; i<=lo; ++i)  tmp1[i] = p1[i];
  
  i=0;
  do{
    tmp1[i]=tmp1[i]/tmp2[0];
    if(i==L+jj)
      break;
    
    k=i;
    ii=imin(o2,L+jj-i);
    j=0;
    for(j=0; j<ii; j++){
      k=k+1;
      tmp1[k]=tmp1[k]-tmp1[i]*tmp2[j+1];
    }
    ++i;
  } while(1);
  
  pshift(tmp1,L+o2+1,-jj,tmp1,L+o2+1);
  
  /* Clear output array */
  memset( (void *) q, 0, lq*DSIZE);
  
  /* Output */
  for(i=0; i<imin(L,lq); ++i) q[i]=tmp1[i];
  
  /* Free space */
  free1double(tmp1);
  free1double(tmp2);
}


/*************************************************************************
recurs - compute recursive polynomials Q_k, P_k
--------------------------------------------------------------------------
Input:
 k	 polynomial order; k>=1
 n	 number of subsurface interfaces in a Goupillaud medium 
	 (horizontal interfaces equally spaced in time);
 r[]	 reflection coefficient series of length n+1;

Output:
 q[]	 the coefficients of a polynomial Q of order k
 p[]	 the coefficients of a polynomial P of order k
--------------------------------------------------------------------------
Note:    see E. Robinson, p.124-125
--------------------------------------------------------------------------
Author:  CWP: Albena Mateeva  2000
------------------------------------------------------------------------*/
void recurs( int k, double q[], double p[], int n, double r[])
{ 
  int i,j;
  double *tq, *qr, *qsh;	
  double *tp, *pr, *psh;
  
  /* Allocate temporary arrays */
  tq = ealloc1double(n+1);
  tp = ealloc1double(n+1);
  qr = ealloc1double(n+1);
  pr = ealloc1double(n+1);
  qsh= ealloc1double(n+1);
  psh= ealloc1double(n+1);
  
  /* Zero temporary arrays */
  memset( (void *) tq, 0, (n+1)*DSIZE);
  memset( (void *) tp, 0, (n+1)*DSIZE);
  memset( (void *) qr, 0, (n+1)*DSIZE);
  memset( (void *) pr, 0, (n+1)*DSIZE);
  memset( (void *) qsh, 0, (n+1)*DSIZE);
  memset( (void *) psh, 0, (n+1)*DSIZE);
  
  /* initial values for the recursion: Q1, P1 */
  tp[0]=1;				
  tq[1]=-r[1];
  
  /* recursive calculation of Qk, Pk */
  for(j=2; j<=k; j++) 
    {
      rev(tq,n+1,qr);
      pshift(qr,n+1,(j-1-porder(tq,n+1)),qr,n+1);	
      /* qr - reverse Q polynomial */
      rev(tp,n+1,pr);
      pshift(pr,n+1,(j-1-porder(tp,n+1)),pr,n+1);	
      /* pr - reverse P polynomial */
      
      pshift(qr,n+1,1,qsh,n+1);
      pshift(pr,n+1,1,psh,n+1);
      
      pcomb(1.,-r[j],tq,psh,n+1,n+1,tq,n+1);	/* now tq = Qj */
      pcomb(1.,-r[j],tp,qsh,n+1,n+1,tp,n+1);	/* now tp = Pj */
      
    }
  
  /* return Qk, Pk */
  for(i=0; i<n+1; ++i)
    {
      q[i]=tq[i];
      p[i]=tp[i];
    }
  
  /* Free space */
  free1double(tq);
  free1double(tp);
  free1double(qr);
  free1double(pr);
  free1double(qsh);
  free1double(psh);
}


/*************************************************************************
upsample - resample at half sampling interval
--------------------------------------------------------------------------
Input:
 p4[]    array of coefficients of a polynomial in Z
 l       length of p4[]
 ll      length of output array p2[]

Output:
 p2[]    input array p4[] padded by a zero at every other sample,
         i.e., the output array p2[] represents a polynomial in z=sqrt(Z)
--------------------------------------------------------------------------
Author:  CWP: Albena Mateeva  2000
------------------------------------------------------------------------*/
void upsample( double p4[], int l, double p2[], int ll )
{
  int i=0, ii;
  double *tmp;
  
  ii=imax(ll,2*l);
  
  tmp = ealloc1double(ii);
  memset( (void *) tmp, 0, ii*DSIZE);
  
  do{
    tmp[2*i]=p4[i];
    ++i;
  } while(i<l);
  
  for(i=0;i<ll;++i) p2[i]=tmp[i];
  
  free1double(tmp);
}