C mpfortran -pesize=8 runge.f
C
C----ROB SALGADO (cps615g)
C    runge.f 10/04/92 (SET 3 [??] PROBLEM 1)
C


C--------------------------------------------------PROGRAM
      PROGRAM runge
      INTEGER i,j,q,hz,hzz
      REAL h,x
      DIMENSION y(4),dydx(4),yout(4)
      COMMON pi,g
      pi=3.1415926
      g=39.478301356


      x=0.0
C------INITIAL DATA
      y(1)=-0.2594970
      y(2)=-0.9809988
      y(3)=5.9726963
      y(4)=-1.6298815
      print* ,y(1),y(2),y(3),y(4)
      CALL DERIVS(x,y,dydx)

      print*,'Step size h through [0,1], h=10^(-q). Input q ='
      read*,q
      
C------STEP THROUGH
      hz=int(10**q)
      hzz=int(hz/10)
      hy=1.0/hz
      print*,q,hz,hzz,hy
      do 200 i=0,hz-1,hzz
C               ----------OKAY, UPDATE THE USER
         print*, 'h=',h
         print*, 'rk4:', yout(1),yout(2),yout(3),yout(4)
C         call ANALYTIC(h,yout)
         do 201 j=0,hzz-1
C               ----------AVOID INTRODUCING PROPAGATING ERRORS
         h=(i+j)*hy
C         print *,i+j
C               ----------CALL RUNGE-KUTTA
         CALL RK4(y,dydx,x,h,yout)
 201     continue
 200     continue
         END

C-------------------------------------DEFINE THE SYSTEM OF ODEs
      SUBROUTINE DERIVS(x,y,dydx)
      DIMENSION y(4),dydx(4)
      REAL x,r
      COMMON pi,g

      r=sqrt(y(1)*y(1)+y(2)*y(2))
C                   -------dydx(1) is (y1)', etc...
      dydx(1)=y(3)
      dydx(2)=y(4)
      dydx(3)=-g*y(1)/(r*r*r)
      dydx(4)=-g*y(2)/(r*r*r)
      RETURN
      END
C--------------------------------------ANALYTIC SOLUTION (WITH ICs!)
      SUBROUTINE ANALYTIC(h,y)
      real h,pi2
      common pi,g
      dimension a(4),y(4)
      pi2=2.0*pi
      a(1)=cos(pi2*h)
      a(2)=sin(pi2*h)
      a(3)=-pi2*sin(pi2*h)
      a(4)=pi2*cos(pi2*h)
      print*,'Ana:',a(1),a(2),a(3),a(4)
      if (a(1)*a(2)*a(3)*a(4).eq.0) goto 400
      print*,'              Errors:',(y(1)-a(1))/a(1),(y(2)-a(2))/a(2)
      print*,'                     ',(y(3)-a(3))/a(3),(y(4)-a(4))/a(4)
 400  RETURN
      END
C-------------------------------------RUNGE-KUTTA
C     y=function(starting point)
C     yt=function(starting point + intermediate steps)
C     yout=function(current point)
C     dym, dyt = (d/dx)(function(intermediate steps))

      SUBROUTINE RK4(y,dydx,x,h,yout)
      DIMENSION y(4),dydx(4),yout(4)
      REAL x,h
      INTEGER i,n
      REAL xh,hh,h6
      DIMENSION dym(4),dyt(4),yt(4)
      n=4
C                    ----DEFINE SOME CONSTANTS, FOR CONVENIENCE
      hh=h*0.5
      h6=h/6.0
      xh=x+hh
C                    ----hh*dydx = (1/2)*h*f(x,y) = (1/2)*k1
      DO 101 i=1,n
         yt(i)=y(i)+hh*dydx(i)
 101     continue
C                    ----evaluate dyt=f(x+(1/2)*h,y+(1/2)*k1)
      CALL DERIVS(xh,yt,dyt)
C                    ----hh*dyt = (1/2)*h*dyt = (1/2)*k2
      DO 102 i=1,n
         yt(i)=y(i)+hh*dyt(i)
 102     continue
C                    ----evaluate dym=f(x+(1/2)*h,y+(1/2)*k2)
      call DERIVS(xh,yt,dym)
C                    ----h*dym  = k3
      DO 103 i=1,n
          yt(i)=y(i)+h*dym(i)
C                    ----h*dym = (k3 + k2)
          dym(i)=dym(i)+dyt(i)
 103      continue
C                    ----evaluate dyt=f(x+h,y+k3)
      call DERIVS(x+h,yt,dyt)
C                    ----yout=y+(h/6)*(k1 + k4 + 2*(k3+k2))
      DO 105 i=1,n
         yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.0*dym(i))
 105     continue
      RETURN
      END








