% Octave (and maybe Matlab compatible) program to simulate
% Junior Solar Sprint races.  This program needs the companion
% f1.m file to make the lsode call work.
%
% Larry Doolittle    <larry@doolittle.boa.org>
%   http://doolittle.faludi.com/~larry/
%
% Adapted from Chuck Wright's fine analysis at
%   http://chuck-wright.com/projects.html
%
% Octave is Free Software that can be downloaded from http://www.octave.org/

% Solar panel characteristics, these are wild guesses that more-or-less
% match Chuck Wright's pictures and his 3.3V and 1.35A numbers.  When our
% parts come in, I'll generate real parameters by measuring the actual
% solar array with some 1,2,3 Ohm resistors.
vt=0.026;    % mV
idss=0.7e-9; % Amps reverse leakage
n=6;         % number of cells in series
Rc=0.2;      % Ohms series resistance
isun=1.35;   % Amps in full sun
vraw=[0:.005:3.5];
i=isun-idss*(exp(vraw/(n*vt))-1);
vreal=vraw-i*Rc;

if (0)
 gset xrange [0:4];   gset xlabel 'Volts'
 gset yrange [0:1.5]; gset ylabel 'Amps'
 plot(vreal,i);
 return
end

% combine with the motor
Kt=0.0054;   % N-m/A or Volt-sec
Rm=0.8;      % Ohms series resistance
If=0.26;     % Amps to overcome friction

w=(vreal-i*Rm)/Kt;
torque=(i-If)*Kt;
rpm=60*w/2/pi;
if (0)
 gset xrange [0:6000]; gset xlabel 'speed (rpm)'
 gset yrange [0:.007]; gset ylabel 'torque (N-m)'
 plot(rpm,torque)
 return
end

% assemble the vehicle
M=0.7;       % kg
R=0.01;      % m effective wheel radius (actual * gear ratio)
accel=torque/R/M;
vel=w*R;
if (0)
 gset xrange [0:6];   gset xlabel 'speed (m/s)'
 gset yrange [0:1.5]; gset ylabel 'accel (m/s^2)'
 plot(vel,accel);
 return
end

ix=find((accel>-0.1));
vel=vel(ix);
accel=accel(ix);

nrange=20;
delv=max(vel)/nrange;
pf=[1:nrange]'*[0:2];
global pf delv
for k= [1:nrange]
   ik=find((vel>(k-1.1)*delv).*(vel<(k+.1)*delv));
   pf(k,:)=polyfit(vel(ik),accel(ik),2)';
end

time=[0:.1:8];
vrace=lsode("f1",[0 0],time);
if (1)
 gset autoscale; gset xlabel 'Time (s)'; gset ylabel ''
 plot(time,vrace(:,1),";speed (m/s);",time,vrace(:,2),";distance (m);");
end