Biscuit B. Collie, 2002?-2017

Old age finally caught up with my Border Collie, and I had to put her down around 1:00 am, Saturday, September 23. She’s been a wonderful companion since I got her in July, 2005, through an informal rescue. I have a bunch of great stories about her life, so I plan on posting them here over the next few days or weeks.

When I was looking for a dog, I saw an ad for Border Collie puppies and an adult dog from someone down in Monument, CO who breeds and rescues Border Collies. I called the number and made arrangements to meet her at an off-leash dog park south of Denver, about halfway between the two of us. Bixby was a two or three year old female, black with white on the front legs, back paws, snout and chest. She greeted me with a low crawl, and was happy to accept petting. There was another medium-sized dog there, and Bixby started to stalk it, lowering her head to shoulder level and staring straight at it. The breeder commanded “leave it,” and Bixby complied.

That was Wednesday; I decided to get her, and made arrangements with the breeder. We met again at the park on Friday (July 15th, 2005) and I picked Bixby up.

When I told Mom about my new dog, she said that Bixby was hard to say, and she’d call the dog Biscuit. I immediately thought, “That’s a great name for a dog!”  So after that, she was Biscuit.

Next post: her first vet visit.


Posted in Biscuit | Tagged | Leave a comment

Unit Testing

Here’s the Keynote file from my talk at the Denver Swift Developers meetup: Unit Testing in Swift Updated

Posted in Uncategorized | Leave a comment


Sunday, April 30th: On our morning walk, we ran into a coyote! After the alarm went off at 6:00 am, I got dressed for our morning walk. We crossed the street to the park, and then crossed the creek and traveled north along it. We hadn’t gotten far before we heard the “ow, ow, owoo” coyote bark. I looked across the creek and saw the coyote. Biscuit saw it, too. Her hackles went up, she raised her tail (usually she holds it down low), and she started jumping around and barking.

I walked a little ways toward the coyote, and Biscuit advanced, barking all the while. I decided it wouldn’t be real smart to get close, so we headed back north along the creek.

Apparently that was the coyote’s way home, because he crossed the creek and tried to circle around us, but I-25 kept him from getting far enough away from us. Biscuit saw him and charged! I tried to call her back, but she kept going, barking and growling. I was a little concerned—coyotes are bigger than Biscuit, and kill things for a living. But the coyote ran off, while Biscuit kept watching him. No predators welcome when there’s a Border Collie on duty!

Posted in Biscuit | Leave a comment

Making a Line Chart for iOS

I wanted to be able to draw a chart for a Swift playground without having to add the full Core Plot or Charts frameworks, and found a tutorial here on Medium. It was not quite two years old (from December 2016), so I had to update it for the latest version of Swift, and then had a bit more work to do to get it functioning in a playground.

Travis Stenerson wrote the code for the class itself, LineChart, and the view controller, MyViewController. I tweaked LineChart a bit so it could handle data with negative values and made a couple of other very minor adjustments. Here’s the complete code of the playground page:

//: A tutorial from 
import UIKit
import PlaygroundSupport
extension String {
    func size(withSystemFontSize pointSize: CGFloat) -> CGSize {
        return (self as NSString).size(withAttributes: [NSAttributedString.Key.font: UIFont.systemFont(ofSize: pointSize)])

extension CGPoint {
    func adding(x: CGFloat) -> CGPoint { return CGPoint(x: self.x + x, y: self.y) }
    func adding(y: CGFloat) -> CGPoint { return CGPoint(x: self.x, y: self.y + y) }

class LineChart: UIView {
    let lineLayer = CAShapeLayer()
    let circlesLayer = CAShapeLayer()

    var chartTransform: CGAffineTransform?

    @IBInspectable var lineColor: UIColor = {
        didSet {
            lineLayer.strokeColor = lineColor.cgColor

    @IBInspectable var lineWidth: CGFloat = 1

    @IBInspectable var circleColor: UIColor = {
        didSet {
            circlesLayer.fillColor = circleColor.cgColor
    @IBInspectable var circleSizeMultiplier: CGFloat = 3
    @IBInspectable var axisColor: UIColor = UIColor.white
    @IBInspectable var showInnerLines: Bool = true
    @IBInspectable var labelFontSize: CGFloat = 10

    var axisLineWidth: CGFloat = 1
    public var deltaX: CGFloat = 10
    public var deltaY: CGFloat = 10
    var xMax: CGFloat = 100
    var yMax: CGFloat = 100
    var xMin: CGFloat = 0
    var yMin: CGFloat = 0
    var data: [CGPoint]?

    override init(frame: CGRect) {
        super.init(frame: frame)

    required init?(coder aDecoder: NSCoder) {
        super.init(coder: aDecoder)
    func combinedInit() {
        lineLayer.fillColor = UIColor.clear.cgColor
        lineLayer.strokeColor = lineColor.cgColor

        circlesLayer.fillColor = circleColor.cgColor

        layer.borderWidth = 1
        layer.borderColor = axisColor.cgColor

    override func layoutSubviews() {
        lineLayer.frame = bounds
        circlesLayer.frame = bounds

        if let d = data{
            setTransform(minX: xMin, maxX: xMax, minY: yMin, maxY: yMax)
    override func draw(_ rect: CGRect) {
        // draw rect comes with a drawing context, so let's grab it.
        // Also, if there is not yet a chart transform, we will bail on performing any other drawing.
        // I like guard statements for this because it's kind of like a bouncer to a bar.
        // If you don't have your transform yet, you can't enter drawAxes.
        guard let context = UIGraphicsGetCurrentContext(), let t = chartTransform else { return }
        drawAxes(in: context, usingTransform: t)

    func drawAxes(in context: CGContext, usingTransform t: CGAffineTransform) {

        // make two paths, one for thick lines, one for thin
        let thickerLines = CGMutablePath()
        let thinnerLines = CGMutablePath()

        // the two line chart axes
        let xAxisPoints = [CGPoint(x: xMin, y: 0), CGPoint(x: xMax, y: 0)]
        let yAxisPoints = [CGPoint(x: 0, y: yMin), CGPoint(x: 0, y: yMax)]
        // add each to thicker lines but apply our transform too.
        thickerLines.addLines(between: xAxisPoints, transform: t)
        thickerLines.addLines(between: yAxisPoints, transform: t)

        // next we go from xMin to xMax by deltaX using stride
        for x in stride(from: xMin, through: xMax, by: deltaX) {

            // tick points are the points for the ticks on each axis
            // we check showInnerLines first to see if we are drawing small ticks or full lines
            // tip for new guys: `let a = someBool ? b : c`  is called a ternary operator
            // in english it means "let a = b if somebool is true, or c if it is false."

            print("drawing x = \(x) gridline")

            let tickPoints = showInnerLines ?
                [CGPoint(x: x, y: yMin).applying(t), CGPoint(x: x, y: yMax).applying(t)] :
                [CGPoint(x: x, y: 0).applying(t), CGPoint(x: x, y: 0).applying(t).adding(y: -5)]

            thinnerLines.addLines(between: tickPoints)

            if x != 0 {  // draw the tick label (it is too busy if you draw it at the origin for both x & y
                let label = "\(Int(x))" as NSString // Int to get rid
                let labelSize = "\(Int(x))".size(withSystemFontSize: labelFontSize)
                let labelDrawPoint = CGPoint(x: x, y: 0).applying(t)
                    .adding(x: -labelSize.width/2)
                    .adding(y: 1)

                label.draw(at: labelDrawPoint,
                    [NSAttributedString.Key.font: UIFont.systemFont(ofSize: labelFontSize),
                     NSAttributedString.Key.foregroundColor: axisColor])
        // repeat for y
        for y in stride(from: yMin, through: yMax, by: deltaY) {

            let tickPoints = showInnerLines ?
                [CGPoint(x: xMin, y: y).applying(t), CGPoint(x: xMax, y: y).applying(t)] :
                [CGPoint(x: 0, y: y).applying(t), CGPoint(x: 0, y: y).applying(t).adding(x: 5)]

            print("drawing y = \(y) gridline")
            thinnerLines.addLines(between: tickPoints)
            if y != 0 {
                let label = "\(Int(y))" as NSString
                let labelSize = "\(Int(y))".size(withSystemFontSize: labelFontSize)
                let labelDrawPoint = CGPoint(x: 0, y: y).applying(t)
                    .adding(x: -labelSize.width - 1)
                    .adding(y: -labelSize.height/2)

                label.draw(at: labelDrawPoint,
                    [NSAttributedString.Key.font: UIFont.systemFont(ofSize: labelFontSize),
                     NSAttributedString.Key.foregroundColor: axisColor])
        // finally set stroke color & line width then stroke thick lines, repeat for thin

        // whenever you change a graphics context you should save it prior and restore it after
        // if we were using a context other than draw(rect) we would have to also end the graphics context
    func plot(_ points: [CGPoint]) {
        lineLayer.path = nil
        circlesLayer.path = nil
        data = nil

        guard !points.isEmpty else { return } = points

        if self.chartTransform == nil {
            setAxisRange(forPoints: points)

        let linePath = CGMutablePath()
        linePath.addLines(between: points, transform: chartTransform!)
        lineLayer.path = linePath

        let showPoints = true

        if showPoints {
            circlesLayer.path = circles(atPoints: points, withTransform: chartTransform!)


    func circles(atPoints points: [CGPoint], withTransform t: CGAffineTransform) -> CGPath {
        let path = CGMutablePath()
        let radius = lineLayer.lineWidth * circleSizeMultiplier/2
        for i in points {
            let p = i.applying(t)
            let rect = CGRect(x: p.x - radius, y: p.y - radius, width: radius * 2, height: radius * 2)
            path.addEllipse(in: rect)


        return path

    func setAxisRange(forPoints points: [CGPoint]) {
        guard !points.isEmpty else { return }


        let xs = { $0.x }
        let ys = { $0.y }

        xMax = ceil(xs.max()! / deltaX) * deltaX
        yMax = ceil(ys.max()! / deltaY) * deltaY
        xMin = floor(xs.min()! / deltaX) * deltaX
        yMin = floor(ys.min()! / deltaY) * deltaY
        print("xMin = \(xMin), xMax = \(xMax), yMin = \(yMin), yMax = \(yMax)")
        setTransform(minX: xMin, maxX: xMax, minY: yMin, maxY: yMax)

    func setAxisRange(xMin: CGFloat, xMax: CGFloat, yMin: CGFloat, yMax: CGFloat) {
        self.xMin = xMin
        self.xMax = xMax
        self.yMin = yMin
        self.yMax = yMax

        setTransform(minX: xMin, maxX: xMax, minY: yMin, maxY: yMax)

    func setTransform(minX: CGFloat, maxX: CGFloat, minY: CGFloat, maxY: CGFloat) {
        let xLabelSize = "\(Int(maxX))".size(withSystemFontSize: labelFontSize)
        let yLabelSize = "\(Int(maxY))".size(withSystemFontSize: labelFontSize)
        let xOffset = xLabelSize.height + 2
        let yOffset = yLabelSize.width + 2

        let xScale = (bounds.width - yOffset - xLabelSize.width/2 - 2)/(maxX - minX)
        let yScale = (bounds.height - xOffset - yLabelSize.height/2 - 2)/(maxY - minY)

        chartTransform = CGAffineTransform(a: xScale, b: 0, c: 0, d: -yScale, tx: yOffset - minX*xScale, ty: bounds.height - xOffset + minY*yScale)


class MyViewController : UIViewController {
    var lineChart: LineChart! // needs to be a forced optional so the class won't need an init() method, much like in an app

    override func loadView() { // note the lack of viewDidLoad()
        let view = UIView()
        view.backgroundColor = .blue
        let frame = CGRect(x: 10, y: 10, width: 355, height: 400)
        let lineChart = LineChart(frame: frame)

        self.view = view

        let f: (CGFloat) -> CGPoint = {
            let noiseY = (CGFloat(arc4random_uniform(2)) * 2 - 1) * CGFloat(arc4random_uniform(4))
            let noiseX = (CGFloat(arc4random_uniform(2)) * 2 - 1) * CGFloat(arc4random_uniform(4))
            let b: CGFloat = -25
            let y = 2 * $0 + b + noiseY
            return CGPoint(x: $0 + noiseX, y: y)

        let xs = [Int](1..<20)

        let points ={f(CGFloat($0 * 10))})

        lineChart.deltaX = 20
        lineChart.deltaY = 40


// Present the view controller in the Live View window
PlaygroundPage.current.liveView = MyViewController()



Posted in Uncategorized | Leave a comment

Moroccan Meskouta Orange Cake

I made this for a cooking meetup’s Moroccan-themed dinner. I had a couple of requests for the recipe, so here are the details.


  • ½ cup/120 ml freshly squeezed orange juice
  • Zest of 1 or 2 oranges (all the oranges you need for the juice)
  • 4 eggs
  • 1½ cups/300g sugar
  • ½ cup/120ml vegetable oil
  • 2 cups/250g flour
  • 4 tsp baking powder
  • ½ tsp salt
  • 1 tsp vanilla


Preheat your oven to 350°, then grease and flour a bundt pan (ideal) or a
loaf pan (will work, just make sure it’s big enough). Zest and then
juice the oranges.

Beat the eggs and sugar thoroughly, then the oil. Now mix in the flour,
baking powder, salt, and then the orange juice. Keep mixing until smooth.
Lastly, mix in the zest and the vanilla.

Pour the batter into the pan and bake for 45 minutes. The original recipe
(here) called for 40 minutes, but it’s taken me 50 minutes for the
loaf pan and 45 for the bundt pan. When it’s done baking, remove it
from the oven and let cool for 7–10 minutes. Remove the cake from the
pan to let it finish cooling.

Posted in Cooking & Baking, Uncategorized | Leave a comment

First Vet Visit: Puppies!

About a week after getting Biscuit, I figured that things were working out well, and she was a keeper, so I’d better take her to the vet and get her spayed before something happened. The vet, Chris Schwarz, spent a lot of time palpating her belly, and asked the tech to go get the ultrasound machine. Yep, she was pregnant already. The ultrasound confirmed his suspicions when he found a heart, too small and in the wrong spot for Biscuit’s. Dr. Schwarz estimated that the puppies were due in a week or two, so I thought I had a little time to get ready.

I called the breeder that I’d gotten her from, and explained situation. She volunteered to take Bixby/Biscuit back until the puppies were weaned (or perhaps transferred to another female who also just had a litter of pups). I was relieved that she was willing; I had planned on suggesting that, and was a bit worried that she might not be willing to help out. (I have no experience with puppies, so figured it would be better to have her handle the situation). We agreed that we’d meet on Friday to transfer Biscuit back to her.

She also told me to watch out for nesting. A few days before a female gives birth, she will gather some objects together to make a nest for the puppies. I told her I hadn’t seen this, and would keep an eye out for it.

The next day (Thursday, July 28th, 2005), I saw her squat several times as if she were trying to take a dump when I was out walking her. I called the vet’s office about that, and they said they’d have a vet tech call back. When they hadn’t called back in an hour, I called again, and got an apology and a promise to have a tech call me back.

Around six o’clock, I saw Biscuit squat in the living room. I thought “Oh, no you don’t,” (she’d had a couple of accidents in the apartment) and grabbed her by the collar to drag her outside. When I got her outside, I heard a splat on the sidewalk. Oh, no, that’s not her water breaking is it? Yep!

I looked back at her to see a puppy hanging halfway out, squirming and waving its little legs. I pulled the cell phone out of my pocket, and hit redial. “Biscuit’s having her puppies, what do I do?” The receptionist laughed and said she’d turn me over to a vet tech for advice.

By this time the puppy had dropped on to the grass. The tech told me to get Biscuit and the puppy back inside and into the kennel. I picked up the pup, who had rolled into the dirt, and carried it into the apartment as Biscuit walked alongside, eyes glued to her puppy, with a look that said, “That’s my puppy, what a beautiful puppy!” Pup did not like this flying through the air thing (isn’t there supposed to be a mama dog here for me?), and squeaked, much louder than I expected. The tech told me not to worry about picking the puppy up, it’s a myth that handling will make the mother reject the pup. That’s good, I told her, because I’ve already picked the puppy up.

I put the puppy inside the kennel, but splat, Biscuit’s water broke again, and she ducked into the closet. There was her nest!

She had gathered some throw pillows into a circle, just big enough for her to curl up in. She squatted, and pup #2 came out before I could get her into the kennel. At this point, the tech was telling me that it’s important to get her into the kennel, because the goop that comes out with the puppy makes a stain that won’t come out of the carpet. Oops. As we were talking, Biscuit ate the amniotic sac off the pup (she must have done that for #1 while I was working the phone), and started licking the puppy. I picked it up and put it in the kennel, and Biscuit followed.

The tech told me the problems that might happen, what to keep an eye on, and when to call for help or drive Biscuit and the pups over to Northglenn Emergency Animal Hospital. As she was telling me this, I saw that Biscuit was lying on her side, and the pups had scooched themselves over and started nursing: mama and puppies knew just what to do! With the advice passed on, the vet tech signed off, and I settled down to watch Biscuit and her pups.

Next call was to Mom. Hey, you’re a great-grandma! She was surprised and happy with the news about the pups. After that I called the breeder to say that the pups had shown up early. She also had some advice on what to keep an eye on, and suggested that I call her before taking Biscuit and the pups to the vet’s or the hospital. We agreed that I’d drive over the next morning with the dogs.

I checked the puppies, and they looked like they were both female. The vet later confirmed this. I was a little worried that I might have gotten it wrong, with the way they squirmed and wiggled when I turned them over to check. Puppies do not like being upside-down.

The evening proceeded without incident. I just watched Biscuit and her pups, petting them occasionally, and enjoying their soft, smooth coats. Biscuit tried to pick them up in her mouth a couple of times to move them, but they always squeaked in protest and she stopped. I moved the pups for her, which they didn’t seem to mind. I brought over her water bowl, since she didn’t want to leave the puppies.

The next morning I fed Biscuit by hand, in the kennel, since she still didn’t want to leave them. I took her and the puppies over to the vet’s office by 8:00 a.m., when they opened. Not having anything else to put the puppies in, I used a t-shirt to line one of Biscuit’s metal food bowls, and put them in head to tail in a circle, a sort of puppy yin-yang symbol.

The vet checked Biscuit and the puppies and reported that they were all doing well. Checking the puppies involved listening to their hearts and lungs with the stethoscope. It turns out they were no wider than the stethoscope itself, so he just draped them over it, with their hind legs dangling. So sweet!

After that, I drove down to Monument. It was a hot summer day, so I turned the air conditioning all the way up. I figured that the puppies might not be able to cope with getting too hot, but if they got too cold they’d just snuggle up to Biscuit to keep warm.

Posted in Biscuit | Tagged | 2 Comments

Determining Photogate Timing Uncertainty

Determining Photogate Timing Uncertainty

In order to do an uncertainty analysis of Front Range’s conservation of momentum experiments, we need to determine the uncertainty of the times reported by the photogates. The key to experiments of this sort is to collect redundant measurements, use least squares techniques to come up with a model of the expected results and then take statistics on the differences between the measurements and the fitted model.

The results, developed below, show that the photogate timing uncertainty is roughly 1.3 ms, and that the resistance force is proportional to the glider’s velocity.

The experiment to determine the photogate timing uncertainties uses four photogates spaced evenly along an air track, connected to a laptop computer via a USB data acquisition system. A glider slides along the air track, and a fin on the glider interrupts a beam of infrared light across the photogate. The data acquisition system records the times that the beam transitions between the broken and clear states. The air track provides the glider with a low-friction environment, but the resistance is not zero, so we will model the glider as having a constant deceleration during each pass along the track.

In the experiment, the glider is given an initial velocity, and allowed to bounce off the ends of the track, reversing direction each time. The data acquisition system collects data for 60 seconds, allowing on the order of 10 round trips along the track for each run.

load run5
n = length(t_meas)
n_trips = fix(n/16)
n = n_trips*16
t_meas = t_meas(1:n);
l_fin = 0.07; % 7 cm fin length
t_est = (t_meas(1:2:end)+t_meas(2:2:end))/2;
v_est = l_fin./(t_meas(2:2:end) - t_meas(1:2:end));
plot(t_est, v_est, '+')
xlabel('Time (s)')
ylabel('Velocity (m/s)')
title('Run 5 velocity history from time measurements')

n =


n_trips =


n =



Estimating Photogate Spacing

The glider goes through photogates 1, 2, 3 and 4, and then bounces off the far end of the track and goes through photogates 4, 3, 2 and 1. It then bounces off the near end of the track and repeats the process. I should be able to estimate the spacing of the photogates from the measurements, but I need to be sure to use the correct times to work out the three distances.

% Spacing between photogates 1 and 2, travelling right
d1_right = (v_est(2:8:end)+v_est(1:8:end))/2.*(t_est(2:8:end)-t_est(1:8:end));
% Spacing between photogates 1 and 2, travelling left
d1_left = (v_est(7:8:end)+v_est(8:8:end))/2.*(t_est(8:8:end)-t_est(7:8:end));
d1_data = [d1_right; d1_left];
d1 = mean(d1_data)
d1_unc = std(d1_data)

% Spacing between photogates 2 and 3, traveling right
d2_right = (v_est(3:8:end)+v_est(2:8:end))/2.*(t_est(3:8:end)-t_est(2:8:end));
% Spacing between photogates 2 and 3, traveling left
d2_left = (v_est(6:8:end)+v_est(7:8:end))/2.*(t_est(7:8:end)-t_est(6:8:end));
d2_data = [d2_right; d2_left];
d2 = mean(d2_data)
d2_unc = std(d2_data)

% Spacing between photogates 3 and 4, traveling right
d3_right = (v_est(4:8:end)+v_est(3:8:end))/2.*(t_est(4:8:end)-t_est(3:8:end));
% Spacing between photogates 3 and 4, traveling left
d3_left = (v_est(5:8:end)+v_est(6:8:end))/2.*(t_est(6:8:end)-t_est(5:8:end));
d3_data = [d3_right; d3_left];
d3 = mean(d3_data)
d3_unc = std(d3_data)

d1 =


d1_unc =


d2 =


d2_unc =


d3 =


d3_unc =


Analyzing the First Rightward Pass

This will take a few steps…

Define the Locations of Block and Clear Events

t_pass = t_meas(1:8) - mean(t_meas(1:8))
x_photogate = d2*[-1 -1 +1 +1]/2 + [-d1 0 0 d3]
A = [1 0 0 0; 1 0 0 0; 0 1 0 0; 0 1 0 0; 0 0 1 0; 0 0 1 0; 0 0 0 1; 0 0 0 1]';
x_event = x_photogate*A + [-1 +1 -1 +1 -1 +1 -1 +1]*l_fin/2

t_pass =


x_photogate =

   -0.3588   -0.1190    0.1190    0.3612

x_event =

   -0.3938   -0.3238   -0.1540   -0.0840    0.0840    0.1540    0.3262    0.3962

Basic Linear Fit

% Find the midpoints of the occluded intervals
t_mid = (t_pass(1:2:end) + t_pass(2:2:end))/2
% Now find estimated velocities
v_est = l_fin./(t_pass(2:2:end) - t_pass(1:2:end))
% Naive fit
[v_fit, v_unc, a, v_0, a_unc, v_0_unc] = uncertainFit(t_mid, v_est);
plot(t_mid, v_est, '+-', t_mid, v_fit, 'x:')
xlabel('Time (s)')
ylabel('Velocity (m/s)')
title('Naive Velocity Fit')

t_mid =


v_est =



Find an initial condition to minimize the squared error

I’ll use the fitted acceleration and velocity to predict the event times based on a range of initial positions, and then choose the initial position that minimizes the sum of the squares of the differences between the estimated and measured event times.

t_pred = t_pass; % create array of correct size, the values will be replaced
x_init = (-5:5)*d2/240;
sse_val = x_init*0; % again, creating a correctly sized array
for j = 1:length(x_init)
  for i = 1:8
    t_pred(i) = 2*(x_event(i)-x_init(j))/(v_0 + sqrt(v_0^2 + 2*a*(x_event(i)-x_init(j))));
  sse_val(j) = sum((t_pred-t_pass).^2);
plot(x_init, sse_val, 'x-')
title('First attempt at finding best initial condition')
xlabel('Initial position (m)')
ylabel('RSS error')

p = polyfit(x_init, sse_val, 2)
x_best = -p(2)/p(1)/2
best_sse_est = polyval(p, x_best)

p =

   15.1811   -0.0542    0.0000

x_best =


ans =


best_sse_est =



Now find an improved estimate of the initial position

for i = 1:8
  t_pred(i) = 2*(x_event(i)-x_best)/(v_0 + sqrt(v_0^2 + 2*a*(x_event(i)-x_best)));

best_sse_calc = sum((t_pred-t_pass).^2)
plot(x_init, sse_val-polyval(p, x_init), 'x-', x_best, best_sse_est, 'o')

x_init = (-5:5)*d2/960+x_best;
sse_val = x_init*0;
for j = 1:length(x_init)
  for i = 1:8
    t_pred(i) = 2*(x_event(i)-x_init(j))/(v_0 + sqrt(v_0^2 + 2*a*(x_event(i)-x_init(j))));
  sse_val(j) = sum((t_pred-t_pass).^2);
plot(x_init, sse_val, 'x-')
title('Refined attempt at finding best initial condition')
xlabel('Initial position (m)')
ylabel('RSS error')

best_sse_calc =



Use a minimization routine to find best x_0, v_0 and a

The calculations above suggest that finding values of initial position, initial velocity and acceleration that minimize the sum of the squares between modeled and measured times is feasible. Finding an analytical expression for the minimizing values would be challenging, but a minimization algorithm can find these values quickly.

v = version; % for checking if it's running in Matlab or Octave

Loop through all passes along track

x_best_list = [];
v_best_list = [];
a_best_list = [];
e_best_list = [];
sse_list = [];
x_fwd = x_event;
x_back = -x_event(8:-1:1);
for i = 1:n_trips
    i_offset = (i-1)*16;
    % forward trip
    t_pass = t_meas(i_offset+1:i_offset+8);
    t_pass = t_pass-mean(t_pass);
    t_mid = (t_pass(1:2:end) + t_pass(2:2:end))/2;
    v_est = l_fin./(t_pass(2:2:end)-t_pass(1:2:end));
    [v_fit, v_unc, a, v_0, a_unc, v_0_unc] = uncertainFit(t_mid, v_est);
    if v(1) == '7' % it's Matlab
        [X_best, sse_best] = fminsearch(@(x) sse(x, x_fwd, t_pass), [0 v_0 a]);
        % Nelder-Mead unconstrained nonlinear minimization
    elseif v(1) == '3' % it's Octave
        [X_best, sse_best, info, iter, nf, lambda] = sqp ([0 v_0 a]', @(x) sse(x, x_fwd, t_pass));
        % sequential quadratic programming algorithm
        warning('Unknown environment')
    x_best_list = [x_best_list X_best(1)];
    v_best_list = [v_best_list X_best(2)];
    a_best_list = [a_best_list X_best(3)];
    e_best_list = [e_best_list sse_best];
    % backward trip
    t_pass = t_meas(i_offset+9:i_offset+16);
    t_pass = t_pass-mean(t_pass);
    t_mid = (t_pass(1:2:end) + t_pass(2:2:end))/2;
    v_est = l_fin./(t_pass(2:2:end)-t_pass(1:2:end));
    [v_fit, v_unc, a, v_0, a_unc, v_0_unc] = uncertainFit(t_mid, v_est);
    if v(1) == '7' % it's Matlab
        [X_best, sse_best] = fminsearch(@(x) sse(x, x_back, t_pass), [0 v_0 a]);
        % Matlab uses Nelder-Mead unconstrained nonlinear minimization
    elseif v(1) == '3' % it's Octave
        [X_best, sse_best, info, iter, nf, lambda] = sqp ([0 v_0 a]', @(x) sse(x, x_back, t_pass));
        % Octave uses a sequential quadratic programming algorithm
        warning('Unknown environment')
    x_best_list = [x_best_list X_best(1)];
    v_best_list = [v_best_list X_best(2)];
    a_best_list = [a_best_list X_best(3)];
    e_best_list = [e_best_list sse_best];
end % of the n_trips/2 round trips along the track

The Results

The main goal was to determine the timing uncertainty of the photogate measurements, but the multiple passes at decreasing speeds allow us to make an estimate of the resistance force as a function of glider velocity. Before doing those calculations, let’s look at the results.

The optimization algorithm provided the resulting sum of the squares of the timing discrepancies, so these provide the information needed to find the timing uncertainty. The number of data points (8 per one-way pass) minus the number of fitted parameters (3: position, velocity and acceleration) gives the number of degrees of freedom to divide the sum by.

t_unc = sqrt(sum(e_best_list)/(5*n_trips*2)) % the uncertainty estimate
m = 0.1478; % kg (mass of glider)
m_unc = 0.00005; % kg (uncertainty of mass measurements)
F = -m*a_best_list;
plot(v_best_list, F, '+-')
title('Analyzing the Resistance Force')
xlabel('Velocity (m/s)')
ylabel('Force (N)')

t_unc =



This is a disturbing result. It looks like the forward and backwards trips have different acceleration vs velocity profiles. The likely cause of this is a tilted track; while we made an effort to level the track, it apparently wasn’t as level as we had hoped. To estimate the acceleration bias, take the average of the acceleration differences between the forward and backward passes:

F_bias = mean(F(1:2:end)-F(2:2:end))/2
F_adj = F - F_bias*[1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];

F_bias =


The results are close enough that we can make a plausible force vs. velocity fit. Theoretically, aerodynamic drag results in a quadratic dependence on velocity, while viscosity in the layer of air between the glider and track would give a linear relationship between force and velocity. Let’s check both.

p_quad = polyfit(v_best_list, F_adj, 2);
v_quad = (0:0.05:1)*max(v_best_list);
F_quad = polyval(p_quad, v_quad);
[F_lin, F_unc, C_v, F_0, C_v_unc, F_0_unc] = uncertainFit(v_best_list, F_adj)

plot(v_best_list, F_adj, '+-', v_quad, F_quad, '--', v_quad, C_v*v_quad+F_0, ':')
legend('Data', 'Quadratic fit', 'Linear fit', 'Location', 'Northwest')
title('Force after removing bias')
xlabel('Velocity (m/s)')
ylabel('Force (N)')

F_lin =

  Columns 1 through 8 

    0.0018    0.0018    0.0017    0.0016    0.0015    0.0015    0.0014    0.0013

  Columns 9 through 16 

    0.0013    0.0012    0.0011    0.0011    0.0010    0.0010    0.0009    0.0009

  Columns 17 through 22 

    0.0008    0.0008    0.0007    0.0007    0.0006    0.0006

F_unc =


C_v =


F_0 =


C_v_unc =


F_0_unc =


ans =



The quadratic fit between velocity and force shows a substantial force at zero velocity, which seems unlikely. The linear fit nearly goes through zero, and in fact the uncertainty in the zero velocity force is greater than the magnitude of the zero velocity force, suggesting that the resistance force is a simple multiple of the velocity.

Published with MATLAB® 7.1

Posted in Uncategorized | Tagged , , | Leave a comment


I’ve figured out how to make movies from my OpenGL animation, and I’ve posted a video on YouTube, here.

I’ll be giving a talk on the work on Monday, January 24th at CU Boulder, at the Astrodynamics and Satellite Navigation Seminar. It’s at 3:00 pm, probably in the Seebass Forum Room (ECAE 153).

Posted in OpenGL, Uncategorized | Leave a comment

Visualizing Rotation of a Rigid Body

My fans (yes, both of you! Or am I down to one? zero?) may be wondering where I’ve gotten to. I’ve been working up a visualization program for a conference in Nanjing, China.

Back in grad school, I was an aerospace engineer, working in spacecraft dynamics and control. Torque-free rotation of a rigid body is one of those basic physics topics that we studied, on the way to figuring out how to control rotational motion of satellites and such.

My PhD advisor has gotten a big award, to be presented at a conference in Nanjing, and they’d like to have his students and colleagues come to present. Not having worked in the aerospace business for a while, I don’t have any research results to present, but I thought it would be fun to use OpenGL to visualize some of the analytical solutions to the rigid body rotation problem.

Recall what a spinning football does. If you get it spinning exactly around the long axis (the “perfect spiral”), it just spins around that axis. More often, it gets that funny wobble. How do we figure that wobble out?

One way to visualize the motion is to think about the angular momentum and angular velocity vectors (not the same, unfortunately, unless you have a spherically-symmetric object, which footballs aren’t). The angular momentum vector is constant in inertial space (what you’d see outside the vehicle, if you could see angular momentum vectors), but isn’t in a body-fixed frame (what you’d see if you were strapped into the vehicle).

What’s the angular velocity vector? If you think of the object as spinning on an imaginary axle, the angular velocity vector points along the axle, and its length is proportional to the speed of rotation. This axle usually has to be imaginary, since its direction can change relative to the body.

So I’ve worked up a simulation of a rotating, rigid body, with angular momentum and angular velocity vectors added.

A rotating body, with its angular momentum shown as a yellow arrow, and its angular velocity shown as a purple arrow.

This is a frame of the animated image. The yellow arrow is the angular momentum vector, and the purple arrow is the angular velocity vector. The red and yellow ice cream cone thing is the “spotlight” that illuminates the scene.

More images and explanations to come as I make progress with the visualizations!

Posted in OpenGL | Leave a comment

Experimenting with Materials

I’ve been having trouble getting the specular highlights to show up on translucent ellipsoids, so I decided to create a program to let me vary the material colors and light strengths with sliders. I partly based it on Spot, but started from scratch because I wanted to control the material properties and lighting intensities in detail. At first I couldn’t get anything to show up in the OpenGL view, it showed up only as a pure white rectangle. I only found the problem after an almost line-by-line comparison of the new code with Spot. At very end of drawRect:, I had forgotten glSwapAPPLE(), so all of that drawing I did never made it to the window.

At first it seemed that specular lighting wasn’t working, but eventually I realized the problem was that the ambient and diffuse components were overwhelming the specular component: apparently the color has maximum level, and once it reaches that, no additional components will make that part of the object brighter.

Demonstrating Materials and Lighting

It was a great help to have the sliders controlling the lighting and material settings in real time, instead of having to go through the edit-compile-link-run cycle. I also had the settings saved in the ID field of the TGA files, so I can (theoretically, at least) refer to them later on. Actually doing it, though, will require me to create a TGA image reader app that will display the ID text along with the image. Yet another app to create, more practice!

Posted in Uncategorized | Leave a comment