Haipeng Li Lab 6 python code (1)

pdf

School

University of California, Irvine *

*We aren’t endorsed by this school

Course

115A

Subject

Electrical Engineering

Date

Jan 9, 2024

Type

pdf

Pages

Uploaded by DeanResolve98810

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 1/9 In [50]: import numpy as np import matplotlib.pyplot as plt from astropy.io import fits image_file = 'corrected_s0311.fits' image_data = fits.getdata(image_file) x1, x2, y1, y2 = 100 , 150 , 100 , 150 star_region = image_data[y1:y2, x1:x2] profile_x = np.sum(star_region, axis = 0 ) profile_y = np.sum(star_region, axis = 1 ) peak_x = np.max(profile_x) half_max_x = peak_x / 2 peak_y = np.max(profile_y) half_max_y = peak_y / 2 fwhm_x = np.sum(profile_x > half_max_x) fwhm_y = np.sum(profile_y > half_max_y) plt.figure(figsize = ( 12 , 6 )) plt.subplot( 1 , 2 , 1 ) plt.plot(profile_x) plt.axhline(y = half_max_x, color = 'r' , linestyle = '--' ) plt.title( f'X Profile (FWHM: {fwhm_x} pixels)' ) plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.subplot( 1 , 2 , 2 ) plt.plot(profile_y) plt.axhline(y = half_max_y, color = 'r' , linestyle = '--' ) plt.title( f'Y Profile (FWHM: {fwhm_y} pixels)' ) plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.tight_layout() plt.show()

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 2/9 In [29]: Model: CompoundModel Inputs: ('x',) Outputs: ('y',) Model set size: 1 Expression: [0] + [1] Components: [0]: <Gaussian1D(amplitude=3.05394402, mean=1.27426551, stddev=0.81712 989)> [1]: <Const1D(amplitude=0.99037165)> Parameters: amplitude_0 mean_0 stddev_0 amplitude_1 ------------------ ---------------- ------------------ --------------- --- 3.0539440238047133 1.27426550645335 0.8171298855086331 0.9903716523469 903 Computed FWHM: 1.9243408803728308 from astropy.modeling import models, fitting np.random.seed( 0 ) xdata = np.linspace( - 5. , 5. , 200 ) ydata = 3 * np.exp( - 0.5 * (xdata - 1.3 ) ** 2 / 0.8 ** 2 ) + 1.0 ydata += np.random.normal( 0. , 0.2 , xdata.shape) g_init = models.Gaussian1D(amplitude = 1. , mean = 0 , stddev = 1. ) + models.Const1 fit_g = fitting.LevMarLSQFitter() g = fit_g(g_init, xdata, ydata) plt.figure(figsize = ( 8 , 5 )) plt.plot(xdata, ydata, 'ko' , label = 'Data' ) plt.plot(xdata, g(xdata), label = 'Gaussian + Constant' ) plt.xlabel( 'Pixel' ) plt.ylabel( 'Counts' ) plt.legend(loc = 2 ) plt.show() print (g) fwhm_computed = 2.355 * g.stddev_0 print ( "Computed FWHM:" , fwhm_computed)

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 3/9 In [30]: x = np.arange( 0 , 50 ) profile_x = np.sum(star_region, axis = 0 ) profile_y = np.sum(star_region, axis = 1 ) fit_g = fitting.LevMarLSQFitter() g_init_x = models.Gaussian1D(amplitude = np.max(profile_x), mean = np.argmax(pro g_x = fit_g(g_init_x, x, profile_x) g_init_y = models.Gaussian1D(amplitude = np.max(profile_y), mean = np.argmax(pro g_y = fit_g(g_init_y, x, profile_y) plt.figure(figsize = ( 12 , 6 )) plt.subplot( 1 , 2 , 1 ) plt.plot(x, profile_x, 'ko' , label = 'X Profile Data' ) plt.plot(x, g_x(x), label = 'Gaussian Fit' , linestyle = '--' ) plt.axhline(y = g_x.amplitude.value / 2 , color = 'r' , linestyle = '--' ) plt.title( f'X Profile and Gaussian Fit\nFWHM: { 2.355 * g_x.stddev.value: .2 f plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.legend() plt.subplot( 1 , 2 , 2 ) plt.plot(x, profile_y, 'ko' , label = 'Y Profile Data' ) plt.plot(x, g_y(x), label = 'Gaussian Fit' , linestyle = '--' ) plt.axhline(y = g_y.amplitude.value / 2 , color = 'r' , linestyle = '--' ) plt.title( f'Y Profile and Gaussian Fit\nFWHM: { 2.355 * g_y.stddev.value: .2 f plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.legend() plt.tight_layout() plt.show()

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 4/9 In [60]: another_image_file = 's0576.fits' another_image_data = fits.getdata(another_image_file) another_x1, another_x2, another_y1, another_y2 = 100 , 150 , 100 , 150 another_star_region = another_image_data[another_y1:another_y2, another_x1: another_profile_x = np.sum(another_star_region, axis = 0 ) another_profile_y = np.sum(another_star_region, axis = 1 ) another_peak_x = np.max(another_profile_x) another_half_max_x = another_peak_x / 2 another_peak_y = np.max(another_profile_y) another_half_max_y = another_peak_y / 2 another_fwhm_x = np.sum(another_profile_x > another_half_max_x) another_fwhm_y = np.sum(another_profile_y > another_half_max_y) plt.figure(figsize = ( 12 , 6 )) plt.subplot( 1 , 2 , 1 ) plt.plot(another_profile_x) plt.axhline(y = another_half_max_x, color = 'r' , linestyle = '--' ) plt.title( f'X Profile for {another_image_file} (FWHM: {another_fwhm_x} pixe plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.subplot( 1 , 2 , 2 ) plt.plot(another_profile_y) plt.axhline(y = another_half_max_y, color = 'r' , linestyle = '--' ) plt.title( f'Y Profile for {another_image_file} (FWHM: {another_fwhm_y} pixe plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.tight_layout() plt.show()

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 5/9 In [61]: another_image_file = 's0586.fits' another_image_data = fits.getdata(another_image_file) another_x1, another_x2, another_y1, another_y2 = 100 , 150 , 100 , 150 another_star_region = another_image_data[another_y1:another_y2, another_x1: another_profile_x = np.sum(another_star_region, axis = 0 ) another_profile_y = np.sum(another_star_region, axis = 1 ) another_peak_x = np.max(another_profile_x) another_half_max_x = another_peak_x / 2 another_peak_y = np.max(another_profile_y) another_half_max_y = another_peak_y / 2 another_fwhm_x = np.sum(another_profile_x > another_half_max_x) another_fwhm_y = np.sum(another_profile_y > another_half_max_y) plt.figure(figsize = ( 12 , 6 )) plt.subplot( 1 , 2 , 1 ) plt.plot(another_profile_x) plt.axhline(y = another_half_max_x, color = 'r' , linestyle = '--' ) plt.title( f'X Profile for {another_image_file} (FWHM: {another_fwhm_x} pixe plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.subplot( 1 , 2 , 2 ) plt.plot(another_profile_y) plt.axhline(y = another_half_max_y, color = 'r' , linestyle = '--' ) plt.title( f'Y Profile for {another_image_file} (FWHM: {another_fwhm_y} pixe plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.tight_layout() plt.show()

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 6/9 In [56]: x = np.arange( 0 , 50 ) another_profile_x = np.sum(another_star_region, axis = 0 ) another_profile_y = np.sum(another_star_region, axis = 1 ) fit_g = fitting.LevMarLSQFitter() g_init_x = models.Gaussian1D(amplitude = np.max(another_profile_x), mean = np.a g_x = fit_g(g_init_x, x, another_profile_x) g_init_y = models.Gaussian1D(amplitude = np.max(another_profile_y), mean = np.a g_y = fit_g(g_init_y, x, another_profile_y) plt.figure(figsize = ( 12 , 6 )) plt.subplot( 1 , 2 , 1 ) plt.plot(x, another_profile_x, 'ko' , label = 'X Profile Data' ) plt.plot(x, g_x(x), label = 'Gaussian Fit' , linestyle = '--' ) plt.axhline(y = g_x.amplitude.value / 2 , color = 'r' , linestyle = '--' ) plt.title( f'X Profile and Gaussian Fit\nFWHM: { 2.355 * g_x.stddev.value: .2 f plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.legend() plt.subplot( 1 , 2 , 2 ) plt.plot(x, another_profile_y, 'ko' , label = 'Y Profile Data' ) plt.plot(x, g_y(x), label = 'Gaussian Fit' , linestyle = '--' ) plt.axhline(y = g_y.amplitude.value / 2 , color = 'r' , linestyle = '--' ) plt.title( f'Y Profile and Gaussian Fit\nFWHM: { 2.355 * g_y.stddev.value: .2 f plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.legend() plt.tight_layout() plt.show()

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version

Access to all documents
Unlimited textbook solutions
24/7 expert homework help

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 7/9 In [57]: final_image_file = 'final_image.fits' final_image_data = fits.getdata(final_image_file) final_star_region = final_image_data[y1:y2, x1:x2] final_profile_x = np.sum(final_star_region, axis = 0 ) final_profile_y = np.sum(final_star_region, axis = 1 ) final_peak_x = np.max(final_profile_x) final_half_max_x = final_peak_x / 2 final_peak_y = np.max(final_profile_y) final_half_max_y = final_peak_y / 2 final_fwhm_x = np.sum(final_profile_x > final_half_max_x) final_fwhm_y = np.sum(final_profile_y > final_half_max_y) plt.figure(figsize = ( 12 , 6 )) plt.subplot( 1 , 2 , 1 ) plt.plot(final_profile_x) plt.axhline(y = final_half_max_x, color = 'r' , linestyle = '--' ) plt.title( f'X Profile for {final_image_file} (FWHM: {final_fwhm_x} pixels)' plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.subplot( 1 , 2 , 2 ) plt.plot(final_profile_y) plt.axhline(y = final_half_max_y, color = 'r' , linestyle = '--' ) plt.title( f'Y Profile for {final_image_file} (FWHM: {final_fwhm_y} pixels)' plt.xlabel( 'Pixel Position' ) plt.ylabel( 'Intensity' ) plt.tight_layout() plt.show()

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 8/9 In [66]: from astropy.io import fits from photutils.aperture import CircularAperture, CircularAnnulus, aperture_ image_file = 'final_image.fits' image_data = fits.getdata(image_file) x_star, y_star = 1000 , 1000 radii = np.arange( 5 , 50 , 5 ) total_fluxes = [] sky_subtracted_fluxes = [] for r in radii: aperture = CircularAperture((x_star, y_star), r) annulus_aperture = CircularAnnulus((x_star, y_star), r_in = r + 5 , r_out = r + flux_table = aperture_photometry(image_data, aperture) annulus_table = aperture_photometry(image_data, annulus_aperture) mean_sky_per_pixel = annulus_table[ 'aperture_sum' ] / annulus_aperture.a total_flux = flux_table[ 'aperture_sum' ] - mean_sky_per_pixel * aperture total_fluxes.append(flux_table[ 'aperture_sum' ][ 0 ]) sky_subtracted_fluxes.append(total_flux[ 0 ]) plt.figure(figsize = ( 10 , 6 )) plt.plot(radii, total_fluxes, label = 'Total Counts' ) plt.plot(radii, sky_subtracted_fluxes, label = 'Sky-Subtracted Counts' ) plt.xlabel( 'Aperture Radius (pixels)' ) plt.ylabel( 'Counts' ) plt.title( 'Curve-of-Growth' ) plt.legend() plt.show()

12/3/23, 1:59 AM Untitled81 - Jupyter Notebook localhost:8888/notebooks/Untitled81.ipynb 9/9 In [70]: In [ ]: Model: CompoundModel Inputs: ('x', 'y') Outputs: ('z',) Model set size: 1 Expression: [0] + [1] Components: [0]: <Gaussian2D(amplitude=-0.40278672, x_mean=2.90962142, y_mean=-1.5 7473144, x_stddev=0., y_stddev=0.87499035, theta=0.)> [1]: <Const2D(amplitude=0.25890895)> Parameters: amplitude_0 x_mean_0 ... theta_0 amplitude_1 -------------------- ------------------ ... ------- ------------------ -0.40278672145522876 2.9096214203514905 ... 0.0 0.2589089536482877 Total counts (per second): 2.603027904683625e-38 image_file = 'final_image.fits' image_data = fits.getdata(image_file) image_data[ ~ np.isfinite(image_data)] = 0 g_init = models.Gaussian2D(amplitude = 1. , x_mean = 0. , y_mean = 0. , x_stddev = 1. , y, x = np.mgrid[:image_data.shape[ 0 ], :image_data.shape[ 1 ]] fit_g = fitting.LevMarLSQFitter() with warnings.catch_warnings(): warnings.simplefilter( 'ignore' ) g = fit_g(g_init, x, y, image_data) print (g) sky_level = g.amplitude_1.value star_data = image_data - sky_level total_counts = g.amplitude_0.value * 2 * np.pi * g.x_stddev_0.value * g.y_st exposure_time = 1 total_counts_per_second = total_counts / exposure_time print ( "Total counts (per second):" , abs (total_counts_per_second))

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version