A Semi-Automatic Method to Measure the Rotation of Sunspots

Sunspots have been observed to undergo rotation about their umbral centre. This is typically a slow rotation, with even the fastest sunspot rotations only reaching angular velocities of a few degrees per hour. This rotation may inject magnetic energy into the Sun’s atmosphere, which can be stored in the coronal magnetic field and later released in eruptive events such as solar flares and coronal mass ejections. To usefully investigate rotating sunspots long periods of data need to be analysed, often of the order of several days, to build up a bulk rotation profile for the sunspot over time. This article outlines a semi-automated approach for analysing series of solar continuum data to extract the rotation profile of a sunspot as it transits across the solar disc. Moving towards an automated approach is vital for generating large, unbiased statistical samples of rotating sunspots in order to understand their contribution to solar activity. Existing methods typically focus on sunspots near disc centre for short time periods, neglecting much of the rotation history of the sunspot. The method is tested on six sunspots observed in continuum data from the Helioseismic and Magnetic Imager (HMI) instrument on board the Solar Dynamics Observatory (SDO). These have been chosen to test the method for a range of different types of sunspots, including well-behaved sunspots, shape-changing sunspots, fast rotators, non-rotators, and interacting sunspots. The rotation profiles are compared by eye to animations of the sunspot from the data and are in acceptable visual agreement with the observed bulk rotation of the sunspot for all of the cases, except for the one which contains two sunspots in a shared penumbra. The method is also tested against sunspot rotations in active region (AR) 11158 that have been reported in the literature. While the results compare to some degree, the method outlined in this article reports lower rotations than those reported in the literature. Some of this discrepancy can be attributed to selection bias by the approaches in the literature, where only features that undergo larger rotation are tracked in sunspots that exhibit non-uniform rotation. The method also provides uncertainties on the calculated rotation profile which can be broken down to allow the principal sources of error to be identified. For the test sunspots in this article, the dominant source of uncertainty is the resolution of the SDO/HMI instrument.


Introduction
Sunspots represent regions where the strongest magnetic field protrudes from the solar interior, where it is generated, into the atmosphere to play a defining role for solar activity. This means that the atmospheric magnetic field is heavily driven by the dynamics of sunspots (Berger et al., 1998). Such dynamics include shear (Kazachenko et al., 2010) and rotation (Evershed, 1910;Stenflo, 1969), and the twist that these introduce can lead to the formation of coronal structures such as sigmoids (Min and Chae, 2009;Alexander, 2006, 2008). Furthermore, the energy injected by sunspot rotation can lead to solar eruptions such as solar flares (Kumar et al., 2013;Vemareddy, Cheng, and Ravindra, 2016;Zhang, Liu, and Zhang, 2008) and coronal mass ejections (CMEs) (Török et al., 2013;Vemareddy, Cheng, and Ravindra, 2016;Yan et al., 2012). Furthermore, modelling of a solar flare from May 2005 has demonstrated that sunspot rotation dramatically increases the energy released by a flare (Kazachenko et al., 2009).
Statistical surveys using data from the Michelson Doppler Imager (MDI) on board the Solar and Heliospheric Observatory (SOHO) have been carried out (Li and Zhang, 2009;Yan, Qu, and Xu, 2008;Yan, Qu, and Kong, 2008), but these suffer from selection bias, i.e. the authors select sunspots that clearly rotate, and the calculation of the rotation is crude. Zhu, Alexander, and Tian (2012) have performed a more robust SOHO/MDI survey using the differential affine velocity estimator (DAVE) method (Schuck, 2006), but focussed on larger active regions near disc centre without linking the rotation to solar activity.
The term "sunspot rotation" is used to describe the rotation of a sunspot about its umbral centre. This is generally taken to be an average or typical rotation or rotation speed at any given time. However, it is known that some sunspots can exhibit non-uniform rotation, and variations in the rotation in the radial (e.g. Brown et al., 2003) and axial (e.g. Jiang et al., 2012) directions are observed. Despite this non-uniform perturbation, the bulk property provided by the average rotation of the sunspot can be a useful descriptor for the broad dynamics of a sunspot when compiling large samples as the less "complicated" quantities enable global inferences to be more easily drawn. For example, when making a statistical comparison between the rotation of sunspots and the activity (such as solar flares or CMEs) of the associated active region, a simple estimate of the energy injected by sunspot rotation (based on the average rotation) against the energy expended by solar activity is needed. Exploration of non-uniform rotations, particularly where the effect is large, may be important for some case studies.
There are various methods for determining sunspot rotation in the literature. The simpler methods just geometrically estimate the location of principal axes or features (such as a light bridge) at different times and take the difference in angle to be the rotation (e.g. Louis et al., 2012;Yan et al., 2012). The method of Gopasyuk and Gopasyuk (2005) uses single frames of spectropolarimetric observations to decompose line-of-sight velocity measurements for sunspots away from disc centre into radial and tangential components, which are used to determine the angular velocity of the sunspot rotation at that time. This method has not been applied to extended time frames to produce rotation profiles of a sunspot. The method of Brown et al. (2003) uses sequences of continuum data and tracks the rotational motion of features in the penumbra (such as fibrils) to determine the sunspot rotation. Several authors follow an approach similar to Brown et al. (2003), but introduce bias to the rotation measurements by manually selecting only a small number of penumbral features for tracking at a limited number of radial positions (e.g. Jiang et al., 2012;Vemareddy, Ambastha, and Maurya, 2012). The DAVE (Schuck, 2006) and differential affine velocity estimator for vector magnetograms (DAVE4VM) (Schuck, 2008) methods use magnetogram and vector magnetogram data to invert the induction equation to determine the photospheric velocity. This approach is good for determining the magnetic energy input, but often relies on complementary methods (such as those described above) to convert the velocities to rotation.
This article describes the implementation of a semi-automatic analysis algorithm, tailored to use SDO/HMI continuum observations. The approach is based on the method of Brown et al. (2003), but adds key new features: i) the approach is optimised to more rapidly handle a large quantity of data; ii) choices of what were previously free parameters in the analysis are now automated (or at least systematically determined and justified), these choices could previously influence the outcome and lead to bias in the rotation profiles; iii) the method dynamically determines and varies penumbral radii to respond to changes in size and shape of the sunspot over time to better capture the average rotation; iv) projection effects for sunspots further away from solar disc centre are removed so they do not contaminate the rotation profile; v) an uncertainty for the rotation measurement at each time is calculated to reflect both errors due to parameters determined by the algorithm and the variation due to nonuniform rotation.
It differs from related approaches in the literature in that it identifies features at all radial positions in the penumbra (and extended features at multiple radial locations) and at as many axial positions as possible to include variations due to a non-uniform rotation. The only inputs now required for the algorithm are a list of SDO/HMI continuum observations containing a sunspot with a reasonable guess for the centre of the sunspot in the first observation (which serves solely as a sunspot identifier). The algorithm then produces a rotation profile for the duration of the observations along with robust uncertainties on the measured rotation, this is done autonomously and reduces the need for manual intervention or parameter choice speeding up the analysis process. The algorithm is tested with six sunspots with different properties, including one that demonstrates when the algorithm breaks down and how this is handled. Development of such an algorithm is an important step for: (i) producing robust, biasfree rotation profiles with well understood uncertainties; (ii) easily capturing the full rotation profile of all sunspots in an active region as it traverses the solar disc (rather than select sunspots at selected times); and (iii) generating large statistical samples of the rotation profiles of sunspots. This is necessary for addressing questions regarding what proportion of sunspots rotate, how much do sunspots typically rotate, why do sunspots rotate, and what (if any) is the link between sunspot rotation and solar activity.

Data
The analysis method described in this article has been developed to work specifically with sequences of continuum observations from the Helioseismic and Magnetic Imager (HMI) instrument (Schou et al., 2012) on NASA's Solar Dynamics Observatory (SDO) (Pesnell, Thompson, and Chamberlin, 2012). In particular, observations are drawn from the hmi.Ic_45s series which produces continuum observations of the surface of the Sun in the spectral region of the 6173 Å, FeI absorption line at a cadence of typically 45 s.
However, as an individual sunspot may be analysed for 8-10 days as it transits the solar disc, only a subset of the available observations are analysed, and observations at a cadence of 3 minutes are selected. As sunspot rotation rates are typically in the range of 0 -3 • h −1 , or up to 0.15 • of rotation in a 3 minute period, this cadence strikes a reasonable balance between capturing the detail of the rotation and data management.
To test the method, six sunspots drawn from five different active regions are analysed. These are selected as they exhibit different properties to fully test the robustness of the methodology. A list of the active regions can be found in Table 1, along with the dates and number of observations analysed.
AR 11147 is a well behaved example. It is a medium sized sunspot that remains circular for the duration of the analysis. It is not close to any other sunspots, though there is some minor merging with and fragmentation of small pores. The sunspot exhibits some rotation in both directions. This sunspot served as the primary test case during the development of the rotation analysis.
AR 11770 is also a well behaved example, though it is slightly smaller than AR 11147. It retains a circular shape throughout the analysis, but does exhibit some shrinking towards the end of it. It was selected as it does not exhibit any rotation, so it will be used to test that no spurious rotation signal is added by the rotation analysis.
AR 12158 is a high rotator, and is used to test whether the rotation analysis will successfully analyse such cases. It is also a large sunspot which exhibits changes in size (due to fragmentation of pores) and shape (becoming noticeably elliptical for a period), so helps to test the robustness of the rotation analysis to these changes.
AR 12420 comprises two sunspots that are close enough together that they may interfere with the analysis of each other. This tests the robustness of the rotation analysis to interference from nearby sunspots. Both sunspots are comparatively small, and exhibit changes of size and shape throughout the analysis.
AR 12422 is a more complex example in which for a period of the analysis two sunspots share the same penumbra. This example is used to demonstrate when the rotation analysis breaks down.
For each sunspot, two further restrictions are placed on the sequence of observations to be analysed. The first is that the sunspot should not be too close to the solar limb. The specific criteria used is that the angular distance from disc centre must be at most 60 • , i.e.
where (x, y) give the two-dimensional projected location of the sunspot in the observation, (x c , y c ) gives the location of the solar disc centre in the observation, and R s is the radius of the Sun in the same units as the (x, y) values. At this angle, the projected area of the sunspot is 50% of what it would be at disc centre, and the sunspot can reasonably be thought of as facing the observer. Beyond this the sunspot can be thought of as facing the limb. This threshold is taken to be a reasonable boundary between the limb and on disc, and provides sufficient coverage of the sunspot to measure any rotation. The second restriction is a minimum size criteria on the sunspot being analysed, as image resolution is lost at small radii. An area limitation on the reprojected area of the sunspot umbra (which will be strictly defined in Section 3.1) is applied based on the radius of an idealised circular sunspot. A circular sunspot with a radius of 7 pixels would have an area of 49π pixels 2 , and this is taken as a limit for the minimum size of the sunspot. The analysed duration of the sunspot is taken to be a continuous run of observations for which the area of the sunspot is above this threshold.
In some cases, the area may dip below this limit for a short period of time (up to a few hours say) before rising above it again. This is particularly true for sunspots whose size is around the limit, or that shrink near the end of their lifetimes. In order to prevent repeated flipping of a sunspot between above and below the limit, a hysteresis zone is created with the previous threshold being the upper soft-limit and a lower hard-limit of 36π pixels 2 (a circle with a 6 pixel radius) applied. If the umbral area dips below the soft-limit into the hysteresis zone for a short period, then the sequence of observations is continued, but if the area falls below the hard limit then the analysis is terminated.
An example of this is shown in Figure 1 for AR 11770. The sunspot rotates into the analysis area for about 19 hours with a reprojected umbral area of around 250 pixel 2 . It stays at an area between 250 -300 pixel 2 until about 145 hours before shrinking. It reaches a brief plateau around 200 pixel 2 with some short drops below the 49π pixel 2 soft limit, before finally dropping below both the hard and soft limits at about 164 hours. In this case, the sunspot has developed a light bridge while shrinking, dividing the sunspot into two sub- size sunspots. Occasional gaps in the light bridge allow the full sunspot to be captured for a few frames later on (e.g. around 190 hours), but analysis would be truncated before this.
These criteria are applied manually before the sunspot rotation code is run to define a list of observations to pass to the code for analysis. Together with an initial guess for the centre of the sunspot in the first image (which is used purely for initial identification of the chosen sunspot in the observation), this represents the semi-part of the semi-automatic method. Once these are selected, the analysis code runs fully automatically. The analysis code will run in violation of these restrictions if desired, but this leads to larger uncertainties on the calculated rotations.

Methodology
The SDO/HMI continuum data is analysed by comparing subsequent observations. The same sunspot is identified in each image and the centre of the sunspot is found. The sunspot is uncurled about its centre from an x-y geometry to an r-θ one so that any rotational motion is decomposed into lateral motion in the θ direction only.
The intensity profile of each radial band in the sunspot is investigated and the turning points (local minima and maxima) are found. In the penumbra, this will pick out the penumbral fibrils and track their lateral (i.e. rotational) motion. Features in the photosphere outside of the sunspot are generally less reliable, partly because they do not usually exhibit a radial structure, and partly because they are not obviously part of the sunspot.
The turning points are matched to equivalent turning points in the previous observation, and the angular discrepancy between the two is used to calculate the amount of rotation that has occurred in the elapsed time between the two observations.
These rotation values can then be combined to produce a cumulative rotation profile for the sunspot along with an uncertainty. This process is outlined in the flowchart in Figure 2.
Inputs to the analysis are a sequence of SDO/HMI continuum observations, and a reasonable approximation to the centre of the sunspot in the first observation of the sequence (the location of an umbral pixel in the image is sufficient).

Defining Umbral and Penumbral Image Pixels
Intensity thresholds for defining umbral and penumbral pixels are required for the analysis. These are obtained by extracting the DATAMEAN value from each image header in the observation sequence and taking the average of these values, i.e.
where N is the number of observations in the sequence. This is multiplied by two scale factors to produce the following pixel thresholds: i) umbral pixel: I pix ≤ 0.6Ī , ii) penumbral pixel: 0.6Ī < I pix ≤ 1.05Ī .
This typically gives an umbral threshold in the range 29 000 -34 000 DN and a penumbral threshold in the range 51 000 -59 000 DN.
A minimum threshold of 1000 DN is applied to ensure that bad pixels, data dropouts, and off-limb pixels are not unintentionally included in the sunspot.

Finding the Centre of the Sunspot
The observations used to locate the centre of the sunspot are first corrected for limb darkening, as outlined in Appendix A, to ensure a consistent sunspot intensity across the solar disc, and to allow fixed intensity thresholds to be used for the definition of umbral and penumbral pixels for a given sunspot.
The centre location of the sunspot from the previous observation is used as an initial estimate for the new observation, with the supplied "guess" being used for the first observation.
The image pixel containing the estimate is tested, if its corrected intensity value is below the umbral threshold then it is used as start point to determine the umbra of the sunspot in Figure 3 Results of the umbral pixel finding and grouping method to an observation of AR 12420. The two sunspots in the active region along with their initial guesses marked with a + sign is shown on the left, with the located grouped umbral pixels for each initial guess on the right. the current observation. If it is not an umbral pixel, then the nearest valid umbral pixel to the estimate is found and used as a start point.
A contiguous group of pixels defining the umbra of the sunspot are found by examining the neighbouring pixels of the start point to see if they are umbral pixels, then testing their neighbours, and so on until no new umbral pixels are found. An example of the results this produces is shown in Figure 3, where a frame from AR 12420 is used. This shows that initial guesses for each of the two sunspots in the active region produces separate umbral groups without contamination from the other sunspot or nearby pores.
The centre of the sunspot is taken to be the mean location of the group of umbral pixels. So if the identified umbral pixels have locations given by (x k , y k ) where k = 1, . . . , N with N being the number of pixels in the group, then the x and y centre position of the sunspot is given by These have standard deviations σ x and σ y , and standard error on the mean of where δ x and δ y are the basic quantised scale error of 1 pixel and are assumed to follow a uniform distribution. These errors are calculated for all observations, however, the errors have a curious property that their dependence on the size of the sunspot is very weak, and there is a stronger dependence on their shape such that the errors change (one increases and the other decreases) if the sunspot is elongated in a particular direction. This is due to the fact that no pixel can be double sampled and there is a direct link between the size of the sunspot and the number of pixels it contains. This is discussed further in Appendix B.
This means that a characteristic error for the location of the centre of the sunspot can be defined, which reasonably approximates (or even overestimates) the individual errors. This can be reasonably defined as S x = S y = 0.6 pixels (see Appendix B), and is useful for estimating the effect of the uncertainty of the centre of the sunspot on the calculated rotation profiles without pre-calculation of these uncertainties.

Figure 4
Illustration of the steps taken to remove projection effects: (a) the sunspot is mapped onto the surface of a sphere the size of the Sun, (b) this is rotated about the north-south axis to transform the sunspot to the central meridian, (c) this is in turn rotated about an east-west axis to transform the sunspot to disc centre.

Uncurling the Sunspot
Once the centre of the sunspot is found, the limb-corrected x-y sunspot image can be transformed to polar r-θ space about the sunspot centre, where r is the radius from the centre of the sunspot and θ is the anti-clockwise angle about the centre from a westward pointing chord.
Projection effects due to the sunspot generally not being near Sun centre are removed by transforming the sunspot to disc centre using a pair of spherical rotations that preserve solar north. This is illustrated in Figure 4, where the sunspot image is mapped to a sphere the size of the Sun, rotated about the north-south axis, then rotated about an east-west axis to move the sunspot to disc centre while preserving solar north.
In practice, it is more useful to perform these steps in reverse, i.e. define a series of r-θ locations about disc centre and convert these to Cartesian coordinates, where (0, 0) is disc centre, using the equations x = r cos (θ ), and y = r sin (θ ).
This is mapped on to the spherical surface of the Sun by defining the z coordinate where R s is the radius of the Sun in the same units as x and y. This can be written as the position vector This defines a Cartesian space where x and y are in the plane of the observation with y pointing to solar north, and z is perpendicular to this plane pointing towards the observer. The spherical rotations are calculated using the position of the centre of the sunspot relative to Sun centre, (x cen , y cen ), via their longitudinal and latitudinal positions, (α, β), given by x cen r 2 s − y 2 cen , and sin (β) = y cen r s .
These are used to define two rotation matrices: A defines a rotation about the y axis to transform the longitudinal position, and B defines a rotation about the x axis to transform the latitudinal position. These are given by In order to preserve the solar north in the transformed image, the B matrix must be applied to the position vector to be transformed first, i.e. v ss = ABv.
This gives the transformed location of the position about the sunspot centre, where (x ss , y ss ) is the location in the observed image. These transformed points are typically not at pixel centres, and for smaller radii, a single observation pixel may be sampled multiple times, so intensity values between neighbouring pixels are calculated using a bilinear interpolation. The transformed location can be written where i and j are the integer positions of the centre of the image pixels in the observation, and 0 ≤ δ i , δ j < 1 are the discrepancies between image pixel centres. The intensity, I r,θ , of the (r, θ ) location is given by where I i,j is the intensity value of the (i, j ) image pixel. Recall that in order to compare across a sequence of observations, sampling is taken from the limb-corrected image. When sampling for the r-θ plots, the resolution is chosen to be 1 pixel in the radial direction (i.e. r is chosen to have the same resolution as x and y in the original observation), and 1 • in the angular direction. The full 360 • is sampled in θ , but a range that defines an annulus is chosen in r. A minimum radius of 5 pixels is chosen as this is generally sufficient to be fully contained within the umbra for most sunspots meeting the selection criteria (Section 2), and below this level the image resolution is not sufficient to obtain useful measurements. The upper radius is chosen so that it completely contains the penumbra of the sunspot. A typical value used is 50 pixels, but for larger sunspots (e.g. AR 12158) a larger upper radius is used. This is refined on an observation by observation basis, by defining an annulus about the penumbra based on the penumbral pixels (defined in Section 3.1) at different radial bands. For each radial location in the r-θ image, the total number of angular positions whose intensities are in the penumbral range are calculated. Working from the inner radius outward, the lower boundary of the refined annulus is defined to occur when the proportion of penumbral pixels in a radial band first goes above a quarter, this is allowed to increase above a proportion of a half and the upper boundary is defined to occur when the proportion next drops below a half.
The different proportions at each radial extent are chosen to reflect the different conditions at each location. On the inner umbral/penumbral boundary, a lower threshold is chosen so that fibril structures and partial light-bridges may be included, making use of the rotation of structures in the umbra. On the outer boundary, the higher threshold is chosen to better eliminate features outside of the sunspot that may not be rotating with it, or do not possess a radial structure suitable for tracking.
The refined annulus is just calculated at this stage, but is utilised at a later step in the analysis.
An example of the sunspot in AR 11147 is shown in Figure 5. It shows the non-limbcorrected image of the sunspot which is located away from Sun centre, the reprojected limbcorrected sunspot, and the uncurled sunspot image. Overlaid on the uncurled image is the calculated refined annulus for this observation, which in this case adequately isolates the penumbra.

Tracking Lateral Motions
The rotational motion of the sunspot is determined by tracking the horizontal motion of features in the sunspot penumbra. As the penumbra has a strong radial structure, characterised by light and dark radial fibrils, this can be done by looking at the locations of intensity peaks and troughs at different radii, and how they move in subsequent observations. This is demonstrated in Figure 6, where seven intensity profiles from subsequent observations at a fixed radius are plotted, their peaks and troughs located and matched.
For a given radial position in an observation, the intensity profile with respect to the angular position is extracted. A boxcar smoothing with a width of 5 pixels is applied to Figure 6 Example of how peaks and troughs are tracked for a selection of observations from AR 11147. Bottom, a time-slice for a fixed radial position (r = 22 pixels) over 200 observations. Penumbral motion is picked out in the bright and dark streaks in the time-slice, corresponding to the radial penumbral fibrils in an HMI image. Top, seven radial intensity profiles in subsequent observations are plotted (with a vertical offset between each profile). The peaks and troughs in these intensity profiles are identified, and where possible, matched to peaks and troughs in the subsequent radial profile. reduce noise, where a periodic boundary is applied and the boxcar is allowed to smooth across the 0 • and 360 • ends.
For each angular location, a 5 pixel window centred on the location is analysed. If this has the maximum/minimum intensity value in the window, then it is designated as a turning point. This location is refined by fitting a quadratic to the intensity profile at the turning point and its immediate neighbouring pixels, and finding the turning point of the quadratic. This adjusts the location of the turning point by less than ±0.5 • (see Appendix C).
Once all of the turning points for each radial location have been found for an observation, these are compared to the turning points from the previous observation, and two turning points are matched if all of the following conditions hold: i) The two turning points are of the same type (i.e. match peaks with peaks and troughs with troughs), ii) the time between the two observations is less than 10 minutes, iii) the two turning points are on radial bands of the same radius, iv) the angular positions of the two turning points are within ±3 • of each other.
Note that it is not a requirement that all turning points are matched, and there are physical reasons why this behaviour is desirable. For example, a penumbral fibril may fade over time, and its related peaks/troughs would merge and disappear.
Once all possible turning points for an observation have been matched, an average rotation profile by radial location is calculated. Consider two turning points with angular positions θ n−1 i,r and θ n i,r that have been matched in consecutive observations, where i ranges over all matched angular positions, r gives the radial location, and n the observation number. The amount the turning point has rotated between observations is given by θ n i,r − θ n−1 i,r , and the average rotation of the radial band between the observations, d n r is given by where N tp is the number of matched turning points for the given radial band. If there are no matched turning points for a given position, then the average rotation is taken to be d n r = 0. This has a standard error given by Here, σ is the standard deviation associated with the average rotation, and the first term is the usual formula for the standard error on the mean (e.g. Barlow, 1989). There are also two scale errors, δ 1 and δ 2 , due to underlying sampling scales in the data, these are assumed to be uniformly distributed so the uncertainty contribution comes from the standard deviation of the uniform distribution. The first scale error is from sampling the data into 360 × 1 • bins. This gives a scale of The second scale error comes from sampling from the original observation. The basic scale here is 1 pixel, but this must be converted to degrees. The conversion factor will vary depending on the radial position, a small radial band will produce a small ring occupying a smaller pixel extent, where a large radial band will produce a large ring occupying a larger pixel extent. The conversion factor can be estimated by dividing the angular extent of the ring (360 • ) by the circumference of the ring in pixels (2πr), giving the scale Note that δ 1 < δ 2 when r < 58 pixels. This is the case for many sunspots in SDO/HMI (for example, the sunspot in Figure 5 has a penumbral range from around 15 -30 pixels), and in these cases the second scale error will dominate the first. This provides justification for not sampling at a higher resolution (e.g. 720 × 0.5 • angular bins) when uncurling the sunspot. Between them, these two scale errors account for the uncertainty introduced from mapping the original observed data to the r-θ image. Analysis of all observations provides a set of average rotations, d n r ± S d n r , at all radial locations at all times. To obtain a characteristic rotation for the sunspot at a given time, a weighted average over the radial locations is applied. However, the refined radial annuli calculated in Section 3.3 are applied to restrict the calculation to using penumbral rotations only. To reduce noise in the size of the annuli, a running average of the penumbral radii from frames taken between ±30 minutes from frame number n is calculated (along with its associated standard deviation). While the size of the annulus can vary substantially over the course of the lifetime of a sunspot, it typically does this on a time scale greater than the ±30 minute window being applied.
The weighted average is then given by where r 0 and r 1 are the inner and outer radii of the average annuli for the frame, and the weighting is given by This has the error Finally, a cumulative rotation profile, n over the observed lifetime of the sunspot can be calculated as with cumulative error The first frame in the sequence (n = 0) is assumed to have no rotation ( 0 = 0 • ) and no cumulative error (σ 0 = 0), although for sunspots that have rotated onto the solar disc this history is actually unknown.

Error Due to the Variation of the Refined Penumbral Annuli
The calculation of the refined penumbral annulus and the associated running averages may introduce additional errors as the portion of the sunspot analysed varies. To test for this, a series of cumulative rotation profiles are generated for annuli perturbed from the running average. The perturbations are defined by the given average radii, r n 0 and r n 1 , and their associated standard deviations, σ r n 0 and σ r n 1 , recalling that these values change over successive frames. The following perturbed annuli are defined: (r n 0 − σ r n 0 , r n 1 − σ r n 1 ), (r n 0 − 2σ r n 0 , r n 1 − 2σ r n 1 ), (r n 0 − σ r n 0 , r n 1 + σ r n 1 ), (r n 0 − 2σ r n 0 , r n 1 + 2σ r n 1 ), (r n 0 + σ r n 0 , r n 1 − σ r n 1 ), (r n 0 + 2σ r n 0 , r n 1 − 2σ r n 1 ), (r n 0 + σ r n 0 , r n 1 + σ r n 1 ), (r n 0 + 2σ r n 0 , r n 1 + 2σ r n 1 ), which are interpreted as being a series of 1σ and 2σ deviations from the calculated penumbral annulus. In cases when r n 0 + σ r n 0 > r n 1 − σ r n 1 or r n 0 + 2σ r n 0 > r n 1 − 2σ r n 1 , where a "negative" annulus would be defined, a 2 pixel wide annulus about the midpoint of the overlap is used instead.
The rotation profiles, n p , are calculated as outlined above except the weighted average rotation which is calculated between these alternative penumbral annuli rather than the annuli calculated in Section 3.3. The residuals between the two are given by where the multiplier m p = 1 for the 1σ annuli and m p = 0.5 for the 2σ annuli, casting all of the residuals to be 1σ equivalents. The weighted average of these residuals are then taken to be a characteristic uncertainty due to the variation of the refined penumbral annuli, i.e.
where the weighting is given by w p = 1 for the 1σ annuli and w p = 0.5 for the 2σ annuli, which gives a preferential weighting to the 1σ annuli over the 2σ annuli.

Error Due to the Uncertainty on the Centre of the Sunspot
The calculation of the centre of the sunspot in each observation contains errors that may introduce additional uncertainty to the cumulative rotation profile. To test this a series of rotation profiles are generated using perturbations from the calculated centre based on Gaussian random walks whose deviation reflects the uncertainty on the centre position. A single random walk, (δ n x , δ n y ), about the centre of the sunspot is generated by application of the formulae δ n x = δ n−1 x + N n x (0, 1) and δ n y = δ n−1 y + N n y (0, 1), where N n x (0, 1) and N n y (0, 1) are random numbers drawn from a normal distribution with mean zero and standard deviation one, and the sequence starts at (δ 1 x , δ 1 y ) = (0, 0). This is rescaled to reflect the uncertainty on the centre of the sunspot by first calculating the mean, (δ x , δ y ), and standard deviation, (σ δx , σ δy ), of the random walk, then applying the formulae where S x = 0.6 pixel and S y = 0.6 pixel are the characteristic errors on the centre of the sunspot taken from Section 3.2. An example of such a random walk is shown in Figure 7 for a sequence of 4000 steps (i.e. a sequence of 4000 observations). A cumulative rotation profile can be generated using the overall method described above, but using the perturbed sunspot centre x n pert = x n cen + n x , and y n pert = y n cen + n y , for reprojecting and uncurling the sunspot. This profile is interpreted as being a 1σ deviation from the unperturbed cumulative rotation profile.

Figure 7
Example of a random walk about the sunspot centre. The asterisk marks the sunspot centre, the diamond the start of the random walk, and the square its end. This particular random walk comprises 4000 steps.
A series of N c different cumulative rotation profiles, n c are calculated from a set of N c generated random walks, and the residuals between these and the unperturbed profile are calculated as n c = n − n c .
The characteristic uncertainty due to the error on the centre of the sunspot is then taken to be the average of the residuals over the different random walks, i.e.
For the sunspots analysed in this article this error is calculated using N c = 25 random walks. The overall error on the unperturbed cumulative rotation profile is taken to be the three error components (σ n , σ p n , σ c n ) added in quadrature, i.e.

Results
This section will look at the analysis of the sunspot in AR 11147 in detail, then pick out the key points for the remaining five sunspots analysed.

Introduction to AR 11147
The  An example image of the sunspot, which has been reprojected and limb-corrected, is shown in Figure 8. Intensity profiles along the horizontal and vertical slices through the centre of the sunspot are also shown, and the umbral and penumbral thresholds are displayed on the profiles by the dotted lines. Dashed lines and dot-dot-dot-dashed lines map the umbral and penumbral regions to the main image, demonstrating that these thresholds do pick out the umbra and penumbra in the image. Figure 9a shows how the reprojected area of the umbra varies over time. In this case the area remains well above the soft limit, so all observations in the sequence can be used for the analysis. Figure 9b shows a histogram of the radii of the calculated penumbral annulus for each observation. In this case, the two tight peaks show that the identified umbral and penumbral radii do not change much during the sequence. These two plots demonstrate that there is only a small variation in the size and shape of this sunspot over the course of the analysis.

Cumulative Rotation Profile
The cumulative rotation profile of the sunspot is shown in Figure 10. The plot starts on 17 January 2011 at 00:00 UT, but the sunspot analysis starts 5.7 hours after this as the sunspot  begins its transit across the solar disc. Initially, the sunspot exhibits anti-clockwise (positive) rotation until it reaches a maximum at 66 hours of 34 ± 6 • . At this point the direction of rotation changes to clockwise (negative) rotation until it reaches a minimum at 124 hours of −18.4 ± 8.3 • (a total of 52.4 ± 10.2 • clockwise rotation). After this, there is a mixture of anti-clockwise and clockwise rotation, but at a smaller amount (of the order of ±10 • ), and the sunspot ends its transit across the solar disc (and so the analysis) at 205 hours with a cumulative rotation of −8.9 ± 11.0 • .
While the cumulative rotation of the sunspot ends up being small in this case, the total absolute rotation (i.e. ignoring rotation direction) is around 100 • , so the sunspot has rotated by an appreciable amount, just in both anti-clockwise and clockwise directions.

Uncertainties on the Cumulative Rotation Profile
A breakdown of the total error into its three components is shown in Figure 11a. The overall error for AR 11147 is dominated by the cumulative error (defined in Section 3.4), with smaller contributions from the error due to uncertainty on the sunspot centre (Section 3.6) and the error due to the variation of the penumbral annuli (Section 3.5). However, as the errors are added in quadrature, the contribution of the latter two is slight.
The 1σ and 2σ cumulative rotation profiles for the perturbed penumbral annuli used to evaluate the uncertainty due to the variation of the refined penumbral annuli (Section 3.5) are shown in Figure 11b. These remain close to the unperturbed cumulative rotation profile and demonstrate why this contribution to the overall error is small. This is due to the small change in size of the sunspot over its analysis as demonstrated in Figure 9. Figure 11c shows a histogram of the uncertainties on the location of the centre of the sunspot over each image in the sequence. The errors on x cen (solid histogram) are in the range 0.35 -0.48 pixel, and the errors on y cen (dashed histogram) are in the range 0.36 -0.51 pixel, so the choice of 0.6 pixel as a characteristic uncertainty is reasonable (at worst it is an over-estimate).
The 25 cumulative rotation profiles from the Gaussian random walk perturbed sunspot centres are shown in grey in Figure 11d along with the unperturbed rotation profile in black. These profiles form a reasonably tight bundle, underlining their relatively small contribution to the uncertainty.
Recall that the cumulative error, the largest error contribution, was comprised of three sub-errors, the standard error and two scale errors (Section 3.4). Typical values for these can be calculated. The first scale error is the easiest and gives δ 2 1 /12 = 1/12 ≈ 0.083 degrees 2 . The second scale error varies with the radial band of the sunspot according to From Figure 9, the range for the radii of the sunspot can be estimated as being between 12 and 33 pixels. For r = 12 pixels this gives δ 2 2 /12 ≈ 1.9 degrees 2 , and for r = 33 pixels this gives δ 2 2 /12 ≈ 0.25 degrees 2 . So the uncertainty contribution from the second scale error ranges between 0.25 -1.9 degrees 2 .
A typical value for the standard error can be estimated by assuming that as the amount an individual tracked feature is allowed to vary is restricted to a maximum of ±3 • , then the standard deviation must be around 1.5 • . The number of values averaged in a given radial band will vary, but counting peaks and troughs for the mid-range radius shown in Figure 6 suggests a value N tp > 30. Lower radial bands will likely have a lower value of N tp due to a lower resolution, so taking the estimate N tp = 20 gives a typical standard error of σ 2 /N tp = 1.5 2 /20 ≈ 0.11 degrees 2 .
From these values it is clear that the second scale error is the dominant source of error. To evaluate its contribution, the second scale error is propagated through the analysis in isolation, i.e. take = 2700 π 2 r 2 .
When averaging over all radial bands this has the weighting w n r = 1 S 2 d n r = π 2 r 2 2700 (32) and the uncertainty Note that this does not vary with n, the observation number.
The error on the cumulative rotation profile is given by as σ 0 = 0. So the error on the cumulative rotation profile due to the second scale error is given by σ n = 0.151 √ n, which at n = 3972 (the final observation) is σ n = 9.5 • , compared to the calculated total error of σ n = 11 • . So the dominant source of error for this sunspot is from the second scale error, which itself derives from the underlying resolution of the HMI observations.

Analysis of Other Sunspots
The rotation profiles for the remaining five sunspots listed in Table 1 are shown in Figure 12.
The cumulative rotation profile of the sunspot in AR 11770 is shown in Figure 12a. The sunspot is smaller than AR 11147 and its size remains stable until the last 20 h of its lifetime when it begins to shrink to below the size limits. This sunspot was chosen for analysis as it does not exhibit any significant rotation, which is evident in the cumulative rotation profile. Initially, the sunspot undergoes clockwise (negative) rotation, until it reaches a maximum of −13.5 • at around 79 h when the rotation becomes anti-clockwise (positive) until it reaches a maximum of 8.7 • at around 183 h. There is a slight reduction in the rotation to 6 • by the end of the period of analysis. However, the associated error is comparatively large (reaching ±16 • by the end of the period of analysis), so it would be difficult to class this measured rotation as being significant. The error on the rotation is dominated by the cumulative error, which is dominated by the second scale error. This has a larger effect than for AR 11147 due to the smaller size of the sunspot. Figure 12b shows the cumulative rotation profile of the sunspot in AR 12158. This is a significantly larger sunspot than the previous cases (with an umbral area 4 -6 times bigger), and also demonstrates a greater change in size as it undergoes fragmentation in the second Figure 13 A sequence of reprojected, limb-corrected images of the sunspot in AR 12422 at different times during its transit across the solar disc. half of its time span. Additionally, there are a number of pores located nearby that potentially could (but actually do not) affect the analysis.
The sunspot exhibits a large anti-clockwise rotation of 205 • over the course of the 220 h. The associated error is ±10 • at this time, and while relatively low compared to the amount of rotation, it is comparable in size to the uncertainties found in the previous two cases. The largest contribution to the uncertainties is the cumulative error, but the error due to variation of the penumbral annuli is also significant.
During the period from 100 -140 h the sunspot becomes non-circular prior to fragmentation and the measured rotation slows to being near standstill. For this case, this is reasonable as the elongated sunspot shape does not exhibit any obvious rotation during the fragmentation that occurs at this time. However, the circular annulus-based analysis will be contaminated with a higher proportion of non-penumbral pixels (both umbral and photospheric), so extra scrutiny is needed. Figure 12c shows the cumulative rotation profile of the sunspot in AR 12422. The sunspot exhibits clockwise rotation, but it is clear that the approach has encountered a problem at around 75 -120 h. The reasons for this can be seen in Figure 13, where the sunspot which is initially reasonably circular (panel a), begins to elongate (panels b-c), before forming a kink and a light bridge (panels d-e). This highlights the process of the initial sunspot fragmenting into two equally sized sunspots, but the two sunspots do not completely separate and continue to share the same penumbra (panel f, the feature to the right of the image forms separately and is not one of the fragmented sunspots). The problems occur in panels d and e, where the calculated centre of the sunspot occurs at the kink and light bridge, causing unreliable calculation of the location of the centre and the penumbral radii. This leads to erratic calculation of the rotation and its associated uncertainty quickly increases.
The measured rotation settles down again at 120 h and the increase in the error becomes low again as the sunspot returns to a more circular shape. While this example is presented as a case where the analysis "does not work", a profile is found and the large increase in the associated uncertainties caused by the problems encountered may be considered a reasonable output that indicates further follow up is required.
The cumulative rotation profile of the two sunspots in AR 12420 are shown in Figure 12d. Both sunspots are small and exhibit brief drops to sizes below the soft size-limit. However, they remain close together and are included to investigate the effect that neighbouring sunspots may have on the analysis. In this case, neither sunspot noticeably effects the analysis of the other.
The leading sunspot (sunspot A) exhibits clockwise rotation from about 40 h and reaches −30 • by 130 h. The trailing sunspot (sunspot B) exhibits predominantly anti-clockwise rotation from 10 h, and has rotated 30 • by 180 h, though a little clockwise rotation between 75 -100 h is detected. Again, the predominant source of error is the cumulative error (through the second scale error), though contributions from the error due to the penumbral radii and the error on the centre of the sunspot are present (due to formation of light bridges changing the size and shape of the sunspots).
For all of the sunspots tested, the characteristic error of 0.6 pixels for the location of the centre of the sunspot remains a good choice. The measured errors (in only one of the coordinates) only creep above this for periods of AR 12158 and AR 12422 when the sunspot becomes highly elliptical prior to fragmentation, representing less than 3% and 6% of observations for each sequence, with outlying errors of 0.7 pixels and 0.64 pixels in each case.

Comparing Measured Rotations to the Data
The rotation profiles generated above show that the methodology described is measuring something, but it has not yet been demonstrated that it is doing a good job of measuring the rotation of the sunspot itself. To do this, animations of each test sunspot are produced. In each, observations of the sunspot are corrected for limb darkening and reprojected to disc centre. An annulus is plotted on top of the observation with radii given by the penumbral annulus radii calculated and used within the analysis. Spokes are added to the annulus to give it orientation, and the annulus is rotated by the cumulative rotation calculated by the analysis.
Viewing the animation allows the measured rotation to be compared by eye to the underlying sunspot. This is carried out for each sunspot analysed. While it is evident that some parts of the penumbra rotate slightly slower/faster than the calculated rotation, the technique outlined in this article does an acceptable job of capturing the bulk rotation of the sunspot, within the calculated uncertainties.
The animations for each sunspot are provided in the supplementary material.

Comparison with Literature Values
It is a useful test to compare results from this approach with measured rotations from the literature. A number of authors have reported on sunspot rotations from AR 11158 which was observed on the Sun from 13 February 2011 till 18 February 2011. The active region was comprised of several sunspots, of which three have rotation measurements reported in the literature, some of which are independently measured by different authors. The method outlined in this article is used to analyse the same sunspots for the same duration as the reference articles.
The comparison results are summarised in Table 2. Some of the sunspots have rotations reported by multiple authors, the number in brackets in the label indicates this. At first it seems that the values calculated by this method are under-reporting the rotation compared to other authors, but further analysis of the approach taken by other authors indicates why this is the case, and why an approach such the one described in this article is required.
For example, for the sunspot studied by Jiang et al. (2012) only three penumbral features are identified and tracked, and all three are located in the same half of the sunspot. A similar approach to uncurling the sunspot is taken, but only one radial location is used to calculate the rotation. The range of angular velocities found suggest a sunspot that may not be rotating uniformly (in this case the variation is in the axial direction), and the larger rotations reported may be the result of selection effects. Li and Liu (2015) analyse one of the sunspots in AR 11158 for 4 days using solid-body cross-correlation on uncurled images. They are able to track the sunspot for longer than the method in this article due to the selection criteria outlined in Section 2. This accounts for around 25 • of the discrepancy. Both methods are in agreement that this sunspot undergoes a large amount of rotation. Vemareddy, Ambastha, and Maurya (2012) analyse each of the sunspots studied by Jiang et al. (2012) and Li and Liu (2015) using an approach similar to the one described in this article, except for each sunspot they track only one penumbral feature at one radial position. It is clear from Figure 3 in their paper that they have selected one of the largest rotating features, and there are clearly other features that rotate significantly less. As such, the rotation measured may be an over-estimate. For example, the rotation of SP1 compares to the largest rotating feature for P2 (the same sunspot) in Jiang et al. (2012), although the time-periods do not perfectly overlap. Wang et al. (2014) take a different approach and use the DAVE method (Schuck, 2006) to calculate the rotation speeds in the umbra of the sunspots analysed. They look at two sunspots in a continuous period of 1 h 30 m which is split into 3 sections, before, during, and after the associated flare. The before and after periods report average rotation speeds, where the during period is a peak rotation speed. For sunspot f, they just report that the typical before and after rotation rates are similar to p. Figure 4 in their article shows the variation of rotation speeds during the analysis period, and exhibits a large degree of fluctuation (of the order of 10 -15 • h −1 in places), which suggests that there may be a large uncertainty on this measurement. This might be expected in small areas, such as a sunspot umbra, where small errors on the calculated velocities and location of the centre of the sunspot may have a disproportionate effect on the rotation rate. The large spike in rotation speed detected during the flare may occur due to the magnetic nature of the flare affecting the magnetogram data used by the DAVE method, but such spikes are not seen in the penumbral measurements calculated here or by Jiang et al. (2012) or Vemareddy, Ambastha, and Maurya (2012).
These examples highlight several problems and why the approach described here is needed. Jiang et al. (2012) and Vemareddy, Ambastha, and Maurya (2012) demonstrate the issue with (i) selecting only a small number of penumbral features at a single radial location, (ii) manual selection of parameters (such as radial location and particular feature), and how such choices can lead to bias in the reported rotation. Additionally, features that are long in the penumbral direction may not respond to sunspot rotations uniformly but may have rotations that vary with radius. By using all radial positions and a large number of penumbral features at all axial locations, this method is not biased towards the clearest and largest rotating features and produces rotations that are balanced by smaller rotating features. Additionally, features that are long in the radial direction will be identified at different radial positions and any variations in rotation are incorporated into the average rotation produced. In obtaining measurement with large fluctuations, Wang et al. (2014) demonstrate the need to evaluate the uncertainty of the measured rotation. One aspect of this is the degree of spread of the measurement that can be estimated using a standard deviation, but the method presented here goes beyond this to include effects which may have a greater contribution such as the resolution of the data and the uncertainty on calculated properties such as the centre of the sunspot.

Discussion and Conclusions
This article has described a method to calculate the rotation undergone by a sunspot about its umbral centre over the course of a transit across the solar disc. This is achieved by tracking penumbral fibrils and other radial structures seen in continuum observations of the photosphere. The definition of the penumbra and identification of features within it for tracking is designed to be both autonomous to avoid bias due to selection effects and wide ranging to ensure that as many features as possible across the full range of the penumbra are selected to get a full picture of the rotation. In addition to calculating a rotation profile for a given sunspot over the course of the transit, uncertainties on the rotation are also calculated and broken down to help identify the key sources of uncertainty.
The method has been tested on six different sunspots to evaluate its performance in different cases. The approach performed well for five of the sunspots, demonstrating robustness to rotation, change in size, interaction with smaller pores, and interference from nearby sunspots. The only case where the method did not "perform well" was for AR 12422 where two similar-sized sunspots share the same penumbra, and deformation to the perceived umbral shape leads to inadequate definition of the sunspot. However, even in this case, the method produced a rotation profile which registered a greatly enhanced uncertainty to the ill-identified sunspot geometry which would in practice lead to further investigation, so labelling this a "failure" would be an overstatement.
The measured rotation profile for each sunspot is compared by eye to observations through the generation of an animation of the sunspot with a wire-frame annulus super-imposed onto the observation which rotates in line with the measured profile. These agree well with the bulk rotation of the sunspots, though it is evident that sunspots do exhibit some degree of non-uniform rotation as well.
The six sunspots analysed in this article were chosen to test the robustness of the approach to different conditions, and are not a representative sample. However, some points can be drawn from these. First, even non-rotating sunspots seem to exhibit some small rotation (though this may be below the measured uncertainty), and it may be more accurate to consider all sunspots as rotating sunspots, but many may not rotate very much. However, a larger sample would be needed to better understand this distribution. Second, the uncertainties associated with the method are of the order 10 -20 • , with smaller sunspots typically having larger uncertainties than larger ones. This is due to the main source of error being the resolution of the instrument, and smaller sunspots are more susceptible to this. Adapting this method to analyse higher-resolution data from, e.g. ground-based observatories, may lead to a lower uncertainty. However, this will be limited as the contribution due to the scale error from image resolution falls below that of other sources of error-such as the deviation of the non-uniform rotation from the bulk rotation.
As a final test, the rotation profiles for sunspots contained in AR 11158 were calculated and compared to literature values, where four different authors have reported sunspot rotations for this region. The calculated rotations are in some agreement with the literature (in that rotation directions and scales are identified), but calculated rotations are noticeably lower than the literature values. Analysis of the literature methods suggests that selection effects could play a large role in the discrepancies.
The method presented here has a number of key advantages over existing approaches: (i) it is robust in its choice of operating parameters, either determining them automatically, or where that is not possible using well investigated choices tested across a range of sunspots for which an uncertainty can be attributed (for example, the choice of umbral/penumbral pixel values feed into the refining of the penumbral annuli which have a specific uncertainty calculation); (ii) it aims to identify and track as many features as possible distributed about the penumbra of a sunspot in order to reduce bias caused by the effects of limited and/or manual selection of features; (iii) it is designed to measure the rotation profile of sunspots for near-complete transits across the solar disc (though may be used for shorter sequences of observations), and it calculates a detailed uncertainty for the profile that can be used to clearly identify the dominant source of error; and (iv) its semi-automatic approach allows for easier extraction of sunspot rotation profiles which enables large statistical samples to be constructed. Such statistical samples will play a key role in addressing questions regarding the nature of sunspot rotation such as: why sunspots rotate and what proportion of sunspots rotate and by how much, and how the sunspot rotation relates to solar activity.
The first planned application of this method will be to build a sample of sunspot rotations in active regions that produce an X-class flare to determine whether such regions exhibit enhanced sunspot rotation and to estimate whether the energy injected into the region by the sunspot rotation is sufficient to provide the energy released by the solar flare. The semiautomated nature of the method and the full transit history it will provide in the run up to the flares will be critical in building and interpreting this sample.

Open Data
The rotation profiles generated for the test active regions in this article are available for further scientific use. FITS files for each active region, along with documentation, can be found under DOI: 10.17030/uclan.data.00000214.
where (x, y) is the two-dimensional location of the image pixel being corrected, and cos (θ ) = 1 − x 2 + y 2 r 2 s .
Again, r s is the radius of the Sun in the same units as x and y. An example of how the limb correction compares to uncorrected data can be seen in Figure 14. This shows a complete observation containing AR 11147 in its uncorrected and corrected form, along with an intensity cut along the solar equator. It can be seen that the limb correction does a good job of removing the large-scale intensity variation due to limb darkening, while retaining the small-scale variation due to photospheric features.

Appendix B: Errors in the Centre of the Sunspot
In order to test the effects of sunspot size and shape on the errors on the centre of a sunspot, a series of idealised elliptical sunspots are considered. Parameters r > 0 and 0 < μ ≤ 1 are used to define the size and degree of eccentricity of the ellipse. The semi-major axis of the ellipse is given by a = r/ √ μ, and the semi-minor axis is b = r √ μ. These are chosen such that the area of the ellipse A = πab = πr 2 is comparable for like values of r, and the eccentricity is given by e = 1 − μ 2 . The ellipse is projected onto a pixellated grid such that the centre of the ellipse is located at the middle of the centre pixel, and a pixel is taken to be part of the idealised sunspot if its centre is contained within the ellipse. The centre of the idealised sunspot and its uncertainties can be found using the method described in Section 3.2. Figure 15a shows the variation of the errors on the centre for varying r, i.e. how the errors change with the size of the sunspot. The solid line shows the

Figure 16
Illustration of how a quadratic fit to the three neighbouring pixels can be used to refine the location of the turning point. error in the direction of the semi-major axis, and the dashed line in the direction of the semiminor axis. The dotted line shows the characteristic error of S x = S y = 0.6. Above r = 6 (the hard-limit for the minimum sunspot size in Section 2) the errors exhibit little variation, and the large sunspots analysed in this article have maximum umbral areas consistent with r = 25. A more notable variation is with parameter μ. This is explored in Figure 15b which shows the variation of the errors with parameter μ. Here, the errors over the range 6 ≤ r ≤ 50 are averaged for both the semi-major and the semi-minor directions, and plotted with a shaded one standard deviation zone. The dashed line shows the characteristic error of S x = S y = 0.6, and the errors can be seen to fall below this for all but the most elliptic sunspots.
This supports the choice of the characteristic error of S x = S y = 0.6 used in this article.

Appendix C: Adjusting the Location of a Turning Point
Once a turning point is located, a quadratic is fitted to the turning point and its two neighbouring pixels to refine the location of the turning point. This is illustrated in Figure 16, where the two large intensity values at locations θ = 0 and θ = 1, in comparison to location θ = −1, suggest that the peak should be closer to the θ = 1 location than the θ = −1 location.
The function fitted is The turning point of the quadratic occurs when its derivative is zero, i.e.
or when This gives the correction to the angular location required to refine the position of the turning point.