Freesteel Blog » Machining

Monday, May 8th, 2017 at 6:49 pm - - Flightlogger, Hang-glide

I’m getting tired of having “learning experiences” when I really want more “fun experiences”. But in the meantime, here goes.

On 2017-04-29 I was at the hang-gliding competition on Camlo hill in Wales, where I didn’t take any thermals over the back — like the folks who knew what they were doing — and I landed back on the top with nil points and little satisfaction.

For this learning experience, I downloaded the tracklog of one of the guys who did make it and snipped out the first 12kms of flight, like so:

darrentrack12

Direction of flight north is to the right (start of flight bottom left near (0,0)), the coordinates are in metres.

Having convinced myself I can extract the SRTM3 terrain, this is the graph of the flight altitude from the hill for those 12km against the terrain altitude:

darrentrackalt

Or you can plot the actual height above ground by subtracting the two:

darrentrackaltg

This looks like he was circling with a pretty consistent height of 350m as the air mass was carried up and down over the terrain, in a “bubble” of rising air.

Except there is no way this could be a bubble, because even if the bubble extended all the way to the ground 350m below, it would have certainly been been expended in less than 6 minutes at a rate of 1 metre per second to support an efficiently flying glider.

After 25 minutes of existence, we can rule out that it has anything to do with any packets of warmed air which may have risen from the ground 12kms to the south.

“Ah,” the usual response goes, “the glider was flying in a chain of thermic bubbles along the track of the flight, each one rising up just in time to pick up from where the last left off.”

This explanation is bogus. I don’t seem to pick up thermic bubbles that easily along my track, while the pilot that stayed with the thermal reported a continuously existing atmospheric formation which strengthened and subsided, but was always there.

I also don’t buy the idea that this structure is somehow kicking off thermic bubbles on the ground 350m below at least five minutes downwind of its track in time to reach his altitude.

No, this must be a self-contained, self-sustaining convection structure that may have been initiated by a thermic bubble, but which has clearly morphed into something altogether different.

(more…)

Monday, April 24th, 2017 at 8:15 pm - - Flightlogger

I had a very agreeable flight on Saturday in Wales, marred by a crappy going down performance on Bradwell on Sunday, just to put me in my place as not having any talent.

Here’s the good place to be, about a mile up from the ground:
wavellan

It was a confused forecast, which meant a lot of people stayed away. Turns out it was because there was wave everywhere. I flew Tim’s U2 (mid-performance glider), but was put off by the fact of everyone pointing out that I was too high on the bar, and so wouldn’t be able to steer it so well. Hang-loops, unlike climbing harness loops, do not come with a buckle, so you can’t do anything about it.

There was a lot of lift straight off the hill, and then Carl said over the radio, “I think I’ve contacted wave, just to the right of the mound in front of the hill.”

There were lenticular (French for “slow”) wave clouds everywhere and no cumulus. I was wearing a skinny pair of gloves and had no base bar mitts (which I do have on my own glider).
(more…)

Wednesday, April 19th, 2017 at 11:41 pm - - Flightlogger

I’m not going to learn too much from myself, so I’ve got the raw data from Tim’s flight as measured by his Kobo/bluefly IGC logger.

timgps

As you can see, he ranged a good deal further than I did from Mam Tor, reaching an altitude of 1500m, and he reported scoring glide angles of 10 to 1 on his glider much-fancier-than-mine-of-which-I-am-now-jealous.

Let’s see if we can corroborate this from his data with the power of pandas and matplotlib.

pIGCtim = fd.LoadIGC("Tim-1492470000.igc")
pQ = pIGCtim.resample("1S").interpolate(method="time")  # fill uneven records

dstep = 30 # rolling windowed measurements in seconds
vario = (pQ.alt.diff(dstep)/dstep)  # m/s averaged over 30seconds
groundvelocity = numpy.sqrt(pQ.x.diff(dstep)**2 + pQ.y.diff(dstep)**2)/dstep

Here is vario:
timsvario

Here is ground velocity:
timsgroundvel

The glide slope (normally called the glide-angle, even though we don’t quote it as an angle) is the groundvelocity divided by the descent rate. When the descent rate is close to zero, you get stupid numbers, so we need to throw those out and only quote for the parts where there is reasonably direct gliding flight.

(more…)

Friday, April 14th, 2017 at 4:32 pm - - Machining

Building on the weighted servo motor (having spent the first 20 minutes of the day trying to find where I left my code) I’ve tried to code it to run on a sine wave as an oscillation.

It’s not easy, of course.
voltssinwave

How did I get this?

From the fixed voltage experiments I computed the function from the voltage and velocity to the acceleration for motor1 as follows:

cv0, cm0, cv1, cm1, cvm, cvc = 53.202, -1036.634, 56.027, -974.375, -7.497, -1521.571
acceleration1 = (velocity1 - (min(volts1-cv0, 0)*cm0 + max(volts1-cv1, 0)*cm1))*cvm + cvc

Therefore, given the velocity, if you want a particular acceleration, you need the following voltage:

tq = velocity1 - (acceleration1 - cvc)/cvm
if tq < 0:
    volts1 = tq/cm1 + cv1
else:
    volts1 = +tq/cm0 + cv0

The above was the result of trying to make the motor oscillate like a pendulum by setting the acceleration to a negative number times the displacement:

acceleration1 = -acc1fac*(position1 - position1c)

Of course, one source of the problem is knowing the velocity. I’ve bodged this by recording a couple of positions up to 0.2seconds in the past and measuring to them.

if currenttime - prevtime > 0.1:
    position1prev2 = position1prev
    position1prev = position1
    prev2time = prevtime
    prevtime = currenttime
velocity1 = (position1 - position1prev2)/(currenttime - prev2time)

Unfortunately, this produces lots of bumps, so I think I’ll have to do this properly by applying, say, an exponential filter to the position value which delays the position, and use the difference between that and the current position as a measure of velocity.

Like all such runtime velocity measurements (unlike the ones done in pandas that have been calculated by picking a position on either side of the current time, like so:

(df.shift(-10).position1 - df.shift(10).position1)/(df.time.shift(-10) - df.time.shift(10))

…it’s also going to be delayed by a certain amount and be invalid under accelerations.

(more…)

Thursday, April 6th, 2017 at 11:22 am - - Flightlogger

On my last sorry flight I was carrying a total of four GPS units. One of them is a pile of poo. It’s the bluefly one I programmed to read at the frequency of 10 times a second, so greedy was I for precision and data.

gpscompare

Why do they have that option if it’s a complete waste of my time? All that work learning how to remotely step up its baudrate so that it could transmit data at this higher rate.

Here’s how the error comparisons line up, including the minus 3second correction on the 6030 vario’s timestamps. The “u” value is the number of milliseconds from midnight according to the GPS satellite timestamp.

gnames = [ "gps6030", "gpsxcsoar", "gpsTopShelf", "gpsBlueFly" ]
print("samples, gpsA, gpsB, stddev")
for i in range(len(g)):
    for j in range(i+1, len(g)):
        qA, qB = g[i], g[j]
        dt = 3000 if i == 0 else 0 # observed timeslip on 6030 unit
        dx = qA.x - numpy.interp(qA.u-dt, qB.u, qB.x)
        dy = qA.y - numpy.interp(qA.u-dt, qB.u, qB.y)
        k = (dx**2 + dy**2)[tf0:tf1]
        print(len(k), gnames[i], gnames[j], math.sqrt(k.mean()))
samples gpsA gpsB stddev
400 gps6030 gpsxcsoar 6.086
400 gps6030 gpsTopShelf 7.914
400 gps6030 gpsBlueFly 17.040
550 gpsxcsoar gpsTopShelf 7.845
550 gpsxcsoar gpsBlueFly 19.095
4120 gpsTopShelf gpsBlueFly 17.451

Now I have to program the BlueFly to sample at 200ms like the spare top-shelf GPS does, and see if that fixes it.

Meantime, has is that the accelerometer doing with this not-crap GPS?
(more…)

Tuesday, April 4th, 2017 at 1:29 pm - - Machining

Suppose we apply a random series of fixed duty cycles to a servo motor, like so:

volts1

A 50% duty cycle means that the volts are applied half one way and half the other at about a frequency of 100kHz, so it’s equivalent to zero volts.

I can plot the position from the encoder at the same time, like so:

volts1pos

There’s about 1000 ticks per revolution of the wheel and it’s wired backwards so a positive voltage makes it spin backwards.

Amazingly, with the control loop written in Python on a beaglebone controlling the H-bridges, it makes 4500 samples per second, which is perfectly adequate and I have not had to resort to any fancy tricks with C++ or the PRUs.

With a pandas timeseries Series object pos, we can calculate the velocity over a window of +-100 samples (about 1/20th of a second) like so:

# I'd like a better way to shift and difference a timeseries index.
vel = (pos.shift(-100) - pos.shift(100)) / \
      (pos.index.to_series().shift(-100) - pos.index.to_series().shift(100))

And then there is the acceleration:

acc = (vel.shift(-100) - vel.shift(100)) / \
      (vel.index.to_series().shift(-100) - vel.index.to_series().shift(100))

We can plot them all up like so:

volts1velacc

Here’s what the first two seconds look like.
(more…)

Friday, March 31st, 2017 at 10:10 am - - Machining

No time for blogging, but here is a video of the servo motor doing some random lifts with the string.

And this is the graph of the data acquired (red is the step-changing voltage applied) including position, velocity and acceleration:
timetrackservo
Back later when I have time to write this up and do some further analysis.

Friday, March 24th, 2017 at 12:19 pm - - Flightlogger

We move on to the temperature sensor work, and the controversial concept that the temperature of the rising air in a thermal is hotter than the non-upwardly-mobile surrounding environmental atmosphere.

I say it’s controversial because the meaning of “the temperature” in relation to a mobile and turbulent airmass whose structure spans hundreds of metres in the vertical dimension and thousands of Pascals in relative pressure is undefined. Remember that the adiabatic temperature differential is about 0.7degrees per 100m change in altitude.

We do, however, have a single point temperature sensor on a mobile glider which is also progressing up and down the air column (depending on the wind currents and skill of the pilot). The location of the glider with the single temperature sensor is imperfectly known (due to bad GPS (see previous post) and inexplicable barometric behavior (see next post)), and the sensor itself has a thermal mass which means its readings have a delay half-life of about 9 seconds in flowing air (see this other post).

I have taken the precaution of jettisoning my slow, accurate and low resolution dallas temperature sensor for two humidity/temperature sensor combos and an infrared thermometer that has quite a good ambient temperature sensor within its metal can.

These sensors tracked one another pretty well, except for one small problem when I spotted a series of spikes in one and then the other humidity/temperature sensor.

temperaturespikes

What is going on?

(more…)

Monday, March 20th, 2017 at 1:45 pm - - Flightlogger 2 Comments »

Last week I finally had my first flight of the year with my newly build flight data logger. I can’t believe the number of issues it’s already thrown up.

At least I may be making quick enough progress to get past the issues (rather than being swamped by them) using this exceptionally efficient Jupyter/Pandas technology.

For example, my code for parsing and loading the IGC file is 15 lines long.

The code for loading in my flight logger data into a timeseries is just as brief, if you consider each data type individually (there are more than 13 of them from humidity sensors to an orientation meter).

The GPS time series from my flight logger (at 10Hz) can be crudely converted it to XYs in metres, like so:

# pQ is the GPS position pandas.DataFrame
earthrad = 6378137
lng0, lat0 = pQ.iloc[0].lng, pq.iloc[0].lat
nyfac = 2*math.pi*earthrad/360
exfac = nyfac*math.cos(math.radians(lat0))
pQ["x"] = (pQ.lng - lng0)*exfac
pQ["y"] = (pQ.lat - lat0)*nyfac
plt.plot(pQ.x, pQ.y)

gpstrack1

Note the suspicious sharp turn near (-1000, -400). Here’s another sharp turn somewhere else in the sequence covering a 1 minute 5 second period using time slicing technology:

t0, t1 = Timestamp("2017-03-09 15:42:55"), Timestamp("2017-03-09 15:44:00")
q = fd.pQ[t0:t1]
plt.plot(q.x, q.y)

gpstrack2

The dot is at the start point, time=t0.

(more…)

Friday, March 3rd, 2017 at 11:04 am - - Flightlogger, Hang-glide, Uncategorized

To be clear, I haven’t got mathematical proofs here (I don’t have the time), but the experimental evidence is quick to get.

Take the differential barometer sensor (used to measure airspeed) of the hang-glider flight logger. The Arduino code which updates the reading every 200ms looks like this:

long lastpx4timestamp; 
void Flylogger::FetchPX4pitot()
{
    long mstamp = millis(); 
    if (mstamp >= lastpx4timestamp + 200) {
        px4pitot->readpitot(); 
        sdlogger->logpitot(px4timestampset, px4pitot-rawpressure, px4pitot->rawtemp); 
        lastpx4timestamp = mstamp; 
    }
}

Why did I choose 200 milliseconds? It sounded like a good number to read it at. This is a quick way to program it to be a regular reading.

A better way is to actually synchronize it with the clock divided rather than simply add 200ms to the next time, like so:

int mstampdivider = 20; 
int prevmstampdivided = 0; 
void loop()
{
    long mstampdivided = millis()/mstampdivider; 
    if (mstampdivided != prevmstampdivided) {
        prevmstampdivided = mstampdivided; 
        P(micros());  P(" ");  P(singlereading());  P("\n"); 
    }
}

Now that code reads at 20ms rather than 200ms, but it prints a load of output which I can cut and paste into a file and read into pandas, like so:

rows = [ (int(s[0]), int(s[1]))  for s in (ln.split()  for ln in open("../logfiles/dmprapidtest.txt").readlines())  if len(s) == 2]
k = pandas.DataFrame.from_records(rows, columns=["t", "d"])

And then we can plot the autocorrelation (the covariance) with itself shifted in time, like so:

d = k.d   # just the measurement Series
dm = d.mean()
ss = [((d - dm)*(d.shift(i) - dm)).mean()  for i in range(400)]

autocov1

Let’s zoom in on the first 50 covariances:
(more…)