Skip to content

Commit

Permalink
Merge pull request #230 from karllark/add_2ndaxis
Browse files Browse the repository at this point in the history 
Add 2nd x-axis for plots with inverse microns as the main x axis
  • Loading branch information
@karllark
karllark authored Aug 16, 2024
2 parents 6601d6a + fade962 commit d8b9e40
Show file tree
Hide file tree
Showing 5 changed files with 257 additions and 46 deletions.
90 changes: 53 additions & 37 deletions docs/dust_extinction/fit_extinction.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,55 +20,63 @@ extinction curve for the LMC outside of the LMC2 supershell region
.. plot::
:include-source:

import warnings
import matplotlib.pyplot as plt
import numpy as np
import warnings
import matplotlib.pyplot as plt
import numpy as np

from astropy.modeling.fitting import LevMarLSQFitter
import astropy.units as u
from astropy.modeling.fitting import LevMarLSQFitter
import astropy.units as u

from dust_extinction.averages import G03_LMCAvg
from dust_extinction.shapes import FM90
from dust_extinction.averages import G03_LMCAvg
from dust_extinction.shapes import FM90

# get an observed extinction curve to fit
g03_model = G03_LMCAvg()
# get an observed extinction curve to fit
g03_model = G03_LMCAvg()

x = g03_model.obsdata_x / u.micron
# convert to E(x-V)/E(B0V)
y = (g03_model.obsdata_axav - 1.0) * g03_model.Rv
# only fit the UV portion (FM90 only valid in UV)
(gindxs,) = np.where(x > 3.125 / u.micron)
x = g03_model.obsdata_x / u.micron
# convert to E(x-V)/E(B0V)
y = (g03_model.obsdata_axav - 1.0) * g03_model.Rv
# only fit the UV portion (FM90 only valid in UV)
(gindxs,) = np.where(x > 3.125 / u.micron)

# initialize the model
fm90_init = FM90()
# initialize the model
fm90_init = FM90()

# pick the fitter
fit = LevMarLSQFitter()
# pick the fitter
fit = LevMarLSQFitter()

# fit the data to the FM90 model using the fitter
# use the initialized model as the starting point
# fit the data to the FM90 model using the fitter
# use the initialized model as the starting point

# ignore some warnings
# UserWarning is to avoid the units of x warning
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=UserWarning)
g03_fit = fit(fm90_init, x[gindxs].value, y[gindxs])
# ignore some warnings
# UserWarning is to avoid the units of x warning
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=UserWarning)
g03_fit = fit(fm90_init, x[gindxs].value, y[gindxs])

# plot the observed data, initial guess, and final fit
fig, ax = plt.subplots()
# plot the observed data, initial guess, and final fit
fig, ax = plt.subplots()

ax.plot(x, y, "ko", label="Observed Curve")
ax.plot(x[gindxs], fm90_init(x[gindxs]), label="Initial guess")
ax.plot(x[gindxs], g03_fit(x[gindxs]), label="Fitted model")
ax.plot(x, y, "ko", label="Observed Curve")
ax.plot(x[gindxs], fm90_init(x[gindxs]), label="Initial guess")
ax.plot(x[gindxs], g03_fit(x[gindxs]), label="Fitted model")

ax.set_xlabel("$x$ [$\mu m^{-1}$]")
ax.set_ylabel("$E(x-V)/E(B-V)$")
ax.set_xlabel("$x$ [$\mu m^{-1}$]")
ax.set_ylabel("$E(x-V)/E(B-V)$")

ax.set_title("Example FM90 Fit to G03_LMCAvg curve")
# for 2nd x-axis with lambda values
axis_xs = np.array([0.12, 0.15, 0.2, 0.3, 0.5, 1.0, 2.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")

ax.legend(loc="best")
plt.tight_layout()
plt.show()
ax.legend(loc="best")
plt.tight_layout()
plt.show()

Example: P92 Fit
================
Expand Down Expand Up @@ -144,7 +152,15 @@ between data points.
ax.set_xlabel('$x$ [$\mu m^{-1}$]')
ax.set_ylabel('$A(x)/A(V)$')

ax.set_title('Example P92 Fit to GCC09_MWAvg average curve')
# for 2nd x-axis with lambda values
axis_xs = np.array([0.1, 0.12, 0.15, 0.2, 0.3, 0.5, 1.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")

ax.legend(loc='best')
plt.tight_layout()
Expand Down
2 changes: 1 addition & 1 deletion docs/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ Interstellar Dust Extinction
extinction curves.

Extinction describes the effects of dust on observations of single star due to
the dust along the line-of-sight to a star removiong flux by absorbing photons
the dust along the line-of-sight to a star removing flux by absorbing photons
and scattering photons out of the line-of-sight. The wavelength dependence of
dust extinction (also know as extinction curves) provides fundamental
information about the size, composition, and shape of interstellar dust grain.
Expand Down
80 changes: 74 additions & 6 deletions dust_extinction/averages.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,7 +260,7 @@ class B92_MWAvg(BaseExtModel):
# generate the curves and plot them
x = np.arange(1.0/ext_model.x_range[1], 1.0/ext_model.x_range[0], 0.1) * u.micron
ax.plot(x,ext_model(x),label='B1992')
ax.plot(x,ext_model(x),label='B92')
ax.plot(1.0/ext_model.obsdata_x, ext_model.obsdata_axav, 'ko',
label='obsdata')
Expand Down Expand Up @@ -295,7 +295,7 @@ class B92_MWAvg(BaseExtModel):

def evaluate(self, in_x):
"""
B1992 function
B92 function
Parameters
----------
Expand Down Expand Up @@ -373,6 +373,16 @@ class G03_SMCBar(BaseExtModel):
ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
ax.set_ylabel(r'$A(x)/A(V)$')
# for 2nd x-axis with lambda values
axis_xs = np.array([0.1, 0.12, 0.15, 0.2, 0.3, 0.5, 1.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")
ax.legend(loc='best')
plt.show()
"""
Expand Down Expand Up @@ -494,6 +504,16 @@ class G03_LMCAvg(BaseExtModel):
ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
ax.set_ylabel(r'$A(x)/A(V)$')
# for 2nd x-axis with lambda values
axis_xs = np.array([0.1, 0.12, 0.15, 0.2, 0.3, 0.5, 1.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")
ax.legend(loc='best')
plt.show()
"""
Expand Down Expand Up @@ -616,6 +636,16 @@ class G03_LMC2(BaseExtModel):
ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
ax.set_ylabel(r'$A(x)/A(V)$')
# for 2nd x-axis with lambda values
axis_xs = np.array([0.1, 0.12, 0.15, 0.2, 0.3, 0.5, 1.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")
ax.legend(loc='best')
plt.show()
"""
Expand Down Expand Up @@ -1048,6 +1078,16 @@ class GCC09_MWAvg(BaseExtModel):
ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
ax.set_ylabel(r'$A(x)/A(V)$')
# for 2nd x-axis with lambda values
axis_xs = np.array([0.09, 0.1, 0.12, 0.15, 0.2, 0.3, 0.5, 1.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")
ax.legend(loc='best')
plt.show()
"""
Expand All @@ -1062,9 +1102,13 @@ def __init__(self, **kwargs):
ref = importlib_resources.files("dust_extinction") / "data"
with importlib_resources.as_file(ref) as data_path:
# GCC09 sigma clipped average of 75 sightlines
a = Table.read(data_path / "GCC09_FUSE.dat", format="ascii.commented_header")
a = Table.read(
data_path / "GCC09_FUSE.dat", format="ascii.commented_header"
)
b = Table.read(data_path / "GCC09_IUE.dat", format="ascii.commented_header")
c = Table.read(data_path / "GCC09_PHOT.dat", format="ascii.commented_header")
c = Table.read(
data_path / "GCC09_PHOT.dat", format="ascii.commented_header"
)

# FUSE range
self.obsdata_x_fuse = a["x"].data
Expand Down Expand Up @@ -1559,6 +1603,16 @@ class G24_SMCAvg(BaseExtModel):
ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
ax.set_ylabel(r'$A(x)/A(V)$')
# for 2nd x-axis with lambda values
axis_xs = np.array([0.1, 0.12, 0.15, 0.2, 0.3, 0.5, 1.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")
ax.legend(loc='best')
plt.show()
"""
Expand All @@ -1573,7 +1627,9 @@ def __init__(self, **kwargs):
ref = importlib_resources.files("dust_extinction") / "data"
with importlib_resources.as_file(ref) as data_path:
# D22 sigma clipped average of 13 diffuse sightlines
a = Table.read(data_path / "G24_SMCAvg.dat", format="ascii.commented_header")
a = Table.read(
data_path / "G24_SMCAvg.dat", format="ascii.commented_header"
)

# data
self.obsdata_x = 1.0 / a["wave"].data
Expand Down Expand Up @@ -1682,6 +1738,16 @@ class G24_SMCBumps(BaseExtModel):
ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
ax.set_ylabel(r'$A(x)/A(V)$')
# for 2nd x-axis with lambda values
axis_xs = np.array([0.1, 0.12, 0.15, 0.2, 0.3, 0.5, 1.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")
ax.legend(loc='best')
plt.show()
"""
Expand All @@ -1696,7 +1762,9 @@ def __init__(self, **kwargs):
ref = importlib_resources.files("dust_extinction") / "data"
with importlib_resources.as_file(ref) as data_path:
# D22 sigma clipped average of 13 diffuse sightlines
a = Table.read(data_path / "G24_SMCBumps.dat", format="ascii.commented_header")
a = Table.read(
data_path / "G24_SMCBumps.dat", format="ascii.commented_header"
)

# data
self.obsdata_x = 1.0 / a["wave"].data
Expand Down
Loading

0 comments on commit d8b9e40

Please sign in to comment.