program slopes

c     Purpose: Print out the slopes of a 2D autonous system of the form
c                 dx/dt = fx(x,y)
c                 dy/dt = fy(x,y)

c     Usage: Edit files fx.f and fy.f and put in the correct functions fx
c            and fy to be plotted.
c            Edit file slopes.ind to set the ranges of x and y to plot,
c            and the number of slopes to plot in each direction.
c            Then execute 'make'.

c     avoid typos
      implicit none

c declarations

c     the d(x)/d(t) and d(y)/d(t) velocity functions
      real fx,fy
c     maximum values
      real fxmax,fymax

c     plot domain size
      real xlo,xhi,ylo,yhi
c     position
      real x,y
c     mesh size
      real dx,dy

c     number of slopes in the x and y-directions
      integer imax,jmax
c     slope indices
      integer i,j

c     whether we scale the slopes according to the local velocity
      integer scaled

c     scale factor
      real scale

c     relative size of the slopes
      real relsiz

c     slope components
      real sx,sy

c     cosine and sine of the arrow angle
      real ca,sa

c     input file line
      integer line

c executable

c     set cosine and sine of the arrow angle
      ca=cos(20./180.*3.141593)
      sa=sin(20./180.*3.141593)

c     read input file
      line=0
      open(1,file='slopes.ind',status='old',err=210)
      line=1
      read(1,*,err=220)xlo
      line=2
      read(1,*,err=220)xhi
      line=3
      read(1,*,err=220)ylo
      line=4
      read(1,*,err=220)yhi
      line=5
      read(1,*,err=230)imax
      line=6
      read(1,*,err=230)jmax
      line=7
      read(1,*,err=230)scaled
      line=8
      read(1,*,err=230)relsiz
      close(1,err=220)
      goto 250
 210  call fatal('Unable to open input file slopes.ind!')
 220  call fatal('Error reading the floating point number in line '//
     & char(line+48)//' of file slopes.ind!')
 230  call fatal('Error reading the integer number in line '//
     & char(line+48)//' of file slopes.ind!')
 240  call fatal('Unable to close file slopes.ind')
 250  continue

c     check data
      if(xlo.ge.xhi)call fatal(
     & 'Lowest x-value is not smaller than the highest!')
      if(ylo.ge.yhi)call fatal(
     & 'Lowest y-value is not smaller than the highest!')
      if(imax.le.0)call fatal(
     & 'Invalid number of slopes in the x-direction!')
      if(jmax.le.0)call fatal(
     & 'Invalid number of slopes in the y-direction!')
      if(scaled.ne.0 .and. scaled .ne.1)call fatal(
     & 'Invalid slope scaling value, must be 0 or 1!')
      if(max(xhi,yhi).gt.999.99 .or.min(xlo,ylo).lt.-99.999)call fatal(
     & 'Plotting ranges are too big, 100 is the limit!')

c     mesh size
      dx=(xhi-xlo)/imax
      dy=(yhi-ylo)/jmax

c     find the maximum values of fx and fy
      fxmax=0.
      fymax=0.
      do 500 i=1,imax
         x=xlo+(i-0.5)*dx
         do 490 j=1,jmax
            y=ylo+(j-0.5)*dy
            fxmax=max(fxmax,abs(fx(x,y)))
            fymax=max(fymax,abs(fy(x,y)))
 490     continue
 500  continue
      if(fxmax+fymax.eq.0.)call fatal(
     & 'Invalid fx and fy, all slopes are zero!')
      if(fxmax.gt.0.)then
         scale=min(dx/fxmax,dy/fymax)
      else
         scale=dy/fymax
      endif
      scale=scale*relsiz

c     plot the slopes
      open(2,file='slopes.plt',err=920)
      write(2,910,err=920)
 910  format('\\c 2')
      goto 930
 920  call fatal(
     & 'Unable to write to file slopes.plt!  Check disk quota!')
 930  do 1000 i=1,imax
         x=xlo+(i-0.5)*dx
         do 990 j=1,jmax
            y=ylo+(j-0.5)*dy
            if(scaled.eq.0)then
               sx=fx(x,y)
               sy=fy(x,y)
               scale=max(abs(sx/dx),abs(sy/dy))
               if(scale.ne.0.)then
                  sx=sx/scale*relsiz
                  sy=sy/scale*relsiz
               endif
            else
               sx=fx(x,y)*scale
               sy=fy(x,y)*scale
            endif
            write(2,950,err=960)
     &       x-.5*sx,y-.5*sy,
     &       x+.5*sx,y+.5*sy,
     &       x+.5*sx-.25*(sx*ca-sy*sa),y+.5*sy-.25*(sx*sa+sy*ca)
 950        format('3'/f9.5,',',f9.5/f9.5,',',f9.5/f9.5,',',f9.5)
            goto 990
 960        call fatal(
     &       'Unable to write to file slopes.plt!  Check disk quota!')
 990     continue
 1000 continue

c     all done
      call exit(0)
      end


      subroutine fatal(msg)

c     Terminate at a fatal error.

c     avoid typos
      implicit none

c input

c     error message
      character*(*) msg

c executable

      print10,char(7),msg
 10   format(a1,'*** ',a)

      call exit(1)
      end