Motorizing a Telescope, Part 2

• Download source code

Introduction

This blog entry is about the software and algorithms I wrote to control the motorized telescope described in Part 1. It would be very helpful to read that post before proceeding with this one.

The Motors

To explain the algorithm, let’s start with the motors, and then work our way to the stars. The motors I am controlling are basically stepper motors, with a built-in intelligent controller and encoder. They are called MDrives. I have already described the motors in detail here in my four-part MDrive series of posts.

If you are working on your own project, I suspect you will not be using MDrives, as they are pricey. Whatever stepper motors you use, you will have to figure out how to control them on your own. You will definitely need to micro-step them. The motors I am using run at 256 micro-steps per step, yielding 51200 (micro-) steps per revolution.

For your project, I doubt you will need an encoder. However, since the MDrive come with a built-in encoder, I decided to use those. Because of the encoder resolution, all my code works as if there are 2048 (encoder-) steps per revolution.

MDrives

I am taking advantage of two distinct modes of movement that the MDrives provide: a constant slew rate, and an absolute movement. A constant slew rate is pretty simple: I specify how many steps per second I want the motor to run at. I can change this rate while the motor is running and it seamlessly transitions to the new rate.

The absolute movement is a bit trickier, where the motor accelerates to a maximum speed and then decelerates to arrive at a final destination. I don’t think you absolutely have to implement this second mode, but it makes for very nice movements from one celestial target to the next. For this project, I am choosing an acceleration rate of 0.5 revolutions per second per second (A=1048) and the same deceleration rate. I am also limiting the maximum velocity to just under 5 rotations per second (VM=10000). Because each rotation of the motor moves the telescope one degree, this system will slew at up to 5 degrees per second, and take 5 seconds to achieve maximum speed.

I am limiting the current to 25% of the rated motor capacity. If the current to a stepper motor is too small, they will lose synchronization. If the current is too large, the motor will vibrate and be noisy.

MDriveViewModel

The interface between a motor and the rest of the software system can be found in MDriveViewModel class. This class exposes just three properties for a motor, Position (read), MoveAbsolute (write), and Slew (write).

The rest of the system can, at any time, read the position (in terms of encoder-steps). Position 0 is defined as where the motor was when the software started. If the motor shaft has rotated 3 times since the software started, the position would be either 6144 or -6144 (depending on direction).

The MoveAbsolute property allows the rest of the system to move the motor shaft to a specific position. To move the telescope 30 degrees from where it was when the software started, setting the MoveAbsolute property to 61440 will do the trick. This would take about 8 seconds with the acceleration and deceleration profiles automatically controlled.

The Slew property allows the rest of the system to move the motor shaft at a constant rate. To constantly move the telescope at 0.5 degrees per second, setting the Slew property to 1024 would work well. Setting the Slew or MoveAbsolute properties preempts any previous motor movements.

There are two instances of this class, one for each motor.

ComViewModel

The class that initializes and maintains communications with a motor is found in ComViewModel. This class runs autonomously, connecting to the MDrive via the serial port. It searches for the appropriate serial port and then configures the motor as needed (enabling the encoder, setting the acceleration rate, etc.).

Once a second, the ComViewModel polls the MDrive to obtain its position and then updates the Position property of the MDriveViewModel class. Whenever the MDriveViewModel class has a new target movement or slew request, the ComViewModel sends the appropriate commands to the MDrive.

There are two instances of this class, one for each motor.

SlewRateCalculatorViewModel

The SlewRateCalculatorViewModel class has a writeable property called CurrentMotorPostion. The SlewRateCalculator also has writeable properties called TargetDegrees and TargetDegreesPerSecond. These are set by higher-level classes to be described later. Given the three input properties, the SlewRateCalculatorViewModel calculates and updates these three output properties: MoveOrSlew, CalculatedMoveAbsolute, and CalculatedSlewRate.

This calculation decides if the current position is close to, or far away from the target. If it is far away, it will call for the motor to be moved to an absolute position. This allows the motor to accelerate, move quickly, and decelerate smoothly. If the target is very close (within a degree), then it will compute the slew rate of the motor. The calculated slew rate is the sum of the requested TargetDegreesPerSecond and a catch up slew rate. The catch up slew rate is the rate at which the motor will achieve the target position in three seconds. In other words, the catch up slew rate is the difference between the CurrentMotorPostion and the TargetDegrees divided by three.

There is an added complication to the SlewRateCalculatorViewModel class having to do with calibration. This should be ignored for now.

There are two instances of this class, one for each motor.

MainViewModel

At this point, I want to point out that the MainViewModel class ties all the other view models together. For example, whenever the Position of the MDriveViewModel changes, the MainViewModel updates the CurrentMotorPostion of the SlewRateCalculatorViewModel class. The MainViewModel also binds the outputs of SlewRateCalculatorViewModels to the inputs of the MDriveViewModels.

Using the MainViewModel to couple the classes together like this allows the MainViewModel to present a user interface where we can connect and disconnect view models at will. This makes, say manually overriding a motor easy to implement.

AltAzimuthCalculatorViewModel

The AltAzimuthCalculatorViewModel class has writeable properties called DeclinationDegrees and RightAscensionHours. These are set by higher-level classes identifying the target to observe. This class also has the writeable properties Latitude and Longitude.

Once a second, the AltAzimuthCalculatorViewModel class calculates and updates these output properties: CalculatedAltitudeDegrees, CalculatedAzimuthDegrees, CalculatedAltitudeDegreesPerSecond, CalculatedAzimuthDegreesPerSecond, and for display purposes, SiderealTime.

To compute SiderealTime, it makes use of the CalculateLocalSiderealTime function in the static class CelestialCalculations.

To compute CalculatedAltitudeDegrees and CalculatedAzimuthDegrees, it makes use of the ComputeAltAlzimuth function in the static class CelestialCalculations.

To calculate CalculatedAzimuthDegreesPerSecond and CalculatedAltitudeDegreesPerSecond, AltAzimuthCalculatorViewModel computes the ComputeAltAlzimuth function twice – once for the current sidereal time and once for a second into the future. The degrees per second is then the difference in these two times.

The MainViewModel class ties the outputs of AltAzimuthCalculatorViewModel to the inputs of the SlewRateCalculatorViewModel classes.

CelestialCalculations

This class is the real brains of the operation. Instead of me transcribing the code to English, here is the code:

static class CelestialCalculations
{
    // http://www.stargazing.net/kepler/altaz.html
    public static void ComputeAltAlzimuth(DateTime localSiderealDateTime, 
	double rightAscentionHours, double declinationDegrees, double latitudeDegrees, 
	out double azimuthDegrees, out double altitudeDegrees)
    {
        var localSiderealDegrees = DegreesFromDateTime(localSiderealDateTime);

        var rightAscentionDegrees = DegreesFromHours(rightAscentionHours);
        var hourAngleDegrees = localSiderealDegrees - rightAscentionDegrees;
        if (hourAngleDegrees < 0)
            hourAngleDegrees += 360;

        var hourAngleRadians = RadiansFromDegrees(hourAngleDegrees);
        var declinationRadians = RadiansFromDegrees(declinationDegrees);
        var latitudeRadians = RadiansFromDegrees(latitudeDegrees);

        // Alt = asin( sin(Dec) * sin(Lat) + Cos(Dec) * Cos(Lat) * Cos(Ha) )
        var altitudeRadians = Math.Asin(Math.Sin(declinationRadians) * 
	Math.Sin(latitudeRadians) + Math.Cos(declinationRadians) * 
	Math.Cos(latitudeRadians) * Math.Cos(hourAngleRadians));

        // Az = acos( sin(Dec) - sin(Alt) * sin(Lat) ) / ( Cos(Alt) * Cos(Lat) )
        var azimuthRadians = Math.Acos((Math.Sin(declinationRadians) - 
	Math.Sin(altitudeRadians) * Math.Sin(latitudeRadians)) / 
	(Math.Cos(altitudeRadians) * Math.Cos(latitudeRadians)));

        altitudeDegrees = DegreesFromRadians(altitudeRadians);
        azimuthDegrees = DegreesFromRadians(azimuthRadians);

        if (Math.Sin(hourAngleRadians) > 0)
        {
            azimuthDegrees = 360 - azimuthDegrees;
        }
    }

    public static double DegreesFromDateTime(DateTime dateTime)
    {
        return dateTime.TimeOfDay.TotalHours * 360 / 24;
    }

    public static double DegreesFromHours(double hours)
    {
        return hours * 15;
    }

    public static double RadiansFromDegrees(double degrees)
    {
        const double toRadians = 2 * Math.PI / 360.0;
        return degrees * toRadians;
    }

    public static double DegreesFromRadians(double radians)
    {
        const double toDegrees = 360.0 / (2 * Math.PI);
        return radians * toDegrees;
    }

    public static void ComputeRightAscentionDeclination
	(double localSiderealDegrees, double azimuthDegrees, double altitudeDegrees, 
	double latitudeDegrees, out double rightAscentionHours, out double declinationDegrees)
    {
        const double toRadians = 2 * Math.PI / 360.0;

        var latitudeRadians = latitudeDegrees * toRadians;
        var altitudeRadians = altitudeDegrees * toRadians;
        var azimuthRadians = azimuthDegrees * toRadians;

        // wolfram: solve a = acos( sin(d) - sin(2) * sin(3) ) / ( Cos(2) * Cos(3) ) for d
        // wolfram solution: d =  -asin( cos( a * cos(2) * cos(3) + sin(2) * sin(3)) + 2 pi n + pi
        var declinationRadians = Math.Asin(Math.Cos(azimuthRadians * 
	Math.Cos(altitudeRadians) * Math.Cos(latitudeRadians) + Math.Sin(altitudeRadians) * 
	Math.Sin(latitudeRadians))) + Math.PI;
        declinationDegrees = declinationRadians / toRadians;

        rightAscentionHours = 0;
    }

    // http://answers.yahoo.com/question/index?qid=20100509150032AAjyxX6
    public static DateTime CalculateLocalSiderealTime(double longitude, DateTime utcNow)
    {
        var epoch = new DateTime(2000, 1, 1, 12, 0, 0, DateTimeKind.Utc);
        var sinceEpoch = utcNow - epoch;
        var hours = (18.697374558 + 24.06570982441908 * sinceEpoch.TotalDays) % 24 + longitude / 15.0;
        while (hours < 0) hours += 24;
        while (hours >= 24) hours -= 24;
        var localSiderealTimeSpan = TimeSpan.FromHours(hours);
        return utcNow.Date + localSiderealTimeSpan;
    }
}