Rotation

Overview

Due to the Earth's rotation as a celestial object, such as the moon or the sun, moves across the sky it appears to rotate. From the perspective of an observer on Earth the rotation is perpendicular to the line of sight. In other words, the same side always faces the observer, but it may appear sideways relative to how it appeared in the past. This page is about such rotation, which is known as "field rotation".

The field rotation depends the latitude of the observer and the amount of time before or after the celestial object's meridian passing. The meridian passing is when the celestial object passes the meridian and is then at its highest point in the sky. For the sun the meridian passing is known as "solar noon". This page will refer to the meridian passing as "MP".

This page uses mathematical brute force with the help of SageMath (a symbolic manipulator) to compute the amount of field rotation as a function of time. The result is a long formula (see Calculations). As can be seen on this RASC page there is a much simpler formula for the rate of field change (degrees / hour):

rate =  K * cos(az) / cos(alt)

But that's hard to integrate since "az" and "alt" are a function of time as well. The RASC page has a lot of other helpful information about field rotation.

Field rotation is avoided with equatorial mounts.

Coordinate System

Humans, cameras on tripods and altazimuth mounted telescopes all observe the same rotation due to the fact that they each have the same two axes, one vertical one one horizontal and perpendicular to the line of sight. To calculate the field rotation consider an observer in the Norther hemisphere facing South. From this point forward the observer is assumed to be an altazimuth mounted camera. The camera is at the origin of a a right handed coordinate system with the following axes:

Axis Description
X East, to the right of the camera.
Y Up, to the zenith of the camera.
Z North, to the rear of the camera.

This means that the camera sees celestial objects rise in the East, move in a positive X direction to the right, and then set in the West.

Variables

By convention Greek letters are often used for some of the variables, but since Greek letters are hard to type (at least for me) and hard to use with mathematics software such as SageMath ordinary ASCII letters will be used. Scalars (single value) variables are lower case. Vectors are italic. Matrices are UPPER case. Suffix "d" is used for right to avoid "or" (a Python keyword). It makes sense in French and Spanish:

Variable Greek Description
a Altitude of target at MP from the Southern horizon
cd Right of the camera
cu Up from the camera
h Hours since (positive) or before (negative) MP
m Minutes since (positive) or before (negative) MP
o Location of the celestial object
l λ Longitude
od To the right of the celestial object
ou Up from the celestial object
p φ Latitude
R General purpose rotation matrix
Ru Rotation matrix that rotates about u
s Seconds since (positive) or before (negative) MP
t θ Rotation about the celestial axis since (positive) or before (negative) MP
u Earth's axis as a unit length rotation vector
w Clockwise field rotation



Calculations

This is an overview of the field rotation calculations. A complete list of SageMath commands, including the examples, can be found here.

When the celestial object is at the meridian passing it will be considered to be right side up and therefore to have no field rotation. Unit vector o points to the celestial object:

oMP  = vector([0, sin(a), -cos(a)])

To track the rotation right pointing vector od will be attached to the celestial object. At the meridian passing the vector is:

odMP = vector([1, 0, 0])

Starting with general purpose rotation matrix R:

R = Matrix([
  [ -ux^2*(cos(t) - 1) + cos(t)   , -ux*uy*(cos(t) - 1) - uz*sin(t), -ux*uz*(cos(t) - 1) + uy*sin(t)],
  [-ux*uy*(cos(t) - 1) + uz*sin(t),  -uy^2*(cos(t) - 1) +    cos(t), -uy*uz*(cos(t) - 1) - ux*sin(t)],
  [-ux*uz*(cos(t) - 1) - uy*sin(t), -uy*uz*(cos(t) - 1) + ux*sin(t), -uz^2 *(cos(t) - 1) +    cos(t)] ])

Rotation matrix Ru is produced by substituting relevant values for unit vector u:

Ru = R.subs({ux:0, uy:-sin(p), uz:-cos(p)})

Now that rotation matrix Ru is known it can be used to determine the od at times other than the MP:

od = Ru * odMP

The cd vector is determined by first projecting vector o into the XZ-plane (discard Y) which gives us a vector that is parallel to the ground and that is pointing toward the front of the camera. Since cd is right angles to that vector it's possible to get it by swapping the X and Z coordinates and making one of them negative. Finally cd should be normalized so it has a length of one:

cd = vector([-o[2], 0, o[0]]).normalized()

Since field rotation is perpendicular to the line of sight (to o) and since vectors od and cd lie in that perpendicular plane the formula for finding the angle between two unit vectors can be used. For two unit vectors this is acos(v1 . v2). In order to determine the sign the field rotation (w) will positive when od is below the horizon (when od[1) is less than 1). This accounts for the unit_step() part. Finally the just less that one "1 - 1e-16" factor is to avoid errors when, due to floating point limitations, the value passed to acos() is out of [-1, 1] range:

w = 2 * (unit_step(-od[1]) - 1/2) * acos((1 - 1e-16) * od.dot_product(cd))

Running the above through SageMath produces:

w = (2*unit_step(cos(p)*sin(t)) - 1)*arccos(1.00000000000000*(cos(p)*sin(a)*sin(t) + cos(a)*sin(p)*sin(t))*sin(p)*sin(t)/sqrt(abs((cos(t) - 1)*cos(p)*sin(a)*sin(p) - ((cos(t) - 1)*cos(p)^2 - cos(t))*cos(a))^2 + abs(cos(p)*sin(a)*sin(t) + cos(a)*sin(p)*sin(t))^2) + 1.00000000000000*((cos(t) - 1)*cos(p)*sin(a)*sin(p) - ((cos(t) - 1)*cos(p)^2 - cos(t))*cos(a))*cos(t)/sqrt(abs((cos(t) - 1)*cos(p)*sin(a)*sin(p) - ((cos(t) - 1)*cos(p)^2 - cos(t))*cos(a))^2 + abs(cos(p)*sin(a)*sin(t) + cos(a)*sin(p)*sin(t))^2))

It's possible to make the above a function of the amount time before (negative) or after (positive) MP by replacing angle t with the correct numbers of hours (h), minutes (m) or seconds (s):

wh = w.subs({t:15*h*(pi/180)})
wm = wh.subs({h:m/60})
ws = wm.subs({m:s/60})

Expanded this looks like:

wh (w as a function of hours (h)):

wh = (2*unit_step(cos(p)*sin(1/12*pi*h)) - 1)*arccos(1.00000000000000*(cos(p)*sin(1/12*pi*h)*sin(a) + cos(a)*sin(1/12*pi*h)*sin(p))*sin(1/12*pi*h)*sin(p)/sqrt(abs((cos(1/12*pi*h) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/12*pi*h) - 1)*cos(p)^2 - cos(1/12*pi*h))*cos(a))^2 + abs(cos(p)*sin(1/12*pi*h)*sin(a) + cos(a)*sin(1/12*pi*h)*sin(p))^2) + 1.00000000000000*((cos(1/12*pi*h) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/12*pi*h) - 1)*cos(p)^2 - cos(1/12*pi*h))*cos(a))*cos(1/12*pi*h)/sqrt(abs((cos(1/12*pi*h) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/12*pi*h) - 1)*cos(p)^2 - cos(1/12*pi*h))*cos(a))^2 + abs(cos(p)*sin(1/12*pi*h)*sin(a) + cos(a)*sin(1/12*pi*h)*sin(p))^2))

wm (w as a function of minutes (m)):

wm = (2*unit_step(cos(p)*sin(1/720*pi*m)) - 1)*arccos(1.00000000000000*(cos(p)*sin(1/720*pi*m)*sin(a) + cos(a)*sin(1/720*pi*m)*sin(p))*sin(1/720*pi*m)*sin(p)/sqrt(abs((cos(1/720*pi*m) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/720*pi*m) - 1)*cos(p)^2 - cos(1/720*pi*m))*cos(a))^2 + abs(cos(p)*sin(1/720*pi*m)*sin(a) + cos(a)*sin(1/720*pi*m)*sin(p))^2) + 1.00000000000000*((cos(1/720*pi*m) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/720*pi*m) - 1)*cos(p)^2 - cos(1/720*pi*m))*cos(a))*cos(1/720*pi*m)/sqrt(abs((cos(1/720*pi*m) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/720*pi*m) - 1)*cos(p)^2 - cos(1/720*pi*m))*cos(a))^2 + abs(cos(p)*sin(1/720*pi*m)*sin(a) + cos(a)*sin(1/720*pi*m)*sin(p))^2))

ws (w as a function of seconds (s)):

ws = (2*unit_step(cos(p)*sin(1/43200*pi*s)) - 1)*arccos(1.00000000000000*(cos(p)*sin(1/43200*pi*s)*sin(a) + cos(a)*sin(1/43200*pi*s)*sin(p))*sin(1/43200*pi*s)*sin(p)/sqrt(abs((cos(1/43200*pi*s) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/43200*pi*s) - 1)*cos(p)^2 - cos(1/43200*pi*s))*cos(a))^2 + abs(cos(p)*sin(1/43200*pi*s)*sin(a) + cos(a)*sin(1/43200*pi*s)*sin(p))^2) + 1.00000000000000*((cos(1/43200*pi*s) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/43200*pi*s) - 1)*cos(p)^2 - cos(1/43200*pi*s))*cos(a))*cos(1/43200*pi*s)/sqrt(abs((cos(1/43200*pi*s) - 1)*cos(p)*sin(a)*sin(p) - ((cos(1/43200*pi*s) - 1)*cos(p)^2 - cos(1/43200*pi*s))*cos(a))^2 + abs(cos(p)*sin(1/43200*pi*s)*sin(a) + cos(a)*sin(1/43200*pi*s)*sin(p))^2))



Test

This section demonstrates that the formula produces reasonable results. To see it applied to actual data skip to the examples section.

It's possible to test the formulas generated with certain values where the amount of field rotation can be determined without calculation. For example, during the equinox the angle the sun makes with the zenith is equal to the latitude (p). In Austin the latitude is 0.5283 radian, so at solar noon (at MP):

w.subs({t:0, p:0.5283, a:pi/2-0.5283}).n()
2.10734242554470e-8

which is very close to the expected value of zero which makes sense since we defined the field rotation to be zero at t=0.

On the equinox the sun rises at solar 6:00 am. In this case we'd expect the field rotation to be determined by the angle the celestial equator makes with the horizon, which happens to be the altitude at MP:

w.subs({t:-pi/2, p:0.5283, a:pi/2-0.5283}).n()
-1.04249632679490

which is close to the expected value of -(pi/2-0.5283) (opposite of a since it's before MP).

The formula is in radians, but it's possible to convert to degrees by multiplying by 180/pi. In this case we find that the it's -59.73°. It's negative because the sun starts rotated counter clockwise at sunrise relative to what it will be at solar noon.

For the same example at sunset we use a positive value for t:

w.subs({t:pi/2, p:0.5283, a:pi/2-0.5283}).n()
1.04249632679490

In Sydney Australia during the equinox a negative value is used for the latitude (p) which is -0.5911. Also, the altitude (a) is measured from the Southern horizon, so we subtract from pi if we're given the altitude from the closest point on the horizon. At solar noon:

w.subs({t:0, p:-0.5911, a:pi-(pi/2-0.5911)}).n()
3.14159263251637

It's rotated 180° since that's the way it's seen at solar noon (at MP) in Australia. Being upright at solar noon when viewed from the North is arbitrary. Australia sees it the other way.

For an object close to the North celestial pole when viewed in Austin Texas we'd expect the field rotation (w) to be almost equal to the rotation about the celestial pole (t). So let's see what values we get for each rotation of pi/2:

w.subs({t:0, p:0.5283, a:pi-(0.5283+0.1*pi/180)}).n()
3.14159263868863 # pi
w.subs({t:-pi/2, p:0.5283, a:pi-(0.5283+0.1*pi/180)}).n()
-1.56977769124850 # close to - pi/2
w.subs({t:pi/2, p:0.5283, a:pi-(0.5283+0.1*pi/180)}).n()
1.56977769124850 # close to pi/2
w.subs({t:pi, p:0.5283, a:pi-(0.5283+0.1*pi/180)}).n()
1.49011611938477e-8 # effectively zero



Examples

This section has examples of the field rotation (w) calculation applied to actual data.

Mathias Kp Sun Spot Video

Update 2018-11-30 - if I were to do this over I would avoid a lot of manual calculations by converting each image of the sun with the sunspot to an equirectangular map using sphere2equirect.py in image-utils. Doing so would reveal the coordinates easily.

YouTuber Mathias Kp created video "The visual size of the Sun during a day". Although the video is about the size of the sun it also recorded a sun spot along with other data that makes it possible to see if the field rotation works in the expected way.

The video's description includes a flickr link to pictures of the sun taken at various times. The sunspot locations will extracted from this image. The download link allows the original image size to be downloaded which results in a 1272x743 image, which is what was used here. Here's the sides of each sun image as well as the coordinates of the sunspot in each sun image based on image coordinates (the "I" prefix):

Hour Sun Left IX Sun Right IX Sun Top IY Sun Bottom IY Sunspot IX Sunspot IY
8 am 170 416 166 382 334 229
10 am 568 814 153 398 744 232
12 pm 992 1238 152 396 1179 249
2 pm 171 415 468 711 363 582
4 pm 568 814 468 711 762 602
5 pm 992 1238 481 699 1186 608

It's possible to scale the sunspot coordinates relative to the sun that contains them so that they fall in a [-1, 1] range. Note that the sign is flipped for the Y coordinate so that Y is up.

Sunspot_X = 2 * (Sunspot_IX - Sun_Left_IX) / (Sun_Right_IX - Sun_Left_IX) - 1
Sunspot_Y = 2 * (Sunspot_IY - Sun_Top_IY)  / (Sun_Bottom_IY - Sun_Top_IY) - 1

The "Hour" column is refers to local clock time, but we need the offset relative to solar noon. Fortunately the image also provides the latitude (56.54N) and longitude (8.32E) coordinates as well as the date (February 13th 2016). With those coordinates it's possible to construct a timeanddate.com URL for the sun at that location during that month. From that page we can see that on the 13th solar noon was at 12:40 pm with an altitude of 20° or 0.3491 radians. This means to get we can subtract 40 minutes from the "Hour" value to get the time relative to solar noon, which will be the "Solar Time" column.

With the above the sunspot data is now:

Hour Solar Time Solar Minute Sunspot X Sunspot Y
8 am -4:40 -280 0.3333 0.4167
10 am -2:40 -160 0.4309 0.3551
12 pm -0:40 -40 0.5203 0.2049
2 pm 1:20 80 0.5738 0.06173
4 pm 3:20 200 0.5772 -0.1029
5 pm 4:20 260 0.5772 -0.1651

To plot the sunspots the following creates a reusable function plot_sunspots() in sage:

# Create a scatter splot of the sunspots given their labels, locations and a
# filename.
def plot_sunspots(labels, locs, radius, fname, **kwargs):
    # A graphics object with all the graphs combined.
    g = Graphics()
 
    # A yellow circle for the sun.
    sun_plot = polar_plot(0, 0, 2*pi, fill=1, fillcolor='yellow')
    g += sun_plot
 
    # A circle with a radius equal to the average radius of the sun spots.
    average_circle = circle((0, 0), average_radius, color='green')
    g += average_circle
 
    # The sunspots as small black circles.
    sunspots_plot = scatter_plot(locs, markersize=18, facecolor='black')
    g += sunspots_plot
 
    # The ranges are needed in order to determine the scale for the labels.
    xmin = kwargs["xmin"] if "xmin" in kwargs else -1.0
    xmax = kwargs["xmax"] if "xmax" in kwargs else  1.0
    ymin = kwargs["ymin"] if "ymin" in kwargs else -1.0
    ymax = kwargs["ymax"] if "ymax" in kwargs else  1.0
    xscale= (xmax - xmin) / 2
    yscale= (ymax - ymin) / 2
 
    # Add the label for each sunspot.
    if labels:
        for i in range(len(locs)):
            loc = locs[i]
            g += text(labels[i], (loc[0] + .13 * xscale, loc[1] + .03 * yscale))
 
    # Save it. Seems to be 470x470 by default.
    g.save(fname, **kwargs)

The plots will use the following variables:

# The sunspot hours taken from the table. This will be used to label the
# sunspots.
hours = ["8 am", "10 am", "12 pm", "2 pm", "4 pm", "5 pm"]
 
# The sunspot locations. This is calculated, but shown here for convenience.
sunspot_locs = [(0.3333, 0.4167),  (0.4309,  0.3551), (0.5203,  0.2049),
                (0.5738, 0.06173), (0.5772, -0.1029), (0.5772, -0.1651)]

A plot for the sunspot locations:

plot_sunspots(hours, sunspot_locs, "sunspots.png")

which produces:

The green circle is the average sunspot radius in order to illustrate that the distance from the center of the sun has changed, which suggests that something other than just field rotation is happening. It is likely that the cause is solar rotation, which be explored later.

The filed rotation is a function of the solar minute was well other variables using the wm formula found above:

solar_minutes = [-280, -160, -40, 80, 200, 260]
field_rotations = [wm.subs({m:sm, p:56.54*(pi/180), a:20*(pi/180)}).n(digits=4)
                   for sm in solar_minutes]
# field_rotations = [-0.5447, -0.3715, -0.1016, 0.1997, 0.4414, 0.5236]

The next step is apply the inverse of the field rotation to each sunspot location in order to remove the field rotation, and see they would've looked like at solar noon. To do this we'll use a rotation matrix where the amount of rotation is w, but in the opposite direction. Fortunately this 2D matrix is much simpler than the 3D one seen previously. The rotation matrix is then applied to each sunspot location.

Rw = Matrix([
  [cos(w), -sin(w)], 
  [sin(w),  cos(w)] ])
 
# The sunspot locations with the effect of field rotation removed.
sunspot_locs_w = []
for i in range(len(sunspot_locs)):
    loc = sunspot_locs[i]
    fr = field_rotations[i]
    wss = Rw.subs({w:fr}) * vector(loc)
    sunspot_locs_w.append(wss)
# sunspot_locs_w = [(0.5010, 0.1837), (0.5304, 0.1745), (0.5384, 0.1511),
#                   (0.5502, 0.1743), (0.5658, 0.1536), (0.5824, 0.1456)]

An updated table with the effect of field rotation removed:

Hour Solar Time Solar Minute Field Rotation Sunspot WX Sunspot WY
8 am -4:40 -280 -0.5447 0.5010 0.1837
10 am -2:40 -160 -0.3715 0.5304 0.1745
12 pm -0:40 -40 -0.1016 0.5384 0.1511
2 pm 1:20 80 0.1997 0.5502 0.1743
4 pm 3:20 200 0.4414 0.5658 0.1536
5 pm 4:20 260 0.5236 0.5824 0.1456

A plot for the sunspot locations from the above table:

plot_sunspots(None, sunspot_locs_w, average_radius, "sunspots-w.png")

which produces:

Solar Rotation

The rest of this example attempts to account for solar rotation, which is a smaller contribution to the location of the sunspot. If you're only interested in field rotation then you can skip ahead to the Conclusion.

It's possible to zoom into the sunspot locations from the previous graph:

plot_sunspots(hours, sunspot_locs_w, average_radius, "sunspots-w-zoom.png",
    xmin=0.5010, xmax=0.5824, ymin=0.1456, ymax=0.1837)

which produces:

The sunspot locations trend mostly to the right and a bit down, which makes sense since the sun spots were recorded after the Winter solstice, so the suns's axis, which is almost at right angles to the orbital plane, should be tilted that way. As to the amount and direction of the sun's rotation we'll simplify things by just considering the first (upper left and 8 am) and last (lower right at 5 pm) sunspot locations. Dividing by a the elapsed time produces sunspot_v, a vector that describes the per hour movement on the XY-plane as shown (not on the curved surface of the sun):

sunspot_v = (vector([0.5824 - 0.5010, 0.1456 - 0.1837]) / ((5 + 12) - (8 + 0))).n(digits=4)
# sunspot_v = (0.03831, -0.01793) # per hour in the XY-plane in the graph.

the above as a speed (a scalar without direction):

sunspot_speed = sunspot_v.norm()
# sunspot_speed = 0.009986

Next step is to figure out what latitude the sunspots are on the sun. First we'll find a unit vector at right angle to sunspot_v which is a meridian on the sun that passes through the origin of the graph:

sunspot_meridian = vector([-sunspot_v[1], sunspot_v[0]]).normalized()
# sunspot_meridian = (0.4239, 0.9057)

The dot product with one of the sunspots gives a vector who's length is equal to the sine of the latitude.

sunspot_lat = asin(sunspot_meridian.dot_product(vector([0.5010, 0.1837])))
# sunspot_lat = 0.3885 # radians

The longitude on the sun of the sunspot is found by taking the component at right angles to sunspot_lat and mapping that onto a circle with radius cos(sunspot_lat) that's at right angles to the viewer:

sunspot_long = asin(sqrt(average_radius^2 - sunspot_lat^2) / cos(sunspot_lat))
# sunspot_long = 0.4662 # radians

The sun's expected rate of rotation in degrees per hour synodic (as seen from Earth):

A = 14.713
B = -2.396
C = -1.787
degrees_per_day_synodic = A + B * sin(sunspot_lat)^2 + \
    C * sin(sunspot_lat)^4 - 360/365.25
# degrees_per_day_synodic = 13.35

To translate this into an expected velocity in the same sense as sunspot_v first degrees must be translated into radians, then the expected velocity must be scaled by the length of the latitude, which is equal to cos(sunspot_lat).

The location on the latitude determines the apparent speed (it appears to be moving less if it's near the edge which means it's moving toward or away from us).

radians_per_hour_synodic  = (pi / 180) * (degrees_per_day_synodic / 24)
sunspot_speed_expected = cos(sunspot_lat) * radians_per_hour_synodic * cos(sunspot_long)
# sunspot_speed_expected = 0.008024

To get the expected direction of the sunspots we'll assume that the sun's axis is perpendicular to the orbital plane (it's actually off by 7°) and therefore that the sun's axis appears to be perfectly horizontal during the solstices. During the equinoxes the suns axis will appear to be tilted by an amount equal to the axial tilt (23.4°). These sunspots were recorded a bit after the Winter solstice, so we'd expect he sun's axis to be tilted clockwise from vertical. It's assumed that the change is sinusoidal as a function of the day of the amount of time that has passed since the Winter solstice. However, since the declination of the sun is also sinusoidal with a minimum at the Winter solstice. The expected direction is right angles to the sun's axis. The first step is to determine the declination of the sun:

axial_tilt = 23.4 * (pi / 180)
sun_declination = (a - (pi/2 - p)).subs({p:56.54*(pi/180), a:20*(pi/180)})
# sun_declination = -0.2349 # -13.46 degrees
observed_axis_tilt = sqrt(axial_tilt^2 - sun_declination^2)
# observed_axis_tilt = 0.3341 # radians

The expected direction is right angles to observed axis tilt:

expected_direction = vector([cos(observed_axis_tilt), -sin(observed_axis_tilt)]).n(digits=4)
# expected_direction = (0.9447, -0.3279)

Now with the expected direction and the expected speed we can calculate the expected velocity:

sunspot_v_expected = vector([i.n(digits=4) for i in
    sunspot_speed_expected * expected_direction])
# sunspot_v_expected = (0.007581, -0.002631)

The previously calculated sunspot locations with the effect of field rotation removed can now be updated so the effect of solar rotation is also removed:

sunspot_locs_w_sr = []
for i in range(len(sunspot_locs)):
    loc = sunspot_locs_w[i]
    min = solar_minutes[i]
    offset = min * (sunspot_v_expected / 60)
    sunspot_locs_w_sr.append(loc - offset)

An updated table with the expected effect of field rotation and solar rotation removed:

Hour Solar Time Solar Minute Field Rotation Sunspot WSRX Sunspot WSRY
8 am -4:40 -280 -0.5447 0.5364 0.1714
10 am -2:40 -160 -0.3715 0.5506 0.1675
12 pm -0:40 -40 -0.1016 0.5435 0.1493
2 pm 1:20 80 0.1997 0.5400 0.1778
4 pm 3:20 200 0.4414 0.5406 0.1624
5 pm 4:20 260 0.5236 0.5496 0.1570

To plot the above:

plot_sunspots(None, sunspot_locs_w_sr, average_radius, "sunspots-w-sr.png")

which produces:

As you can see the sunspot locations are so close now that it almost looks like one sunspot.
To plot it with zoom:

plot_sunspots(hours, sunspot_locs_w_sr, average_radius, "sunspots-w-sr-zoom.png",
    xmin=0.5010, xmax=0.5824, ymin=0.1456, ymax=0.1837)

which produces:

Conclusion

It is possible to calculate the field rotation as a function of time. In the case of sunspots field rotation, and to a lesser extent solar rotation, almost entirely account for the apparent change in a sunspot's location as a day progresses.