%PROGRAM NAME: moisture.m
tau=50; a=1/tau; b=0.1;   %Model parameters
rbar=10; delt=20; %Rainfall parameters
%-----------------
s=tf('s');
G=b/(s+a); %Model transfer function
dt=.01; tmax=300;
tvec=0:dt:tmax-dt;tvec=tvec'; %Time array
ntmax=tmax/dt;
ndelt=delt/dt;
% FORCED RESPONSE:
r=zeros(1,ntmax);
r(1:ndelt)=rbar*ones(1,ndelt); 
figure(1)
plot(tvec,r/rbar,'k','LineWidth',2)
hold on
xr=lsim(G,r,tvec);
plot(tvec,xr,'LineWidth',2)
%INITIAL CONDITION RESPONSE:
x0=1.0;
xic=x0*exp(-a*tvec);
hold on
plot(tvec,xic,'g','LineWidth',2)
x=xr+xic;
plot(tvec,x,'r','LineWidth',2)
legend('rain','forced','initial condition','total')
xlabel('Time (min)')
ylabel('Moisture (units)')
title(['Moisture Response to i.c. + Rain w duration ',num2str(delt),'min.'])
grid

