# SageMath commands to calculate field rotation. See # https://selliott.org/science/rotation # All variables in alphabetical order. var("a cd cdMP cuMP h m o oMP od odMP ou ouMP p R Ru s t u ux uy uz w wh wm ws") # Given vectors and matrices. # Unit normal vectors at the prime meridian crossing ("MP" suffix). cdMP = vector([1, 0, 0]) # not used cuMP = vector([0, cos(a), sin(a)]) # not used oMP = vector([0, sin(a), -cos(a)]) odMP = vector([1, 0, 0]) ouMP = vector([0, cos(a), sin(a)]) # not used # Unit normal rotation vector around Earth's axis. u = vector([0, -sin(p), -cos(p)]) # The rotation matrix for rotating angle t around unit vector u. See # en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle 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)] ]) # Derived vectors and matrices. # Rotation matrix substituting values for ux, uy and uz. Ru = R.subs({ux:0, uy:-sin(p), uz:-cos(p)}) # Unit normal vectors rotated angle t around the celestial pole (no "MP" # suffix). o = Ru * oMP od = Ru * odMP # The camera dosn't rotate around the celestial pole. Instead, it must point to # the celestial object (o). For cr, which is in the X-Z plane, we can come up # with a vector that is right angles to o projected on to the X-Z plane by # switching the coordinates around. It must be normalized since excluding the Y # coordinate makes it shorter than 1. cd = vector([-o[2], 0, o[0]]).normalized() # Since the field rotation is the amount ou has rotated since ouMP and since # that rotation is in a plane perpendicular to the line of sight (perpendicular # to o) the field rotation (w) is the angle between od and cr. The amount of # rotation is positive if od is pointing downward. Value -sgn(od) would work as # the acos() multiplier, but a multiplier of zero to be avoided. 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)) # Calculate w as a function of time (h, m or s) after (positive) or since # (negative) the MP wh = w.subs({t:15*h*(pi/180)}) wm = wh.subs({h:m/60}) ws = wm.subs({m:s/60}) # 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"] # This is the raw data from the first table. # Hour SL IX SR IX ST IY SB IY S IX S IY sunspot_raw = [ (hours[0], 170, 416, 166, 382, 334, 229), (hours[1], 568, 814, 153, 398, 744, 232), (hours[2], 992, 1238, 152, 396, 1179, 249), (hours[3], 171, 415, 468, 711, 363, 582), (hours[4], 568, 814, 468, 711, 762, 602), (hours[5], 992, 1238, 481, 699, 1186, 608)] # Iterate through sunspot_raw to produce sunspot_locs. This conceptually # applies the following formula to each sunspot location: # 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 sunspot_locs = [] for raw in sunspot_raw: x = (2 * ((raw[5] - raw[1]) / (raw[2] - raw[1])) - 1).n(digits=4) y = -(2 * ((raw[6] - raw[3]) / (raw[4] - raw[3])) - 1).n(digits=4) sunspot_locs.append(vector([x, y])) # The following appears in the second table: # 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)] # The radius of each sunspot, and the average. sunspots_radius = [sqrt(l[0]^2 + l[1]^2) for l in sunspot_locs] average_radius = sum(sunspots_radius) / len(sunspots_radius) # 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) # Plot the sunspots as seen. plot_sunspots(hours, sunspot_locs, average_radius, "sunspots.png") # Calculate the field rotations in radians. 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] Rw = Matrix([ [cos(w), -sin(w)], [sin(w), cos(w)] ]) 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)] # Plot the sunspot locations with the effect of field rotation removed. plot_sunspots(None, sunspot_locs_w, average_radius, "sunspots-w.png") # Zoom in on the sunspot locations: plot_sunspots(hours, sunspot_locs_w, average_radius, "sunspots-w-zoom.png", xmin=0.5010, xmax=0.5824, ymin=0.1456, ymax=0.1837) 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. sunspot_speed = sunspot_v.norm() # sunspot_speed = 0.009986 sunspot_meridian = vector([-sunspot_v[1], sunspot_v[0]]).normalized() # sunspot_meridian = (0.4239, 0.9057) sunspot_lat = asin(sunspot_meridian.dot_product(vector([0.5010, 0.1837]))) # sunspot_lat = 0.3885 # radians sunspot_long = asin(sqrt(average_radius^2 - sunspot_lat^2) / cos(sunspot_lat)) # sunspot_long = 0.4662 # radians 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 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 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" means perpendicular to the line of sight, so it can be seen. The # observed part of axial tilt, and the part that can't be seen, the part along # the line of sight (called "hidden" here), are at right angles. Since the # hidden axial tilt is in phase with the sun's declination we expect the # observed axial tilt to be 90 degrees out of phase with the sun's declination. observed_axis_tilt = sqrt(axial_tilt^2 - sun_declination^2) # observed_axis_tilt = 0.3341 # radians # If the obsered_axis_tilt was 0 degrees (Winter Solstice) the expected # direction would be entirely to the right. expected_direction = vector([cos(observed_axis_tilt), -sin(observed_axis_tilt)]).n(digits=4) # expected_direction = (0.9447, -0.3279) sunspot_v_expected = vector([i.n(digits=4) for i in sunspot_speed_expected * expected_direction]) # Using the expected velocity it should be possible to undo the spread in # location due to solar rotation. The "/ 60" is to convert to the velocity to # minutes so that solar minute may be used. 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) plot_sunspots(None, sunspot_locs_w_sr, average_radius, "sunspots-w-sr.png") # A zoomed in zersion of the above with the same scale as the previous zoomed # in graph: 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)