!
! $Id: llst.f,v 1.3 2000/02/14 23:14:12 aai Exp aai $.
!
!{
      program llst
!
!     Llst - Test lls() et. al.
!
!     To use this program, compile and link with libfortex90.a then
!     execute and either type in (X,Y) pairs directly or redirect them
!     as standard input, one (X,Y) pair per line.
!
!     This version differs from the Fortran 77 version by using dynamic
!     memory allocation.
!
!     Phillip A. Cheeseman, Purdue University Computing Center.
!
!     MAXCAP is the initial maximum number of points we can keep in the
!     data arrays.  The working capacity is stored in nmax.  MAXCAP is
!     deliberately set to a small number for testing purposes in case you
!     try this code somewhere where I have not.  If the reallocation code
!     doesn't work, then it should show up on the first try rather than
!     the 1000th.
!
      integer MAXCAP
      parameter (MAXCAP=2)
!
!     Data arrays and pointers for (X,Y) pairs, and allocation loop indicator.
!
      real, allocatable, dimension(:), target ::  x1,x2
      real, allocatable, dimension(:), target ::  y1,y2
      real, pointer, dimension(:) :: x
      real, pointer, dimension(:) :: y
      logical :: odd = .false.
!
!     Least square line parameters and arrays for plotted line.
!
      real :: a
      real :: b
      real :: sigma
      real :: xc(2)
      real :: yc(2)
!
!     Response code returned by lls() and counter for data points.
!
      integer :: iresp
      integer :: n
      integer :: nmax = MAXCAP
!
!     Messages we'll return if we get a negative response from lls().
!     Data statements are used to initialize the strings to avoid clutter
!     in the declarations (I'm not a fan of statement continuations).
!
      character*40, dimension(3) :: problems
!
      data problems(1) / 'No data points.' /
      data problems(2) / 'Calculation of std. dev. not possible.' /
      data problems(3) / 'Least squares line is vertical.' /
!
!     Make initial allocations and set up pointers.
!
      allocate(x1(nmax))
      allocate(y1(nmax))
      x => x1
      y => y1
!
!     Preset count of data points then read (X,Y) pairs until end of file
!     or limit of arrays.
!
      n = 0
00011 read(*,*,end=21) x(n+1),y(n+1)
      n = n+1
      if(n.lt.nmax) go to 11
!
!     We've run out of space.  Allocate a larger block and move data there.
!     Note that we blindly double what's already in use w/o regard to what
!     we've already used.  If we hit a system limit, some type of execution
!     error should result.
!
      write(*,'(a)') '  Llst: Doubling data capacity.'
      odd = .not.odd
!
!     If this is an odd pass through reallocation, allocate x2 and y2.
!
      if(odd) then
         allocate(x2(nmax+nmax))
         allocate(y2(nmax+nmax))
         x2(1:nmax) = x1(1:nmax)
         y2(1:nmax) = y1(1:nmax)
         deallocate(x1)
         deallocate(y1)
         x => x2
         y => y2
      else
         allocate(x1(nmax+nmax))
         allocate(y1(nmax+nmax))
         x1(1:nmax) = x2(1:nmax)
         y1(1:nmax) = y2(1:nmax)
         deallocate(x2)
         deallocate(y2)
         x => x1
         y => y1
      endif
      nmax = nmax+nmax
      go to 11
!
!     We've reached end of file, now calculate the least squares line.
!
00021 call lls(n, x, y, iresp, sigma, a, b)
      xc(1:2) = x(1)
      do  i=2,n
	 xc(1) = max(x(i),xc(1))
	 xc(2) = min(x(i),xc(2))
      enddo
!
!     Here, we issue an error message if iresp is negative.
!     This could be more elegantly done by dropping the trailing blanks
!     but that's more code for messages we shouldn't see too often.
!
      if(iresp.lt.0) then
         write(*,'(a,i2,a)') '  Llst: problem code is (',iresp,').'
         if((iresp.le.-1).and.(iresp.ge.-3)) then
            write(*,'(8x,a)') problems(-iresp)
         endif
      endif
!
!     Write the results.
!
      if((iresp.ge.0).or.(iresp.eq.-2)) then
         write(*,*) ' Llst: A=',a
         write(*,*) '       B=',b
         if(iresp.ne.-2) then
            write(*,*) '       S=',sigma
         endif
      endif
!
!     Calculate the endpoints of the line and plot it with the raw data.
!     The section below requires the PEPL library (-lpepl in the Makefile).
!     If you don't have, or don't care for PEPL, simply strike the lines
!     up to the 'call chtqit'.
!
!
      yc(1) = a+xc(1)*b
      yc(2) = a+xc(2)*b
      call chtini(4)
      call chtfmt(0, 'Least Squares Line', 'Y', 'X')
      call chtcrv(n,  x,  y, 0, 0)
      call chtcrv(2, xc, yc, 0, 0)
      call chtfin
      call chtqit
!
      stop
!
      end
!}