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.   

 
 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)