'''
how to write a code in Python to fit dust SED using a modified blackbody model? 
'''
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Define the modified blackbody function
def modified_blackbody(wavelength, T, beta, nu0, A):
    """
    Modified blackbody function to fit a dust SED.
    Parameters:
    wavelength: array of wavelengths in microns
    T: dust temperature in Kelvin
    beta: dust emissivity index
    nu0: reference frequency in Hz
    A: amplitude of the dust SED
    """
    nu = 3e8 / (wavelength * 1e-6)  # convert wavelength to frequency in Hz
    term1 = (nu / nu0)**(beta+1)
    term2 = np.exp(2.8e-11 * nu0 / T) - 1
    return A * term1 * term2

# Load the dust SED data
wavelength, flux, error = np.loadtxt('dust_sed.txt', unpack=True)

# Plot the dust SED data
plt.errorbar(wavelength, flux, yerr=error, fmt='o', label='data')
plt.xlabel('Wavelength (microns)')
plt.ylabel('Flux (Jy)')
plt.xscale('log')
plt.yscale('log')

# Initial guess for the fit parameters
T_guess = 20  # Kelvin
beta_guess = 1.5
nu0_guess = 2.5e12  # Hz
A_guess = 1e-6  # Jy

# Fit the modified blackbody model to the dust SED data using curve_fit
popt, pcov = curve_fit(modified_blackbody, wavelength, flux, sigma=error,
                       p0=[T_guess, beta_guess, nu0_guess, A_guess])

# Print the best-fit parameters and their errors
print('T = {:.2f} +/- {:.2f} K'.format(popt[0], np.sqrt(pcov[0, 0])))
print('beta = {:.2f} +/- {:.2f}'.format(popt[1], np.sqrt(pcov[1, 1])))
print('nu0 = {:.2e} +/- {:.2e} Hz'.format(popt[2], np.sqrt(pcov[2, 2])))
print('A = {:.2e} +/- {:.2e} Jy'.format(popt[3], np.sqrt(pcov[3, 3])))

# Plot the best-fit model and the data
wavelength_model = np.logspace(np.log10(wavelength.min()), np.log10(wavelength.max()), 100)
flux_model = modified_blackbody(wavelength_model, *popt)
plt.plot(wavelength_model, flux_model, label='model')
plt.legend()

plt.show()