

c  the definition used here for acosh can be found in crc std. math. tab
cles,
c  25th ed., p. 261.  the definition is valid for x >= 1.
c  x:  real number for which acosh is computed.  x should be >= 1.  if n
cot,
c      an error condition results.
c  ipv:  integer indicating whether the principal value of the 
c        function is to be used or not.  if ipv is >= 0, p.v. is compute
cd.
cx is in proper range                         
      function acosh(x, ipv)
      if (x .ge. 1.) then
      sgn = 1.
      if (ipv .lt. 0) sgn = - sgn
      acosh = log(x + (sgn * sqrt((x ** 2) - 1.)))
      else
      write(unit=6, fmt=100) 
      stop 
      end if
  100 format(32h argument of acosh must be >= 1.)
      return 
      end
c
c     this subroutine adds a value to a linked list.
c
      subroutine addlnk(nstart, new, list)
c
      dimension list(1)
      list(new) = 0
      if (nstart .ne. 0) goto 10
      nstart = new
      goto 999
   10 continue
      n = nstart
   20 continue
      if (list(n) .eq. 0) goto 30
      n = list(n)
      goto 20
   30 continue
      list(n) = new
  999 continue
      return 
      end
c************************* arg.for ***********************************
c
c function to determine the phase (argument) of a complex number.
c
c*********************************************************************
      function arg(c)
      complex c
      arg = atan2(aimag(c),real(c))
      return 
      end
c
c*** function atan2pi takes a -pi to pi atan2 and changes to 0 to 2pi
c
      function atan2pi(b, a)
      double precision dpi
c
      data dpi / 3.14159265358979 /
      arg = atan2(b,a)
      if (arg .ge. 0.) then
      atan2pi = arg
      else
      atan2pi = (2. * dpi) + arg
      end if
      return 
      end
c  the definition used here for atanh can be found in crc std. math. tab
cles,
c  25th ed., p. 261.  the definition is valid for x >= 1.
c  x:  real number for which atanh is computed.  x**2 should be < 1.  if
c not,
c      an error condition results.
      function atanh(x)
cx is in proper range                      
      if ((x ** 2) .lt. 1.) then
      atanh = log((1 + x) / (1 - x)) / 2.
      else
      write(unit=6, fmt=100) 
      stop 
      end if
  100 format(47h argument of atanh squared must be less than 1.)
      return 
      end
c
      subroutine bufferin(iunit, k, array, nwrds, ind)
c
      dimension array(1)
c implies the read was succesful         
      ind = 0
      read(unit=iunit, end=100) (array(j),j = 1, nwrds)
      goto 999
  100 continue
cend of file encountered                         
      ind = 1
  999 continue
      return 
      end
c
      subroutine bufferou(iunit, k, array, nwrds, ind)
c
      dimension array(1)
      write(unit=iunit) (array(j),j = 1, nwrds)
      return 
      end
c
      complex function ccdotp(cva, cvb)
c
      complex cva(3), cvb(3)
c
      ccdotp = ((cva(1) * cvb(1)) + (cva(2) * cvb(2))) + (cva(3) * cvb(3
     &))
      return 
      end
      complex function cexpdot(xk, v, w)
c
      dimension v(3), w(3)
      vdotr = dotp(v,w)
      cexpdot = cexp(((0.,1.) * xk) * vdotr)
      return 
      end
      subroutine cntlnk(nstart, link, num)
c
c         this subroutine returns the number of links in the
c         linked list  link using the initial pointer  nstart.
c
      dimension link(1)
      n = nstart
      i = 0
      do 10 while (n .ne. 0)
      i = i + 1
      n = link(n)
   10 continue
      num = i
      return 
      end
c
      complex function crdotp(cva, vb)
      complex cva(3)
c
      dimension vb(3)
      crdotp = ((cva(1) * vb(1)) + (cva(2) * vb(2))) + (cva(3) * vb(3))
      return 
      end
c          kilroy was here
c
c      include 'pscrossp.inc/nolist'
c
      subroutine crossp(x, y, z)
c
      real x(1), y(1), z(1)
      z(1) = (x(2) * y(3)) - (x(3) * y(2))
      z(2) = (x(3) * y(1)) - (x(1) * y(3))
      z(3) = (x(1) * y(2)) - (x(2) * y(1))
      return 
      end
c
      subroutine crveclin(ca, va, cb, vb, cvc)
      complex ca, cb, cvc(3)
c
      dimension va(3), vb(3)
      do k = 1, 3
      cvc(k) = (ca * va(k)) + (cb * vb(k))
      end do
      return 
      end
c
      complex function ctan(x)
c
      complex x
c
      ctan = csin(x) / ccos(x)
      return 
      end
c
      function dblog10(x)
      if (x .le. 0.) then
      dblog10 = -200.
      else
      dblog10 = 10. * alog10(x)
c
      end if
      return 
      end
c
      function dblog20(x)
      if (x .le. 0.) then
      dblog20 = -200.
      else
      dblog20 = 20. * alog10(x)
c
      end if
      return 
      end
c
      function ddotp(x, y)
      double precision x(1), y(1), ddotp
      ddotp = ((x(1) * y(1)) + (x(2) * y(2))) + (x(3) * y(3))
      return 
      end
c
c     this subroutine deletes a value from a linked list.
c
      subroutine dellnk(nstart, nold, list, fail)
      dimension list(1)
c
      logical fail
      fail = .false.
      if (nold .ne. nstart) goto 10
      nstart = list(nold)
      goto 999
   10 continue
      n = nstart
      do 25 while (n .ne. 0)
      if (nold .eq. list(n)) goto 30
      n = list(n)
   25 continue
      fail = .true.
      goto 999
   30 continue
      list(n) = list(nold)
  999 continue
      return 
      end
c
      subroutine dnormze(x)
      double precision xl, x(3)
      xl = ((x(1) ** 2) + (x(2) ** 2)) + (x(3) ** 2)
      xl = dsqrt(xl)
      do 50 j = 1, 3
      x(j) = x(j) / xl
   50 continue
      return 
      end
c
c      include 'psdotp.inc/nolist'
c
      function dotp(x, y)
      real x(3), y(3)
      dotp = ((x(1) * y(1)) + (x(2) * y(2))) + (x(3) * y(3))
      return 
      end
      subroutine fft(l, x, y)
      dimension is(13), iu(13)
      dimension x(256), y(256)
      equivalence (is(2), ibs), (is(3), ics), (is(4), ids), (is(5), ies)
     &, (is(6), ifs), (is(7), igs), (is(8), ihs), (is(9), iis), (is(10)
     &, ijs), (is(11), iks), (is(12), ils), (is(13), ims)
      equivalence (iu(1), ial), (iu(2), ibl), (iu(3), icl), (iu(4), idl)
     &, (iu(5), iel), (iu(6), ifl), (iu(7), igl), (iu(8), ihl), (iu(9), 
     &iil), (iu(10), ijl), (iu(11), ikl), (iu(12), ill), (iu(13), iml)
c
      data twopi / 6.283185 /
      l2n = iabs(l)
      n = 2 ** l2n
      if (l) 6, 6, 8
    6 do 7 j = 1, n
c
    7 y(j) = - y(j)
    8 if (l2n - 1) 15, 15, 9
    9 m = n
      do 12 k = 2, l2n, 2
      m4 = m
      m = m4 / 4
      f = twopi / m4
      do 12 j = 1, m
      theta = f * (j - 1.)
      c1 = cos(theta)
      s1 = sin(theta)
      c2 = (c1 * c1) - (s1 * s1)
      s2 = (2. * c1) * s1
      c3 = (c2 * c1) - (s2 * s1)
      s3 = (c2 * s1) + (s2 * c1)
      jmm4 = j - m4
      do 12 i = m4, n, m4
      j1 = i + jmm4
      j2 = j1 + m
      j3 = j2 + m
      j4 = j3 + m
      r1 = x(j1) + x(j3)
      r2 = x(j1) - x(j3)
      r3 = y(j1) + y(j3)
      r4 = y(j1) - y(j3)
      r5 = x(j2) + x(j4)
      r6 = x(j2) - x(j4)
      r7 = y(j2) + y(j4)
      r8 = y(j2) - y(j4)
      x(j1) = r1 + r5
      y(j1) = r3 + r7
      if (j - 1) 10, 11, 10
   10 tmp1 = r2 + r8
      tmp2 = r4 - r6
      x(j3) = (tmp1 * c1) + (tmp2 * s1)
      y(j3) = (tmp2 * c1) - (tmp1 * s1)
      tmp1 = r1 - r5
      tmp2 = r3 - r7
      x(j2) = (tmp1 * c2) + (tmp2 * s2)
      y(j2) = (tmp2 * c2) - (tmp1 * s2)
      tmp1 = r2 - r8
      tmp2 = r4 + r6
      x(j4) = (tmp1 * c3) + (tmp2 * s3)
      y(j4) = (tmp2 * c3) - (tmp1 * s3)
      goto 12
   11 x(j3) = r2 + r8
      y(j3) = r4 - r6
      x(j2) = r1 - r5
      y(j2) = r3 - r7
      x(j4) = r2 - r8
      y(j4) = r4 + r6
c
   12 continue
   15 if (l2n - (2 * (l2n / 2))) 16, 18, 16
   16 do 17 i = 1, n, 2
      r1 = x(i) + x(i + 1)
      r2 = x(i) - x(i + 1)
      r3 = y(i) + y(i + 1)
      r4 = y(i) - y(i + 1)
      x(i) = r1
      y(i) = r3
      x(i + 1) = r2
   17 y(i + 1) = r4
   18 ims = n / 2
      iml = n
      do 20 j = 12, 2, -1
      is(j) = 1
      if (is(j + 1) - 1) 20, 20, 19
   19 is(j) = is(j + 1) / 2
c
   20 iu(j) = is(j + 1)
      ial = ibs
      jt = 0
      do 22 ia = 1, ial
      do 22 ib = ia, ibl, ibs
      do 22 ic = ib, icl, ics
      do 22 id = ic, idl, ids
      do 22 ie = id, iel, ies
      do 22 ifi = ie, ifl, ifs
      do 22 ig = ifi, igl, igs
      do 22 ih = ig, ihl, ihs
      do 22 ii = ih, iil, iis
      do 22 ij = ii, ijl, ijs
      do 22 ik = ij, ikl, iks
      do 22 il = ik, ill, ils
      do 22 im = il, iml, ims
      jt = jt + 1
      if (jt - im) 22, 22, 21
   21 t = x(jt)
      x(jt) = x(im)
      x(im) = t
      t = y(jt)
      y(jt) = y(im)
      y(im) = t
c
   22 continue
      if (l) 35, 35, 37
   35 do 36 j = 1, n
   36 y(j) = - y(j)
   37 return 
      end
c
      function filenm(infile, dir, pre, ityp)
      character infile*30, dir*30, ext(2)*4, pre*1, filenm*(*), tfile*30
     &, blnk*30
c
      data ext / '.ttj', '.rtj' /
      data blnk / '   ' /
   50 continue
      call fndblk(dir, 30, n)
      call fndblk(infile, 30, k)
      if (ityp .gt. 0) then
      if (n .gt. 1) then
      i = ((n + 1) + k) + 4
      tfile(:i) = ((dir(:n) // pre) // infile(:k)) // ext(ityp)
      else
      i = (1 + k) + 4
      tfile(:i) = (pre // infile(:k)) // ext(ityp)
      end if
      else
      i = (n + k) + 4
      tfile(:i) = (dir(:n) // infile(:k)) // '.dat'
      end if
      filenm = tfile(:i) // blnk(i + 1:30)
      return 
      end
c*********************************************************************
c
c purpose: to determine the position of the 1st blank in a character
c          string.
c
c*********************************************************************
c
      subroutine fndblk(string, maxchr, kend)
c
      character string*(*)
      kend = maxchr
      do 30 k = 1, maxchr
      if (string(k:k) .eq. ' ') then
      kend = k - 1
      goto 50
      end if
c
   30 continue
   50 continue
      return 
      end
      subroutine fndlnk(nstart, link, i, inx)
c
c         nstart - initial pointer into the linked list
c         link   - the linked list
c         i      - logical position of the link to be found
c         inx    - index of ith link
c                - 0 if no ith link is found
c
      dimension link(1)
      n = nstart
      j = 1
      do 10 while (j .ne. i)
      if (n .eq. 0) goto 20
      n = link(n)
      j = j + 1
   10 continue
      inx = n
      goto 30
   20 continue
      inx = 0
   30 continue
      return 
      end
c
c          n  =  the number of weights to be created
c          nrv=  1 or -1; if nrv is -1 the fourier
c                transform of the function will be
c                reordered
c          w  =  the parray where the weights are stored
c
      subroutine hamwgt(n, nrv, w)
      real pi, h1, h2
c
c
c          hamming weight functtion
c               this function can be used for
c               both odd and even values of n
c
      dimension w(256)
      data pi / 3.14159265358979 /
      data h1 / .54 /
      data h2 / .46 /
      no2 = n / 2
c
c       this branch is taken if n
c       is even
c
      s = 1
c
c          this section of the subroutine is
c          for odd values of n
c
      if ((no2 * 2) .eq. n) goto 50
      m = no2 + 1
c
c
      p = (2 * pi) / (n - 1)
      w(m) = 1
      do 40 i = 1, no2
      s = s * nrv
      x = s * (h1 + (h2 * cos(p * i)))
      w(m + i) = x
c
   40 w(m - i) = x
c
c          this section of the subroutine is
c          for even values of n
c
      return 
   50 continue
      p = (2 * pi) / n
      do 60 i = 1, no2
      s = s * nrv
      x = h1 + (h2 * cos(p * i))
      w((no2 + 1) - i) = (nrv * s) * x
   60 w(no2 + i) = s * x
      return 
      end
      subroutine inslnk(nstart, link, i, inx)
c
c         this routine inserts a link of index  inx  into
c         the linked list  link  at logical position  i.
c         nstart - initial pointer into linked list
c         list   - linked list
c         i      - logical position of new link
c         inx    - index of new link
c
      dimension link(1)
      if (i .eq. 1) goto 10
      call fndlnk(nstart, link, i - 1, in)
      link(inx) = link(in)
      link(in) = inx
      goto 20
   10 continue
      link(inx) = nstart
      nstart = inx
   20 continue
      return 
      end
c
c---- subroutine to reads from the terminal an integer
c
      subroutine itmin1(i)
c
      data iu1 / 5 /
      data nu1 / 6 /
      ierr = 0
   50 continue
      read(unit=iu1, fmt=100, err=75) i
      goto 90
   75 continue
      ierr = ierr + 1
      if (ierr .lt. 5) then
      write(unit=nu1, fmt=110) 
      goto 50
      else
      stop 
      end if
   90 continue
  100 format(i5)
  110 format(26h invalid input, try again:,$)
      return 
      end
c
c--- subroutine to read two integer values from the terminal.
c
      subroutine itmin2(i, j)
      data iu1 / 5 /
      data nu1 / 6 /
      ierr = 0
   50 continue
      read(unit=iu1, fmt=100, err=75) i, j
      goto 90
   75 continue
      ierr = ierr + 1
      if (ierr .lt. 5) then
      write(unit=nu1, fmt=110) 
      goto 50
      else
      stop 
      end if
   90 continue
  100 format(2i5)
  110 format(26h invalid input, try again:,$)
      return 
      end
c
c--- subroutine to read 3 integer values from the terminal.
c
      subroutine itmin3(i, j, k)
c
      data iu1 / 5 /
      data nu1 / 6 /
      ierr = 0
   50 continue
      read(unit=iu1, fmt=100, err=75) i, j, k
      goto 90
   75 continue
      ierr = ierr + 1
      if (ierr .lt. 5) then
      write(unit=nu1, fmt=110) 
      goto 50
      else
      stop 
      end if
   90 continue
  100 format(3i5)
  110 format(26h invalid input, try again:,$)
      return 
      end
c
c-- subroutine to read 4 integer values from the terminal.
c
      subroutine itmin4(i, j, k, l)
c
      data iu1 / 5 /
      data nu1 / 6 /
      ierr = 0
   50 continue
      read(unit=iu1, fmt=100, err=75) i, j, k, l
c
      goto 90
   75 continue
      ierr = ierr + 1
      if (ierr .lt. 5) then
      write(unit=nu1, fmt=110) 
      goto 50
      else
      stop 
c
      end if
   90 continue
  100 format(4i5)
  110 format(27h invalid input, try again: ,$)
      return 
      end
c******************* matvec.for **************************************
c
c purpose: to multiply a (3x3) matrix and a vector.
c
c*********************************************************************
c
      subroutine matvec(aa, vecin, vecout)
c
      dimension aa(3, 3), vecin(3), vecout(3)
      do k = 1, 3
      vecout(k) = ((aa(1,k) * vecin(1)) + (aa(2,k) * vecin(2))) + (aa(3,
     &k) * vecin(3))
c
ck=1,3                                                    
      end do
      return 
      end
c
c      include 'psnormze.inc/nolist'
c
      subroutine normze(x)
      real x(3)
      double precision xl, xx(3)
      do i=1,3
        xx(i) = x(i)
      end do

      xl = ((xx(1) ** 2) + (xx(2) ** 2)) + (xx(3) ** 2)
      xl = dsqrt(xl)

      if(xl .ne. 0.) then
        do j = 1, 3
          x(j) = xx(j) / xl
        end do
      else
        do j=1,3
          x(j) = 0.
        end do
      endif

      return 
      end




      subroutine ofrd(fil,ifilnum)
c
c This routine opens a file for reading.  The name of the file to
c be opened is passed in in the parameter 'fil'.  The unit number
c of the opened file is passed back in the parameter 'ifilnum'.  If
c the file 'fil' does not exist, the user will be prompted for the 
c name of an existing file.  If the named file does not exist, a
c directory with the entered mask is displayed.  The user is prompted
c until the name of a valid file is entered.  This "directory masking"
c allows the user to enter names lik "*.dat" to display all files with
c the .dat extension prior to entering a file name.
c
c NOTE: File unit numbers are checked beginning at unit 11.  Thus
c       multiple calls to this routine will produce sequential
c       unit numbers beginning at 11.
c
c UNIX version does not display directory of file names as lib$spawn is not 
c     available
c
      character fil*(*), command*80
      logical ex

10    continue
      ifilnum = 11 + irep
c     Check to see if file exists.
      inquire (FILE=FIL,EXIST=EX)
      if (ex) then
c       If the file exists, we are done.  Open file and return.
        OPEN (IFILNUM, FILE=FIL, STATUS='OLD', READONLY)
      else
c       If file does not exist, we need to prompt for a new file name.
        do while (.not. ex)
           write(6,*) 'Enter name of input data file'
           read (5,'(A30)') fil
c          Check the entered file to see if it exists.
           inquire (FILE=FIL,EXIST=EX)
           if (ex) then
c             If file exists we're done.
              OPEN (IFILNUM, FILE=FIL, STATUS='OLD', READONLY)
c           else
c            If the file does not exist, assume the user entered a
c            directory mask.
c             write(6,*)
c             command = 'dir '
c             write(command(5:80),'(A40)') fil
c             call lib$spawn(command)
c             write(6,*)
           endif
        enddo
      endif
      irep = irep + 1
      GOTO 200

100   irep = irep + 1
c      goto 10
      write(6,*) 'Error opening file'
200   CONTINUE
      RETURN
1400  FORMAT (I5)
1500  FORMAT (2F10.3)
      end
  
      



      subroutine ofwt(fil,ifilnum)
c
c This routine opens a new file for writing.  The name of the file to
c be opened is passed in in the parameter 'fil'.  The unit number
c of the opened file is passed back in the parameter 'ifilnum'.  If
c the file 'fil' exists, the user will be prompted for the 
c name of a new file.  If the entered file name is blank (user entered
c only a <CR>), a directory will be displayed.  If the entered name
c is an existing file, the user is allowed to either over-write the
c file or choose a new file name.  If the user enters a file name that
c cannot be opened, a directory is displayed with the entered name as
c the mask.  This allows the user to enter things like "*.dat" to
c get a listing of all files with the .dat extension prior to entering
c a file name.
c
c NOTE: File unit numbers are checked beginning at unit 60.  Thus
c       multiple calls to this routine will produce sequential
c       unit numbers beginning at 60.
c      
c UNIX version does not display a directory of file names since
c    lib$spawn is not available
c
      character fil*(*), command*80, resp*1
      logical ex

      ifilnum = 60 + irep
      icont = 0
      do while (icont .eq. 0)
        write(6,*) 'Enter name of output data file'
        read (5,'(A30)') fil
        if(fil .eq. ' ') then
c          call lib$spawn('dir')
        else 
          inquire(FILE=fil,EXIST=ex)

          if(ex) then
c  file exists
            write(6,*) ' File already exists. Overwrite [Y]'
            read(5,'(1A)') resp
            if(resp .ne. 'N'  .and.  resp .ne. 'n') then
              icont = 1
              open(ifilnum, FILE=fil, STATUS='new')
            endif
        
          else
c  file does not exist
            open(ifilnum, FILE=fil, STATUS='new', IOSTAT=ios)
            if(ios .eq. 0) then
c  file succesfully opened
              icont = 1
c            else
c  file could not be opened...must want a directory
c              command = 'dir '
c              write(command(5:30),'(a25)') fil
c              call lib$spawn(command)
            endif
          endif
        endif
      enddo
      irep = irep + 1
      RETURN
      END
c*********************************************************************
c
c   purpose:  to read in a file name and open a formatted file 
c             with the given unit #.
c
c            iotyp :: 1 => read
c                     2 => write
c
c            iform :: 1 => formatted
c                     2 => unformatted
c
c*********************************************************************
c
      subroutine openfil(iunit, iotyp, iform, ir1, iw1)
      character fname*30
c
      logical ex
c write                                 
      if (iotyp .eq. 2) then
   10 continue
      write(unit=iw1, fmt=1000) 
      read(unit=ir1, fmt=1010) fname
      inquire(file=fname, exist=ex) 
      if (.not. ex) then
      if (iform .eq. 1) then
      open(unit=iunit, file=fname, status='new', form='formatted') 
      else
      open(unit=iunit, file=fname, status='new', form='unformatted') 
      end if
      else
      write(unit=iw1, fmt=1020) 
      read(unit=ir1, fmt=1030) nansw
      if (nansw .eq. 'y') then
      if (iform .eq. 1) then
      open(unit=iunit, file=fname, status='new', form='formatted') 
      else
      open(unit=iunit, file=fname, status='new', form='unformatted') 
      end if
      else
      goto 10
      end if
      end if
c read                                                 
      else
   20 continue
      write(unit=iw1, fmt=1005) 
      read(unit=ir1, fmt=1010) fname
      inquire(file=fname, exist=ex) 
      if (.not. ex) then
      write(unit=iw1, fmt=1025) 
      goto 20
      end if
      if (iform .eq. 1) then
      open(unit=iunit, file=fname, status='old', form='formatted') 
      else
      open(unit=iunit, file=fname, status='old', form='unformatted') 
      end if
c
      end if
 1000 format(24h enter output file name )
 1005 format(23h enter input file name )
 1010 format(a30)
 1020 format(49h file already exists -- do you want to write over,
     &10h it (y/n) )
 1025 format(34h file does not exist -- try again )
c
 1030 format(a1)
      return 
      end
c
      subroutine rndoff(ymax, ymin)
c
      data up / .99999 /
      data dwn / 0. /
      npowr = alog10(ymax - ymin)
      if (ymin .gt. 0.) then
      ymax = rnd(ymax,npowr,up)
      ymin = rnd(ymin,npowr,dwn)
      else
      if (ymax .lt. 0.) then
      ymax = - rnd(- ymax,npowr,dwn)
      ymin = - rnd(- ymin,npowr,up)
      else
      if (ymin .lt. 0.) ymin = - rnd(- ymin,npowr,up)
      if (ymax .gt. 0.) ymax = rnd(ymax,npowr,up)
      end if
      end if
      return 
      end
c
      function rnd(y, npowr, trnd)
      im = (y / (10. ** npowr)) + trnd
      rnd = im * (10. ** npowr)
      return 
      end
c
c--- subroutine to read a real value from the terminal.
c
      subroutine rtmin1(x)
c
      data iu1 / 5 /
      data nu1 / 6 /
      ierr = 0
   50 continue
      read(unit=iu1, fmt=100, err=75) x
      goto 90
   75 continue
conly allows 5 invalid entries          
      ierr = ierr + 1
      if (ierr .lt. 5) then
      write(unit=nu1, fmt=110) 
      goto 50
      else
      stop 
      end if
   90 continue
  100 format(f10.0)
  110 format(26h invalid input, try again:,$)
      return 
      end
c
c-- subroutine to read 2 real values from the terminal.
c
      subroutine rtmin2(x, y)
c
      data iu1 / 5 /
      data nu1 / 6 /
      ierr = 0
   50 continue
      read(unit=iu1, fmt=100, err=75) x, y
      goto 90
   75 continue
      ierr = ierr + 1
      if (ierr .lt. 5) then
      write(unit=nu1, fmt=110) 
      goto 50
      else
      stop 
      end if
   90 continue
  100 format(2f10.0)
  110 format(26h invalid input, try again:,$)
      return 
      end
c
c-- subroutine to read 3 real values from the terminal.
c
      subroutine rtmin3(x, y, z)
c
      data iu1 / 5 /
      data nu1 / 6 /
      ierr = 0
   50 continue
      read(unit=iu1, fmt=100, err=75) x, y, z
      goto 90
   75 continue
      ierr = ierr + 1
      if (ierr .lt. 5) then
      write(unit=nu1, fmt=110) 
      goto 50
      else
      stop 
      end if
   90 continue
  100 format(3f10.0)
  110 format(26h invalid input, try again:,$)
      return 
      end
      function sinc(x)
      if (abs(x) .lt. .0001) sinc = 1.
      if (abs(x) .lt. .0001) goto 10
      sinc = sin(x) / x
   10 continue
      return 
      end
c
      subroutine trnslt(va, x, y, z, vb)
      dimension va(3), x(3), y(3), z(3), vb(3)
c
      double precision va, vb
      vb(1) = ((va(1) * x(1)) + (va(2) * y(1))) + (va(3) * z(1))
      vb(2) = ((va(1) * x(2)) + (va(2) * y(2))) + (va(3) * z(2))
c
      vb(3) = ((va(1) * x(3)) + (va(2) * y(3))) + (va(3) * z(3))
      return 
      end
c*********************************************************************
c
c  purpose: to test the existance of a file.
c
c*********************************************************************
c
      logical function tstex(infile)
      character infile*30
c
      logical ex
c
      inquire(file=infile, exist=ex) 
      tstex = ex
      return 
      end
c**********************************************************************
c
c   program to take a string and return the same string only in all
c   uppercase letters.
c
c   inputs :  string - character string
c             nchar - # of characters to check
c
c   output : string - character string with all upper case letters.
c
c   using the ascii values for each of the characters determine if the
c   character is a lower case letter (97-122), if it is replace it with
c   the upper case character (65-90)
c
c**********************************************************************
c
      subroutine upcase(string, nchar)
c
      character string*80
      do k = 1, nchar
cget the ascii value of a characte
      nascii = ichar(string(k:k))
      if ((nascii .gt. 96) .and. (nascii .lt. 123)) then
      nascii = nascii - 32
      string(k:k) = char(nascii)
      end if
c
      end do
      return 
      end
c
c--- converts a slant range vector into a aspect angle
c
      subroutine vecang(sr, azim, elev)
      dimension sr(3)
c
      data dtr / .01745329 /
      azim = atan2(- sr(1),- sr(2))
      azim = azim / dtr
      elev = atan2(- sr(3),sqrt((sr(1) ** 2) + (sr(2) ** 2)))
c
      elev = elev / dtr
      return 
      end
c
c---- converts aspect angles into a slant range vector
c
      subroutine angvec(azim, elev, sr)
      dimension sr(3)
c
      data dtr / .01745329 /
      rel = elev * dtr
      raz = azim * dtr
      sr(1) = - (cos(rel) * sin(raz))
      sr(2) = - (cos(rel) * cos(raz))
c
      sr(3) = - sin(rel)
      return 
      end
      subroutine vecdif(u, v, w)
      dimension u(3), v(3), w(3)
      w(1) = u(1) - v(1)
      w(2) = u(2) - v(2)
      w(3) = u(3) - v(3)
      return 
      end
      subroutine veclin(a, u, b, v, w)
      dimension u(3), v(3), w(3)
      w(1) = (a * u(1)) + (b * v(1))
      w(2) = (a * u(2)) + (b * v(2))
      w(3) = (a * u(3)) + (b * v(3))
      return 
      end
      FUNCTION VECMAG(A)
C... function to compute the magnitude of a 3-dimensional vector
      DIMENSION A(3)
      VECMAG = SQRT( A(1)*A(1) + A(2)*A(2) + A(3)*A(3) )
      RETURN
      END
      subroutine vecscal(c, x, y)
      dimension x(3), y(3)
      y(1) = c * x(1)
      y(2) = c * x(2)
      y(3) = c * x(3)
      return 
      end
      subroutine vecsum(u, v, w)
      dimension u(3), v(3), w(3)
      w(1) = u(1) + v(1)
      w(2) = u(2) + v(2)
      w(3) = u(3) + v(3)
      return 
      end
