C+
C	PROGRAM pz-test
C
C	Get amplitdue/phase displacement response for input poles and zeros
C	  entered in sac format
C-
	dimension freqcr(40),ampcr(40),phascr(40)
	character*80 filename, getstring
C
C	VTGDSN.CMN	VTSO version of the DTM SORCM.CMN header
C
	DIMENSION	HEADGD(78)
	INTEGER*4	NP,NZ,NZZ,ISTATB
	COMPLEX		P,Z
	COMMON/GDSNCM/
	1	A0,	!  counts/micrometer factor
	1	NP,	!  number of complex poles in response
	1	P(26),	!  complex poles of response function
	1	NZ,	!  number of complex zeroes in response
	1	Z(9),	!  complex zeroes of response function (first 9)
	1	NZZ,	!  number of complex zeroes which are zero
	1	ISTATB,	!  eleven status bits for record
	1	XLAT,	!  latitude of station
	1	XLON,	!  longitude of station
	1	ELEV	!  elevation of station in meters
	EQUIVALENCE (A0,HEADGD(1))
C
	complex D,S
	twopi = 8.0*atan(1.0)
	filename = getstring('Input p-z filespec')
	open(unit=1,name=filename,status='old')
	open(unit=2,name='pz-test.lst',status='unknown',form='print')
	  write(2,*) ' sac-format polezero input from file ',
     &		filename(1:lenc(filename))
	fzero = rvalue('Starting frequency..[0.01]',0.01)
	freqcr(1) = fzero
	freqcr(11) = 15.0*fzero
	freqcr(21) = 150.0*fzero
	freqcr(31) = 1500.0*fzero
	do j=1,9
	  freqcr(j+1) = fzero*(j+1)
	  freqcr(j+11) = 10.0*freqcr(j+1)
	  freqcr(j+21) = 100.0*freqcr(j+1)
	  freqcr(j+31) = 1000.0*freqcr(j+1)
	end do
	call polezero(1,2)
	    NZERO = MIN0(NZ,9+NZZ)
	    NPOLE = MIN0(NP,26)
	    npz = max0(npole,nzero)
	    nzmnzz = nz - nzz
	    NUMCOR = 40
	    DO J=1,numcor
	      OMEGA = twopi*FREQCR(J)
	      D = CMPLX(A0,0.0)
	      S = CMPLX(0.0,OMEGA)
	      if (npz .gt. 0) then
	        do k=1,npz
	          IF (k .le. nzero) THEN
	            IF (K .LE. nzmnzz) THEN
	              D = D*(S-Z(K))
	            ELSE
	              D = D*S
	            END IF
	          end if
	          if (k .le. npole) D = D/(S-P(K))
	        end do
	      end if
	      AMPCR(J) = CABS(D)
	      PHASCR(J) = ATAN2(-AIMAG(D),REAL(D))
C
C	My Fourier convention has the opposite sign from GDSN
C
	    end do
	    WRITE(2,'(1H0,T6,''Frequency'',T19,''Displ. Amp.'',
	1	T39,''Phase'')')
	    do j=1,numcor
	      if (j .gt. 1) then
	        do while (phascr(j) .lt. phascr(j-1))
	          phascr(j) = phascr(j) + twopi
	        end do
	      end if
	      write(2,'(F12.4,1pg15.4,0p2f15.4))') freqcr(j),ampcr(j),
     1			phascr(j),phascr(j)*57.29578
	    end do
	call exit(0)
	end
C+
	subroutine polezero(nin,nout)
C-
C
C	VTGDSN.CMN	VTSO version of the DTM SORCM.CMN header
C
	DIMENSION	HEADGD(78)
	INTEGER*4	NP,NZ,NZZ,ISTATB
	COMPLEX		P,Z
	COMMON/GDSNCM/
	1	A0,	!  counts/micrometer factor
	1	NP,	!  number of complex poles in response
	1	P(26),	!  complex poles of response function
	1	NZ,	!  number of complex zeroes in response
	1	Z(9),	!  complex zeroes of response function (first 9)
	1	NZZ,	!  number of complex zeroes which are zero
	1	ISTATB,	!  eleven status bits for record
	1	XLAT,	!  latitude of station
	1	XLON,	!  longitude of station
	1	ELEV	!  elevation of station in meters
	EQUIVALENCE (A0,HEADGD(1))
C
	logical more,another
	character*80 fcom
	real*4 twoval(2)
	complex compval
	equivalence(twoval(1),compval)
	character*80 test
	fcom(1:1) = ' '
	nz = 0
	nzz = 0
	np = 0
	a0 = 1.0
	another = .true.
	do while (another)
	  read(nin,'(a)',end = 100) test
	  call uppercase(test)
	  if (test(1:5) .eq. 'ZEROS') then
	    read(test(7:lenc(test)),*) nz
	    nzz = nz
	    more = .true.
	    j = 0
	    do while (more)
	      read(nin,'(a)',end=100) test
	      call uppercase(test)
	      if (test(1:5) .ne. 'POLES' .and. test(1:5) .ne.
     &		'CONST') then
	        read(test,*) twoval
	        j = j + 1
	        nzz = nzz - 1
	        Z(j) = compval
	      else
	        backspace nin
	        more = .false.
	      end if
	    end do
	  else if (test(1:5) .eq. 'POLES') then
	    read(test(7:lenc(test)),*) np
	    j = 0
	    more = .true.
	    do while (more .and. j .lt. 26)
	      read(nin,'(a)',end=100) test
	      call uppercase(test)
	      if (test(1:5) .ne. 'ZEROS' .and. test(1:5) .ne.
     &		'CONST') then
	        read(test,*) twoval
	        j = j + 1
	        P(j) = compval
	      else
	        backspace nin
	        more = .false.
	      end if
	    end do
	  else if (test(1:5) .eq. 'CONST') then
	    read(test(10:lenc(test)),*) a0
	  else
	    fcom = test
	  end if
	end do
100	close(unit=nin)
	if (np .eq. 0 .and. nz .eq. 0) then
	  write(*,*) 'No poles or zero'
	  stop
	end if
	if (nout .le. 0) return
	  if (fcom(1:1) .ne. ' ') write(nout,*) fcom(1:lenc(fcom))
          WRITE(nout,107) A0,NP
	  if (np .gt. 0) write(nout,*) (P(J),J=1,NP)
          WRITE(nout,108) NZ,NZZ
	  if ((nz-nzz) .gt. 0) write(nout,*) (Z(J),J=1,NZ)
107     FORMAT(5X,'A0 = ',1PG11.4,2X,I2,'  RESPONSE POLES: ')
108     FORMAT(5X,I2,' RESPONSE ZEROES  ',I2,' ZERO VALUES:')
	return
	end
	character*(*) function GETSTRING(prompt)

c+
c	character*(*) function GETSTRING(prompt)
c
c  outputs 'prompt' using PRINTX
c  and accepts input character string
c
c  also looks for supplied default
c	which is enclosed in '[]' at end of 'prompt'
c
c				Alan Linde ... Aug 1986
c-
	character*1 tab/9/
	character*(*) prompt
	character*80 temp

c output 'prompt'

	call printx(prompt)

	getstring=' '

c now check for supplied default

	kk=lenc(prompt)

	if (prompt(kk:kk).eq.']') then

	  ll=0
	  do 110 i=kk-1,1,-1
	    if (prompt(i:i).eq.'['.and.ll.eq.0) ll=i+1
110	  continue

	  if (ll.ne.0) getstring=prompt(ll:kk-1)

	end if

c  get the response

	read(*,'(a)') temp
	nout = lenc(temp)
	if (nout .gt. 0) then
	  ntab = index(temp(1:nout),tab)
	  if (ntab .gt. 0) nout = ntab - 1
	  if (nout.gt.0) getstring=temp(1:nout)
	end if

	return
	end
	function lenc(string)
C+
C	function lenc(string)
C
C	Returns length of character variable STRING excluding right-hand
C	  most blanks or nulls
C-
	character*(*) string
	length = len(string)	! total length
	if (length .eq. 0) then
	  lenc = 0
	  return
	end if
	if(ichar(string(length:length)).eq.0)string(length:length) = ' '
	do j=length,1,-1
	  lenc = j
	  if (string(j:j).ne.' ' .and. ichar(string(j:j)).ne.0) return
	end do
	lenc = 0
	return
	end
	SUBROUTINE PRINTX(LINE)
C+
c	SUBROUTINE PRINTX(LINE)
C  OUTPUTS A MESSAGE TO THE TERMINAL
C  PRINTX STARTS WITH A LINE FEED BUT DOES NOT END WITH A CARRIAGE RETURN
C  THE PRINT HEAD REMAINS AT THE END OF THE MESSAGE
C
C  IF THE MESSAGE LENGTH IS LESS THAN 40,
C	DOTS ARE INSERTED UP TO COL. 39
C	AND A COLON IS PUT IN COL. 40.
C
C  USE FOR CONVERSATIONAL INTERACTION
C			Alan Linde ... April 1980.
C	10 Sugust 1985:  Corrected a minor error for  strings > 40 bytes
C	20 June 1986:  Made it compatible with Fortran 77
C-
	character*(*) line
	CHARACTER*60 BUF
	CHARACTER*2 COLON
	CHARACTER*1 DOT,DELIM
	DATA DELIM/'$'/,DOT/'.'/,COLON/': '/
	KK = lenc(LINE)	!  length minus right-hand blanks
	  IF (LINE(KK:KK) .EQ. DELIM) KK = KK - 1
	  IF (KK .GT. 58) KK = 59
	BUF(1:KK) = LINE(1:KK)
	IF (KK .LT. 49) THEN
	  DO J=KK+1,49
	    BUF(J:J) = DOT
	  END DO
	  KK = 49
	END IF
	BUF(KK:KK+1) = COLON
	KK = KK + 2
	WRITE(*,'(A,$)') BUF(1:KK)
	RETURN
	END
C+
	REAL FUNCTION RVALUE(MSG,RDEF)
C
C PURPOSE:
C		THIS FUNCTION ACCEPTS A MESSAGE (ASKING FOR A VALUE)
C		AND RETURNS THE VALUE ENTERED AT THE TERMINAL
C ROUTINES CALLED:
C		PRINTX
C
C USE:
C		ANS=RVALUE('ENTER AN REAL NUMBER',RDEF)
C	If enter a cariage return, rvalue is set to RDEF
C
C AUTHOR:
C			ALAN LINDE ... AUGUST 1980 (for VALUE)
C
C EXTENSIONS:
C	30 JULY 1989:  CAN HANDLE ENTRY FOLLOWED BY A BLANK OR TAB
C-
	CHARACTER*1 TAB,E/'E'/,BLANK/' '/
	parameter (tab=char(9))
	CHARACTER*30 STUFF
	CHARACTER*(*) MSG
C
100	CALL PRINTX(MSG)
	READ(*,'(A)') STUFF
	nin = lenc(stuff)
	IF (NIN .GT. 0) THEN
	  NTAB = INDEX(STUFF(1:NIN),TAB)
	  IF (NTAB .GT. 0) NIN = NTAB - 1
	  NBLANK = INDEX(STUFF(1:NIN),BLANK)
	  IF (NBLANK .GT. 0) NIN = NBLANK - 1
	END IF
	IF (NIN .EQ. 0) THEN
	  RVALUE = RDEF
	ELSE
	  READ(STUFF(1:NIN),*,ERR=100) RVALUE
	END IF
	RETURN
	END