diff --git a/motor_analysis/350kvTP.PNG b/motor_analysis/350kvTP.PNG
new file mode 100644
index 00000000..c0b987d9
Binary files /dev/null and b/motor_analysis/350kvTP.PNG differ
diff --git a/motor_analysis/350kvVelli.PNG b/motor_analysis/350kvVelli.PNG
new file mode 100644
index 00000000..86b63c35
Binary files /dev/null and b/motor_analysis/350kvVelli.PNG differ
diff --git a/motor_analysis/VelliPlot.m b/motor_analysis/VelliPlot.m
new file mode 100644
index 00000000..9d3b69c2
--- /dev/null
+++ b/motor_analysis/VelliPlot.m
@@ -0,0 +1,98 @@
+% limits
+Imax = 64; %A
+Umax = 22; %V
+Irange = 100; %A. Range for plotting current
+omegaMax = 2000; %rad/s mechanical. For plotting voltage ellipses
+omegastep = 200; %rad/s mechanical. For plotting voltage ellipses
+
+
+%%
+%350 kv motor
+lambda = 2.24/1000;
+L = 23e-6;
+R = 32e-3;
+pp = 7; Poles = pp*2;
+Ld = L;
+Lq = L;
+
+% %%
+% % Donkey
+% kv = 820;
+% lambda = 60/(kv*2*pi*pp*sqrt(3));
+% L = 8e-6; %Guess! TODO: measure
+% R = 30e-3; %Guess! TODO: measure
+% pp = 7; Poles = pp*2;
+% Ld = L;
+% Lq = L;
+
+%%
+Istep = Irange/400;
+Idplt = repmat(-Irange:Istep:Irange,801,1);
+Iqplt = repmat((-Irange:Istep:Irange)',1,801);
+UmaxSq = (Umax/sqrt(3))^2;
+
+t = linspace(0,2*pi);
+IdMaxt = Imax*cos(t);
+IqMaxt = Imax*sin(t);
+
+%%
+figure(1)
+plot(IdMaxt, IqMaxt);
+
+hold on;
+
+%Plot torque
+%[c,h] = contour(Idplt,Iqplt, (Poles/2).*(3/2).*(lambda.*Iqplt + (Ld-Lq).*Iqplt.*Idplt), -5:0.25:5);
+%clabel(c,h,'LabelSpacing',500);
+
+%Plot voltage ellipses
+EllRHS = Ld^2.*(lambda/Ld + Idplt).^2 + Lq^2.*Iqplt.^2;
+omegaAtEllipse = sqrt(UmaxSq./EllRHS);
+[c,h] = contour(Idplt,Iqplt, omegaAtEllipse./(Poles/2), 0:omegastep:omegaMax);
+clabel(c,h,'LabelSpacing',500);
+xlabel 'Id (A)'
+ylabel 'Iq (A)'
+colormap(jet)
+c = colorbar;
+%c.Label.String = 'Speed (rad/s)';
+ylabel(c,'Speed (mechanical rad/s)')
+grid on
+axis equal
+
+%plot infinite speed point
+plot(-lambda/Ld, 0, 'r*');
+
+hold off;
+
+%%
+figure(2)
+
+%fw range
+t = linspace(pi/2,pi);
+IdMaxt = Imax*cos(t);
+IqMaxt = Imax*sin(t);
+Tmaxt = (Poles/2).*(3/2).*(lambda.*IqMaxt + (Ld-Lq).*IqMaxt.*IdMaxt);
+EllRHSmaxt = Ld^2.*(lambda/Ld + IdMaxt).^2 + Lq^2.*IqMaxt.^2;
+omegamaxt = sqrt(UmaxSq./EllRHSmaxt);
+
+%MTPA range
+MTPAmaxtId = 0; %TODO make work for salient machines
+MTPAmaxtIq = Imax;
+TmaxtMTPA = (Poles/2).*(3/2).*(lambda.*MTPAmaxtIq + (Ld-Lq).*MTPAmaxtIq.*MTPAmaxtId);
+Tmaxt = [repmat(TmaxtMTPA, 1, 100) Tmaxt];
+omegamaxt = [linspace(0,omegamaxt(1)) omegamaxt];
+
+%Present mechanical speed
+omegamaxt = omegamaxt./(Poles/2);
+
+Pmaxt = Tmaxt.*omegamaxt;
+
+%plotyy(t, Tmaxt, t, omegamaxt);
+%plotyy(t, Tmaxt, t, Pmaxt);
+%plotyy(t, Pmaxt, t, omegamaxt);
+
+h = plotyy(omegamaxt, Tmaxt, omegamaxt, Pmaxt);
+grid on
+xlabel 'Speed (mechanical rad/s)'
+ylabel 'Torque (Nm)'
+ylabel(h(2), 'Power (W)');
\ No newline at end of file
diff --git a/numeric_path_opt/Main.m b/numeric_path_opt/Main.m
new file mode 100644
index 00000000..d41c1dfd
--- /dev/null
+++ b/numeric_path_opt/Main.m
@@ -0,0 +1,43 @@
+%Params
+AccelPerA = 3000;
+phaseR = 0.033;
+Ts = 0.001;
+N = 50;
+thetaFinal = 200;
+x0 = [0;0];
+
+%system definition
+%X = [theta; omega]
+Ac = [0 1;
+     0 0];
+Bc = [0;
+     AccelPerA];
+C = 0;
+D = 0;
+
+SYSC = ss(Ac, Bc, [], []);
+SYSD = c2d(SYSC, Ts, 'zoh');
+[Phi, Gamma] = predictionmatrices(SYSD.a, SYSD.b, SYSD.c, N);
+
+Df = Gamma(end-1:end,:);
+ff = [thetaFinal; 0];
+
+H = 2*eye(N)*(3/2)*phaseR*Ts;
+[u, fval] = quadprog(H,[],[],[],Df,ff);
+
+xv = reshape(Gamma*u, 2,N)';
+x = xv(:,1);
+v = xv(:,2);
+
+kv350_lambda = 2.2e-3;
+power = u.*v.*(3/2)*kv350_lambda;
+
+Vbus = 24;
+Ib = power/Vbus;
+Im = u;
+duty = Ib./Im;
+CapIsqr = (duty.*(Ib-Im)).^2 + ((1-duty).*Ib).^2;
+CapIrms = sqrt(sum(CapIsqr)/N);
+CapR = 0.08;
+Ncaps = 8;
+Cappow = (CapIrms/Ncaps)^2 * CapR
diff --git a/numeric_path_opt/predictionmatrices.m b/numeric_path_opt/predictionmatrices.m
new file mode 100644
index 00000000..929f7464
--- /dev/null
+++ b/numeric_path_opt/predictionmatrices.m
@@ -0,0 +1,14 @@
+function [Phi, Gamma, Lambda] = predictionmatrices(A, B, C, N)
+%UNTITLED2 Summary of this function goes here
+%   Detailed explanation goes here
+
+n = size(A,1);
+Atilde = [A; zeros((N-1)*n, n)];
+k = [zeros(n, N*n); -kron(eye(N-1), A) zeros((N-1)*n,n)] + eye(N*n);
+
+Phi = k\Atilde;
+Gamma = k\kron(eye(N), B);
+Lambda = kron(eye(N), C);
+
+end
+
diff --git a/various/path.c b/various/path.c
new file mode 100644
index 00000000..341e5c4b
--- /dev/null
+++ b/various/path.c
@@ -0,0 +1,47 @@
+#include <stdio.h>
+
+#include <math.h>
+
+void output(float x, float y, float str, float param) {
+  printf("{ %05.1ff, %05.1ff, %1.1ff, %06.1ff },\n", (x+1)*75, (y+1)*75, str, param);
+}
+
+void searchat(float x, float y) {
+  output(x, y, 0.4, 500);
+  //output(x+0.02, y, 0.4, 1000);
+  //output(x, y+0.02, 0.4, 1000);
+  //output(x-0.02, y, 0.4, 1000);
+  //output(x, y-0.02, 0.4, 1000);
+}
+
+void pickupto(float x, float y) {
+  output(x, y, 0.4, 0);
+  //output(x, y, 0.4, 1000);
+  output(x, y, 0, 1000);
+}
+
+void moveto(float x, float y) {
+  output(x, y, 0, 0);
+}
+
+#define POINTS 15
+#define TWOPI (3.1415*2.0)
+
+int main(int argc, char** argv) {
+  moveto(-1,-1);
+  
+  for (int j=0; j<2; j++)
+  for (int i=0; i<POINTS; i++) {
+    float angle = ((float)i/POINTS + 0.5 + 0.125) * TWOPI ;
+    if (!(i%2) != !j)
+      pickupto(cos(angle)*0.8, sin(angle)*0.8);
+    else {
+      moveto(cos(angle)*0.8, sin(angle)*0.8);
+      searchat(cos(angle)*0.8, sin(angle)*0.8);
+    }
+  }
+    
+  
+  moveto(-1,-1);
+  
+}