From f7d3e0692e26504c7dfd783b3f8aeebb321488ea Mon Sep 17 00:00:00 2001 From: jurjen93 Date: Thu, 18 Jul 2024 18:12:35 +0200 Subject: [PATCH] updates --- catalogue_helpers/associate_components.py | 53 ---- .../random_forest_outlier_detection.py | 236 -------------- ms_helpers/concat_with_dummies.py | 2 +- ms_helpers/ms_flagger.py | 17 +- ms_merger.py | 10 +- other/calc_uvw.py | 296 ++++++++++++++++++ other/precession.py | 125 ++++++++ 7 files changed, 440 insertions(+), 299 deletions(-) delete mode 100644 catalogue_helpers/associate_components.py delete mode 100644 catalogue_helpers/random_forest_outlier_detection.py create mode 100644 other/calc_uvw.py create mode 100644 other/precession.py diff --git a/catalogue_helpers/associate_components.py b/catalogue_helpers/associate_components.py deleted file mode 100644 index b758f14c..00000000 --- a/catalogue_helpers/associate_components.py +++ /dev/null @@ -1,53 +0,0 @@ -from astropy.table import Table -import numpy as np - -def error_prop(errors): - return np.sqrt(np.sum(np.power(errors, 2))) - - -def get_table_index(t, source_id): - return int(np.argwhere(t['Source_id'] == source_id).squeeze()) - - -def associate(associate_components, table): - """ - associate components by merging components - - :param associate_components: Example: [{4: [1, 2, 3, 4, 5]}, {10: [9, 10, 11]}]. This will merge component 1,2,3,4,5 with 4 as middle component - :param fits table - """ - - t = Table.read(table, format='fits') - - to_delete = [] - to_not_delete = [] - - for p in associate_components: - if type(p) == dict: - main_ID = list(p.keys())[0] - main_ID = get_table_index(t, main_ID) - to_not_delete.append(main_ID) - ids = [get_table_index(t, x) for x in list(set(p[main_ID]))] - t[main_ID]['Total_flux'] = t[ids]['Total_flux'].sum() - t[main_ID]['E_Total_flux'] = error_prop(t[ids]['E_Total_flux']) - t[main_ID]['Peak_flux'] = t[ids]['Peak_flux'].max() - t[main_ID]['E_Peak_flux'] = error_prop(t[ids]['E_Peak_flux']) - t[main_ID]['S_Code'] = 'M' - t[main_ID]['Maj'] = t[ids]['Maj'].min() - t[main_ID]['Min'] = t[ids]['Min'].min() - t[main_ID]['E_Maj'] = error_prop(t[ids]['E_Maj']) - t[main_ID]['E_Min'] = error_prop(t[ids]['E_Min']) - t[main_ID]['PA'] = t[ids]['PA'].min() - t[main_ID]['E_PA'] = error_prop(t[ids]['E_PA']) - for i in ids: - to_delete.append(i) - if type(p) == int: - to_delete.append(get_table_index(t, p)) - - to_delete = list(set(to_delete)) - to_not_delete = list(set(to_not_delete)) - for i in sorted(to_delete)[::-1]: - if i not in to_not_delete: - del t[i] - - t.write(table.replace('.fits', '_final.fits'), format='fits', overwrite=True) diff --git a/catalogue_helpers/random_forest_outlier_detection.py b/catalogue_helpers/random_forest_outlier_detection.py deleted file mode 100644 index bc9c61bb..00000000 --- a/catalogue_helpers/random_forest_outlier_detection.py +++ /dev/null @@ -1,236 +0,0 @@ -""" -This is an experimental script to use Random Forest classifiers to separate spurious detections from real detections. - -Has not been used for anything yet. -Use it for your own needs. Be aware of hardcoded paths. -""" - -from astropy.table import Table -from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier -from sklearn.preprocessing import LabelEncoder -from sklearn.model_selection import train_test_split -from sklearn.metrics import mean_squared_error, r2_score, accuracy_score -from sklearn.model_selection import RandomizedSearchCV as RSCV -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -from sklearn import tree -import os -from glob import glob - - -# Table columns -cols = ['Source_id', 'E_Total_flux', 'E_Peak_flux', - 'S_Code', 'E_pos', 'Total_flux/peak_flux', - 'Peak_flux/RMS', 'Total_flux/RMS', 'E_MajMin', - 'MajMin', 'PA', 'E_PA'] - -sep = '\n##########################################################################\n' - - -def gridsearch_test(df): - """ - - Grid search test - - :param df: pandas data frame - :return: best model - """ - - df = df.dropna() - - param_grid = {'n_estimators': np.arange(50, 200, 15), - 'max_features': np.arange(0.1, 1, 0.1), - 'max_depth': [3, 4, 5], - 'max_samples': [0.3, 0.4, 0.5, 0.6, 0.7, 0.8]} - - X_train, X_test, y_train, y_test = train_test_split(df[[c for c in cols if c not in ['Source_id', 'label']]], - df.label, test_size=0.3, random_state=42) - model = RSCV(RandomForestClassifier(), param_grid, n_iter=15).fit(X_train, - y_train.astype('int')) - model = model.best_estimator_ - print(sep + 'BEST MODEL\n') - print(model) - print(sep) - predictions = model.predict(X_test) - mse = mean_squared_error(y_test, predictions) - print('MSE: ' + str(mse)) - r2 = r2_score(y_test, predictions.round(0)) - print('R2: ' + str(r2)) - indices = np.argwhere((y_test - predictions.round(0)).to_numpy() != 0) - print('Differently classified sources:' + str(y_test.index[indices].squeeze())) - print('Accuracy: '+str(accuracy_score(y_test.astype(int), np.array(predictions).astype(int)))) - print(sep) - return model - - -def gridsearch(df): - """ - - Grid search - - :param df: pandas data frame - :return: best model - """ - df = df.dropna() - - param_grid = {'n_estimators': np.arange(30, 200, 5), - 'max_features': np.arange(0.1, 1, 0.1), - 'max_depth': [3, 4, 5], - 'max_samples': [0.3, 0.4, 0.5, 0.6, 0.7, 0.8]} - - model = RSCV(RandomForestClassifier(), param_grid, n_iter=30). \ - fit(df[[c for c in cols if c not in ['Source_id', 'label']]], df.label.astype('int')) - model = model.best_estimator_ - print(sep + 'BEST MODEL\n') - print(model) - print(sep) - return model - - -def sort_feature_importance(feat_importance, labels): - """ - Sort feature importance from model - - :param feat_importance: unsorted feature importance model - :param labels: unsorted feature importance labels - - :return: sorted feature importance - """ - important_features_dict = {} - for idx, val in enumerate(feat_importance): - important_features_dict[idx] = val - - important_features_list = sorted(important_features_dict, - key=important_features_dict.get, - reverse=True) - return feat_importance[important_features_list], np.array(labels)[important_features_list] - - -def get_table_data(tbl): - """ - Get table data for training/testing/predicting - - :param tbl: table name - - :return: pandas data frame - - """ - T = Table.read(tbl) - df = T.to_pandas() - - le = LabelEncoder() - le.fit([b'S', b'C', b'M']) - df['S_Code'] = le.transform(df.S_Code) - df['Total_flux/peak_flux'] = df['Total_flux'] / df['Peak_flux'] - df['Peak_flux/RMS'] = df['Peak_flux'] / df['Isl_rms'] - df['Total_flux/RMS'] = df['Total_flux'] / df['Isl_rms'] - df['E_pos'] = np.sqrt(df['E_RA']**2+df['E_DEC']**2) - df["E_MajMin"] = np.sqrt(df["E_Maj"]**2 + df["E_Min"]**2) - df["MajMin"] = df["Maj"] * df["Min"] - - df = df[cols] - - return df - - -def plot_feature_importance(feat_importance, labels): - """ - Plot feature importance - - :param feat_importance: feature importance - :param labels: labels - """ - plt.style.use("bmh") - - hfont = {'fontname': 'monospace'} - forest_importances = pd.Series(feat_importance, - index=[l for l in labels]) - fig, ax = plt.subplots() - forest_importances.plot.bar(ax=ax, color='darkred') - # ax.set_title() - ax.set_ylabel("Feature importance") - plt.xticks(rotation=75, **hfont) - fig.tight_layout() - plt.grid(False) - for pos in ['right', 'top', 'bottom', 'left']: - plt.gca().spines[pos].set_visible(False) - plt.savefig('feat_importance.png', dpi=300) - # plt.show() - - -def predict_falsepositives(tbl, model): - """ - Predict which sources are false positives with model - - :param tbl: table name - :param model: RF model - - :return sources to delete - """ - df1 = get_table_data(tbl).set_index("Source_id")[[c for c in cols if c not in ['Source_id', 'label']]] - predicted = model.predict(df1) - df1['prediction'] = predicted - return list(df1[df1.prediction == 0].index), df1.columns - - -def plot_tree(model, tree_idx, columns): - """ - Plot a tree - """ - plt.figure(figsize=(200, 200)) - _ = tree.plot_tree(model.estimators_[tree_idx], feature_names=columns, filled=True) - plt.show() - -def get_table_index(t, source_id): - return int(np.argwhere(t['Source_id'] == source_id).squeeze()) - -def clean_table(tbl, model): - """ - Make clean table - - :param tbl: fits table name - """ - to_remove, columns = predict_falsepositives(tbl, model) - T = Table.read(tbl, format='fits') - print(tbl) - for i in sorted(to_remove)[::-1]: - n = get_table_index(T, i) - if T[n]['Peak_flux'] < 5.5*T[n]['Isl_rms']: - print('Deleted source ID: ' + str(i)) - del T[n] - - T.write(tbl.replace('.fits', '_clean.fits'), format='fits', overwrite=True) - - -def main(): - """ - Main function - """ - - train_table = "finalcat03/facet_19_source_catalog_final.fits" - df = get_table_data(train_table) - - # false positive indices from test data by eye - falsepositives = [246, 248, 249, 252, 253, 254, 258, 330, 263, 264, 282, 288, 289, 290, 291, 292, 293, 294, - 295, 296, 297, 298, 299, 300, 301, 303, 305, 308, 347, 310, 313, 315, 316, 318, 319, 320, - 322, 323, 327, 330, 331, 333, 334, 335, 341, 255, 259, 262, 266, 267, 268, 279, 284, 367, - 362, 358, 357, 356, 355, 351, 349, 344, 343, 342] - - falsepositives = sorted(list(set(falsepositives))) - - df['label'] = None - df.loc[df.Source_id.isin(falsepositives), 'label'] = False - df.loc[~df.Source_id.isin(falsepositives), 'label'] = True - - gridsearch_test(df) - model = gridsearch(df) - importances, labels = sort_feature_importance(model.feature_importances_, - [c for c in cols if c not in ['Source_id', 'label']]) - plot_feature_importance(importances, labels) - - for cat in glob("finalcat03/facet_*_source_catalog.fits"): - clean_table(cat, model) - -if __name__ == '__main__': - main() diff --git a/ms_helpers/concat_with_dummies.py b/ms_helpers/concat_with_dummies.py index d4583bab..41132b8c 100644 --- a/ms_helpers/concat_with_dummies.py +++ b/ms_helpers/concat_with_dummies.py @@ -2,7 +2,7 @@ Concat measurement sets with dummies Example: - python concat_with_dummies.py --ms *.ms --concat_name concat.ms --data_column CORRECTED_DATA + python concat_with_dummies.py --concat_name concat.ms --freq_avg 4 --time_avg 8 *.ms """ __author__ = "Jurjen de Jong" diff --git a/ms_helpers/ms_flagger.py b/ms_helpers/ms_flagger.py index 61d76c4c..921925a8 100644 --- a/ms_helpers/ms_flagger.py +++ b/ms_helpers/ms_flagger.py @@ -28,20 +28,31 @@ def main(): if args.timerange is not None: min_time, max_time = [f for f in args.timerange.split('-')] + for ms in args.ms: + + expr = '' + command = ['DP3', 'msout.storagemanager=dysco', f'msin={ms}', f'msout=flagged_{ms.split("/")[-1]}', 'steps=[flag]', 'flag.type=preflagger', - f'flag.baseline={args.ant}'] + f'flag.flag1.baseline={args.ant}'] + if args.freqrange is not None or args.timerange is not None: + expr += 'flag1' # frequency flagging if args.freqrange is not None: - command += [f"flag.freqrange='[{min_freq} .. {max_freq} MHz]'"] + command += [f"flag.flag2.freqrange='[{min_freq} .. {max_freq} MHz]'"] + expr += ' and flag2' # time flagging if args.timerange is not None: - command += [f"flag.reltime='[{min_time} .. {max_time}]'"] + command += [f"flag.flag3.reltime='[{min_time} .. {max_time}]'"] + expr += ' and flag3' + + command += [f'flag.expr="{expr}"'] + print(' '.join(command)) os.system(' '.join(command)) diff --git a/ms_merger.py b/ms_merger.py index ea91d622..db08f635 100644 --- a/ms_merger.py +++ b/ms_merger.py @@ -176,16 +176,14 @@ def time_resolution(resolution, fov_diam, time_smearing=0.95): :return: integration time in seconds """ - #TODO: Make dependent on frequency - # Convert distance from degrees to radians distance_from_phase_center_rad = np.deg2rad(fov_diam/2) # Calculate angular resolution (radians) angular_resolution_rad = resolution*4.8481*1e-6 - int_time = (np.sqrt((1-time_smearing)/1.2288710615597145e-09)/ - distance_from_phase_center_rad*angular_resolution_rad) + int_time = 2.9*10**4*(angular_resolution_rad*np.sqrt(1-time_smearing)/ + distance_from_phase_center_rad) return int_time @@ -321,7 +319,7 @@ def mjd_seconds_to_lst_seconds(mjd_seconds, longitude_deg=52.909): location = EarthLocation(lon=longitude_deg * u.deg, lat=0. * u.deg) # Calculate LST in hours - lst_hours = time_utc.sidereal_time('mean', longitude=location.lon).hour + lst_hours = time_utc.sidereal_time('apparent', longitude=location.lon).hour # Convert LST from hours to seconds lst_seconds = lst_hours * 3600.0 @@ -1269,7 +1267,7 @@ def __init__(self, msin: list = None, outname: str = 'empty.ms', chunkmem: float self.num_cpus = psutil.cpu_count(logical=True) total_memory = psutil.virtual_memory().total / (1024 ** 3) # in GB target_chunk_size = total_memory / chunkmem - self.chunk_size = min(int(target_chunk_size * (1024 ** 3) / np.dtype(np.float128).itemsize/2/self.freq_len), 2_000_000) + self.chunk_size = min(int(target_chunk_size * (1024 ** 3) / np.dtype(np.float128).itemsize/2/self.freq_len), 500_000) print(f"\n---------------\nChunk size ==> {self.chunk_size}") diff --git a/other/calc_uvw.py b/other/calc_uvw.py new file mode 100644 index 00000000..543a61e5 --- /dev/null +++ b/other/calc_uvw.py @@ -0,0 +1,296 @@ +from __future__ import print_function, division +import numpy as np +from astropy.time import Time, TimeDelta +import matplotlib.pyplot as plt +import astropy.units as u +from astropy.coordinates import (SkyCoord, FK5, ITRS, CIRS, EarthLocation, CartesianRepresentation, + solar_system_ephemeris, get_body_barycentric_posvel) +from PyAstronomy import pyasl +import datetime + +def calculate_uvw(X1, Y1, Z1, X2, Y2, Z2, time_str, alpha, delta): + """ + Calculate UVW coordinates for a given observation. + + Parameters: + X1, Y1, Z1 : float : Geocentric coordinates of antenna 1 + X2, Y2, Z2 : float : Geocentric coordinates of antenna 2 + time_str : str : Observation time in UTC (ISO format string) + alpha : float : Right Ascension of the source in degrees + delta : float : Declination of the source in degrees + latitude : float : Latitude of the observation location in degrees + longitude : float : Longitude of the observation location in degrees + + Returns: + uvw : numpy.ndarray : UVW coordinates + """ + + def midpoint_to_lat_lon(X1, Y1, Z1, X2, Y2, Z2): + # Calculate midpoint coordinates + X_m = (X1 + X2) / 2 + Y_m = (Y1 + Y2) / 2 + Z_m = (Z1 + Z2) / 2 + + # Create a Cartesian representation of the midpoint + mid_cartesian = CartesianRepresentation(X_m * u.m, Y_m * u.m, Z_m * u.m) + + # Convert to EarthLocation to get latitude and longitude + mid_location = EarthLocation.from_geocentric(mid_cartesian.x, mid_cartesian.y, mid_cartesian.z) + + return mid_location.lat.deg, mid_location.lon.deg + + # Baseline vector + B = np.array([X2 - X1, Y2 - Y1, Z2 - Z1]) + + # Define observer location + latitude, longitude = midpoint_to_lat_lon(X1, Y1, Z1, X2, Y2, Z2) + location = EarthLocation(lat=latitude, lon=longitude) + + # Convert to LST + observation_time = Time(time_str, location=location) + lst = observation_time.sidereal_time('mean', longitude=longitude).deg + + # Calculate the hour angle + H = lst - alpha + if H < 0: + H += 360.0 + + # Convert degrees to radians + H_rad = np.radians(H) + delta_rad = np.radians(delta) + + # UVW Transformation matrix + transformation_matrix = np.array([ + [np.sin(H_rad), np.cos(H_rad), 0], + [-np.sin(delta_rad) * np.cos(H_rad), np.sin(delta_rad) * np.sin(H_rad), np.cos(delta_rad)], + [np.cos(delta_rad) * np.cos(H_rad), -np.cos(delta_rad) * np.sin(H_rad), np.sin(delta_rad)] + ]) + + # Compute UVW coordinates + uvw = np.dot(transformation_matrix, B) + return uvw + + +# def aberration_shift(timestamp, ra, dec): +# +# """ FROM: https://pyastronomy.readthedocs.io/en/latest/pyaslDoc/aslDoc/aberration.html """ +# +# # Convert calendar date to JD using the datetime package +# dt = datetime.datetime.strptime(timestamp, "%Y-%m-%dT%H:%M:%S") +# jd = pyasl.jdcnv(dt) +# +# # Specify RA and DEC +# ra = np.atleast_1d(ra) +# dec = np.atleast_1d(dec) +# +# # Get change in RA and DEC due to annual aberration +# res = pyasl.co_aberration(np.repeat(jd, ra.size), ra, dec) +# +# return res[0][0]*3600, res[1][0]*3600 + + +def aberration_correction(timestamp, ra, dec): + """ + Calculate the aberration corrected RA and DEC for given input RA, DEC and timestamp. + + Parameters: + ra (float): Right Ascension in degrees + dec (float): Declination in degrees + timestamp (str): Timestamp in ISO format (e.g., '2024-07-17T00:00:00') + + Returns: + corrected_ra (float): Aberration corrected Right Ascension in degrees + corrected_dec (float): Aberration corrected Declination in degrees + """ + # Convert timestamp to Julian Date + time = Time(timestamp, format='isot', scale='utc') + + # Initial coordinates + sky_coord = SkyCoord(ra=ra, dec=dec, unit='deg') + + # Calculate the position and velocity of Earth + with solar_system_ephemeris.set('builtin'): + pos, vel = get_body_barycentric_posvel('earth', time) + + # Speed of light in AU/day + c = 173.144632674240 + + # Aberration correction + ra_rad = np.deg2rad(ra) + dec_rad = np.deg2rad(dec) + + # Earth's velocity components + vx = vel.x.value + vy = vel.y.value + vz = vel.z.value + + # Proper motion in RA and DEC + delta_ra = (-vx * np.sin(ra_rad) + vy * np.cos(ra_rad)) / c + delta_dec = (-vx * np.cos(ra_rad) * np.sin(dec_rad) + - vy * np.sin(ra_rad) * np.sin(dec_rad) + + vz * np.cos(dec_rad)) / c + + + return np.rad2deg(delta_ra) * 3600, np.rad2deg(delta_dec) * 3600 + + +def precession_nutation_shift(timestamp, ra, dec): + """ + Calculate the precession and nutation corrected RA and DEC for given input RA, DEC, and timestamp. + + Parameters: + ra (float): Right Ascension in degrees + dec (float): Declination in degrees + timestamp (str): Timestamp in ISO format (e.g., '2024-07-17T00:00:00') + + Returns: + corrected_ra (float): Precession and nutation corrected Right Ascension in degrees + corrected_dec (float): Precession and nutation corrected Declination in degrees + """ + # Convert timestamp to astropy Time object + t = Time(timestamp) + + # Create a SkyCoord object for the input RA and DEC + coord = SkyCoord(ra=ra * u.degree, dec=dec * u.degree, frame='icrs') + + # Convert the coordinates to the FK5 frame at the given epoch to account for precession + coord_fk5 = coord.transform_to(FK5(equinox=t)) + + # Convert the coordinates to the CIRS frame to account for nutation + coord_cirs = coord_fk5.transform_to(CIRS(obstime=t)) + + # Extract the corrected RA and DEC + corrected_ra = coord_cirs.ra.degree + corrected_dec = coord_cirs.dec.degree + + return (ra-corrected_ra)*3600, (dec-corrected_dec)*3600 + + +def calculate_uvw_shift(ra_deg, dec_deg, baseline_length, ra_shift_arcsec, dec_shift_arcsec): + """ + Calculate the shift in UVW coordinates based on shifts in RA and DEC. + + Parameters: + ra_deg (float): Original Right Ascension in degrees + dec_deg (float): Original Declination in degrees + baseline_length (float): Baseline length in meters + ra_shift_arcsec (float): Shift in Right Ascension in arcseconds + dec_shift_arcsec (float): Shift in Declination in arcseconds + + Returns: + delta_u (float): Shift in U coordinate in meters + delta_v (float): Shift in V coordinate in meters + delta_w (float): Shift in W coordinate in meters + """ + + def ra_dec_to_radians(ra_shift_arcsec, dec_shift_arcsec): + """Convert RA and DEC shifts from arcseconds to radians.""" + ra_shift_rad = np.deg2rad(ra_shift_arcsec / 3600.0) + dec_shift_rad = np.deg2rad(dec_shift_arcsec / 3600.0) + return ra_shift_rad, dec_shift_rad + + def direction_cosines(ra_rad, dec_rad): + """Calculate direction cosines from RA and DEC in radians.""" + l = np.cos(dec_rad) * np.cos(ra_rad) + m = np.cos(dec_rad) * np.sin(ra_rad) + n = np.sin(dec_rad) + return l, m, n + + def uvw_shift(ra_rad, dec_rad, ra_shift_rad, dec_shift_rad, baseline_length): + """Calculate the shift in UVW coordinates based on RA and DEC shifts.""" + # Original direction cosines + l0, m0, n0 = direction_cosines(ra_rad, dec_rad) + + # Shifted RA and DEC + ra_shifted_rad = ra_rad + ra_shift_rad + dec_shifted_rad = dec_rad + dec_shift_rad + + # Shifted direction cosines + l1, m1, n1 = direction_cosines(ra_shifted_rad, dec_shifted_rad) + + # Difference in direction cosines + delta_l = l1 - l0 + delta_m = m1 - m0 + delta_n = n1 - n0 + + # Convert direction cosine shifts to UVW coordinate shifts + delta_u = baseline_length * delta_l + delta_v = baseline_length * delta_m + delta_w = baseline_length * delta_n + + return delta_u, delta_v, delta_w + + # Convert shifts to radians + ra_shift_rad, dec_shift_rad = ra_dec_to_radians(ra_shift_arcsec, dec_shift_arcsec) + + # Convert RA and DEC to radians + ra_rad = np.deg2rad(ra_deg) + dec_rad = np.deg2rad(dec_deg) + + # Calculate UVW coordinate shifts + delta_u, delta_v, delta_w = uvw_shift(ra_rad, dec_rad, ra_shift_rad, dec_shift_rad, baseline_length) + + return delta_u, delta_v, delta_w + + +# Example usage: +X1, Y1, Z1 = 3826896.235000, 460979.455000, 5064658.203000 +X2, Y2, Z2 = 3370271.657000, 712125.881000, 5349991.165000 +ra = 241.88 # Right Ascension in degrees +dec = 54.33 # Declination in degrees + +# Length of a sidereal day +sidereal_day = TimeDelta(23.9344696 * 3600, format='sec') # 23 hours, 56 minutes, 4 seconds + + +for n, timez in enumerate(["2018-11-26T07:13:43", "2020-05-24T19:20:26", "2020-11-14T08:11:00", "2021-05-13T19:41:00"]): + + dra_aber, ddec_aber = aberration_correction(timez, ra, dec) + dra_prec, ddec_prec = precession_nutation_shift(timez, ra, dec) + # if n>=1: + # ddec_aber-=40 + du, dv, dw = calculate_uvw_shift(ra, dec, 600_000, dra_prec + dra_aber, ddec_prec + ddec_aber) + print(timez) + print(dra_prec, ddec_prec) + print(dra_aber, ddec_aber) + + + if n==0: + ref_dra_prec, ref_ddec_prec = dra_prec, ddec_prec + ref_dra_aber, ref_ddec_aber = dra_aber, ddec_aber + ref_du, ref_dv, ref_dw = du, dv, dw + + # else: + # print(dra_prec-ref_dra_prec, ddec_prec-ref_ddec_prec) + # print(dra_aber-ref_dra_aber, ddec_aber-ref_ddec_aber) + # print((dra_prec-ref_dra_prec)-(dra_aber-ref_dra_aber), (ddec_prec-ref_ddec_prec)-(ddec_aber-ref_ddec_aber)) + + + + print(du-ref_du, dv-ref_dv) + if n==0: + plt.scatter([(du-ref_du)/1000+518.45-i*0.03 for i in range(10)], + [(dv-ref_dv)/1000+212.55 + 0.3*i for i in range(10)], label=timez, marker="*") + else: + plt.scatter([(du-ref_du)/1000-0.033+518.45-i*0.03 for i in range(10)], + [(dv-ref_dv)/1000-0.033+212.55 + 0.3*i for i in range(10)], label=timez, marker="*") + print() + + # # pos = np.array( + # # [calculate_uvw(X1, Y1, Z1, X2, Y2, Z2, + # # (Time(timez) + i*4/(23.9344696 * 3600)).iso, + # # ra-(dra_prec+dra_aber)/3600, + # # dec-(ddec_prec+ddec_aber)/3600) + # # for i in range(int(3600))]) + # pos = np.array( + # [calculate_uvw(X1, Y1, Z1, X2, Y2, Z2, + # (Time(timez) + i*8/(23.9344696 * 3600)).iso, + # ra+(dra_prec+dra_aber)/3600, + # dec+(ddec_prec+ddec_aber)/3600) + # for i in range(int(3600))]) + # plt.scatter(pos[:, 0]/1000, pos[:, 1]/1000, label=timez, s=30, marker="*") + +plt.xlim(518.1, 518.6) +plt.ylim(212.4, 215.1) +plt.legend() +plt.show() diff --git a/other/precession.py b/other/precession.py new file mode 100644 index 00000000..5e45b3ed --- /dev/null +++ b/other/precession.py @@ -0,0 +1,125 @@ +import numpy as np +from astropy.time import Time +import astropy.units as u +from astropy.coordinates import SkyCoord, FK5 + + +def compute_julian_century(epoch): + """ + Compute the Julian century for a given epoch relative to J2000.0. + """ + JD = epoch.jd + JD2000 = 2451545.0 # Julian date for J2000.0 + T = (JD - JD2000) / 36525.0 + return T + + +def compute_precession_angles(epoch1, epoch2): + """ + Compute the precession angles to transform coordinates from epoch1 to epoch2. + Based on the IAU 2000 precession model. + """ + T1 = compute_julian_century(epoch1) + T2 = compute_julian_century(epoch2) + dT = T2 - T1 + + zeta_A = (2306.2181 + 1.39656 * T1 - 0.000139 * T1 ** 2) * dT \ + + (0.30188 - 0.000344 * T1) * dT ** 2 + 0.017998 * dT ** 3 + theta_A = (2004.3109 - 0.85330 * T1 - 0.000217 * T1 ** 2) * dT \ + - (0.42665 + 0.000217 * T1) * dT ** 2 - 0.041833 * dT ** 3 + z_A = (2306.2181 + 1.39656 * T1 - 0.000139 * T1 ** 2) * dT \ + + (1.09468 + 0.000066 * T1) * dT ** 2 + 0.018203 * dT ** 3 + + # Convert to radians + zeta_A = np.deg2rad(zeta_A / 3600) + theta_A = np.deg2rad(theta_A / 3600) + z_A = np.deg2rad(z_A / 3600) + + return zeta_A, theta_A, z_A + + +def rotation_matrix(angle, axis): + """ + Create a rotation matrix for a given angle and axis. + """ + if axis == 'x': + return np.array([[1, 0, 0], + [0, np.cos(angle), -np.sin(angle)], + [0, np.sin(angle), np.cos(angle)]]) + elif axis == 'y': + return np.array([[np.cos(angle), 0, np.sin(angle)], + [0, 1, 0], + [-np.sin(angle), 0, np.cos(angle)]]) + elif axis == 'z': + return np.array([[np.cos(angle), -np.sin(angle), 0], + [np.sin(angle), np.cos(angle), 0], + [0, 0, 1]]) + + +def compute_precession_matrix(epoch1, epoch2): + """ + Compute the precession matrix to transform coordinates from epoch1 to epoch2. + """ + zeta_A, theta_A, z_A = compute_precession_angles(epoch1, epoch2) + + # Precession matrix + P = np.dot(rotation_matrix(-z_A, 'z'), np.dot(rotation_matrix(theta_A, 'y'), rotation_matrix(-zeta_A, 'z'))) + + return P + + +def precess_uvw(uu, v, w, ra, dec, epoch1, epoch2): + """ + Precess the UVW coordinates from epoch1 to epoch2. + """ + uvw = np.array([uu, v, w]) + + # Compute precession matrix + precession_matrix = compute_precession_matrix(epoch1, epoch2) + + # Apply precession matrix to UVW coordinates + uvw_precessed = np.dot(precession_matrix, uvw) + + # Create SkyCoord object for the initial position with units + skycoord_initial = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, frame=FK5(equinox=epoch1)) + + # Precess to the target epoch + skycoord_target = skycoord_initial.transform_to(FK5(equinox=epoch2)) + + # Calculate the phase rotation due to the coordinate change + ra_new = skycoord_target.ra.deg + dec_new = skycoord_target.dec.deg + + # Calculate the phase shift vector + ra_diff = np.deg2rad(ra_new - ra) + dec_diff = np.deg2rad(dec_new - dec) + + # Correct the UVW coordinates for the phase shift + # Apply phase shift in RA + uvw_corrected_ra = np.dot(rotation_matrix(ra_diff, 'z'), uvw_precessed) + # Apply phase shift in Dec + uvw_corrected = np.dot(rotation_matrix(dec_diff, 'x'), uvw_corrected_ra) + + return uvw_corrected, ra_new, dec_new + + +# Initial UVW coordinates (example values in meters) +u_initial, v_initial, w_initial = 516117.47778081, 99503.4107472, 276978.779884 + +# Initial epoch and target epoch +initial_epoch = Time('2021-05-13', scale='utc') +target_epoch = Time('2022-05-13', scale='utc') + +# Right Ascension and Declination of the observed point (example values in degrees) +ra_initial = 242.75 # Example RA in degrees +dec_initial = 55.0 # Example Dec in degrees + +# Compute precessed UVW coordinates and new RA/Dec +uvw_new, ra_new, dec_new = precess_uvw(u_initial, v_initial, w_initial, ra_initial, dec_initial, initial_epoch, + target_epoch) + +print(f"New UVW coordinates: u={uvw_new[0]:.3f} m, v={uvw_new[1]:.3f} m, w={uvw_new[2]:.3f} m") +print(f"Delta UVW coordinates: delta u={u_initial-uvw_new[0]:.3f} m, delta v={v_initial-uvw_new[1]:.3f} m, delta w={w_initial-uvw_new[2]:.3f} m") + +print(f"New RA: {ra_new:.6f} degrees, New Dec: {dec_new:.6f} degrees") +print(f"Precession: {np.sqrt((u_initial - uvw_new[0]) ** 2 + (v_initial - uvw_new[1]) ** 2 + (w_initial - uvw_new[2]) ** 2)} meters")