c2--567--1--------2--------3---------4---------5---------6---------772
c
      program fcodebungy
c
c     copyright Tatsuki Oda (2009.05.01)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
      parameter (nx=1000)
c
      integer i,j,nstep
      integer method

      double precision,external :: fa,fv

      double precision xk1,xk2,xk3,xk4,xl1,xl2,xl3,xl4,t1,t

      dimension z(nx),v(nx),f(nx), ek(nx), ev(nx), et(nx)
c
c2--567--1--------2--------3---------4---------5---------6---------772
c
c for file-open
      open(unit=18,file='energy.xy')

c    for input
      open(unit=11,file='f-data.in')

c    for output
      open(unit=16,file='bungy.xy')
c
c
c2--567--1--------2--------3---------4---------5---------6---------772
c
c
c for data input
c
c input-data with (11)
c
      read(11,*) nstep     ! number of steps
      read(11,*) dt        ! Delta t
      read(11,*) am        ! mass of particle
      read(11,*) ak0       ! force constant (k) of spring
      read(11,*) alength   ! length of spring
      read(11,*) height    ! height of stage
      read(11,*) gamma     ! friction of resistance
      read(11,*) method    ! algorithm of solving
c                             method =0 velret
c                                    =1 Runge-Kutta (2th)
c                                    =2 Runge-Kutta (4th)
c
      g=9.80d0            ! acceleratio of gravity
c
c
      write(06,*)' Finish to read the data, now !'
c
c2--567--1--------2--------3---------4---------5---------6---------772
c
c
c for initial conditions
c
      z(1)=height
      v(1)=0.0d0           ! free diving
      f(1)=-am*g
c
      hml=height-alength
c
c for solving the equation of motion
c
      if(method.eq.0)then

        z(2)=z(1)+v(1)*dt+f(1)*dt**2*0.5d0/am
        v(2)=v(1)+f(1)*dt/am

         do i=2,nstep
          if(z(i).gt.hml)then
             ak=0.0d0
          else
             ak=ak0
          endif
          f(i)=-am*g-ak*(z(i)-hml)
          eta=gamma*dt/(2.0d0*am)
          etap1=eta+1.0d0
          z(i+1)=2.0d0*z(i)/etap1-
     &           z(i-1)*(1.0d0-eta)/etap1+
     &           f(i)*dt**2/am/etap1
          v(i)=(z(i+1)-z(i-1))/(2.0d0*dt)

          ek(i)=am*v(i)**2*0.5d0
          ev(i)=am*g*z(i)+ak*(z(i)-hml)**2*0.5d0
          et(i)=ek(i)+ev(i)
        enddo


      elseif(method.eq.1) then
           do i=1, nstep
            t = dt*(i-1)
            t1 = t + dt * 0.5d0
            vi=v(i)
            zi=z(i)
            xk1=dt*fa(t,vi,zi,am,g,ak0,hml,gamma)
            xl1=dt*fv(t,vi,zi,am,g,ak0,hml,gamma)
            xk2=dt*fa(t1,vi+xk1*0.5d0,zi+xl1*0.5d0,am,g,ak0,hml,gamma)
            xl2=dt*fv(t1,vi+xk1*0.5d0,zi+xl1*0.5d0,am,g,ak0,hml,gamma)
            xk3=dt*fa(t,vi*xk2*0.5d0,zi+xl2*0.5d0,am,g,ak0,hml,gamma)
            xl3=dt*fv(t,vi*xk2*0.5d0,zi+xl2*0.5d0,am,g,ak0,hml,gamma)
            xk4=dt*fa(t,vi*xk3,zi+xl3,am,g,ak0,hml,gamma)
            xl4=dt*fv(t,vi*xk3,zi+xl3,am,g,ak0,hml,gamma)

            v(i+1)=vi+(xk1+xk2*2d0+xk3*2d0+xk4)/6d0
            z(i+1)=zi+(xl1+xl2*2d0+xl3*2d0+xl4)/6d0

            if(z(i).gt.hml) then
               ak = 0d0
            else
               ak = ak0
            endif

            ek(i)=am*v(i)**2*0.5d0
            ev(i)=am*g*z(i)+ak*(z(i)-hml)**2*0.5d0
            et(i)=ek(i)+ev(i)

         enddo
      endif
c
c
c for output of solutions



       write(16,800)(i,dt*dble(i-1),z(i),v(i),f(i),i=1,nstep)
 800   format(i5,4e22.12)

       write(18,900)(i,dt*dble(i-1),ek(i),ev(i),et(i),i=1,nstep)
 900   format(i5,4e22.12)


       if(abs(gamma).lt.1d-9) then
             etave = 0d0
             etvar = 0d0
             do i=1, nstep
                etave = etave + et(i)
             enddo
              etvar = etvar/nstep
             do i=1, nstep
                etvar = etvar + (et(i)-etave)**2
             enddo
             etvar = etvar/nstep
             write(06,950) etave, etvar
 950         format('energy average : ', e15.6, ' variance :', e15.6)
          endif


 999  continue

      write(06,*)' Finish the execution, now !'
c
      stop
      end

c     ------------- first function program -------------------

       double precision function fa(t,v,z,am,g,ak0,hml,gamma)
c
c     copyright Tatsuki Oda (2009.05.01)
c
       implicit double precision (a-h,o-z)
       implicit integer (i-n)
c
       if(z.gt.hml) then
          ak = 0d0
       else
          ak = ak0
       endif
       fa = -g-(ak*(z-hml)+gamma*v)/am
c
       return
       end

c---------------end of function program --------------------

c---------------second function program --------------------
      double precision function fv(t,v,z,am,g,ak0,hml,gamma)
c
c     copyright Tatsuki Oda (2009.05.01)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
      fv=v
c
      return
      end

c--------------end of second function program --------------