program styrofoam_k
      implicit none
      integer i,i_count,pgopen,ik,k_steps,i_color,i_current,ik_best
      real time_expt(20),position_expt(20),mass,g,t_min,t_max,y_min
      parameter (mass = 2.54, g = -9.8)  
      real time,y,v,dt,k_ini,k_final,dk,F_k,F_k_min,y_max
      real a,a_mid,v_mid,y_new,time_new,k,y_interp,k_best
      

      write(*,*) 'Please input initial k and final k, and steps:'
      read(*,*) k_ini, k_final, k_steps

      open (unit=10,file='styrofoam.txt',status='old',err=100)
      i_count = 0
      write(*,*) 'The experimental time and positions:'
      do i=1,20
        read (10,*,end=200) time_expt(i), position_expt(i)
        position_expt(i) = -1*position_expt(i)
        write(*,*) time_expt(i), position_expt(i)
        i_count = i_count + 1 
      end do
      write(*,*)
 200  close (10)

C Plot experimental curve
      y_min = position_expt(i_count)
      y_max = position_expt(1)

      t_min = time_expt(1)
      t_max = time_expt(i_count)

      if (pgopen('/xwin').le.0) stop
      call pgpap (7.0,0.75)
      call pgenv (t_min,t_max,y_min,y_max,0,0)
      call pglab ('Time','Hight','Determine resistance of air, k')
      call pgpt (i_count,time_expt,position_expt,9)


      dt = 0.0001
c outer-most k-loop
      dk = (k_final-k_ini)/k_steps
      do 1000 ik=1,k_steps
      i_color = mod(ik,15)
      call pgsci(i_color)
      k = k_ini + ik*dk


c Start simulated curve
      time = t_min
      i_current = 1
      v = 0.0
      y = position_expt(1)
      F_k = 0.0
      do while (time .lt. time_expt(i_count))
        a = g + k*v*v 
c       a = g + k*(-1)*v
        v_mid = v + a*(0.5)*dt
        a_mid = g + k*v_mid*v_mid
c       a_mid = g + k*(-1)*v_mid
        v = v + a_mid*dt 
        y_new = y + v_mid*dt
        time_new = time + dt
        if (time_new.gt.time_expt(i_current+1)) then
          i_current = i_current + 1
          call interpolate(time,y,time_new,y_new,time_expt(i_current),
     &         y_interp)
          F_k = F_k + ( y_interp - position_expt(i_current) )**2
        endif
        time = time_new
        y = y_new
        call pgpt (1,time,y,-1)
      end do

      write(*,*) 'F(k) = F(',k,')=',F_k
      if (ik.eq.1) then
         F_k_min = F_k
         k_best = k
         ik_best = ik
      end if
      if (F_k.lt.F_k_min) then
         F_k_min = F_k
         k_best = k
         ik_best = ik
      end if
 1000 continue

      if ( (ik_best.eq.1) .or. (ik_best.eq.k_steps) ) then
        write (*,*)
        write (*,*) 'Warning, no F(k) minimum is found in this k range!'
        write (*,*)
       else
        write (*,*)
        write (*,*) 'This best k is', k_best
        write (*,*)
      end if

      call pgclos
 100  stop
      end

      subroutine interpolate(x1,y1,x2,y2,x,y)
      real x1,y1,x2,y2,x,y, slope
      slope = (y2-y1)/(x2-x1)
      y = slope*(x-x1) + y1
      end