Need help simulating solar limb darkening

Question

edit: At the Solar Dynamics Observatory (SDO) website, I just found the image sdo.gsfc.nasa.gov/assets/img/latest/latest_1024_HMIIC.jpg. The color gradient of the limb darkening seems very similar to the Wikimedia image below. I've discovered that it is called a "colorized intensitygram" and the color gradient is purely artificial - the data is single channel intensity. The limb darkening is certainly real (compare to the artificially "flattented" display!)

From http://www.solarham.net/latest_imagery/hmi1.htm

I would like to try to start to understand what the sun actually "looks like", with human vision, assuming the brightness has been reduced.

Like the gas giant planets, the sun doesn't have an abrupt surface, it just gets denser and hotter and denser and hotter as you go deeper.

One consequence of this is limb-darkening. As the material becomes denser the deeper you go, it becomes more opaque, so - roughly speaking - the light you see, including the color and brightness, is determined by the layers above that point.

If you look at center of the solar disk, you can see deeper and hotter parts. If you look near the edge of the sun, or solar limb, the incidence is oblique and you are not seeing as deep or as hot.

The solar physics and photon transport theory is complicated, but in the visible part of the spectrum it may be OK to think of blue light scattering much more strongly than red light. In the thick, dense solar atmosphere, it scatters enough to attenuate. So in blue light you are seeing even shallower, which is even colder, and therefore dimmer in blue.

Right now I'd just like some good analytical approximation of the wavelength dependent limb darkening of the sun, or some actual linear images (before web-processing) of the sun at various visible wavelengths so I can make my own approximation.

Here is an image from Wikimedia titled 2012_Transit_of_Venus_from_SF which I've separated into RGB channels and plotted brightness across horizontal (solid) and vertical (dashed) diameters, 20 pixels wide. You can see the dramatic difference in behavior of different wavelengths. Since the image is "from the internet" I've no information about the linearity, so it's an illustrative example, but not what you'd think of as data.

import numpy as np
import matplotlib.pyplot as plt

img = plt.imread("sun limb darkening.png") 
# FROM: https://upload.wikimedia.org/wikipedia/commons/thumb/4/4d/2012_Transit_of_Venus_from_SF.jpg/600px-2012_Transit_of_Venus_from_SF.jpg

w, c, hw  = 600, 300, 10   # image happens to be 600x600 pixels square!

hor = img[    c-hw:c+hw].sum(axis=0) / (2*hw)
ver = img[ :, c-hw:c+hw].sum(axis=1) / (2*hw)

plt.figure()

plt.imshow(img)

r, g, b = c - 0.67*c*hor.T

plt.plot(r, '-r')    
plt.plot(g, '-g')    
plt.plot(b, '-b')    

r, g, b = c - 0.67*c*ver.T

plt.plot(r, '--r')    
plt.plot(g, '--g')    
plt.plot(b, '--b')

plt.plot([0, w], [c, c], '-k')
plt.xlim(0, 599)
plt.ylim(599, 0)

plt.text(150, 137, 'R', fontsize=18)
plt.text(156, 192, 'G', fontsize=18)
plt.text(170, 272, 'B', fontsize=18)

plt.show()

Answer 1

I used the paper Wavelength dependency of the Solar limb darkening for solar limb darkening data.

It uses the following model for the normalized brightness distribution across the disk of the Sun:

$$ I(\mu)=1-u(1-\mu^\alpha) $$

Here, $\mu$ is the normalized distance from the limb; when far from the Sun, it can be expressed in terms of the normalized distance from the center of the disk, $r$:

$$ \mu = \sqrt{1-r^2} $$

$u$ and $\alpha$ are parameters. The paper is not very clear about the value of $u$, but it seems like they used $u=0.85$.

This formula tells us how the brightness varies from the center of the solar disk to the limb for a given wavelength. In order to properly represent the color of the Sun, we also need to know the relative brightness of different wavelength bands.

For this we need data on the solar spectrum, which I obtained from here (this data is for "air mass zero," meaning our result will be what the Sun would look like outside the Earth's atmosphere).

I assume that this spectrum is averaged over the whole disk, in which case we need to be able to scale the limb darkening model by average instead of peak (center) brightness. The formula given in the paper is:

$$ I_\text{avg} \propto \int_0^1 I(\mu)\mu\ d\mu \propto \frac{2 + \alpha (1 - u)}{2 + \alpha} $$

The last piece of the puzzle is to convert the spectrum to a color. In this case I first compute the color in the CIE XYZ color space, then transform to linear RGB and then sRGB. (See this answer I posted on photography.SE for more details about computing the XYZ color of a spectrum, and the sRGB Wikipedia page for more detail about the conversion to sRGB.) I used the standard colorimetric observer tables from here.

---

Here is the Mathematica code that I used:

alphaRaw = ImportString[(* copy and paste from the PDF table *), "Table"];
(* remove headers *)
Select[ArrayQ] @ SplitBy[alphaRaw, Head@*First];
(* join corresponding columns *)
Flatten[Partition[%, 2], {{2}, {1, 3}}];
(* resample to 1 nm *)
{alphaPS, alphaNL} = (Interpolation[#] /@ Range[360, 830]) & /@ %;

spectrumRaw = Import["http://rredc.nrel.gov/solar/spectra/am0/E490_00a_AM0.xls"];
(* resample to 1 nm (inputs in um!) *)
spectrum = Interpolation[spectrumRaw[[1, 2 ;;, ;; 2]]] /@ (Range[360, 830] / 1000);

cieRaw = Import["http://www.cis.rit.edu/research/mcsl2/online/CIE/StdObsFuncs.xls"];
(* pair each color function column with the wavelength column *)
Transpose[{cieRaw[[1, 6 ;;, 2 ;; 4]]], #}] & /@ cieRaw[[1, 6 ;;, 2 ;; 4]]];
(* resample to 1 nm *)
cieXYZ = (Interpolation[#] /@ Range[360, 830]) & /@ %;

With[{u = 0.85},
  (* compute color at 101 points evenly spaced in mu *)
  colorXYZ = Transpose @ Table[
               Total[ (* integration -> summation *)
                 Transpose @ cieXYZ
                 * spectrum
                 * (1 - u (1 - mu^alphaPS))
                 / ((2 + alphaPS (1 - u)) / (2 + alphaPS))
               ],
               {mu, 0, 1, 0.01}
             ];
  (* scale so maximum = 1 *)
  colorXYZ /= Max[colorXYZ / {0.97, 1, 0.83}];
  (* dim by 15% *)
  colorXYZ *= 0.85;
]

(I'm not sure why Mathematica renders the XYZ colors with a blue tint.)

colorRGBlinear = {{3.2406, -1.5372, -0.4986}, {-0.9689, 1.8758, 0.0415}, {0.0557, -0.2040, 1.0570}}.colorXYZ;
colorRGBlinear /= Max[colorRGBlinear];
colorRGBlinear *= 0.85;

gamma = With[{a = 0.055}, If[# <= 0.0031308, 12.92 #, (1 + a) #^(1/2.4) - a] &];
colorSRGB = Map[gamma, colorRGBlinear, {2}];

Plot of the color difference (in CIE xy coordinates) from the limb to center (small black arrow). The ticks on the spectral locus are wavelength in nm, and the ticks on the Planckian locus are temperature in K.