11-SourceFunction (template)#
We start by importing the modules
Numpy – operations on arrays and matrixes (and pi)
Matplotlib pyplot – plotting library
Matplotlin patches – module that enables patches in plots
Astropy units – defined quantities with units. We also import the CDS conversions
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib as mpl
%matplotlib inline
from astropy import constants as const
import astropy.units as u
from astropy.units import cds
cds.enable()
<astropy.units.core._UnitContext at 0x12266e910>
0. To execute: Below is the function that creates a gradiant of color between two curves, according to a certain function.#
def rect(ax,x,y,w,h,c):
#ax = plt.gca()
polygon = plt.Rectangle((x,y),w,h,color=c)
ax.add_patch(polygon)
def rainbow_fill_between(ax, X, Y1, Y2, colors=None,
cmap=plt.get_cmap("Reds"),**kwargs):
plt.plot(X,Y1,lw=0) # Plot so the axes scale correctly
dx = X[1]-X[0]
N = X.size
#Pad a float or int to same size as x
if (type(Y2) is float or type(Y2) is int):
Y2 = np.array([Y2]*N)
#No colors -- specify linear
if colors is None:
colors = []
for n in range(N):
colors.append(cmap(n/float(N)))
#Varying only in x
elif len(colors.shape) == 1:
colors = cmap((colors-colors.min())
/(colors.max()-colors.min()))
#Varying only in x and y
else:
cnp = np.array(colors)
colors = np.empty([colors.shape[0],colors.shape[1],4])
for i in range(colors.shape[0]):
for j in range(colors.shape[1]):
colors[i,j,:] = cmap((cnp[i,j]-cnp[:,:].min())
/(cnp[:,:].max()-cnp[:,:].min()))
colors = np.array(colors)
#Create the patch objects
for (color,x,y1,y2) in zip(colors,X,Y1,Y2):
rect(ax,x,y2,dx,y1-y2,color,**kwargs)
1. First, let’s consider a case where there is only emission.#
Note that the lenght of the slab is 4 length units.
a. The properties of the material are constant \(\rho(s)=\rho_o\) and \(j(s)=j_o\).#
The product of the emission coefficient and the density give an intensity of 0.25 per unit of length.
TODO: start from \(dI=...\) find an equation for \(I(s)\) and show your work below. Add a curve for \(I(s)\) in the graph below.
fig, ax = plt.subplots(1,1, figsize=(9,4))
plt.rcParams.update({'font.size': 14})
plt.rcParams.update({'font.family': 'Times New Roman'})
ax.set_ylim(0, 1.2)
ax.set_xlim(-1,5)
# For the first graph, I want a patch with constant color
cmap = mpl.colormaps["Reds"]
rgba = cmap(0.5) # pick the color in the center of the color map
rec = mpatches.Rectangle( (0, 0), 4 , 1.2, fc=rgba) # Create a shaded rectangle
p = ax.add_patch(rec) # add the rectangle to the plot
t = ax.text(5.2, 0.5, 'Observer is here', style='italic',
bbox={'facecolor':'red', 'pad':5})
ax.plot([-1,0],[0,0], lw=4, c='green')
ax.set_xlabel('$s$')
ax.set_ylabel('$I(s)$')
#---------------------------------------
#---------------------------------------
# At home
Text(0, 0.5, '$I(s)$')
TODO: Please write a small paragraph with an interpretation of the result obtained:
b. Now consider the emission inside a slab with density varying such that:#
and for which the emiisivity is still constant \(j(s)=j_o\).
TODO:
a. Find an expression for the intensity \(I(s)\) everywhere in the slab and show your work below
b. Analytically, find what the value of \(\rho_o j_o\) has to be for the final intensity to be the same than that of #1a. Show your work below.
c. Add a curve for \(I(s)\) to the graph below.
fig, ax = plt.subplots(1,1, figsize=(9,4))
plt.rcParams.update({'font.size': 14})
plt.rcParams.update({'font.family': 'Times New Roman'})
ax.set_ylim(0, 1.2)
ax.set_xlim(-1,5)
# Create a patch with a gradient to illustrate
# the change in density
X=np.linspace(0,4,100)
Y1=np.copy(X)*0
Y2=np.copy(X)*0+4
g = 1.0-np.copy(X/4)
rainbow_fill_between(ax,X,Y1,Y2,colors=g)
t = ax.text(5.2, 0.5, 'Observer is here', style='italic',
bbox={'facecolor':'red', 'pad':5})
ax.plot([-1,0],[0,0], lw=4, c='green')
#---------------------------------------
#---------------------------------------
# At home
[<matplotlib.lines.Line2D at 0x1241cccd0>]
TODO: Please write a small paragraph with an interpretation of the result obtained:
2. Now let’s look at the case where there is both absorption and emission.#
In class, we found the ‘formal solution’ of RT to be: $\(I(\tau) ~~=~~ ~ I_o\mathrm{e}^{\tau-\tau_o} ~+~ \int_{\tau'=\tau_o}^{\tau'=\tau} S(\tau')~ \mathrm{e}^{\tau-\tau'}d\tau'\)$
a. Constant source function (\(S(\tau)=So\))#
In this case, the , the formal solution can be written as: $\(I(s) ~~=~~ ~ I_o\mathrm{e}^{\tau(s)-\tau_o} ~+~ S_o\left[1-\mathrm{e}^{\tau(s)-\tau_o}\right]\)$
Here, let’s consider that everything else is also constant: \(\rho(s)=\rho_o\) and \(j(s)=j_o\).
Numerical values are defined in the code below.
TODO:
Find an expression for \(\tau(s)-\tau_o\). Show your work.
In the graph below, add a curve for the first term of the \(I(s)\) equation. Make sure to add a legend.
In the graph below, add a curve for the second term of the \(I(s)\) equation. Make sure to add a legend.
In the graph below, add a curve for \(I(s)\) itself. Make sure to add a legend.
Note: I have shown what the graph should look like in the class slides :)
fig, ax = plt.subplots(1,1, figsize=(9,4))
plt.rcParams.update({'font.size': 14})
plt.rcParams.update({'font.family': 'Times New Roman'})
ax.set_ylim(0, 2.2)
ax.set_xlim(-1,5)
# For the first graph, I want a patch with constant color
cmap = mpl.colormaps["Reds"]
rgba = cmap(0.5) # pick the color in the center of the color map
ax.add_patch(mpatches.Rectangle( (0, 0), 4 , 2.2, fc=rgba)) # add the rectangle to the plot
ax.text(5.2, 1, 'Observer is here', style='italic',
bbox={'facecolor':'red', 'pad':5})
ax.set_xlabel('$s$')
#---------------------------------------------------------
# Here are the values for already defined for this problem
Io = 2
j_rho = 0.5
kappa_rho = 0.5
S = j_rho / kappa_rho
d = 4
ss = np.linspace(0,4, 50)
# Curve for the source function value
ax.plot([-1,0], [Io,Io], c='green', lw=3)
source = np.array( [S]*ss.size )
ax.plot(ss, source, ls='dotted', c='k', lw=2, label='Source function')
#---------------------------------------
#---------------------------------------
# Optical depth
# First term (if there was no emission)
# Second term (if there was no initial intensity)
# Intensity
#-------------
ax.legend(loc=0)
<matplotlib.legend.Legend at 0x122ec54f0>
TODO: Please write a small paragraph with an interpretation of the result obtained:
b. Constant source function (\(S(\tau)=So\)) but linearly decreasing density#
TODO: Re-do what we did in #2a, assuming this time that the density is decreasing linearly such that:
\[ \rho(s) = \rho_o \left( 1- \frac{s}{d} \right), \]where \(d\) is the length of the slab.
Note: You can directly use some of the results from your previous notebook (e.g. for the optical depth).
fig, ax = plt.subplots(1,1, figsize=(9,4))
plt.rcParams.update({'font.size': 14})
plt.rcParams.update({'font.family': 'Times New Roman'})
ax.set_ylim(0, 2.2)
ax.set_xlim(-1,5)
# Create a patch with a gradient to illustrate
# the change in density
X=np.linspace(0,4,100)
Y1=np.copy(X)*0
Y2=np.copy(X)*0+4
g = 1.0-np.copy(X/4)
rainbow_fill_between(ax,X,Y1,Y2,colors=g)
ax.text(5.2, 1, 'Observer is here', style='italic',
bbox={'facecolor':'red', 'pad':5})
ax.set_xlabel('$s$')
#---------------------------------------------------------
# Here are the values for already defined for this problem
Io = 2
j_rho = 0.5
kappa_rho = 0.5
S = j_rho / kappa_rho
d = 4
ss = np.linspace(0,4, 50)
# Curve for the source function value
ax.plot([-1,0], [Io,Io], c='green', lw=3)
source = np.array( [S]*ss.size )
ax.plot(ss, source, ls='dotted', c='k', lw=2,label='Source function')
#---------------------------------------
#---------------------------------------
#---------------------------------------
ax.legend(loc=0)
<matplotlib.legend.Legend at 0x1240d8c70>
TODO: Please write a small paragraph with an interpretation of the result obtained:
3. Case when the source function is not constant.#
Let’s assume that the density and opacity in the slab are constant, such that \(\kappa_o\rho_o=2.0\) per unit length.
The source function is however a function of \(\tau\) such that: $\( S(\tau) = S_0 + S_1\tau\)\( where \)S_0=0.5\( intensity unit, and \)S_1=1.3$ intensity units per optical depth unit.
There is no intial intensity entering the slab so \(I_o = 0\).
Prepare your code such that you can vary the values of the paramters in the equations above.
a. Find an expression for \(I(\tau)\).#
TODO: In this case, you cannot use the expression we found for the constant source function. You will have to start back at the expression for the formal solution itself:
\[I(\tau) ~=~ ~ I_o\mathrm{e}^{\tau-\tau_o} ~+~ \int_{\tau'=\tau}^{\tau'=\tau_o} S(\tau')~ \mathrm{e}^{\tau-\tau'}d\tau'.\]But thankfully, given that \(\rho(s)=\rho_o\) and \(\kappa(s)=\kappa_o\), you can use your previous results for \(\tau(s)\) for constant opacity and density.
Show your work below
b. On your graph add a curve for:#
TODO:
the optical depth \(\tau\) as a function of \(s\).
the source function \(S\) as a function of \(s\).
the intensity \(I\) as a function of \(s\).
Concentrating on the side of your graph where the light exits the slab, estimate the optical depth at which the value of the source function is equal to the value of the intensity emerging from the slab.
fig, ax = plt.subplots(1,1, figsize=(10,5))
plt.rcParams.update({'font.size': 18})
plt.rcParams.update({'font.family': 'Times New Roman'})
ax.set_ylim(0, 8)
ax.set_xlim(-1,5)
# For the first graph, I want a patch with constant color
cmap = mpl.colormaps["Reds"]
rgba = cmap(0.5) # pick the color in the center of the color map
rec = mpatches.Rectangle( (0, 0), 4 , 8, fc=rgba) # Create a shaded rectangle
ax.add_patch(rec) # add the rectangle to the plot
ax.text(5.2, 4, 'Observer is here', style='italic',
bbox={'facecolor':'red', 'pad':5})
ax.set_xlabel('$s$')
# Values to use
d = 4
Io = 0
kappa_rho = 2.0
S0 = 0.5
S1 = 1.3
ss = np.linspace(0,4, 50)
ax.plot([-1,0],[0,0], lw=4, c='green')
#---------------------------------------
#---------------------------------------
# At home
[<matplotlib.lines.Line2D at 0x1242a4b80>]
TODO: Please write a small paragraph with an interpretation of the result obtained (hint: remember what the intensity aspires to become :):