<jupyter_start><jupyter_text>Staying Positive<jupyter_code>%matplotlib inline
%run notebook_setup.py<jupyter_output><empty_output><jupyter_text>Unfortunately, there's no trivial way of setting the spherical harmonic coefficients to ensure the map intensity is non-negative everywhere. One helpful thing to keep in mind is that the average intensity across the surface does not change when you modify the spherical harmonic coefficients. That's because all spherical harmonics other than $Y_{0,0}$ are perfectly anti-symmetric: for every bright region on the surface, there's an equally dark region on the other side that cancels its contribution to the surface-integrated intensity. What this means is that the magnitude of the spherical harmonic coefficients controls the departure of the intensity from this mean value (which is equal to the intensity of the $Y_{0,0}$ harmonic, $\frac{1}{\pi}$). One way to ensure the map is non-negative everywhere is therefore simply to limit the amplitude of all spherical harmonic coefficients to a small value.  Spherical harmonic maps However, in some cases it might be useful to check in the general case whether or not a map is positive semi-definite. This entails running a nonlinear minimization on the intensity across the surface. For spherical harmonic maps, this is implemented in the `minimize()` method. Let's import `starry`, instantiate a map, and take a look.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import starry

starry.config.lazy = False
starry.config.quiet = True<jupyter_output><empty_output><jupyter_text>Let's set the coefficients randomly...<jupyter_code>map = starry.Map(ydeg=5)
np.random.seed(0)
map[1:, :] = 0.1 * np.random.randn(map.Ny - 1)<jupyter_output><empty_output><jupyter_text>... and plot the map on a lat-lon grid:<jupyter_code>image = map.render(projection="rect")
plt.imshow(image, origin="lower", cmap="plasma", extent=(-180, 180, -90, 90))
plt.xlabel("longitude [deg]")
plt.ylabel("latitude [deg]")
plt.colorbar();
# Pre-run this to ensure it's compiled
map.minimize();<jupyter_output><empty_output><jupyter_text>The map clearly goes negative in certain regions. If we call `minimize()`, we can get the latitude, longitude, and value of the intensity at the minimum:<jupyter_code>%%time
lat, lon, value = map.minimize()
print(lat, lon, value)
plt.imshow(image, origin="lower", cmap="plasma", extent=(-180, 180, -90, 90))
plt.xlabel("longitude [deg]")
plt.ylabel("latitude [deg]")
plt.axvline(lon, color="k", ls="--")
plt.axhline(lat, color="k", ls="--")
plt.title("minimum: {0:.3f}".format(value))
plt.colorbar();<jupyter_output><empty_output><jupyter_text>The method is *fairly* fast, so it could be used, for example, as a penalty when doing inference.  Limb-darkened maps For limb-darkened maps, it's a little easier to check whether the map is positive everywhere. Because the limb darkening profile is one-dimensional, we can use [Sturm's theorem](https://en.wikipedia.org/wiki/Sturm%27s_theorem) to verify that   * The intensity is non-negative everywhere   * The intensity is monotonically decreasing toward the limbLimb-darkened maps (or spherical harmonic maps with a limb darkening filter) implement the `limbdark_is_physical` method, which checks whether both points are true. Note, importantly, that the second point is specific to limb *darkening*. In principle the specific intensity could get brighter toward the limb (as is the case at certain wavelengths for the Sun), so you wouldn't want to use in those cases.<jupyter_code>map = starry.Map(udeg=4)<jupyter_output><empty_output><jupyter_text>Let's try this on a few limb-darkened maps:<jupyter_code>map[1:] = [0.5, 0.25, 0.5, 0.25]
map.show()<jupyter_output><empty_output><jupyter_text>Is it physical?<jupyter_code>map.limbdark_is_physical()<jupyter_output><empty_output><jupyter_text>No! Let's plot the intensity as a function of $\mu$ to see why:<jupyter_code>mu = np.linspace(0, 1, 1000)
plt.plot(mu, map.intensity(mu=mu))
plt.axhline(0, color="k", ls="--")
plt.gca().invert_xaxis()
plt.xlabel(r"$\mu$")
plt.ylabel("relative intensity");<jupyter_output><empty_output><jupyter_text>The intensity is negative close to the limb.  Let's try a different coefficient vector:<jupyter_code>map[1:] = [0.1, -2.0, 2.25, 0.5]
map.show()<jupyter_output><empty_output><jupyter_text>Is it physical?<jupyter_code>map.limbdark_is_physical()
mu = np.linspace(0, 1, 1000)
plt.plot(mu, map.intensity(mu=mu))
plt.axhline(0, color="k", ls="--")
plt.gca().invert_xaxis()
plt.xlabel(r"$\mu$")
plt.ylabel("relative intensity");<jupyter_output><empty_output><jupyter_text>Even though it's positive everywhere, it's not monotonic! One last example:<jupyter_code>map[1:] = [0.5, -0.1, 0.25, 0.25]
map.show()<jupyter_output><empty_output><jupyter_text>Is it physical?<jupyter_code>map.limbdark_is_physical()
mu = np.linspace(0, 1, 1000)
plt.plot(mu, map.intensity(mu=mu))
plt.axhline(0, color="k", ls="--")
plt.gca().invert_xaxis()
plt.xlabel(r"$\mu$")
plt.ylabel("relative intensity");<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1226061204343534, 'fraction_numerical': 0.0175715695952616, 'mean_word_length': 6.3314037626628075, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '518178', 'license': 'MIT', 'n_tokens_mistral': 1783, 'n_tokens_neox': 1631, 'n_words': 528, 'path': 'notebooks/StayingPositive.ipynb', 'repo_name': 'fbartolic/starry', 'size': '5065'}