C
C program 16
C
C more illustrations of subroutines and functions
C
C  CWJ SDSU 2/8/2005
C
C  this program computes the first derivative of the gaussian
C  and checks how it changes with dx
C
C  MAIN PROGRAM CALLS:
C 
C  get_gauss_params
C  get_deriv_params
C  countdown
C

      program program16

      implicit none
    
C................parameters of gauss exp(-ax^2)
      real a,x
C................parameters for approximate derivative
      real dxmin,dxmax
      integer nsteps

C...................................................................

      call get_gauss_params(a,x)
      call get_deriv_params(dxmin,dxmax,nsteps)
      call countdown(a,x,dxmin,dxmax,nsteps)

      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine get_gauss_params(a,x)

      implicit none
      real a,x

      print*,' Please enter parameters for gaussian function ',
     & 'exp( - a x^2 ) '
  303 continue
      print*,' please enter a,x '
      read*,a,x
    
C   error trap
       if(a.le.0.0)then
          write(6,*)' a must be > 0 '
          goto 303
       endif

      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine get_deriv_params(dxmin,dxmax,nsteps)

      implicit none

      real dxmin,dxmax
      integer nsteps

      real temp   ! for swapping if needed

  100 continue
      print*,' enter # of steps ' 
      read(*,*)nsteps 
      if(nsteps.eq.0)goto 100

  101 continue
      print*,' enter largest dx '
      read*,dxmax
      print*,' enter smallest dx '
      read*,dxmin

C...........  make sure they are not zero 

      if(dxmax.eq.0.0 .or. dxmin.eq.0.0)then
        print*,' naughty, no naughts! '
        goto 101
      endif
      if(dxmax.eq.dxmin)then
        write(*,*)' must be different '
        goto 101
      endif

C......... make sure they are positive

      dxmax=abs(dxmax)
      dxmin = abs(dxmin)

C.............. make sure dxmax > dxmin

      if(dxmin.gt.dxmax)then
         temp = dxmax
         dxmax = dxmin
         dxmin=temp
      endif

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     subroutine countdown
C
C   changes dx from dxmax down to dxmin
C    in logarithmic steps
C
C   input:
C     a,x   ! parameters of gaussian
C     dxmax,dxmin  !limits of dx
C     n   ! # of steps to take      
C
C  output:
C   NONE (writes to terminal)
C
C  calls: 
C    subroutine dgaussdx
C    function exactderivgauss
C

      subroutine countdown(a,x,dxmin,dxmax,n)

      implicit none

      real a,x  ! gaussian parameters

      real dxmax,dxmin
      integer n      ! parameters for derivatives

      real exactderivgauss    ! exact value of derivative

C................... now to discretize steps in dx
C                    do it logarithmically

      real lndx
      real dx
      real dlndx
      integer i

C.......................
      logical flag

C................
      real fprime
      real error
C------------------- 

       dlndx = (log(dxmax)-log(dxmin))/float(n)

      print*,'   dx      approx    exact      diff  '      
      do i = 0,n
        lndx = log(dxmax)-i*dlndx
        dx = exp(lndx)
        call dgaussdx(a,x,dx,fprime,flag)
         if(flag)then
          print*,' There is a problem with a ',a
          stop
        endif    
        error=(exactderivgauss(a,x)-fprime)
        write(6,101)dx,fprime,exactderivgauss(a,x),error
  101 format(f6.4,3x,3f10.5)           ! we will cover this later
         
      enddo

      return    !actually this is not really needed as we are done
      end

CCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  subroutine dgaussdx
C
C  couputes d/dx exp(-ax^2)
C
C  INPUT:
C      a,x
C      dx : step size for dgauss
C
C  OUTPUT:
C      dgdx = derivate d/dx exp(-ax^2)
C              using symmetric "three-point" formula
C
C      errflag = .true. if a < 0, else = .false
C
C   FUNCTIONS CALLED:
C      gauss
C

      subroutine dgaussdx(a,x,dx,dgdx,errflag)

      implicit none
      real a   ! parameter for gaussian
      real x
      real dx  ! discretization step
      
      real dgdx

      real gauss  ! this user-defined function should be declared 

      logical errflag  

      errflag=.false.
      if(a.le.0.0)then
        errflag = .true.
        return
      endif

C   alternate two-point asymetric function
C      dgdx=(gauss(a,x+dx)-gauss(a,x))/dx

      dgdx = (gauss(a,x+dx) - gauss(a,x-dx))/(2*dx)

      return
      end


CCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C function gauss
C
C gaussian function  - exp(-a*x^2)
C
C  input: a,x
C

      function gauss(a,x)
      implicit none
      real a,x
      real gauss

      gauss=exp(-a*x*x)
      return
      end
CCCCCCCCCCCCCCCCCCCCCCC
C
C  function exactderivgauss
C
C  INPUT: a,x
C
C  FUNCTIONS CALLED:
C   gauss(a,x)
C 
      function exactderivgauss(a,x)

      implicit none

      real exactderivgauss ! note since it was not declared above, 
                                         ! must declare here (since implicit none)

      real a,x
      real gauss	! because this function is called

      exactderivgauss = -2*a*x*gauss(a,x)

      return
      end