anchor.diff:          call pgpt(23, time_stamps, temps_black, 2)
           return
         end
  -
  -      function getR(t_stamps, temps_coffee)
  -      real temps_coffee(23)
  -      getR = 1
  -      end
  
positive.code:       program cooling
        implicit none
        real r,time,t_max,T,Ts,T_ini,delta_t, getR
        integer pgopen,i,n_steps,i_again, addedMilk, zipAble
  
   100  write(*,*) "Cooling rate and surroundind temperature = 17:"
        read(*,*) r  !, Ts
        Ts = 17
        write(*,*)     'Time duration and initial temperature:'
        !read(*,*) t_max, T_ini
        t_max = 20
        T_ini = 90.0
        
        if(pgopen('/xwin').le.0) stop
        call pgpap (5.0,0.75)
  
   102  write (*,*) 'Please type in the small time interval delta_t :'
        delta_t = 0.001 !read (*,*) delta_t
        time = 0
        T = T_ini
        call pgeras
        call pgenv(0.0,t_max,70.0,T_ini+5,0,0)
        call pglab('Time (min)','Temperature (C)','Cooling Time')
        !call readData()  not needed now
        !adding milk first
        zipAble = 0
        n_steps = int(t_max/delta_t)
        call pgsci(2)
  
negative.code:           call pgmove(time,T)
            addedMilk = 1
            T = T-5
            r = r_irish
            call pgsci(1)
            call pgdraw(time, T)
            call pgsci(6)
          end if
          if (zipAble.eq.0.and.T.le.75.0) then
            call pgsci(7)
            !write(*,*) "Milk at end makes coffee sipable after ", time
            time2 = time
            zipAble = 1
          end if
          !write(*,*) "Plot plot"
          call pgpt (1,time,T,-1)
        end do
        call pgsci(1)
        time_diff = time1 - time2
        if(abs(time_diff).le.0.1) then
          write(*,*) "Adding whiskey first better when r_irish=", r_irish
        else
          if(time1.le.time2) then
            r_irish = r_irish - 0.0001
          else
            r_irish = r_irish + 0.0001
          end if
          call sleep_time(2)
          goto 103
        end if
        write (*,*) 'Again ? (1=YES;0=NO)'
        read (*,*) i_again
        if (i_again.eq.1) goto 100
        call pgclos
  
        end program
  
        subroutine Euler(t, T_coffee, T_room, r, delta_t)
          change = -r*(T_coffee - T_room)*delta_t
          T_coffee = T_coffee + change
          t = t + delta_t
        end
  
        subroutine sleep_time(st)
        integer st
        do i = 1, (10000000 * st)
         idummy = i**3
        enddo
        end
  
        subroutine readData()
  