"""Independent SciPy reference for the AC21 finite-volume layer and sky integrals.""" import numpy as np from scipy.integrate import quad from scipy.sparse import diags, lil_matrix from scipy.sparse.linalg import spsolve SIGMA = 5.670374419e-8 STAR_TEMPERATURE = 4520. STAR_RADIUS = .667 * 6.957e8 PLANET_RADIUS = 7e7 ORBIT_DAYS = .81347437 INCLINATION = np.deg2rad(82.27) IRRADIANCE = SIGMA * STAR_TEMPERATURE**4 / 4.88**2 def heating_cells(longitudes, latitudes, albedo): edges = np.linspace(-np.pi, np.pi, longitudes + 1) low = np.maximum(edges[:-1], -np.pi / 2) high = np.minimum(edges[1:], np.pi / 2) longitude_mean = np.where(high > low, (np.sin(high) - np.sin(low)) / np.diff(edges), 0.) latitude_edges = np.linspace(-1, 1, latitudes + 1) primitive = .5 * (latitude_edges * np.sqrt(np.maximum(0, 1 - latitude_edges**2)) + np.arcsin(latitude_edges)) latitude_mean = np.diff(primitive) / np.diff(latitude_edges) return (1 - albedo) * IRRADIANCE * latitude_mean[:, None] * longitude_mean[None, :] + SIGMA * 150.**4 def reference_field(longitudes=128, latitudes=24, capacity=1e7, omega=2000 / 7e7, albedo=.2): source = heating_cells(longitudes, latitudes, albedo) coefficient = capacity * abs(omega) / (2 * np.pi / longitudes) transport = lil_matrix((longitudes, longitudes)) for index in range(longitudes): transport[index, index] = coefficient transport[index, (index - (1 if omega >= 0 else -1)) % longitudes] = -coefficient transport = transport.tocsr() temperature = np.maximum(1000., (source / SIGMA)**.25) for latitude in range(latitudes): trial = temperature[latitude] for _ in range(60): residual = transport.dot(trial) + SIGMA * trial**4 - source[latitude] if np.max(np.abs(residual)) < 1e-8: break correction = spsolve(transport + diags(4 * SIGMA * trial**3), -residual) scale = 1. while np.min(trial + scale * correction) <= 0: scale /= 2 trial += scale * correction else: raise RuntimeError('Independent sparse Newton solve did not converge') return temperature def properties(temperature, albedo=.2): latitudes, longitudes = temperature.shape longitude = -np.pi + (np.arange(longitudes) + .5) * 2 * np.pi / longitudes weight = 4 * np.pi / temperature.size emitted = np.sum(SIGMA * temperature**4) * weight night = np.sum((SIGMA * temperature[:, abs(longitude) > np.pi / 2]**4 - SIGMA * 150.**4)) * weight absorbed = (1 - albedo) * IRRADIANCE * np.pi return dict(minimum_kelvin=float(temperature.min()), maximum_kelvin=float(temperature.max()), emitted_w_per_radius_squared=float(emitted), nightside_net_fraction=float(night / absorbed), relative_energy_residual=float((emitted - absorbed - 4 * np.pi * SIGMA * 150.**4) / absorbed)) def surface_curve(temperature, phases): latitudes, longitudes = temperature.shape latitude_sine = -1 + (np.arange(latitudes) + .5) * 2 / latitudes latitude_cosine = np.sqrt(1 - latitude_sine**2) longitude = -np.pi + (np.arange(longitudes) + .5) * 2 * np.pi / longitudes normalization = (PLANET_RADIUS / STAR_RADIUS)**2 / (SIGMA * STAR_TEMPERATURE**4) * 1e6 values = [] for phase in phases: observer_longitude = np.pi - 2 * np.pi * phase projection = (latitude_cosine[:, None] * np.sin(INCLINATION) * np.cos(longitude[None, :] - observer_longitude) + latitude_sine[:, None] * np.cos(INCLINATION)) values.append(float(np.sum(SIGMA * temperature**4 / np.pi * np.maximum(projection, 0)) * 4 * np.pi / temperature.size * normalization)) return np.array(values) def linear_map_eclipse(phase): separation = 4.88 * STAR_RADIUS / PLANET_RADIUS stellar_x = separation * np.sin(2 * np.pi * phase) stellar_y = separation * np.cos(INCLINATION) * np.cos(2 * np.pi * phase) stellar_radius = STAR_RADIUS / PLANET_RADIUS def blocked_chord(x): star_squared = stellar_radius**2 - (x - stellar_x)**2 if star_squared <= 0: return 0. planet_half = np.sqrt(max(0, 1 - x*x)) star_half = np.sqrt(star_squared) low = max(-planet_half, stellar_y - star_half) high = min(planet_half, stellar_y + star_half) return max(0., high - low) * (1 + .3*x) intersections = [] # Circle intersection x coordinates provide quadrature breakpoints. distance = np.hypot(stellar_x, stellar_y) if abs(stellar_radius - 1) < distance < stellar_radius + 1: along = (1 - stellar_radius**2 + distance**2) / (2 * distance) height = np.sqrt(1 - along**2) intersections = [along * stellar_x / distance + sign * height * stellar_y / distance for sign in [-1, 1]] points = sorted(x for x in intersections if -1 < x < 1) occulted = quad(blocked_chord, -1, 1, points=points, epsabs=1e-10, epsrel=1e-10, limit=150)[0] return dict(phase=phase, star_x_planet_radii=float(stellar_x), star_y_planet_radii=float(stellar_y), unocculted_integral=float(np.pi - occulted), intensity='1 + 0.3 projected-x, positive across the unit disk') def resolution_check(phases, best_node): reference = surface_curve(reference_field(256, 32), phases) coarse = surface_curve(reference_field(128, 24), phases) fine_best = surface_curve(reference_field(256, 48, 1e7, best_node['capacity_omega_w_m2_k']/1e7, best_node['albedo']), phases) coarse_best = surface_curve(reference_field(128, 24, 1e7, best_node['capacity_omega_w_m2_k']/1e7, best_node['albedo']), phases) return dict(default_maximum_128x24_vs_256x32_ppm=float(np.max(np.abs(reference-coarse))), default_rms_difference_ppm=float(np.sqrt(np.mean((reference-coarse)**2))), best_node_maximum_128x24_vs_256x48_ppm=float(np.max(np.abs(fine_best-coarse_best))), interpretation='Forward-model resolution differences, compared with the specified synthetic noise of twelve ppm; not observational calibration uncertainty.') def global_east_map_eclipse(phase): """Independent vertical-chord integral for fixed intensity 1 + 0.3 east-normal.""" separation = 4.88 * STAR_RADIUS / PLANET_RADIUS star_x = separation*np.sin(2*np.pi*phase) star_y = separation*np.cos(INCLINATION)*np.cos(2*np.pi*phase) star_radius = STAR_RADIUS/PLANET_RADIUS observer_longitude = np.pi-2*np.pi*phase observer_latitude = np.pi/2-INCLINATION ax = -.3*np.cos(observer_longitude) ay = -.3*np.sin(observer_latitude)*np.sin(observer_longitude) az = .3*np.cos(observer_latitude)*np.sin(observer_longitude) unocculted = np.pi + az*2*np.pi/3 if np.cos(2*np.pi*phase) >= 0: return float(unocculted) def blocked(x): planet_squared = max(0.,1-x*x) star_squared = star_radius**2-(x-star_x)**2 if star_squared <= 0: return 0. half = np.sqrt(planet_squared) star_half = np.sqrt(star_squared) low=max(-half,star_y-star_half) high=min(half,star_y+star_half) if high <= low: return 0. def z_primitive(y): if half == 0: return 0. return .5*(y*np.sqrt(max(0.,planet_squared-y*y)) + planet_squared*np.arcsin(np.clip(y/half,-1,1))) return ((1+ax*x)*(high-low)+ay*(high*high-low*low)/2 +az*(z_primitive(high)-z_primitive(low))) distance=np.hypot(star_x,star_y) points=[] if abs(star_radius-1) < distance < star_radius+1: along=(1-star_radius**2+distance**2)/(2*distance) height=np.sqrt(1-along**2) points=[along*star_x/distance+sign*height*star_y/distance for sign in [-1,1]] hidden=quad(blocked,-1,1,points=sorted(x for x in points if -1