diff --git a/docs/dust_extinction/fit_extinction.rst b/docs/dust_extinction/fit_extinction.rst index 35afb25..5bc6751 100644 --- a/docs/dust_extinction/fit_extinction.rst +++ b/docs/dust_extinction/fit_extinction.rst @@ -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 ================ @@ -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() diff --git a/docs/index.rst b/docs/index.rst index 6f332d7..d077803 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -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. diff --git a/dust_extinction/averages.py b/dust_extinction/averages.py index 197aa12..0135b8c 100644 --- a/dust_extinction/averages.py +++ b/dust_extinction/averages.py @@ -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') @@ -295,7 +295,7 @@ class B92_MWAvg(BaseExtModel): def evaluate(self, in_x): """ - B1992 function + B92 function Parameters ---------- @@ -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() """ @@ -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() """ @@ -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() """ @@ -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() """ @@ -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 @@ -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() """ @@ -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 @@ -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() """ @@ -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 diff --git a/dust_extinction/parameter_averages.py b/dust_extinction/parameter_averages.py index 7434ab8..6b1023e 100644 --- a/dust_extinction/parameter_averages.py +++ b/dust_extinction/parameter_averages.py @@ -70,6 +70,16 @@ class CCM89(BaseExtRvModel): 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() """ @@ -196,6 +206,16 @@ class O94(BaseExtRvModel): 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() """ @@ -325,6 +345,16 @@ class F99(BaseExtRvModel): 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() """ @@ -464,6 +494,16 @@ class F04(BaseExtRvModel): 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() """ @@ -599,6 +639,16 @@ class VCG04(BaseExtRvModel): 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.12, 0.15, 0.2, 0.3]) + 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() """ @@ -709,6 +759,16 @@ class GCC09(BaseExtRvModel): 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]) + 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() """ @@ -834,6 +894,16 @@ class M14(BaseExtRvModel): ext_model = M14(Rv=cur_Rv) ax.plot(x,ext_model(x),label='R(V) = ' + str(cur_Rv)) + # for 2nd x-axis with lambda values + axis_xs = np.array([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.set_xlabel(r'$x$ [$\mu m^{-1}$]') ax.set_ylabel(r'$A(x)/A(V)$') @@ -1029,6 +1099,16 @@ class G16(BaseExtRvAfAModel): 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', title=r'$f_A = 1.0$') plt.show() @@ -1058,6 +1138,16 @@ class G16(BaseExtRvAfAModel): 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', title=r'$R_A(V) = 3.1$') plt.show() """ @@ -1165,6 +1255,16 @@ class F19(BaseExtRvModel): 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() """ diff --git a/dust_extinction/shapes.py b/dust_extinction/shapes.py index 33f3cee..6a7a74c 100644 --- a/dust_extinction/shapes.py +++ b/dust_extinction/shapes.py @@ -227,14 +227,27 @@ class FM90(Fittable1DModel): ax.set_xlabel(r'$x$ [$\mu m^{-1}$]') ax.set_ylabel(r'$E(\lambda - V)/E(B - V)$') + # for 2nd x-axis with lambda values + axis_xs = np.array([0.12, 0.15, 0.2, 0.3]) + 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() """ + n_inputs = 1 n_outputs = 1 # bounds based on Gordon et al. (2024) results - C1 = Parameter(description="linear term: y-intercept", default=0.10, bounds=(-10.0, 5.0)) + C1 = Parameter( + description="linear term: y-intercept", default=0.10, bounds=(-10.0, 5.0) + ) C2 = Parameter(description="linear term: slope", default=0.70, bounds=(-0.1, 5.0)) C3 = Parameter(description="bump: amplitude", default=3.23, bounds=(-1.0, 6.0)) C4 = Parameter(description="FUV rise: amplitude", default=0.41, bounds=(-0.5, 1.5)) @@ -395,14 +408,27 @@ class FM90_B3(Fittable1DModel): ax.set_xlabel(r'$x$ [$\mu m^{-1}$]') ax.set_ylabel(r'$E(\lambda - V)/E(B - V)$') + # for 2nd x-axis with lambda values + axis_xs = np.array([0.12, 0.15, 0.2, 0.3]) + 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() """ + n_inputs = 1 n_outputs = 1 # bounds based on Gordon et al. (2024) results - C1 = Parameter(description="linear term: y-intercept", default=0.10, bounds=(-10.0, 5.0)) + C1 = Parameter( + description="linear term: y-intercept", default=0.10, bounds=(-10.0, 5.0) + ) C2 = Parameter(description="linear term: slope", default=0.70, bounds=(-0.1, 5.0)) B3 = Parameter(description="bump: amplitude", default=3.23, bounds=(-1.0, 6.0)) C4 = Parameter(description="FUV rise: amplitude", default=0.41, bounds=(-0.5, 1.5)) @@ -592,6 +618,7 @@ class P92(Fittable1DModel): ax.legend(loc='best') plt.show() """ + n_inputs = 1 n_outputs = 1