What could micro-propulsion do?

I have been following the NASA developments on the EM Drive with great interest, and so I was very interested to see that Paul March (a NASA Eagleworks engineer) today stated that he is measuring a force of over 100µN for 80W of electrical power.

I was curious to put this into context. A quick google reveals that you can purchase solar panels for spacecraft that produce 300W/kg.

Let’s say that you can scale up the solar panels such that the mass of the engine is negligible compared to the solar panel mass. (A big assumption)

Taking his Paul March’s figure of 100µN/80W at face value, 300W per kg solar panels would produce a force per kg of = 0.0004 Newtons/kg, so an acceleration of:

a = F/m = 0.0004 m / s^2

This sounds extremely small, but we can accelerate constantly over long periods of time.

After 1 year of this acceleration, we would have a delta-velocity, dV, of:

dV = at = 12 km/s  (43,000 km/hour)

A delta-V of 12km/s isn’t shabby at all. Here’s a delta-V map of our solar system:

DeltaV map of our solar system

So from Low Earth orbit, a dV of 12km/s gets us to the moon in very low orbit, and coming back again! (2.44+0.68+0.14+0.68+1.73 = 5.67)

(of course you can’t actually land and take off since the force is less than that of gravity at the surface, but you could do a very low fly-by).

So we could have a shuttle going to the moon and back every year. Not sure if that’s useful, but I’m sure someone could think of a use!

What else could we do with it?

From that delta-V map, we could get to a low orbit around Venus in about 6 months (2.44+0.68+0.09 + 0.28+0.36+2.94 = 6.79 km/s).

We could get to pretty much any planet in up to a year. (e.g. low orbit Neptune is 2.44+0.68+0.09+0.39+2.7+0.99+0.69+0.27+0.35+6.75 = 15.35 km/s )

Imagine what we could do with a spacecraft capable of visiting the furthest planets and returning in just over a couple years.

And that’s all without any sort of gravitational slingshots!

Effect of moving away from the Sun

A reader pointed out that the solar energy collected would of course decrease as you moved away from the Sun.

We can do a quick calculation – consider a round trip to Mars and back.

Distance from Sun to Mars: 1.38 au

(1/1.38)^2 = 52%

So the power generated by the solar panels would be reduced by 50%. That means, roughly, on average you’d get 75% of the acceleration that I mentioned, and so would increase the round trip from lower earth orbit to lower Mars orbit from 0.95 years to 1.3 years.

(I’ve obviously linearised where I shouldn’t have, the effect being an underestimation of the time, but you can get a crude back-of-the-envelope feeling for this.)

Solar Panels or Nuclear Power?

Nasa’s “Safe Affordable Fission Engine” can generate 100kW with 500Kg mass.  So 200W/kg.  Less than the 300W/kg that the solar panels can produce near earth, and on par with the average for a trip to Mars.  For a trip much further than the Mars then it starts to make sense.

Batteries?

Batteries store around 250Wh/kg. So 1kg can generate 250 watts for one hour… which is what a solar panel can do continuously.

Theoretical maximums

If it could be scaled to produce 1g acceleration, a human could explore this galaxy and even a neighbouring galaxy within their lifetime: http://math.ucr.edu/home/baez/physics/Relativity/SR/Rocket/rocket.html

The problem however is that at high velocities, space dust becomes more and more dangerous. At even high velocities, the very bumps in the spacetime from gravity waves becomes dangerous.

USB i2c convertor chip. A lesson in USB Latency.

I was interested in transmitting data from an i2c sensor to an ARM board running FreeRTOS along a long wire, and decided to use a i2c to i2c converter. Specifically a chip called MCP2221.

A single sensor reading from the i2c sensor is 2 bytes. I wanted to read 128 samples per second, so I need 256 bytes per second of data.

I figured we’d be able to do easily, what with i2c at 400khz and USB at 12 M/s. There is considerable overhead – each sensor reading needs four USB packets:

  1. Send USB packet from arm board to usb chip telling it to read from i2c sensor
  2. Receive ACK from USB chip
  3. Send USB packet from arm board to usb chip telling it to send us the data that it just read
  4. Receive the data read

Each USB packet is 64 bytes, so we have an overhead of 256 bytes for every 2 bytes of data that we want to read. This means we need a throughput of (256 bytes per sample) * (128 samples per second) = 32kb/sec.

Trying to get 32kb/s out of a 1.5MB/s USB link should be easy, I naively thought.

I wrote a driver for the chip for FreeRTOS, and tested it out.

I got a miserable 40 samples per second. Far out from the target of 128 samples per second.

So what was going wrong? Obviously I suspected a mistake with my driver. So first step is to try to reproduce on linux. I downloaded the MCP2221 Linux driver, and tested the speed of it, like so:

#!/usr/bin/env python
import smbus
import time
import sys

CONTINUOUS_MODE = True
I2C_ADDRESS = 0x48

if (len(sys.argv) == 1):
print("Please specify the bus, like:\n sudo i2c_sensor_read.py $(sudo i2cdetect -l | sed -ne 's/^i2c-\([0-9]\+\).*i2c-mcp2221.*$/\1/p')");
exit(1);

bus = smbus.SMBus(int(sys.argv[1]))

def getContinuousValue():
values = bus.read_i2c_block_data(I2C_ADDRESS, 0x0);
value = (values[0] << 8) + values[1]
return value;

last_epoch_time = 0

write_bytes([0x01, 0x4, 0x83]);
while True:
epoch_time = int(time.time()*1000);
epoch_diff = epoch_time - last_epoch_time;
last_epoch_time = epoch_time;
hz = 0;
if (epoch_diff!=0):
hz = int(1000/epoch_diff);
print str(getContinuousValue()) + " (took " + str(int(epoch_diff)) + " ms - " + str(hz) + " per second)";

This produced a miserable 62.5 samples per second. i.e. 16ms. I even tested on windows and got the same result.

So what was going on?

I tested with a oscilloscope. It was taking 1ms to send a USB packet to receive the reply, compared to the 8ms for a roundtrip that we see in the driver.

This 8ms for a roundtrip is 4ms average for sending and receiving a packet. The code in the driver (after I modified it) looks like:

static int mcp2221_ll_cmd(struct i2c_mcp2221 *dev)
{
       int rv;

       /* tell everybody to leave the URB alone */
       dev-&gt;ongoing_usb_ll_op = 1;

       /* submit the interrupt out ep packet */
       if (usb_submit_urb(dev-&gt;interrupt_out_urb, GFP_KERNEL)) {
               dev_err(&dev-&gt;interface-&gt;dev,
                               "mcp2221(ll): usb_submit_urb intr out failed\n");
               dev-&gt;ongoing_usb_ll_op = 0;
               return -EIO;
       }

       /* wait for its completion */
       rv = wait_event_interruptible(dev-&gt;usb_urb_completion_wait,
                       (!dev-&gt;ongoing_usb_ll_op));
       if (rv interface-&gt;dev, "mcp2221(ll): wait interrupted\n");
               goto ll_exit_clear_flag;
       }

       /* tell everybody to leave the URB alone */
       dev-&gt;ongoing_usb_ll_op = 1;

       /* submit the interrupt in ep packet */
       if (usb_submit_urb(dev-&gt;interrupt_in_urb, GFP_KERNEL)) {
               dev_err(&dev-&gt;interface-&gt;dev, "mcp2221(ll):
usb_submit_urb intr in failed\n");
               dev-&gt;ongoing_usb_ll_op = 0;
               return -EIO;
       }

       /* wait for its completion */
       rv = wait_event_interruptible(dev-&gt;usb_urb_completion_wait,
                       (!dev-&gt;ongoing_usb_ll_op));
       if (rv interface-&gt;dev, "mcp2221(ll): wait interrupted\n");
               goto ll_exit_clear_flag;
       }

ll_exit_clear_flag:
       dev-&gt;ongoing_usb_ll_op = 0;
       return rv;
}

I spoke to the Alan Stern and Greg KH, the Linux kernel USB maintainers, who were not surprised by this 8ms delay at all, and told me:

Why are we seeing such a large latency?
I don’t see anything wrong here, USB isn’t guaranteed to have any specific latency, what you are doing here is the worst-possible-case for a USB system (i.e. send a packet, wait for it to complete, send another one, etc.) There are lots of ways to get much better throughput if that is what you are wanting to do.

Can we make this faster?

For a horrid protocol like this, no, there isn’t any way to make it go faster, sorry. That is designed in such a way to make it the worst possible thing for a USB system, go kick the person who designed such a thing (hint, they have no idea how USB works…)

I strongly suggest going and getting a different sensor chip, don’t encourage such behaviour by actually buying their hardware.

Where is the delay coming from?

Multiple places: time to submit the request, time to reserve bandwidth for the previously unused interrupt endpoint, time to complete the transfer, all multiplied by 2.

You can get more information from usbmon (see Documentation/usb/usbmon.txt in the kernel source). But Greg is right; the protocol you described is terrible. There’s no need for a multiple ping-pong interchange like that; all you should need to do is wait for the device to send the next bit (or whatever) of data as soon as it becomes available.

Balancing on a stick

Consider a pole of length h , with a motor and flywheel attached to one end.  The motor can spin up the flywheel.  The other end is on a non-slip floor.

wheel

If the stick is at angle \theta , and the total mass of motor M_m plus mass of flywheel, M_f is M_m + M_f = M_t , the torque due to gravity, \tau_g is:

\tau_g = M_t gh \cdot \sin(\theta)

For a solid and hollow flywheel, the moment of inertia is:

I_s = \dfrac{M_f r^2}{2} and I_h = M_f r^2

If the motor applies torque \tau_m, then the angular acceleration of the flywheel is:

\alpha = \dfrac{d\omega(t)}{dt} = \dfrac{\tau_m}{I}

And this will require a power at time t, P(t) of:

P(t) = \tau_m \cdot \omega(t) = \tau_m \cdot \int{\dfrac{\tau_m}{I}}dt

If we start the flywheel from rest, and keep the torque constant:

P(t) = \dfrac{\tau_m^2}{I}\cdot t

(In units, we’ve got on the left: J/s = Nm/s and on the right (Nm)^2 (kg^{-1} m^{-2}) s = Nm (kg m^2 s^{-2}) (kg^{-1} m^{-2}) s = Nm/s, thus the units match and check our work)

So the power required to keep a given torque increases linearly with time.

Notice that I \propto M_f and the torque required to counter gravity is \tau_m \propto M_t , so if the mass is dominated by the flywheel then P \propto M_t .

So we want to minimize the mass of the flywheel (while keeping the assumption that the flywheel mass dominates the total mass) while maximizing the moment of inertia.

That means that we want to use a hollow flywheel, since this has twice the moment of inertial for a given mass.

Combining the equation P(t) with the equation for \tau_m, we get

P(t) = \dfrac{\tau_m^2}{I}\cdot t = \dfrac{(M_t gh \cdot \sin(\theta))^2}{I}\cdot t = \dfrac{M_t^2 g^2h^2 \cdot \sin^2(\theta)}{I}\cdot t

And substituting in the equation we had for $I_h = M_f r^2$:

P(t) = \dfrac{M_t^2 g^2h^2 \cdot \sin^2(\theta)}{M_f r^2}\cdot t

If we want to minimize the power, then we want to minimize M_t^2/M^f = (M_m + M_f)^2/M_f. Differentiating gives (M_f^2 - M_m^2)/M_f^2 which has a zero at M_f = M_m .

So the minimum power required is when M_f = M_m !

Three flywheels

If we consider the case of having two or three flywheels on the box, perpendicular to each other to give us more control, then we can consider each flywheel individually, and consider M_m to include the mass of the other flywheels. For minimum power usage, this would mean that we’d want to minimize (M_m + 3M_f)^2/M_f which is minimum at $M_f = M_m/3$.

So if we have 1kg of motors, payload, etc, then we’d want the flywheels to be 1/3 kg each.

Radius of flywheel

Since P \propto 1/r^2 we want to maximize the radius of the flywheel.

Example

A random example EDF motor has weight 148 grams (which we’ll round up to 0.2kg) and produces 450 Watts of power. We’ll consider this on a stick of length 1 meter and with a hollow flywheel.

If we begin at an angle of 10 degrees (0.17 radians) and desire a torque sufficient to balance against gravity:

\tau_g = (M_f + M_m) gh \cdot \sin(\theta) = (M_f + M_m) \cdot (1m * 10m/s^2) \cdot 0.17 = (M_f + M_m) \cdot (1.7m^2/s^2)

So:

P(t) = \dfrac{\tau_m^2}{I}\cdot t = \dfrac{(M_f + M_m)^2 \cdot (2.89 m^4/s^4)}{M_f r^2} \cdot t

Let’s make the radius 0.1m – i.e. 10% of the length of the stick. We’ll consider the case of using a single flywheel and minimizing the power usage.

So this gives a power of P(t) = (2M_f)^2/M_f \cdot (0.0289 m^2/s^4) \cdot t = 4 M_f \cdot (0.0289 m^2/s^4) \cdot t = M_f \cdot (0.1156 m^2/s^4) \cdot t

Using the EDF motor mass of 0.2kg:

P(t) = (0.023 kg m^2/s^4) \cdot t = (0.023 J/s^2) \cdot t

The motor can provide a maximum of 450 Watts, meaning that it could maintain a 10 degree angle for 450/0.023 s = 19565 seconds = 5.4 hours.

Lie Detector

I put together a ‘lie detector’ testing for galvanic skin response.

I didn’t work.  I tested various obvious different ideas, such as looking at variances from a sliding window mean etc, but didn’t get anywhere.

It does light up if I hyperventilate though, which is er, useful?  It is also, interestingly, lights up when my daughter uses it and shouts out swear words, which she finds highly amusing.  I think it’s just detecting her jerking her hand around though.  She doesn’t exactly sit still.

IMG_20150627_020447

I used an op-amp to magnify the signal 200 times analogue-y (er, as opposed to digitally..), but I now wonder if I could use a high resolution (e.g. 40bit) digital to analog converter directly, and do more processing in software.

Perhaps even use a neural network to train to train it to detect lies.

I did order a 40bit AD evaluation board from TI, but haven’t had the chance to actually use it yet.

Minimum time to run a test case

A friend told me of a non-deterministic bug that he had found that only crashed occasionally.  He wrote a test case, ran it once, and the bug appeared after 3.5 hours of testing.  He sent the bug report and test case to code owners, who ran it for 4 hours but couldn’t reproduce the bug.  Intuitively we can see that they might have just been unlucky and thus not seen the crash.

So it got me wondering what the minimum amount of time that they need to run the test for, all other factors being equal, to be 95%  likely to reproduce the bug again. (i.e. 2 sigma)

Assuming that the bug appears randomly, with an equal probability of x per hour, then after 3.5 hours of running, there’s a probability p of the bug appearing of:

(1+x)^{3.5} - 1 = p

Rearranging:

(1+x)^{3.5} = 1+ p

1+x = (1+p)^{1/3.5}

x = (1+p)^{1/3.5} - 1

To have the tester reproduce the results in the end with a 95% chance, we need to 90% confidence in the probability of us producing the bug in the first place, and a 90% confidence of the reproducer reproducing the bug, so that combined we get a sqrt(0.90) = 0.95 confidence in overall result.  So the percentage chance, p , of the observed outcome of seeing the crash after 3.5 hours is between 10% to 90%.  Setting p to 0.1 and 0.9 in the above equation, we get that the probability of the bug appearing per hour is between 2.8%  to 20%.

Taking the lowest probability, to get a total of 95% chance of reproducing the bug (and thus 90% chance of reproducing the bug GIVEN the probability of the bug appearing for us being 90%) we need to run for:

(1+0.028)^y = 1.9

y = \log(1.9)/\log(1.028)

y = 23 hours.

So to be 95% confident that the bug does not appear on the reproducer’s system, they would need to run the test case for 23 hours, assuming similar hardware etc.

This can be drastically brought down if the initial tester does a second run, to increase the value for p .

Rerunning the test

The original tester reran the test himself and reproduced the bug 4 times and found that it always appeared in less than 6 hours.

To be 95% confidence in the overall results, and thus 90% confidence for just p,  we want the probability of the occurring 4 times in 4 runs in less than 6 hours to be 10%.  Thus the probability of it occurring in 1 run in less than 6 hours is simply:

1 - (1 - 0.1)^4 = 0.34

So setting p = 34%  and using the equation above but with the number of hours set to 6, we get:

x = (1+0.34)^{1/6} - 1 = 0.50

Which gives x = 5.0%. This means that with 90% confidence, the bug appears with a minimum probability of 5% per hour.  So the amount of time, y, that the reproducer needs to run to be 95% likely to reproduce the bug is:

y = \log(1.90)/\log(1.05)

y = 13 hours

D3.js in QML

I wanted to do a physics-based layout in QML, to have bouncing bubbles etc.

D3.js  has everything that I needed, but it needed to be tweaked to get it to run:

  • Create some dummy functions to stub the webbrowser javascript functions clearTimeout() and setTimeout()  that we don’t have
  • Create a QML Timer to replace them and manually call force.tick()
  • Use dummy objects for D3 to manipulate, and then set the position of our real graphics to the determined position of the dummy objects.  This is needed because D3 takes x,y coordinates to be the center of the object, whereas qml uses x,y to mean the top left.

It worked pretty well, and now I can do things like this in QML:

output


diff --git a/d3.js b/d3.js
index 8868e42..7aac164 100644
--- a/d3.js
+++ b/d3.js
@@ -1,4 +1,9 @@
-!function() {
-  var d3 = {
+function clearTimeout() {
+};
+function setTimeout() {
+};
+var d3;
+!function(){
+  d3 = {
     version: "3.5.5"
   };
@@ -9503,2 +9508,1 @@
-  this.d3 = d3;
-}();
\ No newline at end of file
+}();

(The ds.min.js file is also patched in the same way)

And in QML:


import QtQuick 2.2
import "."
import QtSensors 5.0
import "d3.js" as D3

Rectangle {
    id: bubbleContainer
    width: 1000
    height: 1000
    //clip: true

    property int numbubbles: 200;
    property real bubbleScaleFactor: 0.1
    property real maxBubbleRadius: width*0.15*bubbleScaleFactor + width/20*bubbleScaleFactor
    property variant bubblesprites: []

    function createBubbles() {
        var bubblecomponent = Qt.createComponent("bubble.qml");
        var radius = Math.random()*width*0.15*bubbleScaleFactor + width/20*bubbleScaleFactor
        bubblesprites.push(bubblecomponent.createObject(bubbleContainer, {"x": width/2-radius , "y": height/2 - radius, "radius":radius, "color":  "#ec7e78" } ));
        var xLeft = -radius
        var xRight = radius
        for (var i=1; i < numbubbles; ++i) {             var radius2 = Math.random()*width*0.15*bubbleScaleFactor + width/20*bubbleScaleFactor             var y = height*(0.5 + (Math.random()-0.5))-radius2             if (i % 2 === 0) {                 bubblesprites.push(bubblecomponent.createObject(bubbleContainer, {"x": xLeft, "y": y, "radius":radius2 } ));                 xLeft += radius2*2             } else {                 xRight -= radius2*2                 bubblesprites.push(bubblecomponent.createObject(bubbleContainer, {"x": xRight, "y": y, "radius":radius2 } ));             }         }     }     function boundParticle(b) {         if (b.y > height - b.radius)
            b.y = height - b.radius
        if (b.y < b.radius)             b.y = b.radius         if (b.x > width*2.5)
            b.x = width*2.5
        if (b.x < -width*2)
            b.x = -width*2
    }

    property real padding: width/100*bubbleScaleFactor;
    function collide(node) {
        var r = node.radius + maxBubbleRadius+10,
              nx1 = node.x - r,
              nx2 = node.x + r,
              ny1 = node.y - r,
              ny2 = node.y + r;
        boundParticle(node)
        return function(quad, x1, y1, x2, y2) {
            if (quad.point && (quad.point !== node)) {
              var x = node.x - quad.point.x,
                  y = node.y - quad.point.y,
                  l = Math.sqrt(x * x + y * y),
                  r = node.radius + quad.point.radius + padding;
              if (l < r) {                 l = (l - r) / l * .5;                 node.x -= x *= l;                 node.y -= y *= l;                 quad.point.x += x;                 quad.point.y += y;               }             }             return x1 > nx2 || x2 < nx1 || y1 > ny2 || y2 < ny1;
        };
    }

    property var nodes;
    property var force;

    Component.onCompleted: {
        initializeTimer.start()
    }

    Timer {
        id: initializeTimer
        interval: 20;
        running: false;
        repeat: false;
        onTriggered: {
            createBubbles();
            nodes = D3.d3.range(numbubbles).map(function() { return {radius: bubblesprites[this.index]}; });
            for(var i = 0; i < numbubbles; ++i) {
                nodes[i].radius = bubblesprites[i].radius;
                nodes[i].px = nodes[i].x = bubblesprites[i].x+bubblesprites[i].radius
                nodes[i].py = nodes[i].y = bubblesprites[i].y+bubblesprites[i].radius
            }
            nodes[0].fixed = true;

            force = D3.d3.layout.force().gravity(0.05).charge(function(d, i) { return 0; }).nodes(nodes).size([width,height])
            force.start()
            nodes[0].px = width/2
            nodes[0].py = height/2
            force.on("tick", function(e) {
                var q = D3.d3.geom.quadtree(nodes),i = 0,
                        n = nodes.length;
                while (++i < n) q.visit(collide(nodes[i]));
                for(var i = 0; i < numbubbles; ++i) {
                    bubblesprites[i].x = nodes[i].x - nodes[i].radius;
                    bubblesprites[i].y = nodes[i].y - nodes[i].radius;
                }
            });
            timer.start();
        }
    }

    Timer {
        id: timer
        interval: 26;
        running: false;
        repeat: true;
        onTriggered: {
            force.resume();
            force.tick();
        }
    }
}

And bubble.qml is trivial:



import QtQuick 2.0
Rectangle {
    color: "#f3bab3"
    radius: 5
    width: radius*2
    height: width
}

The full code is here

Health badge

Electronics

I bought a 128×128 SPI Nokia 5110 display (The ILI9163)  to hook up to an Arduino.  I used an Arduino Mini which runs at 3.3v, since the display uses 3.3v logic.

Then wired it up like:

- LED       -->     +3.3V  -  (5V would work too) 
- SCK       -->     Pin 13 -  Sclk,  (3.3V level)
- SDA       -->     Pin 11 -  Mosi,  (3.3V level)
- A0        -->     Pin 10 -  DC or RS pin (3.3V level)
- RST       -->     Pin 9  -  Or just connect to +3.3V
- CS        -->     Pin 8  -  CS pin (3.3V level)
- Gnd       -->     Gnd
- Vcc       -->     +3.3V

(Note that SDA and SCK pins aren’t optional, because we’re using the arduino SPI hardware to drive this. The DC and CS pins are changeable if you change the code)

layout1_bb

And found a library to use this:

cd /libraries/
git clone https://github.com/adafruit/Adafruit-GFX-Library
git clone https://github.com/sumotoy/TFT_ILI9163C

I had to then edit TFT_ILI9163C/TFT_ILI9163C.h to comment out the #define __144_RED_PCB__ , and uncomment the #define __144_BLACK_PCB__ since my PCB was black. Without this change I got corruption on the top 1/4 of the screen.

Now to get something up and running, using the above pin definition:

#include <SPI.h>
#include <Adafruit_GFX.h>
#include <TFT_ILI9163C.h>

/* SPI hardware runs on:
 MOSI:  11
 SCK:   13
 */
#define __CS 8
#define __DC 10
#define __RESET 9

TFT_ILI9163C tft = TFT_ILI9163C(__CS, __DC, __RESET);

const uint16_t WIDTH = 128;
const uint16_t HEIGHT = 128;

void setup(void) {
  tft.begin();
  tft.drawFastHLine(0, 0, WIDTH, 0xFFE0);
}

void loop() {
}

This now gives us a black sceen with a yellow line. Woohoo!

Now, I have a bunch of images that I got an artist to draw, that look like:

face10

I have 10 of these, and their raw size would be 128 * 128 * 4 * 10 = 655KB. Far too large for our 30KB of available flash memory on the nano!

But given the simplicity of the images, we can use a primitive compression algorithm and instead store these as Run Length Encoded. E.g. store as them as pixel index color followed by the number of pixels of that color in a line. We turn to trusty python to do this:

import os,sys
import itertools
from PIL import Image

rleBytes = 0
colorbits = 3

if len(sys.argv) != 3:
    print sys.argv[0], " <input file.png> <output file.h>"
    sys.exit()

def writeHeader(f):
    f.write("/* A generated file for a run length encoded image of " + sys.argv[1] + " */\n\n")
    f.write("#ifndef PIXEL_VALUE_TO_COLOR_INDEX\n")
    f.write("#define PIXEL_VALUE_TO_COLOR_INDEX(x) ((x) >> " + str(8-colorbits) + ")\n")
    f.write("#define PIXEL_VALUE_TO_RUN_LENGTH(x) ((x) & 0b" + ("1" *(8-colorbits)) + ")\n")
    f.write("#endif\n\n")
    f.write("namespace " + os.path.splitext(sys.argv[2])[0] + " {\n")
def writeFooter(f):
    f.write("}")

def writePalette(f, image):
    palette = image.im.getpalette()
    palette = map(ord, palette[:3*2**colorbits]) # This is now an array of bytes of the rgb565 values of the palette
    f.write("/* The RGB565 value for color 0 is in palette[0] and so on */\n")
    # convert from RGB888 to RGB565
    palette565 = [((palette[x] * 2**5 / 256) << 11) + ((palette[x+1] * 2**6 / 256) << 5) + ((palette[x+2] * 2**5/256)) for x in xrange(0, len(palette), 3)]
    f.write("const uint16_t palette[] PROGMEM = {" + ','.join(hex(x) for x in palette565) + "};\n\n")

def writePixel(pixelIndexColor, count):
    global rleBytes
    f.write(hex((pixelIndexColor << (8 - colorbits)) + count) + "," )
    rleBytes += 1

def writePixels(f):
    f.write("/* The first (MSB) " + str(colorbits) + "bits of each byte represent the color index, and the last " + str(8-colorbits) + " bits represent the number of pixels of that color */\n")
    f.write("const unsigned char pixels[] PROGMEM = {")
    (width, height) = image.size

    lastColors = []
    lastColor = -1
    lastColorCount = 0
    for y in xrange(height):
        for x in xrange(width):
            color = pix[x,y]
            if lastColorCount +1 < 2**(8 - colorbits) and (lastColorCount == 0 or color == lastColor):
                lastColorCount = lastColorCount + 1
            else:
                writePixel(lastColor, lastColorCount)
                lastColorCount = 1
            lastColor = color
        # To make the decoder easier, truncate at the end of a row
        # This adds about 5% to the total file size though
        writePixel(lastColor, lastColorCount)
        lastColorCount = 0
    f.write("};\n\n");


image = Image.open(sys.argv[1])
image = image.convert("P", palette=Image.ADAPTIVE, colors=2**colorbits)
#image.save("tmp.png")
pix = image.load()

f = open(sys.argv[2], 'w')
writeHeader(f)
writePalette(f, image)
writePixels(f)
writeFooter(f)
f.close()

print "Summary: Total Bytes: ", rleBytes

The result is that the 10 pictures take up 14kb. This doesn’t produce an optimal image, but it’s a decent starting point. This generates header files like:

/* A generated file for a run length encoded image of face1.png */

#ifndef PIXEL_VALUE_TO_COLOR_INDEX
#define PIXEL_VALUE_TO_COLOR_INDEX(x) ((x) >> 5)
#define PIXEL_VALUE_TO_RUN_LENGTH(x) ((x) & 0b11111)
#endif

namespace face1 {
/* The RGB565 value for color 0 is in palette[0] and so on */
const uint16_t palette[] PROGMEM = {0xef0b,0xf4a3,0x0,0xdb63,0xffff,0xa302,0xe54a,0x5181};

/* The first (MSB) 3 bits of each byte represent the color index, and the last 5 bits represent the number of pixels of that color */
const unsigned char pixels[] PROGMEM = {0x5f,0x56,0xe2...}
}

Which represents the image:

tmp

Which isn’t too bad for a first go. (This was generated by adding a image.save('tmp.png')" in the python code after converting the image).

Now we need an arduino sketch to pull this all in:

#include <SPI.h>
#include <Adafruit_GFX.h>
#include <TFT_ILI9163C.h>
#include "face1.h"
#include "face2.h"
#include "face3.h"
#include "face4.h"
#include "face5.h"
#include "face6.h"
#include "face7.h"
#include "face8.h"
#include "face9.h"
#include "face10.h"

const unsigned char *pixels[] = { face1::pixels, face2::pixels, face3::pixels, face4::pixels, face5::pixels, face6::pixels, face7::pixels, face8::pixels, face9::pixels, face10::pixels } ;
const unsigned char numpixels[] = { sizeof(face1::pixels), sizeof(face2::pixels), sizeof(face3::pixels), sizeof(face4::pixels), sizeof(face5::pixels), sizeof(face6::pixels), sizeof(face7::pixels), sizeof(face8::pixels), sizeof(face9::pixels), sizeof(face10::pixels) };
const uint16_t *palette[] = { face1::palette, face2::palette, face3::palette, face4::palette, face5::palette, face6::palette, face7::palette, face8::palette, face9::palette, face10::palette } ;


/*
 We are using 4 wire SPI here, so:
 MOSI:  11
 SCK:   13
 the rest of pin below:
 */
#define __CS 8
#define __DC 10
#define __RESET 9
TFT_ILI9163C tft = TFT_ILI9163C(__CS, __DC, __RESET);

const uint16_t width = 128;
const uint16_t height = 128;

void drawFace(int i);

void setup(void) {
  tft.begin();
}


void drawImage(const unsigned char *pixels, const unsigned char numPixels, const uint16_t *palette) {
  int16_t x = 0;
  int16_t y = 0;
  for( int i = 0; i < sizeof(face1::pixels); ++i) {
    unsigned char pixel = pgm_read_byte_near(pixels + i);
    /* The first (MSB) 3 of each byte represent the color index, and the last 5 bits represent the number of pixels of that color */
    unsigned char colorIndex = PIXEL_VALUE_TO_COLOR_INDEX(pixel);
    int16_t count = PIXEL_VALUE_TO_RUN_LENGTH(pixel);
    uint16_t color = pgm_read_word_near(palette + colorIndex);
    tft.drawFastHLine(x,y,count+1,color);
    x += count;
    if(x >= width) {
      x = 0;
      y+=1;
    }
  }
}

void loop() {
  for(int i = 1; i <= 10; ++i)
    drawFace(i);
}

void drawFace(int i) {
  drawImage(pixels[i-1], numpixels[i-1], palette[i-1]);
}

Pretty straight forward.

Using 4 bits per color uses up slightly too much space: 34,104 bytes (111%) of program storage.
Using 3 bits per color introduces slight artifcats, but presumably removable if we tweak the images, but it fits with plenty of room to space:

Sketch uses 20,796 bytes (67%) of program storage space. Maximum is 30,720 bytes.
Global variables use 134 bytes (6%) of dynamic memory, leaving 1,914 bytes for local variables. Maximum is 2,048 bytes.

The result is this:

Note that this is actually updating at maximum speed. There is room for improvement in the adafruit graphics code however, so I may have a go at improving this if I need faster updates. It might be nice to faster updates anyway in order to minimize battery use.

Android Phone App

I wanted to keep the Android app side of it as clean and simple as possible. There’s not much to document regarding the GUI side of things – it’s a pretty straightforward Android app.

The idea here is to connect via bluetooth to both the arduino and to a Healthband. It measures the heart rate and, more importantly, the very slightly changes in time between the heart pulses. By the knowing the users age, sex and weight, we can calculate the amount of stress that they are currently experiencing. We can then display that in the app, and send that information via Bluetooth Low Energy to our arduino heart badge.

Screenshot_2015-05-29-23-49-29

Screenshot_2015-05-29-23-51-15

It perhaps needs some tweaking – I’m not really that angry, honest!  Lots of polishing is still needed, but the barebones are starting to come together!

Woodworking – Make a bridge

I wanted to make a bridge out of wooden bricks, with no glue or connectors, held up by gravity only.

The bridge is an inverted catenary:

Which follows the equation:

y = a \cdot \cosh(x/a)

Now, turning this into IPython code:

L = 300.0 # 300mm wide bridge.  All lengths in millimetres
A = 100.0 # Scaling factor

def y(x):
    return -A*( math.cosh(x/A) - math.cosh((L/2)/A))

X = linspace(-L/2, L/2)
Y = [y(x) for x in X]

fig, ax = subplots(1,1)
fig.set_size_inches(14,7)
ax.plot(X,Y,'r')
xlabel('x')
ylabel('height')
title('Bridge curve')
show()

Which produces:

bridge

This is looking good.

Now, I want to build this out of wood, so we won’t actually have this curve, but a ‘blocky’ version of this.

class Block:
    '''A wooden block'''
    def __init__(self, start_x, end_x):
        self.start_x = start_x
        self.start_y = y(start_x)
        self.end_x = end_x
        self.end_y = y(end_x)
        self.length_x = end_x - start_x
        self.midpoint_x = start_x + self.length_x/2.0
        self.tangent = (self.end_y - self.start_y)  / self.length_x        
        self.midpoint_y = self.start_y + (self.end_y - self.start_y)/2.0
        self.coords = [] # filled in later
    def angle(self):
        return math.atan(self.tangent)

step = 60.0

def plotTangent(block):
    xstart = block.midpoint_x-block.length_x/2
    xend = block.midpoint_x+block.length_x/2
    ystart = block.midpoint_y + block.tangent * (xstart - block.midpoint_x)
    yend = block.midpoint_y + block.tangent * (xend - block.midpoint_x)
    ax.plot([block.start_x, block.end_x], [block.start_y, block.end_y], 'k-', lw=1)

blocks = []

for start_x in linspace(-L/2, L/2, L/step, endpoint=False):
    block = Block(start_x, start_x+step)
    blocks.append(block)
    plotTangent(block)
fig

Which produces:

bridge2
Now we have to decide how to cut the wood. To produce the arch, the angle of each piece of wood (the tangent lines), relative to the horizontal, is just: \tan(\theta) = tangent. So the angle between two adjacent pieces of wood is the difference between their angles.

The wood will be parallel with those tangent lines, so we can draw on the normal lines for those tangent lines, of the thickness of the wood:

def dy(x):
    '''return an approximation of dy/dx'''
    return (y(x+0.05) - y(x-0.05)) / 0.1

def normalAngle(x):
    if x == -L/2:
        return 0
    elif x == L/2:
        return math.pi
    return math.atan2(-1.0, dy(x));

WoodThickness = 50.0 # in millimeters

def getNormalCoords(x, ymidpoint, lineLength):
    tangent = dy(x)
    if tangent != 0:
        if x == -L/2 or x == L/2:
            normalSlope = 0
        else:
            normalSlope = -1 / dy(x)
        xlength = lineLength * math.cos(normalAngle(x))
        xstart = x-xlength/2
        xend = x+xlength/2
        ystart = ymidpoint + normalSlope * (xstart - x)
        yend = ymidpoint + normalSlope * (xend - x)
    else:
        xstart = x
        xend = x
        ystart = ymidpoint + lineLength/2
        yend = ymidpoint - lineLength/2
    return (xstart, ystart, xend, yend)

for block in blocks:
    (xstart, ystart, xend, yend) = getNormalCoords(block.midpoint_x, block.midpoint_y, WoodThickness)
    ax.plot([xstart, xend], [ystart, yend], 'b-', lw=1)
fig

bridge3

Now, the red curve is going to represent the line of force. My reasoning is that if there was any force normal to the red curve, then the string/bridge/weight would thus accelerate and change position.

So to minimize the lateral forces where the wood meets, we want the wood cuts to be normal to the red line.  We add that, then draw the whole block:

def getWoodCorners(block, x):
    ''' Return the top and bottom of the piece of wood whose cut passes through x,y(x), and whose middle is at xmid '''
    adjusted_thickness = WoodThickness / math.cos(normalAngle(x) - normalAngle(block.midpoint_x))
    return getNormalCoords(x, y(x), adjusted_thickness)

def drawPolygon(coords, linewidth):
    xcoords,ycoords = zip(*coords)
    xcoords = xcoords + (xcoords[0],)
    ycoords = ycoords + (ycoords[0],)
    ax.plot(xcoords, ycoords, 'k-', lw=linewidth)

#Draw all the blocks
for block in blocks:
    (xstart0, ystart0, xend0, yend0) = getWoodCorners(block, block.start_x) # Left side of block
    (xstart1, ystart1, xend1, yend1) = getWoodCorners(block, block.end_x) # Right side of block
    block.coords = [ (xstart0, ystart0), (xstart1, ystart1), (xend1, yend1), (xend0, yend0) ]
    drawPolygon(block.coords, 2) # Draw block

fig

bridge4

There’s a few things to note here:

  • The code appears to do some pretty redundant calculations when calculating the angles – but much of this is to make the signs of the angles correct.  It’s easier and nicer to let atan2 handle the signs for the quadrants, than to try to simplify the code and handle this ourselves.
  • The blocks aren’t actually exactly the same size at the points where they meet.  The differences are on the order of 1%, and not noticeable in these drawings however.

Now we need to draw and label these blocks in a way that makes it easiest to cut:

fig, ax = subplots(1,1)
fig.set_size_inches(17,7)
axis('off')

def rotatePolygon(polygon,theta):    
    """Rotates the given polygon which consists of corners represented as (x,y), around the ORIGIN, clock-wise, theta radians"""
    return [(x*math.cos(theta)-y*math.sin(theta) , x*math.sin(theta)+y*math.cos(theta)) for (x,y) in polygon]

def translatePolygon(polygon, xshift,yshift):
    """Shifts the polygon by the given amount"""
    return [ (x+xshift, y+yshift) for (x,y) in polygon]

sawThickness = 3.0 # add 3 mm gap between blocks for saw thickess
#Draw all the blocks
xshift = 10.0 # Start at 10mm to make it easier to cut by hand
topCutCoords = []
bottomCutCoords = []
for block in blocks:
    coords = translatePolygon(block.coords, -block.midpoint_x, -block.midpoint_y)
    coords = rotatePolygon(coords, -block.angle())
    coords = translatePolygon(coords, xshift - coords[0][0], -coords[3][1])
    xshift = coords[1][0] + sawThickness
    drawPolygon(coords,1)
    (topLeft, topRight, bottomRight, bottomLeft) = coords
    itopLeft = int(round(topLeft[0]))
    itopRight = int(round(topRight[0]))
    ibottomLeft = int(round(bottomLeft[0]))
    ibottomRight = int(round(bottomRight[0]))
    topCutCoords.append(itopLeft)
    topCutCoords.append(itopRight)
    bottomCutCoords.append(ibottomLeft)
    bottomCutCoords.append(ibottomRight)
    ax.text(topLeft[0], topLeft[1], itopLeft)
    ax.text(topRight[0], topRight[1], itopRight, horizontalalignment='right')
    ax.text(bottomLeft[0], bottomLeft[1], ibottomLeft, verticalalignment='top', horizontalalignment='center')
    ax.text(bottomRight[0], bottomRight[1], ibottomRight, verticalalignment='top', horizontalalignment='center')

print "Top coordinates:", topCutCoords
print "Bottom Coordinates:", bottomCutCoords

Which produces:

Top coordinates: [10, 141, 144, 228, 231, 307, 310, 394, 397, 528]
Bottom Coordinates: [43, 131, 156, 214, 247, 291, 324, 382, 407, 495]

bridge5

Now to cut this!

Flight aerodynamics simulator

I wanted to have a modern aerodynamics simulator, to test out my flight control hardware and software.

Apologies in advance for a very terse post.  I spent a lot of hours in a very short timespan to do this as a quick experiment.

So I used Unreal Engine 4 for the graphics engine, and built on top of that.

The main forces for a low speed aircraft that I want to model:

  1. The lift force due to the wing
  2. The drag force due to the wind
  3. Gravity
  4. Any horizontal forces from propellers in aircraft-style propulsion
  5. Any vertical forces from propellers in quadcopter-style propulsion

Lift force due to the wing

The equation for the lift force is:

L = C_l \cdot A \cdot .5 \cdot r \cdot V^2

I created a function for the lift coefficient, C_l, based on the angle, by calculating it theoretically.  To get proper results, I would need to actually measure this in a wind tunnel, but this is good enough for a first approximation:

micron lift coefficients

The horizontal axis is the angle of attack, in degrees.  When the aircraft is flying “straight”, the angle of attack is not usually 0, but around 5 to 10 degrees, thus still providing an upward force.

I repeat this for each force in turn.

3D Model

To visualize it nicely, I modelled the craft in blender, manually set the texture space, and painted the texture in gimp.  As you can tell from the texture space, there are several horrible problems with the geometry ‘loops’.  But it took a whole day to get the top and bottom looking decent, and it was close enough for my purposes.

 

I imported the model into Unreal Engine 4, and used a high-resolution render for the version in the top right, and used a low-resolution version for the game. screenshot2

textured_micron_screenshot

Next, here’s the underneath view.  You can see jagged edges in the model here, because at the time I didn’t understand how normal smoothing worked.  After a quick google, I selected those vertexes in blender and enabled normal smoothing on them and fixed that.

ue4_2(1)

and then finally, testing it out on the real thing:

Screenshot_20180626_232400

Screenshot_20180626_232350

Architecture

The architecture is a fairly standard Hardware-in-loop system.

The key modules are:

simulator

  • The Flight Controller which controls the craft.  This is not used if we are connected to external real hardware, such as the my controller.
  • The Communication Module to the real hardware, receiving information about the desired thrust of the engines and the Radio Control inputs from the user, and sending information about the current simulated position and speed.
  • The Physics Simulator which calculates the physical forces on the craft and applies them.
  • The User Interface which displays a lot of information about the craft as well as the internal controller.

The low level flight controller and network communication code is written in C++.  Much of the high level logic is written in a visual language called ‘Blueprint’.

User Interface

From the GUI User Interface, you can control:

  • The aerodynamic forces on the wings
  • The height above ground to air pressure curve
  • The wing span and aerofoil chord length
  • The moment of inertia
  • The thrust of the turbines
  • The placement of the turbines
  • The PID values for the internal controller

Result

It worked pretty well.

It uses Hardware-In-Loop to allow the real hardware to control this, including a RC-Transmitter.
In this video I am allowing the PID algorithms to control the roll, pitch, yaw and height, which I then add to in order to control it.

 

 

PID Tuning

I implemented an auto-tuner for the PID algorithm, which you can tune and trigger from the GUI.

 Hardware

I used two arduinos to control the system:

overview

And set up the Arduinos as so:  (Btw, I used the Fritzing software for this – it’s pretty cool).

neva arduino_bb(1)

 

And putting together the hardware for testing (sorry for the mess).  (I’m using QGroundControl to test it out).

And then mounting the hardware on a bit of wood as a base to keep it all together and make it more tidy:

20140822_202401.jpg

I will hopefully later make a post about the software controlling this.

ESC Motor Controller Delay

I was particularly worried about the delay that the controller introduces.

I modified the program and used a basic UFO style quadcopter, then added in a 50ms buffer, to test the reaction.

For reference, here’s a photo of the real thing that I work on:

The quadcopter is programmed to try to hover over the chair.

 

I also tested with different latencies:

 

50ms really is the bare minimum that you can get away with.

These are manually tuned PIDs.

I did also simulate this system in 1D in Matlab’s Simulink:

simulink_roll

A graph of the amplitude:

matlab estimated transfer function sin

And finally, various bode plots etc.  Just click on any for a larger image.  Again, apologies for the awful terseness.

Electronics for children

I’ve been trying to introduce my daughter, 6 years old, to electronics.  I started her on the Velleman MK100 Electronic Christmas Tree which went fairly smoothly.

Then we progressed to a Velleman MK127 Microbug which was a lot more difficult.

The problem with these kits is that beginners will make mistakes that will be unable to recover from.  Not only that, but they don’t explain how they work so you can’t debug them or understand what you are doing.

On top of that, the bots can be outright hackish – for example, requiring you to scrap the edge of the motors, and then solder the chasis of the motor directly to the board at an angle!   And this is even more difficult than it sounds, since the motor is highly magnetic and so requires a steady strong hand to keep the soldering iron still.

This was the result:

Notice the four wires from the motor. This was because we put the motors in the wrong way around, despite following the instructions. So I had to reverse the motors. This sort of repair is probably beyond a child. But worse:

Notice the yellow and white wires? This is another repair job after finding that the light sensitive resistor had broken the bottom of the board and was not making a proper connection. Pretty much impossible to detect this for a novice or if you didn’t have a voltmeter.

I’m not sure what the solution is. In the past I’ve used various non-soldering electronics kits, but they all have their own disadvantages.