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
{
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);
var altitudeRadians = Math.Asin(Math.Sin(declinationRadians) *
Math.Sin(latitudeRadians) + Math.Cos(declinationRadians) *
Math.Cos(latitudeRadians) * Math.Cos(hourAngleRadians));
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;
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;
}
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;
}
}