#! /usr/bin/octave -qf

1;

#old_r = 100;
#old_t = 0;

function ris = r(t,M,D)
  ris = r1(t);
  #ris = r1bis(t);
  %ris = r1ter(t);
  #ris = r2(t,D);
  #ris = r3(t,M,D);
endfunction

function ris = r1(t)
  if (t<25)
	ris = 100;
  else
	ris = 0;
  endif
endfunction

function ris = r1bis(t)
  diff = t-floor(t);
  if ((t<5)&& (
      ((0<diff)&&(diff<0.05))
        || (0.2<diff)&&(diff<0.25)
        || (0.4<diff)&&(diff<0.45)
        || (0.6<diff)&&(diff<0.65)
        || (0.8<diff)&&(diff<0.85)
      ))
	ris = 1000000;
  else
	ris = 0;
  endif
endfunction

function ris = r1ter(t,M,D)
  global old_t;
  global old_r;
  if (t-old_t>0.0001) 
    old_t = old_t + 0.0001;
    if (t<5)
	    ris = 100;
    else
	    ris = 0;
    endif
    old_r = ris;
  else
    ris = old_r;
  endif
endfunction


function ris = r2(t,D)
  global old_t;
  global old_r;
  if (t-old_t>0.1) 
    old_t = old_t + 0.1;
    if (D<0.5)
	    ris = 100;
    else
	    ris = 0;
    endif
    old_r = ris;
  else
    ris = old_r;
  endif
endfunction

function ris = r3(t,M,D)
  global old_t;
  global old_r;
  if (t-old_t>0.0001) 
    old_t = old_t + 0.0001;
    if (D<0.5)
	    ris = 100;
    else
	    ris = 0;
    endif
    old_r = ris;
  else
    ris = old_r;
  endif
endfunction


function ris = stimolo(r,M,D,soglia)
  if (r-M >= soglia)
    ris = 1;
  else 
	if (r-M >= -30 && r-M<soglia)
		ris = 0;
	else
		ris = -M*D/2;
	endif
  endif
endfunction

function ris = ricordo(r,M)
  if (r>M)
    ris = (r-M)/2;
  else
    ris = 0;
  endif;
endfunction

function xdot = f(x,t)

soglia = 80;
k = 0.2;

D = x(1);
M = x(2);

xdot(1) = (0.2)*(-D + stimolo(r(t,M,D),M,D,soglia) + k);
xdot(2) = (0.2)*(-M +ricordo(r(t,M,D),M));

endfunction


function model2()

# INITIAL CONDITIONS
x0 = [0.2; 0];


simtime = 50;


t = linspace(0,simtime,200);


y = lsode("f",x0,t);


plot (t,y(:,1),"-");

pause();


plot (t,y(:,2),":");

%pause();


endfunction

model2();