import numpy as np
from scipy.optimize import curve_fit
from astropy import units as u
from planck import planck_function

def bb_flux(wavelength, temperature, angular_radius):
    bb_flux = (np.pi) * planck_function(wavelength, temperature) * (angular_radius)**2

    return bb_flux

def bb_flux_nounits(wavelength, temperature, angular_radius):
    flux = bb_flux(wavelength, temperature, angular_radius)

    return flux.value

def calculate_chisq(y_data, y_data_uncertainties, x_data, func, parameters):
    chisq = np.sum(((y_data - func(x_data, *parameters))/y_data_uncertainties)**2)
    return chisq
    
def bb_fit_parameters(wavelengths, fluxes, flux_uncertainties):
    popt, pcov = curve_fit(bb_flux_nounits, wavelengths, fluxes, p0=[5000, 1.0e-10])
    temperature = popt[0]
    angular_radius = popt[1]
    perr = np.sqrt(np.diag(pcov))
    

    chisq = calculate_chisq(fluxes, flux_uncertainties, wavelengths, bb_flux_nounits, popt)

    return temperature, angular_radius, perr

import numpy as np
from scipy.optimize import curve_fit
from astropy import units as u
from planck import planck_function, dplanck_dT

def bb_flux(wavelength, temperature, angular_radius):
    bb_flux = (np.pi) * planck_function(wavelength, temperature) * (angular_radius)**2

    return bb_flux

def bb_flux_nounits(wavelength, temperature, angular_radius):
    flux = bb_flux(wavelength, temperature, angular_radius)

    return flux.value

def dBB_dT(wavelength, temperature, angular_radius):
    dBB_dT = (np.pi) * dplanck_dT(wavelength, temperature) * (angular_radius)**2

    return dBB_dT

def dBB_dT_nounits(wavelength, temperature, angular_radius):
    dBB_dT_nounits = dBB_dT(wavelength, temperature, angular_radius)

    return dBB_dT_nounits

def calculate_chisq(y_data, y_data_uncertainties, x_data, func, parameters):
    chisq = np.sum(((y_data - func(x_data, *parameters))/y_data_uncertainties)**2)
    return chisq
    
def bb_fit_parameters(wavelengths, fluxes, flux_uncertainties):
    popt, pcov = curve_fit(bb_flux_nounits, wavelengths, fluxes, p0=[5000, 1.0e-10])
    temperature = popt[0]
    angular_radius = popt[1]
    perr = np.sqrt(np.diag(pcov))
    
    return temperature, angular_radius, perr

Add derivative of BB flux with temperature
Add derivative of BB flux with temperature

Added functions which calculate the derivative of the blackbody flux
with temperature (with units and without units)

  import numpy as np
  from scipy.optimize import curve_fit
  from astropy import units as u
- from planck import planck_function
+ from planck import planck_function, dplanck_dT
  
  def bb_flux(wavelength, temperature, angular_radius):
      bb_flux = (np.pi) * planck_function(wavelength, temperature) * (angular_radius)**2
  
      return bb_flux
  
  def bb_flux_nounits(wavelength, temperature, angular_radius):
      flux = bb_flux(wavelength, temperature, angular_radius)
  
      return flux.value
  
+ def dBB_dT(wavelength, temperature, angular_radius):
+     dBB_dT = (np.pi) * dplanck_dT(wavelength, temperature) * (angular_radius)**2
+ 
+     return dBB_dT
+ 
+ def dBB_dT_nounits(wavelength, temperature, angular_radius):
+     dBB_dT_nounits = dBB_dT(wavelength, temperature, angular_radius)
+ 
+     return dBB_dT_nounits
+ 
  def calculate_chisq(y_data, y_data_uncertainties, x_data, func, parameters):
      chisq = np.sum(((y_data - func(x_data, *parameters))/y_data_uncertainties)**2)
      return chisq
      
  def bb_fit_parameters(wavelengths, fluxes, flux_uncertainties):
      popt, pcov = curve_fit(bb_flux_nounits, wavelengths, fluxes, p0=[5000, 1.0e-10])
      temperature = popt[0]
      angular_radius = popt[1]
      perr = np.sqrt(np.diag(pcov))
      
- 
-     chisq = calculate_chisq(fluxes, flux_uncertainties, wavelengths, bb_flux_nounits, popt)
- 
      return temperature, angular_radius, perr
  
Add derivative of BB flux with temperature
## Code Before:
import numpy as np
from scipy.optimize import curve_fit
from astropy import units as u
from planck import planck_function

def bb_flux(wavelength, temperature, angular_radius):
    bb_flux = (np.pi) * planck_function(wavelength, temperature) * (angular_radius)**2

    return bb_flux

def bb_flux_nounits(wavelength, temperature, angular_radius):
    flux = bb_flux(wavelength, temperature, angular_radius)

    return flux.value

def calculate_chisq(y_data, y_data_uncertainties, x_data, func, parameters):
    chisq = np.sum(((y_data - func(x_data, *parameters))/y_data_uncertainties)**2)
    return chisq
    
def bb_fit_parameters(wavelengths, fluxes, flux_uncertainties):
    popt, pcov = curve_fit(bb_flux_nounits, wavelengths, fluxes, p0=[5000, 1.0e-10])
    temperature = popt[0]
    angular_radius = popt[1]
    perr = np.sqrt(np.diag(pcov))
    

    chisq = calculate_chisq(fluxes, flux_uncertainties, wavelengths, bb_flux_nounits, popt)

    return temperature, angular_radius, perr

## Instruction:
Add derivative of BB flux with temperature
## Code After:
import numpy as np
from scipy.optimize import curve_fit
from astropy import units as u
from planck import planck_function, dplanck_dT

def bb_flux(wavelength, temperature, angular_radius):
    bb_flux = (np.pi) * planck_function(wavelength, temperature) * (angular_radius)**2

    return bb_flux

def bb_flux_nounits(wavelength, temperature, angular_radius):
    flux = bb_flux(wavelength, temperature, angular_radius)

    return flux.value

def dBB_dT(wavelength, temperature, angular_radius):
    dBB_dT = (np.pi) * dplanck_dT(wavelength, temperature) * (angular_radius)**2

    return dBB_dT

def dBB_dT_nounits(wavelength, temperature, angular_radius):
    dBB_dT_nounits = dBB_dT(wavelength, temperature, angular_radius)

    return dBB_dT_nounits

def calculate_chisq(y_data, y_data_uncertainties, x_data, func, parameters):
    chisq = np.sum(((y_data - func(x_data, *parameters))/y_data_uncertainties)**2)
    return chisq
    
def bb_fit_parameters(wavelengths, fluxes, flux_uncertainties):
    popt, pcov = curve_fit(bb_flux_nounits, wavelengths, fluxes, p0=[5000, 1.0e-10])
    temperature = popt[0]
    angular_radius = popt[1]
    perr = np.sqrt(np.diag(pcov))
    
    return temperature, angular_radius, perr
