function deltah=simulate2017(ip,q,deltat);

// transient simulator
// ip and q are row vectors of the same size.
// q is infiltration in meters/day

K=1; W=1000; l=1000; Sy=0.05; T=60;
// basic hydrolgic data
// T is thickness of aquifer

[ig,N]=size(ip);
// input heads as an array 1 X N

for i=2:N-1 // note that ip(1) and ip(N) are treated as constants

// hsallow aquifer assumption
//ll=(ip(i)-ip(i-1))*W*min(ip(i),ip(i-1))*K/l;
//rr=(ip(i)-ip(i+1))*W*min(ip(i),ip(i+1))*K/l; 

// thick aquifer assumption
ll=(ip(i)-ip(i-1))*W*T*K/l;
rr=(ip(i)-ip(i+1))*W*T*K/l;

deltah(1,i)=(q(i)*W*l-ll-rr)/(W*l*Sy)*deltat;
// the basic transient calculation
end;

deltah(1,1)=0;
deltah(1,N)=0;

return;

