Smoky Cogs

To find the position of the Sun, the first thing we need to do is set the epoch. The epoch is very important since a lot of the values used in astronomical calculations changes over time due to things like precession. The other reason for the epoch is to know the starting point of our calculation.

The trick to finding the Sun’s position is to know the position a a particular epoch, and then, based on the time interval between the epoch and now, determine how far the Sun has moved and thus determine where it is.

The things we need to be able to calculate the sun’s position are the Sun’s mean ecliptic longitude, perigee (point where the earth is closest to the Sun), eccentricity of the Earth’s orbit. All these values are adjusted to be to the values at the time of the epoch set.

We then calculate the eccentric anomaly which I covered in another tutorial, and then find the true anomaly from this, and at the end we get a value for the ecliptic longitude.

Since the Sun travels along the ecliptic, we can set the ecliptic latitude to 0.

The last step in the calculation is to convert the values to equatorial coordinates.

		public static void CalcSunPos(DateTime dDate, DateTime dEpoch, ref double fRA, ref double fDecl)
		{
			double fD;
			double fN;
			double fM;
			double fE;
			double fLambda;
			double fBeta;
			double fSolarMEL;
			double fSolarPL;
			double fSunEarthEcc;
			double fOblique;

            fBeta = 0;
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fSolarMEL = GetSolarMEL(dEpoch, true);
			fSolarPL = GetSolarPerigeeLong(dEpoch, true);
			fSunEarthEcc = GetSunEarthEcc(dEpoch, true);
			fOblique = GetEarthObliquity(dEpoch, true);
			fN = (360.0 / 365.242191) * fD;
			fN = Trig.PutIn360Deg(fN);
			fM = fN + fSolarMEL - fSolarPL;
			fM = Trig.PutIn360Deg(fM);
			fM = Trig.DegToRad(fM);
   			fE = CalcEccentricAnomaly(fM, fM, fSunEarthEcc, 0.0000001);
			fTanV2 = Math.Sqrt((1.0 + fSunEarthEcc) / (1.0 - fSunEarthEcc)) * Math.Tan(fE / 2.0);
			fV = Math.Atan(fTanV2) * 2.0;

            fLambda = fN + fSolarMEL + fE;
			fLambda = Trig.PutIn360Deg(fLambda);

			UraniaCoord.ConvEclToEqu(fOblique, fLambda, fBeta, ref fRA, ref fDecl);
		}

For an orbiting body, the anomaly is the angle between the the semi-major axis and the position of the body in the orbit. Since the vast majority of orbits are not quite true circles, but are instead elliptical, we have a mean anomaly and an eccentric anomaly.

The mean anomaly is the anomaly if the object in question was orbiting in a perfect circle, and can be found using Kepler’s equation

MA = E – esin(E).

The eccentric anomaly is the anomaly of the elliptical orbit, and gives a more realistic value.

Finding the eccentric anomaly is based on finding the mean anomaly, and then correcting for the eccentricity. The function is a recursive function, that takes a guess at the value, and then refines the value, recursively, until the required accuracy is achieved.

		public static double CalcEccentricAnomaly(double fEGuess, double fMA, double fEcc, double fAcc)
		{
			//Calc Ecctrentric Anomaly to specified accuracy
			double fDelta;
			double fDeltaE;
			double fE;
			double fETmp;
			double fEG;

			fEG = fEGuess;

			fDelta = fEG - (fEcc * Math.Sin(fEG)) - fMA;
			if (Math.Abs(fDelta) > fAcc)
			{
				fDeltaE = (fDelta / (1.0 - (fEcc * Math.Cos(fEG))));
				fETmp = fEG - fDeltaE;
				fE = CalcEccentricAnomaly(fETmp, fMA, fEcc, fAcc);
			}
			else
			{
				fE = fEGuess;
			}
			return fE;
		}

Airy discs are a troublesome phenomenon encountered with every telescope. If you magnify an object with a high enough magnification, in every telescope, you reach a stage where the object appears like a disc. This disc has nothing to do with seeing the object itself – like you would see a disc when viewing Jupiter, for example), but rather is side-effect of diffraction of the light coming from the object. This means that it is impossible to get a completely sharp image of a star through a telescope. The size of the airy disc is one of the main determinants of the resolution of a telescope.

The size of the airy disc is dependent on the size of the aperture of the telescope, and the wavelength of the light coming from the object. The larger the wavelength (more red) of the light, the larger the airy disc. This also means for radio telescopes, that “see” in relatively large wavelengths compared to optical telescopes, the airy disc is very large, and can sometimes be half a degree wide or more, severely limiting the resolution of these telescopes.

The formula for calculating the airy disc is:
sine(Airy disc) = 1.22 * Wavelength / Aperture
but since the angle of the airy disc is very small we can approximate by:
Airy disc = 1.22 * Wavelength / Aperture

		public static double CalcWavelengthAiry(double pdAiryDiam, double pdApertureAiry)
		{
			return pdAiryDiam * pdApertureAiry / 1.22;
		}

		public static double CalcApertureAiry(double pdAiryDiam, double pdWavelengthAiry)
		{
			return pdWavelengthAiry * 1.22 / pdAiryDiam;
		}

		public static double CalcAiryDiam(double pdApertureAiry, double pdWavelengthAiry)
		{
			return pdWavelengthAiry * 1.22 / pdApertureAiry;
		}

The focal ratio, or f-ratio as it is also known, is a very important number when we start talking about astro-photography.

The f-ratio is defined as the focal length of the telescope divided by the aperture of the telescope, and is known as the speed of the telescope. that is, a telescope with a f-ratio of f/4 is “faster” than a telescope with a f-ratio of f/8.

The “faster” the telescope, the brighter the image will be, which means that the length of exposures for photography don’t need to be as long to capture the same amount of light. The advantage of a “slower” telescope though is that images will have greater depth, meaning that sharpness over the whole image is better, so there is a little bit of a trade-off depending on what you are trying to do.

		public static double CalcFRatio(double FocalLen, double Aperture)
		{
			return FocalLen / Aperture;
		}

		public static double CalcFocalLen(double FRatio, double Aperture)
		{
			return FRatio * Aperture;
		}

		public static double CalcAperture(double FRatio, double FocalLen)
		{
			return FocalLen / FRatio;
		}

One of the most misunderstood concepts in astronomy for a person who has never owned a telescope is how the magnification is calculated.

Calculating the magnification is dependent on the focal length of the telescope itself and the focal length of the eyepiece.

The two main types of optical telescopes are reflectors and refractors. A reflector uses a large curved mirror which focuses the incoming light, which eventually passes through an eyepiece, which diverges the light, so that it can form an image on your retina. A refractor has a lens which focuses the light instead of a mirror, and the light then passes through the eyepiece to diverge the light just as in with a reflector.

The focal length is the distance from the mirror or lens to the point where the light focuses.

For a particular telescope, the focal length of the telescope remains the same. Different magnifications are found by using eyepieces with different focal lengths.

The formula for working out magnification is:
Magnification = Telescope Focal Length / Eyepiece Focal Length

So, for example, if you have a telescope with a focal length of 900mm, and use an eyepiece with a focal length of 9mm, you will see objects magnified by 100x.

The inverse calculations are also easy, so you can find which eyepiece you need to use to magnify by, say, 50x very easily.

		public static double CalcTeleFocalLenMag(double pdEPFocalLen, double pdMag)
		{
			return pdEPFocalLen * pdMag;
		}

		public static double CalcEPFocalLenMag(double pdTeleFocalLen, double pdMag)
		{
			return pdTeleFocalLen / pdMag;
		}

		public static double CalcMagnification(double pdTeleFocalLen, double pdEPFocalLen)
		{
			return pdTeleFocalLen / pdEPFocalLen;
		}

It is often useful to know in which constellation an object appears, since it is more intuitive to find something in the sky when you know which constellation you are looking for the object in.

Finding the constellation for objects which travel along near the ecliptic this is very easy. The ecliptic is the path along which the Moon, Sun and planets move, so this is particularly applicable to these objects.

The twelve constellations on the ecliptic form the zodiac, and to find the constellation which the object is in, all you need is the ecliptic longitude of the object, since the boundaries of the constellations are fixed. By passing the right ascension and declination of the object, and then using the conversion to ecliptic coordinates, we just need to check which constellation is within the ecliptic longitude value.

		public static string GetCurrentConst(double fOblique, double fRA, double fDecl)
		{

			double fELong = 0;
			double fELat = 0;
			string sConst;

			ConvEquToEcl(fOblique, fRA, fDecl, ref fELong, ref fELat);
            
			if (fELong <  33.18)
			{
				sConst = "Pisces";
			}
			else if (fELong <  51.16)
			{
				sConst = "Aries";
			}
			else if (fELong <  93.44)
			{
				sConst = "Taurus";
			}
			else if (fELong <  119.48)
			{
				sConst = "Gemini";
			}
			else if (fELong <  135.30)
			{
				sConst = "Cancer";
			}			
			else if (fELong <  173.34)
			{
				sConst = "Leo";
			}
			else if (fELong <  224.17)
			{
				sConst = "Virgo";
			}
			else if (fELong <  242.57)
			{
				sConst = "Libra";
			}
			else if (fELong <  271.26)
			{
				sConst = "Scorpius";
			}			
			else if (fELong <  302.49)
			{
				sConst = "Sagittarius";
			}
			else if (fELong <  311.72)
			{
				sConst = "Capricorn";
			}
			else if (fELong <  348.58)
			{
				sConst = "Aquarius";
			}
			else
			{
				sConst = "Pisces";
			}
			return sConst;
		}

I have already addressed equatorial coordinates, so will now focus on galactic coordinates.

The galactic coordinate system is the coordinate system based on the plane of the Milky Way galaxy, which, on a clear night can be seen arching across the sky.

Galactic coordinates are particularly useful when studying objects pertaining to the Milky Way galaxy itself as opposed to our local galactic neighbourhood, and would include objects such as globular clusters or supernovae.

The galactic longitude is denoted by ℓ and the galactic latitude by b. The galactic longitude is measured from the galactic centre, and is in the range of 0°-360°, using the Sun as the centre of the celestial sphere. The galactic latitude is in the range -90° – 90° with the galactic poles being perpentidular to the galactic plane.

To convert between the equatorial and galactic systems, we need the right ascension and declination of the galactic north pole, as well as the ascending node of the milky way, which is the point (in degrees) where the galactic plane crosses the equatorial plane.

For the Milky Way, the right ascension is 192.25° and the declination is 27.4°, with the ascending node being 33°.

		public static void ConvEquToGal(double fGalNPRA, double fGalNPDecl, double fGalPlaneAscNode, double fRA, double fDecl, ref double fL, ref double fB)
		{
			double fSinB;
			double fX;
			double fY;
			double fRADiff;

			fGalNPRA = fGalNPRA * 15;
			fRA = fRA * 15;

			fGalNPRA = Trig.DegToRad(fGalNPRA);
			fGalNPDecl = Trig.DegToRad(fGalNPDecl);
			fRA = Trig.DegToRad(fRA);
			fDecl = Trig.DegToRad(fDecl);
			fRADiff = fRA - fGalNPRA;
			fSinB = (Math.Cos(fDecl) * Math.Cos(fGalNPDecl) * Math.Cos(fRADiff)) + (Math.Sin(fDecl) * Math.Sin(fGalNPDecl));
			fB = Math.Asin(fSinB);
			fY = (Math.Sin(fDecl)) - (Math.Sin(fB) * Math.Sin(fGalNPDecl));
			fX = (Math.Cos(fDecl) * Math.Sin(fRADiff) * Math.Cos(fGalNPDecl));
			fL = Math.Atan(fY / fX);
			fL = Trig.RadToDeg(fL);
			fB = Trig.RadToDeg(fB);
			fL = Trig.TanQuadrant(fX, fY, fL);
			fL = fL + fGalPlaneAscNode;
			fL = Trig.PutIn360Deg(fL);
		}

		public static void ConvGalToEqu(double fGalNPRA, double fGalNPDecl, double fGalPlaneAscNode, double fL, double fB, ref double fRA, ref double fDecl)
		{
			double fSinDecl;
			double fX;
			double fY;
			double fLDiff;

			fGalNPRA = fGalNPRA * 15;

			fGalPlaneAscNode = Trig.DegToRad(fGalPlaneAscNode);
			fGalNPDecl = Trig.DegToRad(fGalNPDecl);
			fL = Trig.DegToRad(fL);
			fB = Trig.DegToRad(fB);

			fLDiff = fL - fGalPlaneAscNode;
			fSinDecl = (Math.Cos(fB) * Math.Cos(fGalNPDecl) * Math.Sin(fLDiff)) + (Math.Sin(fB) * Math.Sin(fGalNPDecl));
			fDecl = Math.Asin(fSinDecl);
			fY = (Math.Cos(fB) * Math.Cos(fLDiff));
			fX = (Math.Sin(fB) * Math.Cos(fGalNPDecl)) - (Math.Cos(fB) * Math.Sin(fLDiff) * Math.Sin(fGalNPDecl));
			fRA = Math.Atan(fY / fX);
			fRA = Trig.RadToDeg(fRA);
			fDecl = Trig.RadToDeg(fDecl);
			fRA = Trig.TanQuadrant(fX, fY, fRA);
			fRA = fRA + fGalNPRA;
			fRA = Trig.PutIn360Deg(fRA);

			fRA = fRA / 15.0;
		}