*********************************************************************** 
*     Program diffeq.f
*     Author: Jesse Welch, for Jim Mahoney's Programming Workshop SP06
*     Date: 4/05/06
*     Version: 2.0
*
*     This program calculates the values of a function y(x) over the 
*     interval [x0,deltx*nstep], each of which the user is given full
*     choice.  The function which is being solved is y'(x) = f(x, y(x))
*     With f(x, y(x)) = const1*((x1)**(const2))*((x2)**(const3))
*     The user also inputs the initial value of y(x) at x = x0.
*     In this way, the user can use this program to calculate any 
*     ordinary differential equation with initial values and then
*     using the outputted file, plot the values of y(x) over a range
*     supplied by the user.
*
*     Version History:
*     1.0 - Before the program was at all utilizable i went through
*               the stages of creating first the subroutine which would
*               carry out the calculation, followed by the deriv funct-
*               ion.  Once this was done, i followed up with a front
*               end so to speak with the user controls + file access.
*     2.0   Once i had a working system, i tightened up all the
*               loose ends, added in the ability for user controlled 
*               derivative function.  
*     3.0   Cleaned up the notation and variable names throughout.  
*               I also added in the beginnings of an adaptive step size
*     3.5       Though its not a perfect step size algorithm, there is
*               something to it now, it does work.
*
*     Bugs: None known of at this point, many testing cycles checked
*               Though I consider my step size algorithm to be bugged
*               There is currently no way to go back from a decreased
*               step size, and the only way it does it is cutting it in
*               half.  These could be made better, though the program 
*               does run without error.
*
*     Future: I am going to follow up with a batch file which runs and 
*               plots the outputted file. I also hope to extend my 
*               program to include an extension to systems of ordinary
*               differential equations linked to a partial differential
*               equation.
*     
*     Usage:  The program runs with a very explicit nature, just answer
*               the questions as they are asked.
*
***********************************************************************

***********************************************************************
*     MAIN PROGRAM
***********************************************************************

      program diffeq
      
      real*8 deltx, xtemp1, xtemp2, xtemp3, ytemp1, ytemp2, ytemp3
      real*8 a, b, c, error, x0
      integer j, nstep


***********************************************************************
*     Debug help, easy solution to dy/dx = y(x) which is known to me,
*     just used to speed things up while testing
***********************************************************************
*     nstep = 100
*     deltx = .1
*     error = .1    
*     x0 = 0
*     ytemp1 = 1
*     a = 1
*     b = 0
*     c = 1
***********************************************************************

      write(*,*) 'How many steps do you want to calculate? '
      read(*,*) nstep
      write(*,*) 'What increment do you want to use? '
      read(*,*) deltx
      write(*,*) ' What is the percentage error limit you are '
      write(*,*) ' willing to accept? '
      read(*,*) error
      write(*,*) 'What is the initial value of x, x0? '
      read(*,*) x0
      write(*,*) 'What is the initial value of y(x) at x = x0? '
      read(*,*) ytemp1
      write(*,*) 'What is the form of the function of x and y(x) '
      write(*,*) 'f(x,(y(x)) which equals the derivative of y(x)? '
      write(*,*) 'In the form dy(x)/dx = f(x,y(x)), with '
      write(*,*) ''
      write(*,*) 'f(x,y(x)) = a * x^b * y(x)^c '
      read(*,*) a,b,c
      open(1, file='runge1.gnu')
      write(1,*) x0,ytemp1
      open(2, file='runge2.gnu')
      write(2,*) x0,ytemp1
      open(3, file='runge3.gnu')
      write(3,*) x0,ytemp1  

      xtemp2 = x0 + deltx/2.0
      xtemp3 = x0 + deltx

***********************************************************************
*     Loop which creates file data points, calling on the sub runge
***********************************************************************

      ytemp2 = ytemp1

      do 60 j = 1, nstep
         xtemp1 = x0 + j*deltx

***********************************************************************
*     This is the first time runge is called, calculates the next value
*     of the solution function for the given specifications
***********************************************************************
         call runge(xtemp1, ytemp1, deltx, a, b, c)
         write(1,*) xtemp1 , ytemp1           


***********************************************************************
*     This is the second time runge is called, calculating the next 
*     value of the function of half a step, then proceeds to calculate
*     the next step to complete the full step of the original
***********************************************************************
         call runge(xtemp2, ytemp2, deltx/2.0, a, b, c)
         write(2,*) xtemp2 , ytemp2           
         ytemp3 = ytemp2
         call runge(xtemp3, ytemp2, deltx/2.0, a,b,c)
         write(2,*) xtemp3, ytemp2

***********************************************************************
*     This is the adaptive step size part, if your calculation is with
*     in your error percentage, the calculation carries on as normal,
*     however, if the difference is greater than your error, the 
*     computer will cut your step size in half and carry on with twice
*     as many total steps.
***********************************************************************
         if(Abs(((ytemp2 - ytemp1)/ytemp1)) .lt. error) then
            write(3,*) xtemp1, ytemp1
         else
            write(3,*) xtemp2, ytemp3
            write(3,*) xtemp3, ytemp2
            nstep = 2 * nstep
         endif
         
         xtemp2 = xtemp2 + deltx
         xtemp3 = xtemp3 + deltx
 60   continue

      close(1)
      stop
      end

***********************************************************************
*     This is the subroutine which calculates the value of y(x) for a 
*     given indepentant and dependand variable, as well as three 
*     constants which define the derivative function.
***********************************************************************

      subroutine runge(indvar, depvar, deltax, constant, xpower, ypower)
      
      real*8 deriv, h, indvar, deltax, depvar,k1, k2, k3, k4, constant
      real*8 xpower, ypower


      h = deltax / 2.0

      k1 = deltax * deriv(indvar,depvar,constant,xpower,ypower)
    
      k2 = deltax *
     *     deriv(indvar + h, depvar+k1 / 2.0, constant, xpower, ypower)
      
      k3 = deltax *
     *     deriv(indvar + h, depvar+k2 / 2.0, constant, xpower, ypower)
      
      k4 = deltax *
     *     deriv(indvar + deltax, depvar + k3, constant, xpower, ypower)
      
      depvar = depvar + (k1 + (2.0 * (k2 + k3)) + k4) / 6.0
      
      return
      end


***********************************************************************
*     This is the function which calculates the value of the user
*     defined derivative function:  x2'(x1) = f(x1,x2(x1)) with 
*     f(x1, x2(x1)) = a * x**b * y**c
***********************************************************************

      function deriv(indvar,depvar,constant,xpower,ypower)
      
      real*8 deriv, indvar, depvar, constant, xpower, ypower


      deriv = constant*((indvar)**(xpower))*((depvar)**(ypower))

      return
      end
