gplib.f

linux环境下利用fortran语言开发

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

* These routines are from the gplib library as revised by K. Araya
* and Sebastien Geoltrain
*
* We went in and quickly trapped a bunch of potentially illegal values,
* however, this stuff should really be re-coded to be passing in
* array lengths instead of hard-wiring work areas, etc.  In the
* pre-release version, a bad value of nsweep in klaud86 could send
* the computations into limbo.  Credit to Darrell Kramer of Mobil on
* this and some related problems.  Jack Cohen, 08/04/89
*
*-------------------------------------------------------------------
* This subroutine performs the one-sided cross correlation of
* two functions by FFT, Bartlet truncation:
*
* Usage:
*     call corr86(x,lx,y,ly,z,lz)
*
*     x          the first array
*     lx         the length of the array x
*     y          the second array to move (*)
*     ly         the length of the array y
*     z          the output array
*     lz         length of the output array z
*
* Note:
*     lx + lz - 1 must be less than or equal to 8192.
*-----------------------------------------------------------------------
*
* Author:
* Source: /src/public/gplib
* Revision: Sebastien Geoltrain Jan.22,1989
*----------------------------------------------------------------------

      subroutine corr86(x,lx,y,ly,z,lz)

      real x(1),y(1),z(1),times
      integer lx,ly,lz,n,i,leng
      complex    w1(8192),w2(8192),w3(8192)

*    	check length parameters

	if(lx .LE. 0 .OR. ly .LE. 0) then
		stop ' corr86: illegal lx or ly'
	endif

*    	computes length of output as closest power of 2, not to 
*	exceed 2**13 (i.e. 8192).

	leng=lx+ly-1

        n=2

        do 10 i=1,13
           if(leng.le.n)then
              leng = n
              go to 1200
           else
              n = n * 2
           endif
10      continue

1200    continue

        if (lx .gt. 8192 .or. ly .gt. 8192 .or. leng .gt. 8192) then
	   stop ' corr86: illegal lx, ly or leng'
	endif

*	clears work arrays

	do 3000 i=1,leng
	   w1(i)=cmplx(0.0,0.0)
	   w2(i)=cmplx(0.0,0.0)
	   w3(i)=cmplx(0.0,0.0)
3000	continue

*    	moves input to work arrays

	do 4000 i=1,lx
	   w1(i)=cmplx(x(i),0.0)
4000	continue

	do 4100 i=1,ly
	   w2(i)=cmplx(y(i),0.0)
4100	continue

*	Fourier transforms both inputs

	call fft85(w1,leng,-1)
	call fft85(w2,leng,-1)

*     	takes complex conjugate of second input

	do 4200 i=1,leng
       	   w2(i)=conjg(w2(i))
4200	continue

*	multiplies first input and conjugate of second input 
*	in frequency domain

	do 4300	i=1,leng
	   w3(i)=w1(i)*w2(i)
4300	continue

*    	performs inverse FFT of resulting product

	call fft85(w3,leng,1)

*    	outputs result as real array of length lz

	times=sqrt(float(leng))

	do 4400	i=1,lz
	   z(i)=times*real(w3(i))
4400	continue

        return
	end
*--------------------------------------------------------------
* This subroutine computes the FFT of a function:
*
* Usage:
*     call fft85(cx,lx,sign)
*
*     cx         the input array (complex)
*                after the calling it is the output.
*     lx         the length of the input array
*     sign       (-1) forward transform
*                (1)  inverse transform
*-----------------------------------------------------------------
*
* Author: J. Claerbout
* Source: /src/public/gplib
* Revision: Sebastien Geoltrain Jan.22,1989
*-----------------------------------------------------------------
      subroutine fft85(cx,lx,sign)
      integer sign,lx,i,j,l,istep,m
      real sc
      complex  cx(lx),carg,cw,ctemp

      j=1
      sc=sqrt(1.e0/lx)

      do 30 i=1,lx

          if(i.gt.j) goto 10

          ctemp=cx(j)*sc
          cx(j)=cx(i)*sc
          cx(i)=ctemp

10        m=lx/2

20        if(j.le.m) goto 30

          j=j-m
          m=m/2
          if(m.ge.1) goto 20

30    j=j+m

      l=1

40    istep=2*l

      do 50 m=1,l

          carg=(0.e0,1.e0)*(3.14159265e0*sign*(m-1))/l
          cw=cexp(carg)

          do 50 i=m,lx,istep

              ctemp=cw*cx(i+l)
              cx(i+l)=cx(i)-ctemp

50    cx(i)=cx(i)+ctemp

      l=istep

      if(l.lt.lx) goto 40

      return
      end
*---------------------------------------------------------------------
* This subroutine generates a two-sided Klauder wavelet
* by a specified frequency array:
*
* Usage:
*     call klaud86(wave,leng,del,nsweep,fl,fh,ialfa)
*
*     wave       array
*     leng       the length of the generated Klauder
*                wavelet (should be an odd number)
*                hardwired at 255 in trisubs.getarray
*     del        sampling rate in seconds
*     nsweep     length of original vibroseis sweep (le. 8192)
*     fl         lower cut frequency
*     fh         higher cut frequency
*     ialfa      taper length (GE.2.and.le.500)
*-----------------------------------------------------------------------
*
* Author:
* Source: /src/public/gplib
* Revision: Sebastien Geoltrain Jan.22,1989
*----------------------------------------------------------------------

      subroutine klaud86(wave,leng,del,nsweep,fl,fh,ialfa)
      real wave(1),sweep(8192),w1(1024),fl,fh,del
      integer i,leng,lengs,lengh,lengh1,nsweep,ialfa,ifirst

*	parameter checks

*	if length is too big or too small, exit

 	if(leng .LE. 2 .OR. leng .GT. 1024) then
 		stop ' klauder86: leng illegal'
 	endif
 
*	if length is even, resets to closest smaller odd number

 	if(mod(leng,2).eq.0) then
 	      leng=leng-1
         endif
 
* 	if sampling rate is negative or null, exit

 	if(del .LE. 0.e0) then
 		stop ' klauder86: del illegal'
 	endif

*	if sweep is too long or too short, exit

	if(nsweep .LE. 0 .OR. nsweep .GT. 8192) then
		stop ' klauder86: nsweep illegal'
	endif

*	if taper is too small or too large, exit

	if(ialfa .LT. 2 .OR. ialfa .GT. 500) then
		stop ' klauder86: ialfa illegal'
	endif 

*	generates vibroseis sweep

	lengs = nsweep

        call vibro86(sweep,lengs,del,fl,fh,ialfa)

*    	crosscorrelates the sweep

	lengh=leng/2+1

	call corr86(sweep,lengs,sweep,lengs,w1,lengh)

	ifirst=2

	call norm85(w1,lengh,1,lengh,ifirst)

*	makes output into a two-sided function

        lengh1=lengh-1

        do 100 i=1,lengh1
              wave(lengh+i)=w1(i+1)
              wave(lengh-i)=w1(i+1)
100     continue

        wave(lengh)=1.e0

	return
	end
*--------------------------------------------------------------------------
* This subroutine normalizes an array by its RMS energy, or the
* first value in the array or by the largest magnitude in the 
* Array:
*
* Usage:
*     call norm85(x,lx,is,ie,in)
*
*     x		the input array
*     lx	the length of array
*     is        the initial point to generate factor
*     ie        the last point to generate factor
*     in        operation flag
*               = 1 normalizes by the RMS energy
*               = 2 normalizes by the first value
*               = 3 normalizes by the largest magnitude
*---------------------------------------------------------------------------
*
* Author:
* Source: /src/public/gplib
* Revision: Sebastien Geoltrain Jan.22,1989
*---------------------------------------------------------------------------

      subroutine norm85(x,lx,is,ie,in)

      integer lx,is,ie,in,i
      real x(lx),b,x1,p

*     test input parameters

      if(lx .LT. 1 .OR. is .GE. ie .OR. is .LE. 0 .OR. ie .GT. lx) then
	  stop ' norm85: illegal parameters'
      endif

      if(in .EQ. 1) then

               p=0.e0

               do 150 i=is,ie
                      p=p+x(i)*x(i)
150            continue

               p=sqrt(p)

               do 160 i=1,lx
                      x(i)=x(i)/p
160            continue

       else if(in .EQ. 2) then

               x1=x(is)

               if(x1 .NE. 0.e0) then

                      do 250 i=1,lx
                             x(i)=x(i)/x1
250                   continue

               endif

        else if(in .EQ. 3) then

               b=0.e0

               do 350 i=is,ie
                      b=amax1(abs(x(i)),b)
350            continue

               do 360 i=1,lx
                      x(i)=x(i)/b
360            continue

	endif

      return
      end
c------------------------------------------------------
	subroutine fasthilbert(p,q,ntrace)
	integer n,ntrace,i2,l1,l2
	real p(ntrace),q(ntrace),h(32),a,b

	call hilbertcoef(h)

	do5 i2=1,ntrace
5	q(i2) = 0.0
	do20 n=1,ntrace
	do30 i2=1,32,2
	l1=n-i2
	l2=n+i2
	if(l1.le.0.or.l1.gt.ntrace) then
	a = 0.0
	else
	a = p(l1)
	endif
	if(l2.le.0.or.l2.gt.ntrace) then
	b = 0.0
	else
	b = p(l2)
	endif
30	q(n) = q(n) + h(i2)*(a-b)
20	continue
	return
	end
c
	subroutine hilbertcoef(h)
	real h(32),opi,pi
	integer i2
	pi = 3.1415927
	opi = 2./pi
	do10 i2 = 1,32,2
10  	h(i2) = -opi/i2*(0.42+0.52*cos(pi*i2/32.)+.08*cos(pi*i2/16.))
	return
	end
c
*-----------------------------------------------------------
* This subroutine generates a vibroseis linear sweep:
*
* Usage:
*     call vibro86(sweep,leng,del,fa,fb,ltaper)
*
*     sweep       the generated sweep output array
*     leng        length of sweep: samples
*     del         sampling rate, in seconds
*     fa          starting cut-off frequency (no restriction)
*                 in hertz
*     fb          ending cut-off frequency (no restriction)
*                 in hertz
*     ltaper      sample length of 2 quadrant cosine taper
*                 (0 = no taper)
*--------------------------------------------------------------
*
* Author:
* Source: /src/public/gplib
* Revision: Sebastien Geoltrain Jan.22,1989
*--------------------------------------------------------------

      subroutine vibro86(sweep,leng,del,fa,fb,ltaper)

      integer i,leng,ltaper,l
      real sweep(leng),fa,fb,t,f,durat,fbad,pi,pi2,weight,del

      pi = 3.1415927e0

      pi2=pi*2.e0

*     tests input

      if(del .LT. 0.e0)then 
         stop ' vibro86: illegal del'
      endif

*     computes sweep

      durat=(leng-1.e0)*del
      fbad=(fb-fa)/(durat*2.e0)

      do 10 i=1,leng
         t = (i-1)*del
         f = fa+fbad*t
         sweep(i) = cos(pi2*f*t)
10    continue

*     applies taper

      do 20 i=1,ltaper-1
         l=leng+1-i
         weight=.5e0-.5e0*cos(pi*(i-1)/(ltaper-1) )
         sweep(i)=sweep(i)*weight
         sweep(l)=sweep(l)*weight
20    continue

      return
      end