> 
> 
> ########## - Chaos Model by Wyatt Houtz-1998 - ##########
> 
> #Initial Conditions
> 
> restart;
> with(plots):
> 
> n:=4000:                #Number of steps
> q:=9:                   #Heat constant
> he:=14:                 #Heat of Engine
> hr:=14:                 #Heat of Radiator
> tau:=3:                 #Delay Steps
> c:=10:                  #Constant
> dt:=0.025:
> el:=floor(tau/dt):
> slow_omega:=1:         #1 if use omega, 0 if use omega n+1
> 
> 
> 
> 
> #########FUNCTIONS#######
> 
> Rbottom:=8:
> Rslope:=5:
> Rbottom_s:=convert(Rbottom, string):
> Rslope_s:=convert(Rslope, string):
> R:=x -> Rslope*(x-Rbottom):
> Rtop:=Rbottom+1/Rslope:
> 
> 
> #Left Function with upper and lower bounds at 0 and 1
> Lbottom:=5:
> Lslope:=5:
> Lbottom_s:=convert(Lbottom, string):
> Lslope_s:=convert(Lslope, string):
> L:=x -> Lslope*(x-Lbottom):
> Ltop:=Lbottom+1/Lslope:
> 
> 
> Fr:=x -> piecewise(x>=Rbottom and x<=Rtop,R(x),x>Rtop,1,x Fl:=x -> piecewise(x>=Lbottom and x<=Ltop,L(x),x>Ltop,1,x 
> 
> 
> 
> 
> 
> #########ARRAYS##########
> 
> Time:=array(-el..n):
> theta:=array(-el..n):
> Te:=array(-el..n):
> Tr:=array(-el..n):
> omega:=array(-el..n):
> 
> ######Initial Array Conditions########
> 
> for i from -el to 0 do
> theta[i]:=0:
> Te[i]:=0:
> Tr[i]:=0:
> omega[i]:=0:
> od:
> 
> for i from -el to n do
> Time[i]:=i*dt:
> od:
> 
> 
> ########Calculations with delay calculations########
> 
> 
> for i from 1 to n do
> 
> 
> theta[i]:= (1 - he*dt)*theta[i-1] + he*dt*Te[i-1]:
> 
> 
> if (theta[i] = theta[i-1]) then
>    omega[i]:=omega[i-1]
> fi;
> 
> if (theta[i] > theta[i-1]) then
>    if (Fr(theta[i]) >= omega[i-1]) then
>          omega[i]:= Fr(theta[i])
>    else
>          omega[i]:=omega[i-1]
>    fi
> fi;
> 
> if (theta[i] < theta[i-1]) then
>    if (Fl(theta[i]) <= omega[i-1]) then
>          omega[i]:= Fl(theta[i])
>    else
>          omega[i]:=omega[i-1]
>    fi
> fi;
> 
> Te[i]:=(1 - c*dt)*Te[i-1] + c*dt*omega[i-slow_omega]*(Tr[i-el] - Te[i-el]) + c*dt*Te[i-el] + dt*q:
> 
> Tr[i]:=(1 - (omega[i-slow_omega] + hr)*dt)*Tr[i-1] + omega[i-slow_omega]*dt*Te[i-el]:
> od:
> 
> 
> 
> 
> #####Graph/Data Review#######
> 
> Tmin:=0:
> Tmax:=n*dt:
> time_int:=[floor(Tmin/dt),ceil(Tmax/dt)]:
> step_size:=convert(dt,string):
> key_parameter:=convert(Rbottom, string):
> 
> 
> ####Functions##########
> #RTPlot:= Rtop + 5:
> #RBPlot:= Lbottom - 5:
> 
> #range_low:=min(seq(Fl(i),i=time_int[1]..time_int[2])):
> #range_high:=max(seq(Fl(i),i=time_int[1]..time_int[2])):
> #n_high:=RTPlot:
> #header:=cat(`TL = `,Lbottom_s,`, TR = `,Rbottom,`, Lm = `,Lslope_s,`, Rm = `,Rslope_s):
> #Fr_plot:=[seq([i, Fr(i)], i=2..n)]:
> #Fl_plot:=[seq([i, Fl(i)], i=2..n)]:
> #plot([Fr_plot,Fl_plot], view=[RBPlot..RTPlot, 0..1], thickness=1, color=orange, title=header);
> 
> 
> #####Review Interval of Theta#######
> 
> range_low:=min(seq(theta[i],i=time_int[1]..time_int[2])):
> range_high:=max(seq(theta[i],i=time_int[1]..time_int[2])):
> time_low:=convert(Tmin,string):
> time_high:=convert(Tmax,string):
> header:=cat(`Theta for dt =`,step_size,`, `, `Parameter = `,key_parameter):
> theta_plot:=[seq([Time[i], theta[i]], i=-el..n)]:
> plot(theta_plot, view=[Tmin..Tmax,range_low..range_high], color=black, thickness=1, title=header);
> 
> 
> #####Review Interval of Te #######
> 
> range_low:=min(seq(Te[i],i=time_int[1]..time_int[2])):
> range_high:=max(seq(Te[i],i=time_int[1]..time_int[2])):
> time_low:=convert(Tmin,string):
> time_high:=convert(Tmax,string):
> header:=cat(`Te for dt =`,step_size,`, `, `Parameter = `,key_parameter):
> Te_plot:=[seq([Time[i],Te[i]], i=-el..n)]:
> plot(Te_plot, view=[Tmin..Tmax, range_low..range_high], color=black, thickness=1, title=header);
> 
> 
> #####Review Interval of Tr #######
> 
> range_low:=min(seq(Tr[i],i=time_int[1]..time_int[2])):
> range_high:=max(seq(Tr[i],i=time_int[1]..time_int[2])):
> time_low:=convert(Tmin,string):
> time_high:=convert(Tmax,string):
> header:=cat(`Tr for dt =`,step_size,`, `, `Parameter = `,key_parameter):
> Tr_plot:=[seq([Time[i], Tr[i]], i=-el..n)]:
> plot(Tr_plot, view=[Tmin..Tmax, range_low..range_high], color=black, thickness=1, title=header);
> 
> 
> #####Review Interval of Omega #######
> 
> range_low:=min(seq(omega[i],i=time_int[1]..time_int[2])):
> range_high:=max(seq(omega[i],i=time_int[1]..time_int[2])):
> time_low:=convert(Tmin,string):
> time_high:=convert(Tmax,string):
> header:=cat(`Omega for dt =`,step_size,`, `, `Parameter = `,key_parameter):
> omega_plot:=[seq([Time[i], omega[i]], i=-el..n)]:
> plot(omega_plot, view=[Tmin..Tmax, 0..1.1], color=black,thickness=1, title=header);
> 
> 
> #####Check Numerical Values#####
> 
> #for i from n-100 to n do
> #[i,theta[i],omega[i],Te[i],Tr[i]];
> #od:
> 
> 
>