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,nstep
      integer method 
      dimension z(nx),v(nx),f(nx)
c
c
c2--567--1--------2--------3---------4---------5---------6---------772
c
c for file-open

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)
        enddo

      else
        write(*,*)' Runge-Kutta (2nd or 4th) not installed yet ! '
        stop
      endif
c
c
c for output of solutions
c
      write(16,800)(i,dt*dble(i-1),z(i),v(i),f(i),i=1,nstep)
 800  format(i5,4e22.12)


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