program cooling_r_n_expt
c This program will find cooling constant r by matching with experimental data
      implicit none
      integer pgopen, i, j, j_max, n_steps, i_color
      parameter (j_max=24)
      real r, time, time_max, T, time_expt(j_max),T_expt(j_max), delta_t
      real T_s, T_ini, dummy

c Ask various input data
      write(*,*) 'Please type in how long (in min) you will wait :'
      read (*,*) time_max
      write(*,*) 'Please type in room temperature :'
      read (*,*) T_s
      write(*,*) 'Please type in the initial termperature :'
      read (*,*) T_ini

c Open PGPLOT devices and set environments
      if(pgopen('/xwin').le.0) stop
      call pgpap (5.0,0.75)
      call pgenv (0.0,time_max,0.0,100.0,0,0)
      call pglab ('Time (min)', 'Temperature (C)', 'Coffee Cooling')

c Then open/read external file draw the experimental value (black coffee only)
      open (unit=10, file='coffee_expt.dat')
      do j=1,j_max
        read (10,*) time_expt(j), dummy, T_expt(j)
      end do
      close (10)

c Draw the black coffee values
      i_color = 1
      call pgsci (i_color)
      call pgpt (j_max, time_expt,T_expt,-1)

c Ask r and Delta_t
      write(*,*) 'Please type in cooling constant r :'
      read (*,*) r
      write(*,*) 'Please type in small time-step delta_t :'
      read (*,*) delta_t

c Run cooling simulation and draw thw results
 100  n_steps = int(time_max/delta_t)
      time = 0.0
      T = T_ini
      i_color = i_color + 1
      call pgsci (i_color)
      call pgpt (1, time, T, -1)
      do i=1,n_steps
        time = time + delta_t
        T = T - r*(T-T_s)*delta_t
        call pgpt (1, time, T, -1)
      end do

c Do a new r draw again, stop if r is smaller or equal zero.

      write(*,*) 'Now you can give next r, or let r smaller than zero to
     & stop'
      read (*,*) r
      if (r.gt.0.0) goto 100

c Stopping
      call pgclos

      end