From d04a665c37fbaa74a772d8b0771a8d9dc917b5a0 Mon Sep 17 00:00:00 2001 From: EyanWeissbluth <133696906+EyanWeissbluth@users.noreply.github.com> Date: Wed, 19 Jul 2023 22:15:35 -0700 Subject: [PATCH 1/4] Update analysis.py --- src/chart/analysis.py | 155 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 155 insertions(+) diff --git a/src/chart/analysis.py b/src/chart/analysis.py index 558c4d3..4f414d8 100644 --- a/src/chart/analysis.py +++ b/src/chart/analysis.py @@ -66,3 +66,158 @@ def read_run(directory=None): def concat(data_list): raise NotImplementedError() + + +def LSR_shift(lon, lat, ele, time, RA, Dec): + """ + Identifies the exact postion at which the observations were taken and corrects for the Local Standard of Rest. + + :param lon: longitude + :param lat: latitude + :param ele: elevation + :param time: date and time + :param RA: Right Ascension + :param Dec: Declination + """ + + location = EarthLocation.from_geodetic(lon, lat, ele*u.m) #Lon, Lat, elevation + location = location.get_itrs(obstime=Time(time)) #To ITRS frame, makes Earth stationary with Sun + pointing_45deg = SkyCoord(RA,Dec, frame='icrs') #Center of CHART pointing + frequency = SpectralCoord(1.420405751768e9 * u.Hz, observer=location, target=pointing_45deg) #Shift expected from just local motion + f0_shifted = frequency.with_observer_stationary_relative_to('lsrk') #correct for kinematic local standard of rest + f0_shifted = f0_shifted.to(u.GHz) + v = doppler(f0_shifted,f0) + v_adjustment = v.to(u.km/u.second) + print(v_adjustment) + return v_adjustment + + +def average_overlapping(x1, y1, x2, y2): + """ + Averages the y values where the x values are shared between two arrays and keeps y values for x values that are not shared. + + :param x1: First x array + :param y1: First y array + :param x2: Second x array + :param y2: Second y array + :return: Tuple of combined x values and averaged/kept y values + """ + # Find the unique x values in both arrays + unique_x = np.union1d(x1, x2) + + # Create an array to store the averaged/kept y values + avg_y = np.zeros_like(unique_x) + + # Iterate over the unique x values + for i in range(len(unique_x)): + # Find the indices of the current x value in the two x arrays + ind1 = np.where(x1 == unique_x[i])[0] + ind2 = np.where(x2 == unique_x[i])[0] + + # If the current x value is in both arrays + if len(ind1) > 0 and len(ind2) > 0: + # Compute the average of the two corresponding y values + avg_y[i] = (y1[ind1[0]] + y2[ind2[0]]) / 2 + # If the current x value is only in the first array + elif len(ind1) > 0: + # Keep the corresponding y value from the first array + avg_y[i] = y1[ind1[0]] + # If the current x value is only in the second array + elif len(ind2) > 0: + # Keep the corresponding y value from the second array + avg_y[i] = y2[ind2[0]] + + return unique_x, avg_y + + +def interactive_plot(unique_x): + """ + Creates a plot that can be modified with sliders. + + param unique_x: x values of overlapping CHART data + """ + x = unique_x + e = 2.71828 + fig = plt.figure() + ax = fig.add_subplot(1, 1, 1) + a1 = a2 = a3 = a4 = b1 = b2 = b3 = b4 = c1 = c2 = c3 = c4 = 1 + + line1, = ax.plot(x, (a1*(e ** -(((x-b1)**2) / (2*(c1**2)))))) + line2, = ax.plot(x, (a2*(e ** -(((x-b2)**2) / (2*(c2**2)))))) + line3, = ax.plot(x, (a3*(e ** -(((x-b3)**2) / (2*(c3**2)))))) + line4, = ax.plot(x, (a4*(e ** -(((x-b4)**2) / (2*(c4**2)))))) + line5, = ax.plot(x, ((a1*(e ** -(((x-b1)**2) / (2*(c1**2))))) + + (a2*(e ** -(((x-b2)**2) / (2*(c2**2))))) + + (a3*(e ** -(((x-b3)**2) / (2*(c3**2))))) + + (a4*(e ** -(((x-b4)**2) / (2*(c4**2))))))) + + a1_slider = FloatSlider(min=-100, max=100, step=1, value=1) + a2_slider = FloatSlider(min=-100, max=100, step=1, value=1) + a3_slider = FloatSlider(min=-100, max=100, step=1, value=1) + a4_slider = FloatSlider(min=-100, max=100, step=1, value=1) + b1_slider = FloatSlider(min=-100, max=100, step=1, value=1) + b2_slider = FloatSlider(min=-100, max=100, step=1, value=1) + b3_slider = FloatSlider(min=-100, max=100, step=1, value=1) + b4_slider = FloatSlider(min=-100, max=100, step=1, value=1) + c1_slider = FloatSlider(min=-100, max=100, step=1, value=1) + c2_slider = FloatSlider(min=-100, max=100, step=1, value=1) + c3_slider = FloatSlider(min=-100, max=100, step=1, value=1) + c4_slider = FloatSlider(min=-100, max=100, step=1, value=1) + + color1_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='black') + color2_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='black') + color3_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='black') + color4_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='black') + color5_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='red') + + def update(a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, a3=1, b3=1, c3=1, a4=1, b4=1, c4=1, + color1='black', color2='black', color3='black', color4='black', color5='red'): + line1.set_ydata((a1*(e ** -(((x-b1)**2) / (2*(c1**2)))))) + line2.set_ydata((a2*(e ** -(((x-b2)**2) / (2*(c2**2)))))) + line3.set_ydata((a3*(e ** -(((x-b3)**2) / (2*(c3**2)))))) + line4.set_ydata((a4*(e ** -(((x-b4)**2) / (2*(c4**2)))))) + line5.set_ydata(((a1*(e ** -(((x-b1)**2) / (2*(c1**2))))) + + (a2*(e ** -(((x-b2)**2) / (2*(c2**2))))) + + (a3*(e ** -(((x-b3)**2) / (2*(c3**2))))) + + (a4*(e ** -(((x-b4)**2) / (2*(c4**2))))))) + line1.set_color(color1) + line2.set_color(color2) + line3.set_color(color3) + line4.set_color(color4) + line5.set_color(color5) + fig.canvas.draw_idle() + + interact(update, a1=a1_slider, b1=b1_slider, c1=c1_slider, a2=a2_slider, b2=b2_slider, c2=c2_slider, + a3=a3_slider, b3=b3_slider, c3=c3_slider, a4=a4_slider, b4=b4_slider, c4=c4_slider, color1=color1_dropdown, + color2=color2_dropdown, color3=color3_dropdown, color4=color4_dropdown, color5=color5_dropdown); + + return ax + + +def goodness_of_fit(unique_x, combined_gauss, avg_y): + """Performs a chi-squared goodness of fit test between the CHART data and the user created combined Gaussian curve. + + param unique_x: x values of overlapping CHART data + param combined_gauss: y values of combined Gaissian curve + param avg: y values of overlapping CHART data + """ + x_gauss = np.array(unique_x) + y_gauss = np.array(combined_gauss) + x_data = np.array(unique_x) + y_data = np.array(avg_y) + + gauss = np.concatenate((x_gauss.reshape(-1,1), y_gauss.reshape(-1,1)), axis=1) + data = np.concatenate((x_data.reshape(-1,1), y_data.reshape(-1,1)), axis=1) + + + mask = (gauss[:,0] >= -100) & (gauss[:,0] <= 100) + gauss_masked = gauss[mask] + data_masked = data[mask] + + chi_squared_statistic_x = np.sum((gauss_masked[:,0] - data_masked[:,0])**2 / data_masked[:,0]) + chi_squared_statistic_y = np.sum((y_gauss - y_data)**2 / y_data) + + p_value_x = chi2.sf(chi_squared_statistic_x, len(gauss_masked[:,0]) - 1) + p_value_y = chi2.sf(chi_squared_statistic_y, len(y_gauss) - 1) + + print("P-value for y: ", p_value_y) From e4a9ffd58a9f3daea48510189f514fac4ce06cd7 Mon Sep 17 00:00:00 2001 From: EyanWeissbluth <133696906+EyanWeissbluth@users.noreply.github.com> Date: Tue, 25 Jul 2023 13:58:08 -0700 Subject: [PATCH 2/4] Update analysis.py --- src/chart/analysis.py | 89 +++++++++++++++++++------------------------ 1 file changed, 39 insertions(+), 50 deletions(-) diff --git a/src/chart/analysis.py b/src/chart/analysis.py index 4f414d8..d6f3807 100644 --- a/src/chart/analysis.py +++ b/src/chart/analysis.py @@ -70,7 +70,8 @@ def concat(data_list): def LSR_shift(lon, lat, ele, time, RA, Dec): """ - Identifies the exact postion at which the observations were taken and corrects for the Local Standard of Rest. + Identifies the exact postion at which the observations were taken and corrects for the Local Standard of Rest. + Along with this it also converts the celestial coordinates to galactic coordinates. :param lon: longitude :param lat: latitude @@ -89,9 +90,15 @@ def LSR_shift(lon, lat, ele, time, RA, Dec): v = doppler(f0_shifted,f0) v_adjustment = v.to(u.km/u.second) print(v_adjustment) + + + coor=SkyCoord(RA, Dec, frame='icrs') + l=coor.galactic.l.degree + b=coor.galactic.b.degree + print(l,b) + return v_adjustment - def average_overlapping(x1, y1, x2, y2): """ Averages the y values where the x values are shared between two arrays and keeps y values for x values that are not shared. @@ -140,60 +147,39 @@ def interactive_plot(unique_x): e = 2.71828 fig = plt.figure() ax = fig.add_subplot(1, 1, 1) - a1 = a2 = a3 = a4 = b1 = b2 = b3 = b4 = c1 = c2 = c3 = c4 = 1 - - line1, = ax.plot(x, (a1*(e ** -(((x-b1)**2) / (2*(c1**2)))))) - line2, = ax.plot(x, (a2*(e ** -(((x-b2)**2) / (2*(c2**2)))))) - line3, = ax.plot(x, (a3*(e ** -(((x-b3)**2) / (2*(c3**2)))))) - line4, = ax.plot(x, (a4*(e ** -(((x-b4)**2) / (2*(c4**2)))))) - line5, = ax.plot(x, ((a1*(e ** -(((x-b1)**2) / (2*(c1**2))))) + - (a2*(e ** -(((x-b2)**2) / (2*(c2**2))))) + - (a3*(e ** -(((x-b3)**2) / (2*(c3**2))))) + - (a4*(e ** -(((x-b4)**2) / (2*(c4**2))))))) - - a1_slider = FloatSlider(min=-100, max=100, step=1, value=1) - a2_slider = FloatSlider(min=-100, max=100, step=1, value=1) - a3_slider = FloatSlider(min=-100, max=100, step=1, value=1) - a4_slider = FloatSlider(min=-100, max=100, step=1, value=1) - b1_slider = FloatSlider(min=-100, max=100, step=1, value=1) - b2_slider = FloatSlider(min=-100, max=100, step=1, value=1) - b3_slider = FloatSlider(min=-100, max=100, step=1, value=1) - b4_slider = FloatSlider(min=-100, max=100, step=1, value=1) - c1_slider = FloatSlider(min=-100, max=100, step=1, value=1) - c2_slider = FloatSlider(min=-100, max=100, step=1, value=1) - c3_slider = FloatSlider(min=-100, max=100, step=1, value=1) - c4_slider = FloatSlider(min=-100, max=100, step=1, value=1) - - color1_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='black') - color2_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='black') - color3_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='black') - color4_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='black') - color5_dropdown = Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value='red') - - def update(a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, a3=1, b3=1, c3=1, a4=1, b4=1, c4=1, + a = b = c = [1]*4 + + lines = [ax.plot(x, (a[i]*(e ** -(((x-b[i])**2) / (2*(c[i]**2))))))[0] for i in range(4)] + lines.append(ax.plot(x, sum([a[i]*(e ** -(((x-b[i])**2) / (2*(c[i]**2)))) for i in range(4)]))[0]) + + sliders = [FloatSlider(min=-100, max=100, step=1, value=1) for _ in range(12)] + colors = ['black']*4 + ['red'] + color_dropdowns = [Dropdown(options=['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black'], value=colors[i]) for i in range(5)] + + def update(a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, a3=1, b3=1, c3=1, a4=1, b4=1, c4=1, color1='black', color2='black', color3='black', color4='black', color5='red'): - line1.set_ydata((a1*(e ** -(((x-b1)**2) / (2*(c1**2)))))) - line2.set_ydata((a2*(e ** -(((x-b2)**2) / (2*(c2**2)))))) - line3.set_ydata((a3*(e ** -(((x-b3)**2) / (2*(c3**2)))))) - line4.set_ydata((a4*(e ** -(((x-b4)**2) / (2*(c4**2)))))) - line5.set_ydata(((a1*(e ** -(((x-b1)**2) / (2*(c1**2))))) + - (a2*(e ** -(((x-b2)**2) / (2*(c2**2))))) + - (a3*(e ** -(((x-b3)**2) / (2*(c3**2))))) + - (a4*(e ** -(((x-b4)**2) / (2*(c4**2))))))) - line1.set_color(color1) - line2.set_color(color2) - line3.set_color(color3) - line4.set_color(color4) - line5.set_color(color5) + a = [a1,a2,a3,a4] + b = [b1,b2,b3,b4] + c = [c1,c2,c3,c4] + for i in range(4): + lines[i].set_ydata((a[i]*(e ** -(((x-b[i])**2) / (2*(c[i]**2)))))) + lines[4].set_ydata(sum([a[i]*(e ** -(((x-b[i])**2) / (2*(c[i]**2)))) for i in range(4)])) + for i in range(5): + lines[i].set_color([color1,color2,color3,color4,color5][i]) fig.canvas.draw_idle() - interact(update, a1=a1_slider, b1=b1_slider, c1=c1_slider, a2=a2_slider, b2=b2_slider, c2=c2_slider, - a3=a3_slider, b3=b3_slider, c3=c3_slider, a4=a4_slider, b4=b4_slider, c4=c4_slider, color1=color1_dropdown, - color2=color2_dropdown, color3=color3_dropdown, color4=color4_dropdown, color5=color5_dropdown); + interact(update, + a1=sliders[0], b1=sliders[1], c1=sliders[2], + a2=sliders[3], b2=sliders[4], c2=sliders[5], + a3=sliders[6], b3=sliders[7], c3=sliders[8], + a4=sliders[9], b4=sliders[10], c4=sliders[11], + color1=color_dropdowns[0], color2=color_dropdowns[1], + color3=color_dropdowns[2], color4=color_dropdowns[3], + color5=color_dropdowns[4]) return ax - + def goodness_of_fit(unique_x, combined_gauss, avg_y): """Performs a chi-squared goodness of fit test between the CHART data and the user created combined Gaussian curve. @@ -220,4 +206,7 @@ def goodness_of_fit(unique_x, combined_gauss, avg_y): p_value_x = chi2.sf(chi_squared_statistic_x, len(gauss_masked[:,0]) - 1) p_value_y = chi2.sf(chi_squared_statistic_y, len(y_gauss) - 1) + #print("Chi-squared Statistic for x: ", chi_squared_statistic_x) + #print("P-value for x: ", p_value_x) + #print("Chi-squared Statistic for y: ", chi_squared_statistic_y) print("P-value for y: ", p_value_y) From f8f9a786fd6dd788bdac7b4da5692df57174a055 Mon Sep 17 00:00:00 2001 From: EyanWeissbluth <133696906+EyanWeissbluth@users.noreply.github.com> Date: Mon, 31 Jul 2023 15:20:08 -0700 Subject: [PATCH 3/4] Update analysis.py Please add comments to this. --- src/chart/analysis.py | 75 ++++++++++++++++++++++--------------------- 1 file changed, 38 insertions(+), 37 deletions(-) diff --git a/src/chart/analysis.py b/src/chart/analysis.py index d6f3807..1be87f3 100644 --- a/src/chart/analysis.py +++ b/src/chart/analysis.py @@ -68,36 +68,37 @@ def concat(data_list): raise NotImplementedError() -def LSR_shift(lon, lat, ele, time, RA, Dec): +def LSR_shift(longitude, latitude, elevation, time, altitude, azimuth): """ Identifies the exact postion at which the observations were taken and corrects for the Local Standard of Rest. - Along with this it also converts the celestial coordinates to galactic coordinates. + Along with this it also converts location, altitude, and azimuth to galactic coordinates. - :param lon: longitude - :param lat: latitude - :param ele: elevation - :param time: date and time - :param RA: Right Ascension - :param Dec: Declination + :param latitude: latitude in degrees + :param longitude: longitude in degrees + :param elevation: elevation in meters + :param time: observation time in UTC format string + :param altitude: altitude in degrees + :param azimuth: azimuth in degrees """ - location = EarthLocation.from_geodetic(lon, lat, ele*u.m) #Lon, Lat, elevation + loc = EarthLocation(lat=latitude*u.deg, lon=longitude*u.deg, height=elevation*u.m) + altaz = AltAz(obstime=Time(time), location=loc, alt=altitude*u.deg, az=azimuth*u.deg) + skycoord = SkyCoord(altaz.transform_to(ICRS)) + location = EarthLocation.from_geodetic(longitude, latitude, elevation*u.m) #Lon, Lat, elevation location = location.get_itrs(obstime=Time(time)) #To ITRS frame, makes Earth stationary with Sun - pointing_45deg = SkyCoord(RA,Dec, frame='icrs') #Center of CHART pointing + pointing_45deg = SkyCoord(altaz.transform_to(ICRS)) #Center of CHART pointing frequency = SpectralCoord(1.420405751768e9 * u.Hz, observer=location, target=pointing_45deg) #Shift expected from just local motion f0_shifted = frequency.with_observer_stationary_relative_to('lsrk') #correct for kinematic local standard of rest f0_shifted = f0_shifted.to(u.GHz) v = doppler(f0_shifted,f0) v_adjustment = v.to(u.km/u.second) - print(v_adjustment) - - - coor=SkyCoord(RA, Dec, frame='icrs') - l=coor.galactic.l.degree - b=coor.galactic.b.degree - print(l,b) - - return v_adjustment + return v_adjustment, skycoord.galactic + +def find_array_with_number(freqs, j, number): + for k_index, k in enumerate(freqs[j]): + if numpy.any((k[:-1] <= number) & (number <= k[1:])): + return k_index, k + return None, None def average_overlapping(x1, y1, x2, y2): """ @@ -179,7 +180,6 @@ def update(a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, a3=1, b3=1, c3=1, a4=1, b4=1, c4= return ax - def goodness_of_fit(unique_x, combined_gauss, avg_y): """Performs a chi-squared goodness of fit test between the CHART data and the user created combined Gaussian curve. @@ -187,26 +187,27 @@ def goodness_of_fit(unique_x, combined_gauss, avg_y): param combined_gauss: y values of combined Gaissian curve param avg: y values of overlapping CHART data """ - x_gauss = np.array(unique_x) - y_gauss = np.array(combined_gauss) - x_data = np.array(unique_x) - y_data = np.array(avg_y) + x_observed = np.array(unique_x) + y_observed = np.array(combined_gauss) + x_expected = np.array(unique_x) * 2 + y_expected = np.array(avg_y) * 2 - gauss = np.concatenate((x_gauss.reshape(-1,1), y_gauss.reshape(-1,1)), axis=1) - data = np.concatenate((x_data.reshape(-1,1), y_data.reshape(-1,1)), axis=1) + observed = np.concatenate((x_observed.reshape(-1,1), y_observed.reshape(-1,1)), axis=1) + expected = np.concatenate((x_expected.reshape(-1,1), y_expected.reshape(-1,1)), axis=1) + mask = (observed[:,0] >= -100) & (observed[:,0] <= 100) + observed_masked = observed[mask] + expected_masked = expected[mask] - mask = (gauss[:,0] >= -100) & (gauss[:,0] <= 100) - gauss_masked = gauss[mask] - data_masked = data[mask] + chi_squared_statistic_x = np.sum((observed_masked[:,0] - expected_masked[:,0])**2 / expected_masked[:,0]) + chi_squared_statistic_y = np.sum((y_observed - y_expected)**2 / y_expected) - chi_squared_statistic_x = np.sum((gauss_masked[:,0] - data_masked[:,0])**2 / data_masked[:,0]) - chi_squared_statistic_y = np.sum((y_gauss - y_data)**2 / y_data) + p_value_x = chi2.sf(chi_squared_statistic_x, len(observed_masked[:,0]) - 1) + p_value_y = chi2.sf(chi_squared_statistic_y, len(y_observed) - 1) - p_value_x = chi2.sf(chi_squared_statistic_x, len(gauss_masked[:,0]) - 1) - p_value_y = chi2.sf(chi_squared_statistic_y, len(y_gauss) - 1) + chi_squared_statistic_xy = chi_squared_statistic_x + chi_squared_statistic_y + degrees_of_freedom_xy = len(observed) - 2 + reduced_chi_squared_statistic_xy = chi_squared_statistic_xy / degrees_of_freedom_xy + p_value_xy = chi2.sf(reduced_chi_squared_statistic_xy, degrees_of_freedom_xy) + return reduced_chi_squared_statistic_xy, p_value_xy - #print("Chi-squared Statistic for x: ", chi_squared_statistic_x) - #print("P-value for x: ", p_value_x) - #print("Chi-squared Statistic for y: ", chi_squared_statistic_y) - print("P-value for y: ", p_value_y) From de9233e2e1479ce84d1bf034745b61e8525d2edf Mon Sep 17 00:00:00 2001 From: EyanWeissbluth <133696906+EyanWeissbluth@users.noreply.github.com> Date: Mon, 31 Jul 2023 15:24:41 -0700 Subject: [PATCH 4/4] Update analysis.py Please add comments to this one. --- src/chart/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/chart/analysis.py b/src/chart/analysis.py index 1be87f3..cb77afd 100644 --- a/src/chart/analysis.py +++ b/src/chart/analysis.py @@ -180,6 +180,7 @@ def update(a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, a3=1, b3=1, c3=1, a4=1, b4=1, c4= return ax + def goodness_of_fit(unique_x, combined_gauss, avg_y): """Performs a chi-squared goodness of fit test between the CHART data and the user created combined Gaussian curve. @@ -210,4 +211,3 @@ def goodness_of_fit(unique_x, combined_gauss, avg_y): reduced_chi_squared_statistic_xy = chi_squared_statistic_xy / degrees_of_freedom_xy p_value_xy = chi2.sf(reduced_chi_squared_statistic_xy, degrees_of_freedom_xy) return reduced_chi_squared_statistic_xy, p_value_xy -