Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change variable names in the gaussian fittings #219

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Change variable names in the gaussian fittings #219

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9cbf240
Update qlook.py (issue218)
ShinjiFujita Oct 8, 2024
90d17e4
Update qlook.py (issue218)
ShinjiFujita Oct 8, 2024
f415716
Update qlook.py (fixed typos "::" -> ":")
ShinjiFujita Dec 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 56 additions & 14 deletions decode/qlook.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -247,8 +247,15 @@ def daisy(
)
perr = np.sqrt(np.diag(pcov))
data_fitted = gaussian_2d((x, y), *popt).reshape(x.shape)
chi2, reduced_chi2 = calc_chi2(
data, data_fitted, sigma=1.0, num_params=len(initial_guess)
)
fit_res_params_dict = make_fit_res_params_dict(
popt, perr, chi2, reduced_chi2
)
is_gaussfit_successful = True
except:
except Exception as e:
print(f"An error occurred: {e}")
is_gaussfit_successful = False

# save result
Expand Down Expand Up @@ -280,12 +287,12 @@ def daisy(
linestyles="-",
)
ax.set_title(
f"Peak = {popt[0]:+.2e} [{cont.units}], "
f"dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], "
f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}],\n"
f"FWHM_x = {popt[3]*popt[4]*2.354820:+.2f} [{skycoord_units}], "
f"FWHM_y = {popt[4]*2.354820:+.2f} [{skycoord_units}], "
f"P.A. = {-np.rad2deg(popt[5]+np.pi/2.0):+.1f} [deg],\n"
f"Peak = {fit_res_params_dict['peak']:+.2e} [{cont.units}], "
f"dAz = {fit_res_params_dict['offset_az']:+.2f} [{cont.lon.attrs['units']}], "
f"dEl = {fit_res_params_dict['offset_el']:+.2f} [{cont.lat.attrs['units']}],\n"
f"FWHM_maj = {fit_res_params_dict['hpbw_major']:+.2f} [{skycoord_units}], "
f"FWHM_min = {fit_res_params_dict['hpbw_minor']:+.2f} [{skycoord_units}], "
f"P.A. = {fit_res_params_dict['position_angle']:+.1f} [deg],\n"
f"min_frequency = {min_frequency}, "
f"max_frequency = {max_frequency}",
fontsize=10,
Expand Down Expand Up @@ -518,8 +525,15 @@ def raster(
)
perr = np.sqrt(np.diag(pcov))
data_fitted = gaussian_2d((x, y), *popt).reshape(x.shape)
chi2, reduced_chi2 = calc_chi2(
data, data_fitted, sigma=1.0, num_params=len(initial_guess)
)
fit_res_params_dict = make_fit_res_params_dict(
popt, perr, chi2, reduced_chi2
)
is_gaussfit_successful = True
except:
except Exception as e:
print(f"An error occurred: {e}")
is_gaussfit_successful = False

# save result
Expand Down Expand Up @@ -551,12 +565,12 @@ def raster(
linestyles="-",
)
ax.set_title(
f"Peak = {popt[0]:+.2e} [{cont.units}], "
f"dAz = {popt[1]:+.2f} [{cont.lon.attrs['units']}], "
f"dEl = {popt[2]:+.2f} [{cont.lat.attrs['units']}],\n"
f"FWHM_x = {popt[3]*popt[4]*2.354820:+.2f} [{skycoord_units}], "
f"FWHM_y = {popt[4]*2.354820:+.2f} [{skycoord_units}], "
f"P.A. = {-np.rad2deg(popt[5]+np.pi/2.0):+.1f} [deg],\n"
f"Peak = {fit_res_params_dict['peak']:+.2e} [{cont.units}], "
f"dAz = {fit_res_params_dict['offset_az']:+.2f} [{cont.lon.attrs['units']}], "
f"dEl = {fit_res_params_dict['offset_el']:+.2f} [{cont.lat.attrs['units']}],\n"
f"FWHM_maj = {fit_res_params_dict['hpbw_major']:+.2f} [{skycoord_units}], "
f"FWHM_min = {fit_res_params_dict['hpbw_minor']:+.2f} [{skycoord_units}], "
f"P.A. = {fit_res_params_dict['position_angle']:+.1f} [deg],\n"
f"min_frequency = {min_frequency}, "
f"max_frequency = {max_frequency}",
fontsize=10,
Expand Down Expand Up @@ -1369,6 +1383,34 @@ def gaussian_2d(xy, amp, x0, y0, sigma_x_over_y, sigma_y, theta, offset):
return g.ravel()


def calc_chi2(data_obs, data_fit, sigma, num_params):
chi2 = np.sum(((data_obs - data_fit) / sigma) ** 2)
dof = len(data_obs) - num_params
reduced_chi2 = chi2 / dof
return chi2, reduced_chi2


def make_fit_res_params_dict(popt, perr, chi2, reduced_chi2):
res = {}
res["offset_az"] = popt[1]
res["offset_el"] = popt[2]
res["offset_az_error"] = perr[1]
res["offset_el_error"] = perr[2]
res["hpbw_major"] = popt[3] * popt[4] * 2.354820
res["hpbw_minor"] = popt[4] * 2.354820
res["hpbw_major_error"] = popt[3] * perr[4] * 2.354820
res["hpbw_minor_error"] = perr[4] * 2.354820
res["position_angle"] = -np.rad2deg(popt[5] + np.pi / 2.0)
res["position_angle_error"] = np.rad2deg(perr[5])
res["peak"] = popt[0]
res["peak_error"] = perr[0]
res["floor"] = popt[6]
res["floor_error"] = perr[6]
res["chi2"] = chi2
res["reduced_chi2"] = reduced_chi2
return res


def main() -> None:
"""Entry point of the decode-qlook command."""

Expand Down