## 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;
}
}
```