My graph is not picking up my time increments, its suppose to have a nice curve, but its coming out as a straight line.

clear
Ts=30; % ms
Td=60; % ms
Chs=.001; % L/mmhg
Chd=.015; % L/mmhg
N=800; % number of elements
t=0:1:800; % all time values together
for i=0:N
    if i<=99
        Ch(i+1)=Chd; % starts at .015 L/mmhg in chd phase
    elseif i>99 & i<=349 
        Ch(i+1)=(Chd-Chs)*exp((-t(i))/Ts)+Chs; % switches to chs
    elseif i>349
        Ch(i+1)=(Chs-Chd)*exp((-t(i))/Ts)+Chd; % switches back to chd
    end
end
timebase=0:.01:8; % putting back in s
plot(timebase,Ch)
grid on
ylim([0 .02])
xlabel('time(s)')
ylabel('Ch (L/mmHg')

It's suppose to look like this.

Answer 1

[William Rose]

@BAILEY MCMASTER,

clear
Ts=30; % ms
Td=60; % ms
Chs=.001; % L/mmhg
Chd=.015; % L/mmhg
N=800; % number of elements
t=0:1:800; % all time values together
for i=0:N
    if i<=99
        Ch(i+1)=Chd; % starts at .015 L/mmhg in chd phase
    elseif i>99 & i<=349 
        Ch(i+1)=(Chd-Chs)*exp(-(t(i)-99)/Ts)+Chs; % switches to chs
    elseif i>349
        Ch(i+1)=(Chs-Chd)*exp(-(t(i)-349)/Ts)+Chd; % switches back to chd
    end
end
timebase=0:.01:8; % putting back in s
plot(timebase,Ch)
grid on
ylim([0 .02])
xlabel('time(s)')
ylabel('Ch (L/mmHg')

Good luck.

[William Rose]

@BAILEY MCMASTER, note the lines I changed are

Ch(i+1)=(Chd-Chs)*exp(-(t(i)-99)/Ts)+Chs; % switches to chs

and

Ch(i+1)=(Chs-Chd)*exp(-(t(i)-349)/Ts)+Chd; % switches back to chd

Good luck with your modeling.

[BAILEY MCMASTER]

Thank so much, the time change make sense because its the intial change for the start of the new phase. Thank you!

[William Rose]

@BAILEY MCMASTER, I basically copied your equations when I adjusted your script, but I added "-99" and "-349". I suspect you intended to have different time constants, Ts and Td, for the fall and the rise. Fix if you wish.

[William Rose]

@BAILEY MCMASTER, you are welcome. The answer from @Star Strider is more elegant than my answer. @Star Strider's answer modifies your code in more significant ways. You and I can both learn from it.

[William Rose]

@BAILEY MCMASTER,

Your equation for the second phase (rising phase) does something which I suspect is not what you want: It starts the rising phase at Ch()=Chs. This is not a problem with the time constants and start times in this particular case, since the signal Ch() is very close to Chs before the rising phase begins. But what if the timing was such that Ch() had not yet decayed to Chs? With your current method, Ch() would instantly jump down to Chs, to begin its rising phase. That would be weird. See below, where I have changed Ts from 30 to 150 ms to illustrate my point.

clear
Ts=150; % 150 ms, not 30
Td=60; % ms
Chs=.001; % L/mmhg
Chd=.015; % L/mmhg
N=800; % number of elements
t=0:1:800; % all time values together
for i=0:N
if i<=99
    Ch(i+1)=Chd; % starts at .015 L/mmhg in chd phase
elseif i>99 & i<=349 
    Ch(i+1)=(Chd-Chs)*exp(-(t(i)-99)/Ts)+Chs; % switches to chs
elseif i>349
    Ch(i+1)=(Chs-Chd)*exp(-(t(i)-349)/Ts)+Chd; % switches back to chd
end
end
timebase=0:.01:8; % putting back in s
plot(timebase,Ch)
grid on
ylim([0 .02])
xlabel('time(s)')
ylabel('Ch (L/mmHg')

A better approach is to have Ch() begin its rising phase from whatever value it currently has.

Answer 2

[Star Strider]

Vectorised version —

clear
Ts=30; % ms
Td=60; % ms
Chs=.001; % L/mmhg
Chd=.015; % L/mmhg
N=800; % number of elements
t=0:1:800; % all time values together

Lv1 = t<=99;
Ch(Lv1) = Chd;
Lv2 = t>99 & t<=349;
Ch(Lv2) = (Chd-Chs)*exp(-t(1:nnz(Lv2))/Ts)+Chs; % switches to chs
Lv3 = t > 349;
Ch(Lv3) = (Chs-Chd)*exp(-t(1:nnz(Lv3))/Ts)+Chd; % switches back to chd

timebase=0:.01:8; % putting back in s
plot(timebase,Ch)
grid on
ylim([0 .02])
xlabel('time(s)')
ylabel('Ch (L/mmHg')

.