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:       time = 0
        zipAble = 0
        n_steps = int(t_max/delta_t)
        T = T_ini
        call pgsci(6)
        do i=1,n_steps
  
negative.code:       program cooling
        implicit none
        real r,time,t_max,T,Ts,T_ini,delta_t, getR
        real r_irish, r_black, time1, time2
        real time_diff
        integer pgopen,i,n_steps,i_again, addedMilk, zipAble
        integer black_greater
        if(pgopen('/xwin').le.0) stop
        call pgpap (8.0,0.8)
        call pgask(.false.)
   100  write(*,*) "Enter cooling rate coffee, temperature = 17:"
        read(*,*) r_black  !, Ts
        write(*,*) "Enter irish cooling rate guess for it to be" 
        write(*,*) "more effective to add it first"
        read(*,*) r_irish
        Ts = 17
        write(*,*)     'Time duration and initial temperature:'
        !read(*,*) t_max, T_ini
        t_max = 20
        T_ini = 90.0
   102  write (*,*) 'Please type in the small time interval delta_t :'
        delta_t = 0.001 !read (*,*) delta_t
   103  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')
  
        black_greater = 1
        if(r_black.le.r_irish) black_greater = -1
        !call readData()  not needed now
        !adding milk first
        r=r_black
        zipAble = 0
        n_steps = int(t_max/delta_t)
        call pgsci(2)
  