import java.io.*; class MinorPlanet extends SolarSystemBody { /* DECLARE INSTANCE VARIABLES */ /* Auxiliary Variables: */ /* A list (mp_list) is maintained of the minor planets for whom a data file (observations and orbit) exists. */ public int number_of_minor_planets = 0; public String[] mp_list = new String[1001]; /* The current minor planet on that list */ public String current_minor_planet = " "; /* In method "probability_of_collision", it is necessary to determine the state transition matrix, which essentially represents the sensitivity of the current state vector to small changes in the epoch state vector. This requires integrating forward from epoch along eight different variant paths; rather than start from scratch for each collision/near miss event, these instance variables will be used to keep a "running count" of the variant trajectories. */ public double[][] collision_sensitivity_r_variant = new double[9][4]; public double[][] collision_sensitivity_rprime_variant = new double[9][4]; public double collision_sensitivity_time = 0; public double[][][] collision_state_matrix = new double[2000][9][9]; public double[][] collision_nominal_state = new double[2000][9]; /* nominal_observable is an array containing the nominal computed value of each ra, dec, delay, or doppler observation, as determined by method get_residuals. Note that the observables treat a ra and a dec as separate, so the index will not correspond to the observation number. */ public double[] nominal_observable = new double[40001]; /* In response to feedback from Dr. Gerhard Hahn of DLR, I am adding a GUI-selectable non-database variable representing the threshold approach distance defining a near-miss. The default value of 0.05 AU matches the MPC definition. */ public double near_miss_threshold = 0.05; /* In searching for virtual impactors, users have encountered collisions on variant trajectories; to capitalize on these fortuitous accidents, I'm adding three flags that will control whether variant trajectories are checked for collisions. "VIflag1" will indicate that a search for VIs is underway. "VIflag2" is set by detect_collision to prevent the variant trajectories from being tested for collision; the absence of this flag is interpreted to mean that variant trajectories should be tested for collision. "VIflag3" is set for variant trajectories by probability_of_collision to prevent a recursive call to probability_of_collision from within detect_collision. "VIflag4" allows detect_collision while suppressing calls to probability_of_collision, thus avoiding the use of the six variant trajectories. "VIflag5" is used within detect_collision to prevent bisection calls to update from resetting the VI_trajectory_r and _rprime vectors. */ public boolean VIflag1 = false, VIflag2 = false, VIflag3 = false, VIflag4 = false, VIflag5 = false; /* In searching for virtual impactors, we often want the variant trajectories to be centered on a non-nominal trajectory. These arrays will be used to store the state vector passed to update by the VI search algorithm; if the VIflag1 flag is set, probability_of_collision will center its variant trajectories around this state vector, rather than the nominal epoch state vector. */ public double[] VI_trajectory_r = new double[4]; public double[] VI_trajectory_rprime = new double[4]; /* In method "least_squares_solution", the user can exclude observations from the final solution on the basis of excessive chi or residuals. The "testing" flags and thresholds are set in the GUI (although default values are defined here). */ public double chi_threshold = 2.8; public double optical_residual_threshold = 2.3*Math.PI/(180.0*3600.0); public double delay_residual_threshold = 30.0E-6; public double doppler_residual_threshold = 30.0; public boolean chi_testing = true; public boolean residual_testing = false; public int[] excluded_observations = new int[20001]; /* The GUI provides the user with the ability to evaluate collision or near-miss events with all bodies, or Earth-only. To save processing time, the earth_only boolean will tell detect_collision which planets to test. */ boolean earth_only = true; /* Variant collision analysis processing time is unreasonably long when the highest integration accuracy is used. To speed processing, I will allow the user to select between high (6e-9) and medium (1e-5) integration accuracy; this selection will also change the precision with which near-miss and collision events are localized. */ public double integration_error = 6e-9; public double event_localization = 1e-5; /* Since the dopri8 integration coefficients are reinitialized for every integration step, I'm going to try making them class variables; they will be initialized only once. */ double[][] dopri8_a = new double[14][14]; double[] dopri8_c = new double[14]; double[] dopri8_b = new double[14]; double[] dopri8_bhat = new double[14]; /* A colleague once said, "Show me a core, and I'll show you a matrix that fills it." Since the number of observations can be quite large, we must divide matrix_T in the least-squares algorithm into several subsections, each of which can be manipulated separately. */ public double[][] matrix_Ti = new double[501][18001]; public double[][] matrix_Tk = new double[501][18001]; /* Instance Variables Comprising Database for Each Minor Planet */ /* comet_identifier is a boolean that indicates whether the minor planet is a comet */ public boolean comet_identifier = false; /* The following variables and arrays contain the orbital solution (state vector and covariances at epoch) for the minor planet */ /* epoch_time is the TDB jultime to which the calculated or a'priori minor planet state vector applies */ public double epoch_time = 0; /* epoch_r and epoch_rprime are the calculated or a'priori minor planet state vector */ public double[] epoch_r = new double[4]; public double[] epoch_rprime = new double[4]; /* A1_A2_DT is an array containing parameters that model a comet's outgassing thrust. A1 and A2 are the radial and transverse components respectively, and are on the order of 10^(-9) AU/(day)^2. DT represents an offset from perihelion (in days, < 0 if before, > 0 if after) in the peak of a comet's light (and water-vaporization) curve. It is small, on the order of 0-80 days, and is not always required. */ public double[] A1_A2_DT = new double[4]; /* big_p is the epoch error covariance matrix, either calculated using least squares, or constructed using assumed state vector covariances. Diagonal entries 1-3 represent squared covariances in epoch position, entries 4-6 represent squared covariances in epoch velocity, and entries 7 and 8 represent squared covariances in A1 and A2. */ public double[][] big_p = new double[9][9]; /* absolute_magnitude_H is the mean absolute visual magnitude of the asteroid. */ public double absolute_magnitude_H = 0; /* slope_parameter_G is used to calculate the visual magnitude of the asteroid. */ public double slope_parameter_G = 0; /* comet_absolute_magnitude is the mean absolute visual magnitude of the comet */ public double comet_absolute_magnitude = 0; /* comet_slope_parameter is used to calculate the visual magnitude of the comet */ public double comet_slope_parameter = 0; /* The following variables and arrays contain the actual optical, radar delay and radar doppler observations of the minor planet */ /* Note that one optical observation consists of both a ra and a dec */ public int number_of_optical_observations = 0; public int number_of_delay_observations = 0; public int number_of_doppler_observations = 0; /* observation type is an integer array containing the type of each observation: = 1 for optical; = 2 for delay; = 3 for doppler; */ public int[] observation_type = new int[20001]; /* incremental_state_vector_epoch is an array containing the TDB times of each observation: = TDB time of optical observation, or = TDB receiver time for radar observations. Note: incremental_state_vector_epoch[0] = epoch_time; */ public double[] incremental_state_vector_epoch = new double[20001]; /* In method "least_squares_solution", a call to "get_residuals" requires calculation of the nominal state vector and observables at the TDB time of each observation. These state vectors are retained for use in "least_squares" and elsewhere, so as to avoid recalculation. */ /* incremental_state_vector_r is an array containing the state vector positions at the TDB times of each observation */ public double[][] incremental_state_vector_r = new double[20001][4]; /* incremental_state_vector_rprime is an array containing the state vector velocities at the TDB times of each observation */ public double[][] incremental_state_vector_rprime = new double[20001][4]; /* OPTICAL OBSERVATIONS */ /* ra and dec are in radians */ public double[] optical_ra = new double[20001]; public double[] ra_dev = new double[20001]; public double[] optical_dec = new double[20001]; public double[] dec_dev = new double[20001]; /* TDT times of optical observations */ public double[] optical_time = new double[20001]; /* visual_magnitude[i] is the observed visual-range brightness of the minor planet, made as part of optical observation i. A value of "-99" means that no visual-range brightness measurement was made. (Note that we will only solve for brightness parameters in the case of non-thrusting minor planets [asteroids]; thus, "-99" magnitude is unlikely to occur naturally.) An uncertainty of 0.1 will be assumed for all magnitude observations. */ public double[] visual_magnitude = new double[20001]; public double[] visual_magnitude_dev = new double[20001]; /* observation_geocentric is a boolean array; entry i = true iff optical observation i was geocentric */ public boolean[] observation_geocentric = new boolean[20001]; /* observation_latitude contains the geodetic latitude in radians from which each optical observation was made */ public double[] observation_latitude = new double[20001]; /* observation_longitude contains the east longitude in radians from which each optical observation was made */ public double[] observation_longitude = new double[20001]; /* observation_altitude contains the altitude in meters from which each optical observation was made */ public double[] observation_altitude = new double[20001]; /* observation_EOP[i][j] contains the Earth Orientation Parameters for optical observation i: - observation_EOP[i][1] contains the value of xp in radians at the time of observation i; - observation_EOP[i][2] contains the value of yp in radians at the time of observation i; - observation_EOP[i][3] contains the value of (UT1 - TAI) in seconds at the time of observation i; - observation_EOP[i][4] contains the value of delta psi in radians at the time of observation i; - observation_EOP[i][5] contains the value of delta epsilon in radians at the time of observation i; - observation_EOP[i][6] contains the value of mean obliquity in radians at the time of observation i; - observation_EOP[i][7] contains the value of Earth's angular velocity in radians per second at the time of observation i; */ public double[][] observation_EOP = new double[20001][8]; /* observation_r contains the barycentric rectangular equatorial position (in AUs) of the site from which each optical observation was made */ public double[][] observation_r = new double[20001][4]; /* observation_rprime contains the barycentric rectangular equatorial velocity (in AU/day) of the site from which each optical observation was made */ public double[][] observation_rprime = new double[20001][4]; /* DELAY OBSERVATIONS */ /* delay is in seconds */ public double[] radar_delay = new double[20001]; public double[] delay_dev = new double[20001]; /* UTC Julian dates of radar delay receiver observations */ public double[] radar_delay_receiver_time = new double[20001]; /* radar_delay_receiver_latitude contains the geodetic latitude in radians at which each radar delay receiver is located */ public double[] radar_delay_receiver_latitude = new double[20001]; /* radar_delay_receiver_longitude contains the east longitude in radians at which each radar delay receiver is located */ public double[] radar_delay_receiver_longitude = new double[20001]; /* radar_delay_receiver_altitude contains the altitude in meters at which each radar delay receiver is located */ public double[] radar_delay_receiver_altitude = new double[20001]; /* radar_delay_receiver_EOP[i][j] contains the Earth Orientation Parameters for the receiver in radar delay observation i, in format identical to observation_EOP */ public double[][] radar_delay_receiver_EOP = new double[20001][8]; /* radar_delay_receiver_tdb_minus_utc contains the value of (TDB-UTC) in seconds for the receiver from which each radar delay observation was made */ public double[] radar_delay_receiver_tdb_minus_utc = new double[20001]; /* radar_delay_receiver_r contains the barycentric rectangular equatorial position (in AUs) of the receiver site from which each radar delay observation was made */ public double[][] radar_delay_receiver_r = new double[20001][4]; /* radar_delay_receiver_rprime contains the barycentric rectangular equatorial velocity (in AU/day) of the receiver site from which each radar delay observation was made */ public double[][] radar_delay_receiver_rprime = new double[20001][4]; /* radar_delay_transmitter_latitude contains the geodetic latitude in radians at which each radar delay transmitter is located */ public double[] radar_delay_transmitter_latitude = new double[20001]; /* radar_delay_transmitter_longitude contains the east longitude in radians at which each radar delay transmitter is located */ public double[] radar_delay_transmitter_longitude = new double[20001]; /* radar_delay_transmitter_altitude contains the altitude in meters at which each radar delay transmitter is located */ public double[] radar_delay_transmitter_altitude = new double[20001]; /* radar_delay_transmitter_EOP[i][j] contains the Earth Orientation Parameters for the transmitter in radar delay observation i, in format identical to observation_EOP */ public double[][] radar_delay_transmitter_EOP = new double[20001][8]; /* radar_delay_transmitter_tdb_minus_utc contains the value of (TDB-UTC) in seconds for the transmitter from which each radar delay observation was made */ public double[] radar_delay_transmitter_tdb_minus_utc = new double[20001]; /* radar_delay_transmitter_frequency contains the frequency in Hertz for the transmitter from which each radar delay observation was made */ public double[] radar_delay_transmitter_frequency = new double[20001]; /* DOPPLER OBSERVATIONS */ /* doppler is in Hz */ public double[] radar_doppler = new double[20001]; public double[] doppler_dev = new double[20001]; /* UTC Julian dates of radar doppler receiver observations */ public double[] radar_doppler_receiver_time = new double[20001]; /* radar_doppler_receiver_latitude contains the geodetic latitude in radians at which each radar doppler receiver is located */ public double[] radar_doppler_receiver_latitude = new double[20001]; /* radar_doppler_receiver_longitude contains the east longitude in radians at which each radar doppler receiver is located */ public double[] radar_doppler_receiver_longitude = new double[20001]; /* radar_doppler_receiver_altitude contains the altitude in meters at which each radar doppler receiver is located */ public double[] radar_doppler_receiver_altitude = new double[20001]; /* radar_doppler_receiver_EOP[i][j] contains the Earth Orientation Parameters for the receiver in radar doppler observation i, in format identical to observation_EOP */ public double[][] radar_doppler_receiver_EOP = new double[20001][8]; /* radar_doppler_receiver_tdb_minus_utc contains the value of (TDB-UTC) in seconds for the receiver from which each radar doppler observation was made */ public double[] radar_doppler_receiver_tdb_minus_utc = new double[20001]; /* radar_doppler_receiver_r contains the barycentric rectangular equatorial position (in AUs) of the receiver site from which each radar doppler observation was made */ public double[][] radar_doppler_receiver_r = new double[20001][4]; /* radar_doppler_receiver_rprime contains the barycentric rectangular equatorial velocity (in AU/day) of the receiver site from which each radar doppler observation was made */ public double[][] radar_doppler_receiver_rprime = new double[20001][4]; /* radar_doppler_transmitter_latitude contains the geodetic latitude in radians at which each radar doppler transmitter is located */ public double[] radar_doppler_transmitter_latitude = new double[20001]; /* radar_doppler_transmitter_longitude contains the east longitude in radians at which each radar doppler transmitter is located */ public double[] radar_doppler_transmitter_longitude = new double[20001]; /* radar_doppler_transmitter_altitude contains the altitude in meters at which each radar doppler transmitter is located */ public double[] radar_doppler_transmitter_altitude = new double[20001]; /* radar_doppler_transmitter_EOP[i][j] contains the Earth Orientation Parameters for the transmitter in radar doppler observation i, in format identical to observation_EOP */ public double[][] radar_doppler_transmitter_EOP = new double[20001][8]; /* radar_doppler_transmitter_tdb_minus_utc contains the value of (TDB-UTC) in seconds for the transmitter from which each radar doppler observation was made */ public double[] radar_doppler_transmitter_tdb_minus_utc = new double[20001]; /* radar_doppler_transmitter_frequency contains the frequency in Hertz for the transmitter from which each radar doppler observation was made */ public double[] radar_doppler_transmitter_frequency = new double[20001]; /* minor_planet_radius is an estimate of the radius of the minor planet, in AU; it is used in radar calculations to estimate the bounce point */ public double minor_planet_radius = 0; public static void main(String args[]) { MinorPlanet testBody = new MinorPlanet(); } public void least_squares_solution(double[] residuals, double[] rms_residuals) { /* Method to calculate the least squares best-fit orbit of a minor planet, given optical, radar delay, and or radar doppler observations. The process is repeated until the magnitude of the corrections to the epoch state vector is less than some percentage of their corresponding covariances. The method presumes that an a priori state vector has been selected, and resides in epoch_r, epoch_rprime, and epoch_time. The method also assumes that the observations have been ordered from earliest to most recent. The output will be the epoch state vector, and the epoch error covariance matrix (all of which are instance variables of the MinorPlanet class). We will assume a maximum of 20000 observations. Provision is made for excluding observations whose chi or residuals exceed user-defined thresholds; although these observations will not be deleted, they will not be used in the final iterations of the algorithm. Adapted from: D.K. Yeomans, S.J. Ostro, and P.W. Chodas (1987). "Radar Astrometry of Near-Earth Asteroids". Astronomical Journal 94, 189-200. Modified to the "Square Root Information Filter" algorithm from Gerald J. Bierman (1974), "Sequential Square Root Filtering and Smoothing of Discrete Linear Systems", Automatica 10, 147-158. NOTE: In case of divergence, epoch_r[0] will be returned as 6.0. */ int i = 0, j = 0, k = 0, index = 0, retained_index = 0, counter = 0; boolean convergence = false, divergence = false, edit = false; double rms_optical = 0, rms_delay = 0, rms_doppler = 0; double old_rms_optical = 0, old_rms_delay = 0, old_rms_doppler = 0; double convergence_factor = 0.10; double delta = 0, delay = 0, doppler = 0; double observation_residual = 0, observation_chi = 0; int number_of_observations = number_of_optical_observations + number_of_delay_observations + number_of_doppler_observations; double[] nominal_r = new double[4]; double[] nominal_rprime = new double[4]; double[] initial_r = new double[4]; double[] initial_rprime = new double[4]; double[][] r_variant = new double[7][4]; double[][] rprime_variant = new double[7][4]; double[] this_r = new double[4]; double[] this_rprime = new double[4]; double[] next_r = new double[4]; double[] next_rprime = new double[4]; double[] variant_delta = new double[7]; double[] dev = new double[40002]; double[][] matrix_a = new double[40002][7]; double[][] collision_variant = new double[2000][9]; double[][] phi = new double[7][7]; double[] inter = new double[7]; double[] sensitivity_ra = new double[7]; double[] sensitivity_dec = new double[7]; double[] sensitivity_delay = new double[7]; double[] sensitivity_doppler = new double[7]; double[] radec = new double[3]; double[] row = new double[7]; double[] corrections = new double[7]; double[] observer_r = new double[4]; double[] observer_rprime = new double[4]; double[] delay_doppler = new double[3]; double[] receiver_r = new double[4]; double[] receiver_rprime = new double[4]; double[] EOP = new double[8]; double[][] matrix_R = new double[7][7]; double[] matrix_z = new double[7]; /* How to initialize? */ double[][] matrix_temp = new double[40002][7]; double[] retained_residuals = new double[40002]; double[] initial_residuals = new double[40002]; double[] initial_rms_residuals = new double[4]; double[][] initial_big_p = new double[9][9]; /* Initialize all observations as non-excluded for get_residuals */ for (i=1;i<=20000;i++) excluded_observations[i] = 0; /* Using the a priori state vector, determine initial residuals, nominal observables, and the state vector at the tdb time of each observation. Save the a priori state vector in case of divergence. */ get_residuals(residuals,rms_residuals); rms_optical = rms_residuals[1]; rms_delay = rms_residuals[2]; rms_doppler = rms_residuals[3]; for (i=1;i<=3;i++) { initial_r[i] = epoch_r[i]; initial_rprime[i] = epoch_rprime[i]; initial_rms_residuals[i] = rms_residuals[i]; } for (i=1;i<=8;i++) { for (j=1;j<=8;j++) initial_big_p[i][j] = big_p[i][j]; } index = 2*number_of_optical_observations + number_of_delay_observations + number_of_doppler_observations; for (i=1;i<=index;i++) { initial_residuals[i] = residuals[i]; } /* The convergence condition is: The RMS residuals on the current loop must differ by no more than some percentage from the RMS residuals on the preceding loop. (Note that this condition was not my original choice. I tried requiring that the corrections to the epoch state vector must be less than some percentage of their corresponding covariances. However, real-world data is not sufficiently smooth. This new condition uses the observed fact that the RMS residuals decrease with successive loops until the next loop results in almost no change. At this point, the solution has converged.) Clearly, this condition will not be satisfied initially, so we will go through the loop at least once. NOTE: MUST ACCOUNT FOR DIVERGENCE, AS MEASURED BY INCREASING RMS RESIDUALS. If each of the RMS residuals increase from one loop to the next, the divergence flag will be set, and the original solution returned. */ /* Note: For the asteroid ephemeris application, we're starting with a good orbit provided by the MPC. We will therefore begin editing errant observations immediately, to promote convergence. */ edit = true; while ((convergence == false) && (divergence == false)) { /* Clear the counters */ index = 0; retained_index = 0; counter ++; /* Use the observations to generate rows of matrix_a */ for (i=1;i<=number_of_observations;i++) { /* Assume that the observation will be excluded; if not, the flag will be reset below. */ excluded_observations[i] = 1; /* Generate the state transition matrix corresponding to this observation time */ for (j=1;j<=3;j++) { nominal_r[j] = incremental_state_vector_r[i][j]; nominal_rprime[j] = incremental_state_vector_rprime[i][j]; } /* Initialize the variant trajectories */ if (i == 1) { for (j=1;j<=6;j++) { for (k=1;k<=3;k++) { r_variant[j][k] = epoch_r[k]; rprime_variant[j][k] = epoch_rprime[k]; } } for (j=1;j<=3;j++) { variant_delta[j] = 0.000001*epoch_r[j]; r_variant[j][j] = r_variant[j][j] + variant_delta[j]; variant_delta[j+3] = 0.000001*epoch_rprime[j]; rprime_variant[j+3][j] = rprime_variant[j+3][j] + variant_delta[j+3]; } } /* Propogate the variant trajectories, and calculate the transition matrix entries */ for (j=1;j<=6;j++) { for (k=1;k<=3;k++) { this_r[k] = r_variant[j][k]; this_rprime[k] = rprime_variant[j][k]; } update(this_r,this_rprime,incremental_state_vector_epoch[i-1],next_r,next_rprime,incremental_state_vector_epoch[i],false,collision_variant); for (k=1;k<=3;k++) { r_variant[j][k] = next_r[k]; rprime_variant[j][k] = next_rprime[k]; } for (k=1;k<=3;k++) { phi[k][j] = (r_variant[j][k] - nominal_r[k])/variant_delta[j]; phi[k+3][j] = (rprime_variant[j][k] - nominal_rprime[k])/variant_delta[j]; } } /* Generate the sensitivity matrix for this observation */ if (observation_type[i] == 1) { /* Optical observation */ for (j=1;j<=3;j++) { observer_r[j] = observation_r[i][j]; observer_rprime[j] = observation_rprime[i][j]; } for (j=1;j<=3;j++) { delta = 0.000001*nominal_r[j]; nominal_r[j] = nominal_r[j] + delta; get_observed_radec(nominal_r,nominal_rprime,incremental_state_vector_epoch[i],observation_geocentric[i],observation_latitude[i],observation_longitude[i],observation_altitude[i],observation_EOP[i][3],observer_r,observer_rprime,optical_time[i],radec); sensitivity_ra[j] = radec[1] - nominal_observable[index+1]; while (Math.abs(sensitivity_ra[j]) > Math.abs(sensitivity_ra[j] - 2.0*Math.PI)) sensitivity_ra[j] = sensitivity_ra[j] - 2.0*Math.PI; while (Math.abs(sensitivity_ra[j]) > Math.abs(sensitivity_ra[j] + 2.0*Math.PI)) sensitivity_ra[j] = sensitivity_ra[j] + 2.0*Math.PI; sensitivity_ra[j] = sensitivity_ra[j]/delta; sensitivity_dec[j] = (radec[2] - nominal_observable[index+2])/delta; nominal_r[j] = nominal_r[j] - delta; delta = 0.000001*nominal_rprime[j]; nominal_rprime[j] = nominal_rprime[j] + delta; get_observed_radec(nominal_r,nominal_rprime,incremental_state_vector_epoch[i],observation_geocentric[i],observation_latitude[i],observation_longitude[i],observation_altitude[i],observation_EOP[i][3],observer_r,observer_rprime,optical_time[i],radec); sensitivity_ra[j+3] = radec[1] - nominal_observable[index+1]; while (Math.abs(sensitivity_ra[j+3]) > Math.abs(sensitivity_ra[j+3] - 2.0*Math.PI)) sensitivity_ra[j+3] = sensitivity_ra[j+3] - 2.0*Math.PI; while (Math.abs(sensitivity_ra[j+3]) > Math.abs(sensitivity_ra[j+3] + 2.0*Math.PI)) sensitivity_ra[j+3] = sensitivity_ra[j+3] + 2.0*Math.PI; sensitivity_ra[j+3] = sensitivity_ra[j+3]/delta; sensitivity_dec[j+3] = (radec[2] - nominal_observable[index+2])/delta; nominal_rprime[j] = nominal_rprime[j] - delta; } /* Calculate the residuals and chi; if observations meet the applicable thresholds, calculate the resulting rows of matrix_a */ observation_residual = Math.acos(Math.sin(nominal_observable[index+2])*Math.sin(optical_dec[i]) + Math.cos(nominal_observable[index+2])*Math.cos(optical_dec[i])*Math.cos(nominal_observable[index+1] - optical_ra[i])); if ( (ra_dev[i] > 0.0) && (dec_dev[i] > 0.0) ) observation_chi = Math.sqrt(Math.pow(residuals[index+1]*Math.cos(optical_dec[i])/ra_dev[i],2.0) + Math.pow(residuals[index+2]/dec_dev[i],2.0)); else observation_chi = 0.0; if (((edit) && (chi_testing) && (observation_chi < chi_threshold) && (ra_dev[i] < 5.0*Math.PI/(3600.0*180.0)) && (dec_dev[i] < 5.0*Math.PI/(3600.0*180.0))) || ((edit) && (residual_testing) && (Math.abs(observation_residual) < optical_residual_threshold) && (ra_dev[i] < 5.0*Math.PI/(3600.0*180.0)) && (dec_dev[i] < 5.0*Math.PI/(3600.0*180.0))) || (!edit)) { /* Observation is not excluded; reset flag */ excluded_observations[i] = 0; for (j=1;j<=6;j++) { row[j] = 0.0; for (k=1;k<=6;k++) row[j] = row[j] + sensitivity_ra[k]*phi[k][j]; matrix_a[retained_index+1][j] = row[j]; } for (j=1;j<=6;j++) { row[j] = 0.0; for (k=1;k<=6;k++) row[j] = row[j] + sensitivity_dec[k]*phi[k][j]; matrix_a[retained_index+2][j] = row[j]; } /* Fill the dev array with standard deviations */ dev[retained_index+1] = ra_dev[i]; dev[retained_index+2] = dec_dev[i]; /* Add these residuals to the list of those being retained */ retained_residuals[retained_index+1] = residuals[index+1]; retained_residuals[retained_index+2] = residuals[index+2]; /* Increment the retained observation counter */ retained_index = retained_index + 2; } /* Increment the index */ index = index + 2; } else if (observation_type[i] == 2) { /* Delay observation */ for (j=1;j<=3;j++) { receiver_r[j] = radar_delay_receiver_r[i][j]; receiver_rprime[j] = radar_delay_receiver_rprime[i][j]; } for (j=1;j<=7;j++) EOP[j] = radar_delay_transmitter_EOP[i][j]; for (j=1;j<=3;j++) { delta = 0.000001*nominal_r[j]; nominal_r[j] = nominal_r[j] + delta; get_observed_radar_delay_doppler(radar_delay_receiver_latitude[i],radar_delay_receiver_longitude[i],radar_delay_receiver_altitude[i],receiver_r,receiver_rprime,radar_delay_receiver_time[i],radar_delay_transmitter_latitude[i],radar_delay_transmitter_longitude[i],radar_delay_transmitter_altitude[i],EOP,radar_delay_transmitter_tdb_minus_utc[i],radar_delay_receiver_tdb_minus_utc[i],radar_delay_transmitter_frequency[i],minor_planet_radius,nominal_r,nominal_rprime,incremental_state_vector_epoch[i],delay_doppler); delay = delay_doppler[1]; sensitivity_delay[j] = (delay - nominal_observable[index+1])/delta; nominal_r[j] = nominal_r[j] - delta; delta = 0.000001*nominal_rprime[j]; nominal_rprime[j] = nominal_rprime[j] + delta; get_observed_radar_delay_doppler(radar_delay_receiver_latitude[i],radar_delay_receiver_longitude[i],radar_delay_receiver_altitude[i],receiver_r,receiver_rprime,radar_delay_receiver_time[i],radar_delay_transmitter_latitude[i],radar_delay_transmitter_longitude[i],radar_delay_transmitter_altitude[i],EOP,radar_delay_transmitter_tdb_minus_utc[i],radar_delay_receiver_tdb_minus_utc[i],radar_delay_transmitter_frequency[i],minor_planet_radius,nominal_r,nominal_rprime,incremental_state_vector_epoch[i],delay_doppler); delay = delay_doppler[1]; sensitivity_delay[j+3] = (delay - nominal_observable[index+1])/delta; nominal_rprime[j] = nominal_rprime[j] - delta; } /* Calculate the residual and chi; if observations meet the applicable threshold, calculate the resulting row of matrix_a */ observation_residual = Math.abs(residuals[index+1]); if (delay_dev[i] != 0.0) observation_chi = Math.abs(residuals[index+1]/delay_dev[i]); else observation_chi = 0.0; if (((edit) && (chi_testing) && (observation_chi < chi_threshold)) || ((edit) && (residual_testing) && (Math.abs(observation_residual) < delay_residual_threshold)) || (!edit)) { /* Observation is not excluded; reset flag */ excluded_observations[i] = 0; for (j=1;j<=6;j++) { row[j] = 0.0; for (k=1;k<=6;k++) row[j] = row[j] + sensitivity_delay[k]*phi[k][j]; matrix_a[retained_index+1][j] = row[j]; } /* Fill the dev array with standard deviations */ dev[retained_index+1] = delay_dev[i]; /* Add this residual to the list of those being retained */ retained_residuals[retained_index+1] = residuals[index+1]; /* Increment the retained observation counter */ retained_index = retained_index + 1; } /* Increment the index */ index = index + 1; } else if (observation_type[i] == 3) { /* Doppler observation */ for (j=1;j<=3;j++) { receiver_r[j] = radar_doppler_receiver_r[i][j]; receiver_rprime[j] = radar_doppler_receiver_rprime[i][j]; } for (j=1;j<=7;j++) { EOP[j] = radar_doppler_transmitter_EOP[i][j]; } for (j=1;j<=3;j++) { delta = 0.000001*nominal_r[j]; nominal_r[j] = nominal_r[j] + delta; get_observed_radar_delay_doppler(radar_doppler_receiver_latitude[i],radar_doppler_receiver_longitude[i],radar_doppler_receiver_altitude[i],receiver_r,receiver_rprime,radar_doppler_receiver_time[i],radar_doppler_transmitter_latitude[i],radar_doppler_transmitter_longitude[i],radar_doppler_transmitter_altitude[i],EOP,radar_doppler_transmitter_tdb_minus_utc[i],radar_doppler_receiver_tdb_minus_utc[i],radar_doppler_transmitter_frequency[i],minor_planet_radius,nominal_r,nominal_rprime,incremental_state_vector_epoch[i],delay_doppler); doppler = delay_doppler[2]; sensitivity_doppler[j] = (doppler - nominal_observable[index+1])/delta; nominal_r[j] = nominal_r[j] - delta; delta = 0.000001*nominal_rprime[j]; nominal_rprime[j] = nominal_rprime[j] + delta; get_observed_radar_delay_doppler(radar_doppler_receiver_latitude[i],radar_doppler_receiver_longitude[i],radar_doppler_receiver_altitude[i],receiver_r,receiver_rprime,radar_doppler_receiver_time[i],radar_doppler_transmitter_latitude[i],radar_doppler_transmitter_longitude[i],radar_doppler_transmitter_altitude[i],EOP,radar_doppler_transmitter_tdb_minus_utc[i],radar_doppler_receiver_tdb_minus_utc[i],radar_doppler_transmitter_frequency[i],minor_planet_radius,nominal_r,nominal_rprime,incremental_state_vector_epoch[i],delay_doppler); doppler = delay_doppler[2]; sensitivity_doppler[j+3] = (doppler - nominal_observable[index+1])/delta; nominal_rprime[j] = nominal_rprime[j] - delta; } /* Calculate the residual and chi; if observations meet the applicable thresholds, calculate the resulting row of matrix_a */ observation_residual = Math.abs(residuals[index+1]); if (doppler_dev[i] != 0.0) observation_chi = Math.abs(residuals[index+1]/doppler_dev[i]); else observation_chi = 0.0; if (((edit) && (chi_testing) && (observation_chi < chi_threshold)) || ((edit) && (residual_testing) && (Math.abs(observation_residual) < doppler_residual_threshold)) || (!edit)) { /* Observation is not excluded; reset flag */ excluded_observations[i] = 0; for (j=1;j<=6;j++) { row[j] = 0.0; for (k=1;k<=6;k++) row[j] = row[j] + sensitivity_doppler[k]*phi[k][j]; matrix_a[retained_index+1][j] = row[j]; } /* Fill the dev array with standard deviations */ dev[retained_index+1] = doppler_dev[i]; /* Add this residual to the list of those being retained */ retained_residuals[retained_index+1] = residuals[index+1]; /* Increment the retained observation counter */ retained_index = retained_index + 1; } /* Increment the index */ index = index + 1; } } /* At this point, matrix A is complete. We proceed to convert to a unity covariance form */ /* Should be commented out */ for (i=1;i<=retained_index;i++) { for (j=1;j<=6;j++) matrix_temp[i][j] = matrix_a[i][j]/dev[i]; retained_residuals[i] = retained_residuals[i]/dev[i]; } for (i=1;i<=retained_index;i++) { for (j=1;j<=6;j++) matrix_a[i][j] = matrix_temp[i][j]; } /* Construct an orthogonal transform T using Givens rotations, and left multiply. matrix_a will now be upper triangular; isolate the upper 6x6 as matrix_R, and isolate the top 6x1 of T x residuals as matrix_z. */ givens_qr(retained_index, 6, matrix_a); read_matrix_Tk(0); for (i=1;i<=6;i++) { matrix_z[i] = 0; for (j=1;j<=retained_index;j++) matrix_z[i] = matrix_z[i] + matrix_Tk[i][j]*retained_residuals[j]; } for (i=1;i<=6;i++) { for (j=1;j<=6;j++) matrix_R[i][j] = matrix_a[i][j]; } /* Calculate the state vector corrections and covariances. */ invert_nxn(6,matrix_R); for (i=1;i<=6;i++) { corrections[i] = 0; for (k=1;k<=6;k++) /* Square Root Information Filtering assumes that we are in unit covariance form. Since we're often unable to invert matrix_R in that form, we must make the conversion to unit covariance form now, so the results are valid. In this particular calculation, we would multiply matrix_R[i][k] by dev[k] (since matrix_a has been inverted), and we would divide matrix_z by dev[k]; since these precisely cancel, no action is really needed here. */ corrections[i] = corrections[i] + matrix_R[i][k]*matrix_z[k]; } for (i=1;i<=6;i++) { for (j=1;j<=6;j++) { big_p[i][j] = 0; for (k=1;k<=6;k++) /* Square Root Information Filtering assumes that we are in unit covariance form. Since we're often unable to invert matrix_R in that form, we must make the conversion to unit covariance form now, so the results are valid. In this particular calculation, we would multiply matrix_R[i][k] by dev[k] (since matrix_a has been inverted), and we would multiply matrix_R[j][k] by dev[k] (for the same reason). */ /* big_p[i][j] = big_p[i][j] + matrix_R[i][k]*matrix_R[j][k]*(dev[k]*dev[k]); */ big_p[i][j] = big_p[i][j] + matrix_R[i][k]*matrix_R[j][k]; } } /* Update the epoch state vector */ for (i=1;i<=3;i++) { epoch_r[i] = epoch_r[i] + corrections[i]; epoch_rprime[i] = epoch_rprime[i] + corrections[i+3]; } /* Recalculate the residuals using the updated epoch state vector */ get_residuals(residuals,rms_residuals); old_rms_optical = rms_optical; old_rms_delay = rms_delay; old_rms_doppler = rms_doppler; rms_optical = rms_residuals[1]; rms_delay = rms_residuals[2]; rms_doppler = rms_residuals[3]; System.out.println("old_rms_optical = " + old_rms_optical); System.out.println("old_rms_delay = " + old_rms_delay); System.out.println("old_rms_doppler = " + old_rms_doppler); /* Calculate the asteroid's absolute magnitude, slope parameter, and radius. This last is especially important when radar observations are used, to account for the "bounce point". */ solve_for_absolute_magnitude(); /* Based on the RMS residuals, evaluate the convergence and divergence conditions */ convergence = true; if (Math.abs(old_rms_optical - rms_optical) > (convergence_factor*old_rms_optical)) convergence = false; if (Math.abs(old_rms_delay - rms_delay) > (convergence_factor*old_rms_delay)) convergence = false; if (Math.abs(old_rms_doppler - rms_doppler) > (convergence_factor*old_rms_doppler)) convergence = false; if ((counter > 1) && (rms_optical >= 1.2*old_rms_optical) && (rms_delay >= 1.2*old_rms_delay) && (rms_doppler >= 1.2*old_rms_doppler)) divergence = true; if ((convergence) && (!edit) && ((chi_testing) || (residual_testing))) { edit = true; convergence = false; } if ((counter >= 1) && (divergence == false)) convergence = true; /* We may not have yet converged. However, if the residuals have been reduced, shouldn't we keep this state vector and covariance matrix as the initial data in case the next iteration diverges? */ if ((counter == 1) || ((rms_optical <= old_rms_optical) && (rms_delay <= old_rms_delay) && (rms_doppler <= old_rms_doppler))) { for (i=1;i<=3;i++) { initial_r[i] = epoch_r[i]; initial_rprime[i] = epoch_rprime[i]; initial_rms_residuals[i] = rms_residuals[i]; } for (i=1;i<=8;i++) { for (j=1;j<=8;j++) initial_big_p[i][j] = big_p[i][j]; } index = 2*number_of_optical_observations + number_of_delay_observations + number_of_doppler_observations; for (i=1;i<=index;i++) initial_residuals[i] = residuals[i]; } System.out.println("convergence = " + convergence); System.out.println("divergence = " + divergence); /* The loop is complete. If the convergence or divergence condition is satisfied, we will exit; if not, we will go around again */ } if (divergence == true) { /* Return initial orbit and residuals, and set epoch_r[0] = 6.0. */ for (i=1;i<=3;i++) { epoch_r[i] = initial_r[i]; epoch_rprime[i] = initial_rprime[i]; rms_residuals[i] = initial_rms_residuals[i]; } epoch_r[0] = 6.0; for (i=1;i<=6;i++) { for (j=1;j<=6;j++) big_p[i][j] = initial_big_p[i][j]; } index = 2*number_of_optical_observations + number_of_delay_observations + number_of_doppler_observations; for (i=1;i<=index;i++) residuals[i] = initial_residuals[i]; } else if (divergence == false) epoch_r[0] = 0.0; } public void get_residuals(double residual[], double rms_residual[]) { /* Method to calculate the (observed - computed) and RMS residuals in the optical, delay, and doppler observations of a minor planet, given the state vector (an instance variable in the MinorPlanet class). Note that, to improve efficiency, it is assumed that the observations are ordered from earliest to most recent. The algorithm will therefore integrate the epoch state vector successively from observation i to i+1, saving the state vector at observation i+1 in the instance variables "incremental_state_vector_r[i+1][j]" and "incremental_state_vector_rprime[i+1][j]". This incremental state vector can then be used in method "least_squares" to represent the nominal state vector at the TDB time of observation i+1 (incremental_state_vector_epoch[i+1]). */ double delay = 0, doppler = 0, rms_optical = 0, rms_delay = 0, rms_doppler = 0; double[] radec = new double[3]; double[] observer_r = new double[4]; double[] observer_rprime = new double[4]; double[] delay_doppler = new double[3]; double[] receiver_r = new double[4]; double[] receiver_rprime = new double[4]; double[] EOP = new double[8]; double[] state_vector_r = new double[4]; double[] state_vector_rprime = new double[4]; double[] new_state_vector_r = new double[4]; double[] new_state_vector_rprime = new double[4]; double[][] collision_variant = new double[2000][9]; int i = 0, j = 0, k = 0, index = 0, number_of_excluded_optical_observations = 0, number_of_excluded_delay_observations = 0, number_of_excluded_doppler_observations = 0; int number_of_observations = number_of_optical_observations + number_of_delay_observations + number_of_doppler_observations; /* Initialize the required variables */ for (i=1;i<=3;i++) { incremental_state_vector_r[0][i] = epoch_r[i]; incremental_state_vector_rprime[0][i] = epoch_rprime[i]; state_vector_r[i] = epoch_r[i]; state_vector_rprime[i] = epoch_rprime[i]; } incremental_state_vector_epoch[0] = epoch_time; rms_optical = 0.0; rms_delay = 0.0; rms_doppler = 0.0; for (i=1;i<=number_of_observations;i++) { /* For each observation, determine the nominal state vector */ update(state_vector_r,state_vector_rprime,incremental_state_vector_epoch[i-1],new_state_vector_r,new_state_vector_rprime,incremental_state_vector_epoch[i],false,collision_variant); for (j=1;j<=3;j++) { incremental_state_vector_r[i][j] = new_state_vector_r[j]; incremental_state_vector_rprime[i][j] = new_state_vector_rprime[j]; state_vector_r[j] = new_state_vector_r[j]; state_vector_rprime[j] = new_state_vector_rprime[j]; } /* Depending on the observation type, determine the nominal observables, and the residuals */ if (observation_type[i] == 1) { /* Optical observation */ /* Note that method get_observed_radec is being sent the observer posvel vectors that were determined when the observation was entered. That method will see these non-zero vectors (possibly created using precise EOP values), and use them instead of calling method get_observer (which would use estimates of the EOPs to calculate observer posvel). */ for (j=1;j<=3;j++) { observer_r[j] = observation_r[i][j]; observer_rprime[j] = observation_rprime[i][j]; } get_observed_radec(new_state_vector_r,new_state_vector_rprime,incremental_state_vector_epoch[i],observation_geocentric[i],observation_latitude[i],observation_longitude[i],observation_altitude[i],observation_EOP[i][3],observer_r,observer_rprime,optical_time[i],radec); nominal_observable[index+1] = radec[1]; nominal_observable[index+2] = radec[2]; residual[index+1] = optical_ra[i] - radec[1]; while (Math.abs(residual[index+1]) > Math.abs(residual[index+1] - 2.0*Math.PI)) residual[index+1] = residual[index+1] - 2.0*Math.PI; while (Math.abs(residual[index+1]) > Math.abs(residual[index+1] + 2.0*Math.PI)) residual[index+1] = residual[index+1] + 2.0*Math.PI; /* System.out.println("residual[" + (index+1) + "] = " + residual[index+1]); */ residual[index+2] = optical_dec[i] - radec[2]; /* System.out.println("residual[" + (index+2) + "] = " + residual[index+2]); */ /* Don't add these residuals to the rms total if the observation is excluded */ if (excluded_observations[i] == 0) { rms_optical = rms_optical + Math.pow(residual[index+1],2); rms_optical = rms_optical + Math.pow(residual[index+2],2); } else { number_of_excluded_optical_observations = number_of_excluded_optical_observations + 1; } /* Increment index to reflect two observables (ra and dec) */ index = index + 2; } else if (observation_type[i] == 2) { /* Delay observation */ for (j=1;j<=3;j++) { receiver_r[j] = radar_delay_receiver_r[i][j]; receiver_rprime[j] = radar_delay_receiver_rprime[i][j]; } for (j=1;j<=7;j++) EOP[j] = radar_delay_transmitter_EOP[i][j]; get_observed_radar_delay_doppler(radar_delay_receiver_latitude[i],radar_delay_receiver_longitude[i],radar_delay_receiver_altitude[i],receiver_r,receiver_rprime,radar_delay_receiver_time[i],radar_delay_transmitter_latitude[i],radar_delay_transmitter_longitude[i],radar_delay_transmitter_altitude[i],EOP,radar_delay_transmitter_tdb_minus_utc[i],radar_delay_receiver_tdb_minus_utc[i],radar_delay_transmitter_frequency[i],minor_planet_radius,new_state_vector_r,new_state_vector_rprime,incremental_state_vector_epoch[i],delay_doppler); delay = delay_doppler[1]; nominal_observable[index+1] = delay; residual[index+1] = radar_delay[i] - delay; /* System.out.println("residual[" + (index+1) + "] = " + residual[index+1]); */ /* Don't add this residual to the rms total if the observation is excluded */ if (excluded_observations[i] == 0) { rms_delay = rms_delay + Math.pow(residual[index+1],2); } else { number_of_excluded_delay_observations = number_of_excluded_delay_observations + 1; } /* Increment index to relect one observable (delay) */ index = index + 1; } else if (observation_type[i] == 3) { /* Doppler observation */ for (j=1;j<=3;j++) { receiver_r[j] = radar_doppler_receiver_r[i][j]; receiver_rprime[j] = radar_doppler_receiver_rprime[i][j]; } for (j=1;j<=7;j++) EOP[j] = radar_doppler_transmitter_EOP[i][j]; get_observed_radar_delay_doppler(radar_doppler_receiver_latitude[i],radar_doppler_receiver_longitude[i],radar_doppler_receiver_altitude[i],receiver_r,receiver_rprime,radar_doppler_receiver_time[i],radar_doppler_transmitter_latitude[i],radar_doppler_transmitter_longitude[i],radar_doppler_transmitter_altitude[i],EOP,radar_doppler_transmitter_tdb_minus_utc[i],radar_doppler_receiver_tdb_minus_utc[i],radar_doppler_transmitter_frequency[i],minor_planet_radius,new_state_vector_r,new_state_vector_rprime,incremental_state_vector_epoch[i],delay_doppler); doppler = delay_doppler[2]; nominal_observable[index+1] = doppler; residual[index+1] = radar_doppler[i] - doppler; /* System.out.println("residual[" + (index+1) + "] = " + residual[index+1]); */ /* Don't add this residual to the rms total if the observation is excluded */ if (excluded_observations[i] == 0) { rms_doppler = rms_doppler + Math.pow(residual[index+1],2); } else { number_of_excluded_doppler_observations = number_of_excluded_doppler_observations + 1; } /* Increment index to relect one observable (doppler) */ index = index + 1; } } if ((number_of_optical_observations - number_of_excluded_optical_observations) > 0) rms_optical = Math.sqrt(rms_optical/(2.0*(number_of_optical_observations - number_of_excluded_optical_observations))); else rms_optical = 0; System.out.println("rms_optical = " + rms_optical); rms_residual[1] = rms_optical; if ((number_of_delay_observations - number_of_excluded_delay_observations) > 0) rms_delay = Math.sqrt(rms_delay/(number_of_delay_observations - number_of_excluded_delay_observations)); else rms_delay = 0; System.out.println("rms_delay = " + rms_delay); rms_residual[2] = rms_delay; if ((number_of_doppler_observations - number_of_excluded_doppler_observations) > 0) rms_doppler = Math.sqrt(rms_doppler/(number_of_doppler_observations - number_of_excluded_doppler_observations)); else rms_doppler = 0; System.out.println("rms_doppler = " + rms_doppler); rms_residual[3] = rms_doppler; } void get_observed_radar_delay_doppler(double receiver_latitude, double receiver_longitude, double receiver_altitude, double receiver_r[], double receiver_rprime[], double receiver_time, double transmitter_latitude, double transmitter_longitude, double transmitter_altitude, double transmitter_EOP[], double transmitter_tdb_minus_utc, double receiver_tdb_minus_utc, double transmitter_frequency, double asteroid_radius, double asteroid_r[], double asteroid_rprime[], double asteroid_epoch_time, double delay_doppler[]) { /* Algorithm to determine the expected radar delay time and doppler shift for an astrometric radar observation. Inputs: - receiver_latitude, transmitter_latitude are in radians; - receiver_longitude, transmitter_longitude, measured east, are in radians; - receiver_altitude, transmitter_altitude are in meters; - receiver_r[] is a barycentric equatorial position in AU - receiver_rprime[] is a barycentric equatorial velocity in AU/day; - receiver_time is a UTC Julian date; - transmitter_EOP[] are the Earth Observation Parameters at transmitter_time; - transmitter_tdb_minus_utc, receiver_tdb_minus_utc are in seconds; - transmitter_frequency is in Hz; - asteroid_radius is in AU; - asteroid_r, asteroid_rprime, asteroid_epoch_time is the barycentric equatorial state vector of the asteroid, given in AUs, AUs/day, and TDB Julian date respectively. Output: - delay_doppler[1] = the predicted UTC delay from transmitter to asteroid to receiver, in seconds; - delay_doppler[2] = the predicted Doppler frequency shift, in Hz; This algorithm follows the outline provided in "Asteroid and Comet Orbits Using Radar Data", Yeomans, Chodas, Keesey, and Ostro, in The Astronomical Journal, Vol. 103, number 1 (January 1992), pp. 303-317. */ double[] EOP = new double[8]; double[][] collision = new double[2000][9]; double[] r_at_bounce = new double[4]; double[] rprime_at_bounce = new double[4]; double[] rho_receiver = new double[4]; double[] rho_dot_receiver = new double[4]; double[] rho_transmitter = new double[4]; double[] rho_dot_transmitter = new double[4]; double[] dummy_rho = new double[4]; double[] sun_r = new double[4]; double[] sun_rprime = new double[4]; double[] transmitter_r = new double[4]; double[] transmitter_rprime = new double[4]; double receiver_tdb_time = 0, bounce_time = 0, tau_receiver_delay = 0, tau_transmitter_delay = 0, dist_receiver = 0, old_tau_delay = 0, tau_diff = 1, transmitter_time = 0, dist_transmitter = 0, dot_one = 0, dot_two = 0, heliodist1 = 0, heliodist2 = 0, speed_transmitter = 0, speed_receiver = 0, rho_dot_t = 0, rho_dot_r = 0, bigT = 0; int i = 0; /* Note: I feel obligated to provide default values for TDB-UTC. However, since the difference TDB-TDT is being calculated so cavalierly, I SINCERELY HOPE THIS IS NEVER USED!!! Nevertheless, it is more accurate than a default of zero. */ if (receiver_tdb_minus_utc == 0.0) { receiver_tdb_minus_utc = (tdt_to_tdb(receiver_time,receiver_latitude,receiver_longitude,receiver_altitude,transmitter_EOP[3]) - receiver_time) + (get_leap_second(receiver_time) + 32.184); } if (transmitter_tdb_minus_utc == 0.0) { transmitter_tdb_minus_utc = (tdt_to_tdb(receiver_time,receiver_latitude,receiver_longitude,receiver_altitude,transmitter_EOP[3]) - receiver_time) + (get_leap_second(receiver_time) + 32.184); } /* First, get the barycentric position and velocity of the receiver at receiver_time. */ /* The receiver coordinates determined when the observation was entered may have been passed as an argument; these may have used precise values of the EOPs, so we don't want to recalculate. On the other hand, the receiver coordinates may be zero; this method may be creating an ephemeris, so precise values of the EOPs may not be available. We'll test and act accordingly. */ /* Convert the UTC receiver time to TDB */ receiver_tdb_time = receiver_time + receiver_tdb_minus_utc/86400.0; if (receiver_r[1] == 0.0) { /* Receiver coordinates are unavailable; calculate them. */ /* Set EOPs to zero so that method get_radar_position will use estimates and defaults. */ for (i=1;i<=7;i++) EOP[i] = 0.0; /* get receiver barycentric equatorial position and velocity at TDB receiver_time */ get_observer_position(receiver_latitude,receiver_longitude,receiver_altitude,receiver_tdb_time,EOP,receiver_r,receiver_rprime); } /* If receiver_r[1] is not equal to 0.0, receiver coordinates are available; use them. */ /* The downleg iteration - determine bounce_time, and the asteroid's position */ /* As a first approximation, determine the asteroid's position at receiver_time, and estimate bounce_time */ update(asteroid_r,asteroid_rprime,asteroid_epoch_time,r_at_bounce,rprime_at_bounce,receiver_tdb_time,false,collision); dist_receiver = Math.sqrt(Math.pow((r_at_bounce[1] - receiver_r[1]),2) + Math.pow((r_at_bounce[2] - receiver_r[2]),2) + Math.pow((r_at_bounce[3] - receiver_r[3]),2)); tau_receiver_delay = (dist_receiver - asteroid_radius)/(speed_of_light/86400.0); /* Iteratively improve the estimate of bounce_time; repeat until successive values of tau_delay differ by no more than 0.05 microseconds. */ while (tau_diff > 5E-08) { old_tau_delay = tau_receiver_delay; bounce_time = receiver_tdb_time - tau_receiver_delay/86400.0; update(asteroid_r,asteroid_rprime,asteroid_epoch_time,r_at_bounce,rprime_at_bounce,bounce_time,false,collision); for (i=1;i<=3;i++) { rho_receiver[i] = r_at_bounce[i] - receiver_r[i]; rho_dot_receiver[i] = rprime_at_bounce[i] - receiver_rprime[i]; } dist_receiver = Math.sqrt(Math.pow(rho_receiver[1],2) + Math.pow(rho_receiver[2],2) + Math.pow(rho_receiver[3],2)); tau_receiver_delay = (dist_receiver - asteroid_radius)/(speed_of_light/86400.0) + radar_corrections(receiver_tdb_time,r_at_bounce,receiver_r,transmitter_frequency); rho_dot_r = (rho_receiver[1]*rho_dot_receiver[1] + rho_receiver[2]*rho_dot_receiver[2] + rho_receiver[3]*rho_dot_receiver[3])/dist_receiver; tau_diff = Math.abs(tau_receiver_delay - old_tau_delay); } /* The up-leg iteration - determine transmitter_time, */ /* Using an initial approximation of transmitter_time and transmitter_r and _rprime, find the distance from the transmitter at transmitter_time to the asteroid at bounce_time, and iteratively refine transmitter_time */ transmitter_time = bounce_time - tau_receiver_delay/86400.0; /* get transmitter barycentric equatorial position and velocity at TDB transmitter_time */ get_observer_position(transmitter_latitude,transmitter_longitude,transmitter_altitude,transmitter_time,transmitter_EOP,transmitter_r,transmitter_rprime); dist_transmitter = Math.sqrt(Math.pow((r_at_bounce[1] - transmitter_r[1]),2) + Math.pow((r_at_bounce[2] - transmitter_r[2]),2) + Math.pow((r_at_bounce[3] - transmitter_r[3]),2)); tau_transmitter_delay = (dist_transmitter - asteroid_radius)/(speed_of_light/86400.0); /* Iteratively improve the estimate of transmitter_time; repeat until successive values of tau_delay differ by no more than 0.05 microseconds. */ tau_diff = 1.0; while (tau_diff > 5E-08) { old_tau_delay = tau_transmitter_delay; transmitter_time = bounce_time - tau_transmitter_delay/86400.0; get_observer_position(transmitter_latitude,transmitter_longitude,transmitter_altitude,transmitter_time,transmitter_EOP,transmitter_r,transmitter_rprime); for (i=1;i<=3;i++) { rho_transmitter[i] = r_at_bounce[i] - transmitter_r[i]; rho_dot_transmitter[i] = rprime_at_bounce[i] - transmitter_rprime[i]; } dist_transmitter = Math.sqrt(Math.pow(rho_transmitter[1],2) + Math.pow(rho_transmitter[2],2) + Math.pow(rho_transmitter[3],2)); tau_transmitter_delay = (dist_transmitter - asteroid_radius)/(speed_of_light/86400.0) + radar_corrections(transmitter_time,r_at_bounce,transmitter_r,transmitter_frequency); rho_dot_t = (rho_transmitter[1]*rho_dot_transmitter[1] + rho_transmitter[2]*rho_dot_transmitter[2] + rho_transmitter[3]*rho_dot_transmitter[3])/dist_transmitter; tau_diff = Math.abs(tau_transmitter_delay - old_tau_delay); } /* Calculate the predicted delay time, and Doppler shift */ delay_doppler[1] = (tau_receiver_delay + tau_transmitter_delay) + transmitter_tdb_minus_utc - receiver_tdb_minus_utc; dot_one = rho_transmitter[1]*transmitter_rprime[1] + rho_transmitter[2]*transmitter_rprime[2] + rho_transmitter[3]*transmitter_rprime[3]; dot_two = rho_receiver[1]*rprime_at_bounce[1] + rho_receiver[2]*rprime_at_bounce[2] + rho_receiver[3]*rprime_at_bounce[3]; get_planet_posvel(transmitter_time,11,sun_r,sun_rprime); heliodist1 = Math.sqrt(Math.pow((transmitter_r[1] - sun_r[1]),2) + Math.pow((transmitter_r[2] - sun_r[2]),2) + Math.pow((transmitter_r[3] - sun_r[3]),2)); get_planet_posvel(receiver_tdb_time,11,sun_r,sun_rprime); heliodist2 = Math.sqrt(Math.pow((receiver_r[1] - sun_r[1]),2) + Math.pow((receiver_r[2] - sun_r[2]),2) + Math.pow((receiver_r[3] - sun_r[3]),2)); speed_transmitter = Math.sqrt(Math.pow(transmitter_rprime[1],2) + Math.pow(transmitter_rprime[2],2) + Math.pow(transmitter_rprime[3],2)); speed_receiver = Math.sqrt(Math.pow(receiver_rprime[1],2) + Math.pow(receiver_rprime[2],2) + Math.pow(receiver_rprime[3],2)); delay_doppler[2] = -transmitter_frequency*(rho_dot_t + rho_dot_r)/speed_of_light - (transmitter_frequency/Math.pow(speed_of_light,2))*(rho_dot_t*dot_one/dist_transmitter - rho_dot_r*dot_two/dist_receiver - rho_dot_t*rho_dot_r + Math.pow(kappa,2)*(1.0/heliodist1 - 1.0/heliodist2) + 0.5*(Math.pow(speed_transmitter,2) - Math.pow(speed_receiver,2))); } double radar_corrections(double jultime,double asteroid_r[],double site_r[],double frequency) { /* This method computes corrections to radar delay times due to general relativity, the solar corona, and the Earth's troposphere and ionosphere. Inputs: - jultime is a TDB Julian date corresponding to the epoch of site_r[]; - asteroid_r[] and site_r[] are barycentric equatorial positions in AU; asteroid_r[] is evaluated at bounce_time, while site_r[] is evaluated at jultime (either transmitter_time or receiver_time); - frequency is in Hz; Outputs: - delay in seconds. Algorithm based on P. Kenneth Seidelmann, "Explanatory Supplement to the Astronomical Almanac", 1992, pp. 294-296. */ double e = 0, p = 0, q = 0, delay1 = 0, delay2 = 0, delay3 = 0, delay = 0, s = 0, density = 0, bigA = 1.22E+08, r = 0, a = 0.44E+06, b = 0.44E+06, beta = 0, z = 0, mag_zenith = 0, mag_rho = 0; double[] rho = new double[4]; double[] zenith = new double[4]; int i = 0, j = 0; planetary_ephemeris(jultime,rho,0); /* The time delay due to general relativity */ e = Math.sqrt(Math.pow((site_r[1] - planet_r[0][1]),2) + Math.pow((site_r[2] - planet_r[0][2]),2) + Math.pow((site_r[3] - planet_r[0][3]),2)); p = Math.sqrt(Math.pow((asteroid_r[1] - planet_r[0][1]),2) + Math.pow((asteroid_r[2] - planet_r[0][2]),2) + Math.pow((asteroid_r[3] - planet_r[0][3]),2)); q = Math.sqrt(Math.pow((asteroid_r[1] - site_r[1]),2) + Math.pow((asteroid_r[2] - site_r[2]),2) + Math.pow((asteroid_r[3] - site_r[3]),2)); delay1 = 86400.0*(2.0*Math.pow(kappa,2)/Math.pow(speed_of_light,3))*Math.log(Math.abs((e+p+q)/(e+p-q))); /* The time delay due to the solar corona */ /* We will break this integral up into 100 slices from site to asteroid, and evaluate the density at each. */ /* slice length s is in cm */ s = (q/100.0)*au*100000.0; for (i=1;i<=100;i++) { for (j=1;j<=3;j++) rho[j] = (i/100.0)*(asteroid_r[j] - site_r[j]) + (site_r[j] - planet_r[0][j]); beta = Math.asin(rho[3]/Math.sqrt(Math.pow(rho[1],2) + Math.pow(rho[2],2) + Math.pow(rho[3],2))); r = Math.sqrt(Math.pow(rho[1],2) + Math.pow(rho[2],2) + Math.pow(rho[3],2))/planet_radius_0; density = bigA/Math.pow(r,6) + (a*b/Math.sqrt(Math.pow((a*Math.sin(beta)),2) + Math.pow((b*Math.cos(beta)),2)))/Math.pow(r,2); delay2 = delay2 + (40.3/((speed_of_light*au*100000.0/86400.0)*Math.pow(frequency,2)))*density*s; } /* The time delay due to the Earth's troposphere and ionosphere */ for (i=1;i<=3;i++) { zenith[i] = site_r[i] - planet_r[3][i]; rho[i] = asteroid_r[i] - site_r[i]; } mag_zenith = Math.sqrt(Math.pow(zenith[1],2) + Math.pow(zenith[2],2) + Math.pow(zenith[3],2)); mag_rho = Math.sqrt(Math.pow(rho[1],2) + Math.pow(rho[2],2) + Math.pow(rho[3],2)); z = Math.acos((zenith[1]*rho[1] + zenith[2]*rho[2] + zenith[3]*rho[3])/(mag_zenith*mag_rho)); delay3 = 7E-09/(Math.cos(z) + 0.0014/(0.045 + 1.0/Math.tan(z))); /* The output is the sum of the three delays */ delay = delay1 + delay2 + delay3; return delay; } void get_state_transition_matrix(double rnominal[], double rprimenominal[], double time, double phi[][]) { /* Method to calculate the state transition matrix at time, given the epoch state vector (available as an instance variable of the MinorPlanet object class) and the nominal state vector at time. */ double[] epoch_rvariant = new double[4]; double[] epoch_rprimevariant = new double[4]; double[] rvariant = new double[4]; double[] rprimevariant = new double[4]; double[][] collision_variant = new double[2000][9]; double delta = 0; int i = 0, j = 0; /* Initialize the variant state vector */ for (i=1;i<=3;i++) { epoch_rvariant[i] = epoch_r[i]; epoch_rprimevariant[i] = epoch_rprime[i]; } /* Calculate the state transition matrix phi at time */ for (i=1;i<=3;i++) { /* Vary epoch_r[i], and determine the change in state vector at time */ delta = 0.000001*epoch_r[i]; epoch_rvariant[i] = epoch_r[i] + delta; update(epoch_rvariant, epoch_rprimevariant, epoch_time, rvariant, rprimevariant, time, false, collision_variant); for (j=1;j<=3;j++) { phi[j][i] = (rvariant[j] - rnominal[j])/delta; phi[j+3][i] = (rprimevariant[j] - rprimenominal[j])/delta; } epoch_rvariant[i] = epoch_r[i]; /* Vary epoch_rprime[i], and determine the change in state vector at time */ delta = 0.000001*epoch_rprime[i]; epoch_rprimevariant[i] = epoch_rprime[i] + delta; update(epoch_rvariant, epoch_rprimevariant, epoch_time, rvariant, rprimevariant, time, false, collision_variant); for (j=1;j<=3;j++) { phi[j][i+3] = (rvariant[j] - rnominal[j])/delta; phi[j+3][i+3] = (rprimevariant[j] - rprimenominal[j])/delta; } epoch_rprimevariant[i] = epoch_rprime[i]; /* If non-gravitational thrust parameters A1_A2_DT are being utilized, vary A1_A2_DT[i], and determine the change in the state vector at time */ if ((Math.abs(A1_A2_DT[i]) > 0.0) && (i < 3)) { delta = 0.1*A1_A2_DT[i]; A1_A2_DT[i] = A1_A2_DT[i] + delta; update(epoch_rvariant, epoch_rprimevariant, epoch_time, rvariant, rprimevariant, time, false, collision_variant); for (j=1;j<=3;j++) { phi[j][i+6] = (rvariant[j] - rnominal[j])/delta; phi[j+3][i+6] = (rprimevariant[j] - rprimenominal[j])/delta; } phi[i+6][i+6] = 1.0; A1_A2_DT[i] = A1_A2_DT[i] - delta; } } } public void advance_epoch(double newtime) { /* Method to advance the epoch state vector and error covariance matrix (accessible as instance variables of the class MinorPlanet) from epoch_time to another specified date (newtime). Note that the epoch state vector and error covariance matrix are actually changed here; thus, they are not sent back as arguments. */ int i = 0, j = 0, k = 0; double[] newr = new double[4]; double[] newrprime = new double[4]; double[][] new_big_p = new double[9][9]; double[][] phi = new double[9][9]; double[][] inter = new double[9][9]; double[][] collision_variant = new double[2000][9]; /* Advance the epoch state vector */ update(epoch_r,epoch_rprime,epoch_time,newr,newrprime,newtime,false,collision_variant); /* Advance the epoch error covariance matrix */ get_state_transition_matrix(newr,newrprime,newtime,phi); for (i=1;i<=8;i++) { for (j=1;j<=8;j++) { inter[i][j] = 0; for (k=1;k<=8;k++) inter[i][j] = inter[i][j] + big_p[i][k]*phi[j][k]; } } for (i=1;i<=8;i++) { for (j=1;j<=8;j++) { new_big_p[i][j] = 0; for (k=1;k<=8;k++) new_big_p[i][j] = new_big_p[i][j] + phi[i][k]*inter[k][j]; } } /* Update the state vector and error covariance matrix */ epoch_time = newtime; for (i=1;i<=3;i++) { epoch_r[i] = newr[i]; epoch_rprime[i] = newrprime[i]; } for (i=1;i<=8;i++) { for (j=1;j<=8;j++) big_p[i][j] = new_big_p[i][j]; } } double probability_of_collision(int planet, double r_at_time[], double rprime_at_time[], double time, double covariance[][], int collision_counter) { /* Method to determine the probability of collision between an object and a planet, given the posvel of the body at the time of closest approach. The method also uses the epoch time and posvel of the body, and its epoch error covariance matrix, which are assumed to be instance variables. The method also permits use of the non-gravitational thrust parameters. Inputs: time = jultime of closest approach collision_counter = event counter for this call to update r_at_time = object barycentric equatorial position at time of closest approach rprime_at_time = object barycentric equatorial velocity at time of closest approach planet = number of the planet with which the object may collide (Sun=0,Mercury=1,...,Pluto=9,Moon=10) Outputs: probability = probability of collision covariance = propogated covariance matrix at time References: D.K. Yeomans, S.J. Ostro, and P.W. Chodas (1987). "Radar Astrometry of Near-Earth Asteroids". Astronomical Journal 94, 189-200. Donald K. Yeomans and Paul W. Chodas (1995). "Near Earth Objects: Orbit Determination and Impact Prediction". JPL Technical Report 95-0032. Howard Anton (1984). "Elementary Linear Algebra". pp. 204-205. Mary L. Boas (1983). "Mathematical Methods in the Physical Sciences". pp. 723-726. */ double[] rnominal = new double[4]; double[] rprimenominal = new double[4]; double[] rvariant = new double[4]; double[] rprimevariant = new double[4]; double[] epoch_rvariant = new double[4]; double[] epoch_rprimevariant = new double[4]; double[][] collision_variant = new double[2000][9]; double[][] phi = new double[9][9]; double[][] inter = new double[9][9]; double[][] transition_matrix = new double[4][4]; double[][] abc_covariance = new double[4][4]; double[] rho = new double[4]; double[] normal = new double[4]; double[] d = new double[4]; double[] r_project = new double[4]; double[] vector_a = new double[4]; double[] vector_b = new double[4]; double[] vector_c = new double[4]; double distance = 0, delta = 0, angle_cos = 0, angle_sin = 0, angle = 0, mag_normal = 0, mag_d = 0, mag_a = 0, probability = 0, probability1 = 0, probability2 = 0, probability3 = 0, start_time = 0; int i = 0, j = 0, k = 0; /* Note: epoch_time, epoch_r, and epoch_rprime are attributes of object (instance variables), hence available to all methods */ /* IMPORTANT: update must have a boolean that controls whether we test for a collision. Otherwise, we risk infinite regression */ /* Initialize variant arrays and nominal state vector */ for (i=1;i<=3;i++) { if (!VIflag1) { /* This is a standard search; variant trajectories should be centered on epoch state vector */ epoch_rvariant[i] = epoch_r[i]; epoch_rprimevariant[i] = epoch_rprime[i]; } else { /* This is a VI search; variant trajectories should be centered on Monte Carlo state vector input to update */ epoch_rvariant[i] = VI_trajectory_r[i]; epoch_rprimevariant[i] = VI_trajectory_rprime[i]; } rnominal[i] = r_at_time[i]; rprimenominal[i] = rprime_at_time[i]; collision_nominal_state[collision_counter][i] = rnominal[i]; collision_nominal_state[collision_counter][i+3] = rprimenominal[i]; } /* Calculate the position of the target planet at time, for use in determining delta and the target plane */ planetary_ephemeris(time,rho,0); distance = Math.sqrt(Math.pow((r_at_time[1] - planet_r[planet][1]), 2.0) + Math.pow((r_at_time[2] - planet_r[planet][2]), 2.0) + Math.pow((r_at_time[3] - planet_r[planet][3]), 2.0)); /* Calculate the state transition matrix phi at time */ for (i=1;i<=3;i++) { /* Skip if a VI search has resulted in a variant collision */ if (probability != -999.9) { /* If there is a previous collision/near miss event, we can start integrating where the variant trajectories left off. Otherwise, we will establish new variant trajectories */ /* delta = 0.000001*epoch_r[i]; */ /* delta = 0.000001*(distance/0.05); */ delta = 0.0000001; if (collision_sensitivity_time == 0.0) { if (!VIflag1) { /* This is a standard search; variant trajectories should be centered on epoch state vector */ epoch_rvariant[i] = epoch_r[i] + delta; start_time = epoch_time; } else { /* This is a VI search; variant trajectories should be centered on Monte Carlo state vector input to update */ epoch_rvariant[i] = VI_trajectory_r[i] + delta; start_time = epoch_time; } } else { for (j=1;j<=3;j++) { epoch_rvariant[j] = collision_sensitivity_r_variant[2*i-1][j]; epoch_rprimevariant[j] = collision_sensitivity_rprime_variant[2*i-1][j]; } start_time = collision_sensitivity_time; } if (VIflag2) { /* This is either a standard search, or a VI search in which the nominal trajectory is a VI; suppress collision detection for variant trajectories. */ update(epoch_rvariant, epoch_rprimevariant, start_time, rvariant, rprimevariant, time, false, collision_variant); } else { /* This is a VI search in which the nominal trajectory was a near-miss; allow collision detection of the variant trajectories. Set VIflag3 temporarily to prevent a recursive call back to probability_of_collision. */ VIflag3 = true; update(epoch_rvariant, epoch_rprimevariant, start_time, rvariant, rprimevariant, time, true, collision_variant); VIflag3 = false; /* If the variant trajectory results in a collision, overwrite the current Monte Carlo trajectory with this variant VI and return probability = -999.9 (which will result in no further processing in this method). After this method ends, will return to detect_collision which will set collision[1][1] to 1.0 and immediately end, returning finally to update. Update will see the collision[1][1] and probability = -999.9, overwrite the current Monte Carlo trajectory with this variant VI, and end, returning to the VI search algorithm. */ if (collision_variant[(int)collision_variant[0][0]][1] == 1.0) { /* Virtual Impactor */ probability = -999.9; VI_trajectory_r[i] = VI_trajectory_r[i] + delta; } } for (j=1;j<=3;j++) { phi[j][i] = (rvariant[j] - rnominal[j])/delta; phi[j+3][i] = (rprimevariant[j] - rprimenominal[j])/delta; collision_sensitivity_r_variant[2*i-1][j] = rvariant[j]; collision_sensitivity_rprime_variant[2*i-1][j] = rprimevariant[j]; collision_state_matrix[collision_counter][i][j] = (rvariant[j] - rnominal[j])/delta; collision_state_matrix[collision_counter][i][j+3] = (rprimevariant[j] - rprimenominal[j])/delta; } if (collision_sensitivity_time == 0.0) { if (!VIflag1) /* This is a standard search; variant trajectories should be centered on epoch state vector */ epoch_rvariant[i] = epoch_r[i]; else /* This is a VI search; variant trajectories should be centered on Monte Carlo state vector input to update */ epoch_rvariant[i] = VI_trajectory_r[i]; } else { /* No action required; the arrays will be switched to a new set before the next calculation */ } } /* Skip if a VI search has resulted in a variant collision */ if (probability != -999.9) { /* delta = 0.000001*epoch_rprime[i]; */ /* delta = 0.0000001*(distance/0.05); */ delta = 0.000000001; if (collision_sensitivity_time == 0.0) { if (!VIflag1) { /* This is a standard search; variant trajectories should be centered on epoch state vector */ epoch_rprimevariant[i] = epoch_rprime[i] + delta; } else { /* This is a VI search; variant trajectories should be centered on Monte Carlo state vector input to update */ epoch_rprimevariant[i] = VI_trajectory_rprime[i] + delta; } } else { for (j=1;j<=3;j++) { epoch_rvariant[j] = collision_sensitivity_r_variant[2*i][j]; epoch_rprimevariant[j] = collision_sensitivity_rprime_variant[2*i][j]; } } if (VIflag2) { /* This is either a standard search, or a VI search in which the nominal trajectory is a VI; suppress collision detection for variant trajectories. */ update(epoch_rvariant, epoch_rprimevariant, start_time, rvariant, rprimevariant, time, false, collision_variant); } else { /* This is a VI search in which the nominal trajectory was a near-miss; allow collision detection of the variant trajectories. Set VIflag3 temporarily to prevent a recursive call back to probability_of_collision. */ VIflag3 = true; update(epoch_rvariant, epoch_rprimevariant, start_time, rvariant, rprimevariant, time, true, collision_variant); VIflag3 = false; /* If the variant trajectory results in a collision, overwrite the current Monte Carlo trajectory with this variant VI and return probability = -999.9 (which will result in no further processing in this method). After this method ends, will return to detect_collision which will set collision[1][1] to 1.0 and immediately end, returning finally to update. Update will see the collision[1][1] and probability = -999.9, overwrite the current Monte Carlo trajectory with this variant VI, and end, returning to the VI search algorithm. */ if (collision_variant[(int)collision_variant[0][0]][1] == 1.0) { /* Virtual Impactor */ probability = -999.9; VI_trajectory_rprime[i] = VI_trajectory_rprime[i] + delta; } } for (j=1;j<=3;j++) { phi[j][i+3] = (rvariant[j] - rnominal[j])/delta; phi[j+3][i+3] = (rprimevariant[j] - rprimenominal[j])/delta; collision_sensitivity_r_variant[2*i][j] = rvariant[j]; collision_sensitivity_rprime_variant[2*i][j] = rprimevariant[j]; collision_state_matrix[collision_counter][i+3][j] = (rvariant[j] - rnominal[j])/delta; collision_state_matrix[collision_counter][i+3][j+3] = (rprimevariant[j] - rprimenominal[j])/delta; } if (collision_sensitivity_time == 0.0) { if (!VIflag1) /* This is a standard search; variant trajectories should be centered on epoch state vector */ epoch_rprimevariant[i] = epoch_rprime[i]; else /* This is a VI search; variant trajectories should be centered on Monte Carlo state vector input to update */ epoch_rprimevariant[i] = VI_trajectory_rprime[i]; } else { /* No action required; the arrays will be switched to a new set before the next calculation */ } } /* Skip if a VI search has resulted in a variant collision */ if (probability != -999.9) { /* If non-gravitational thrust parameters A1_A2_DT are being utilized, vary A1_A2_DT[i], and determine the change in the state vector at time */ if ((Math.abs(A1_A2_DT[i]) > 0.0) && (i < 3)) { delta = 0.1*A1_A2_DT[i]; A1_A2_DT[i] = A1_A2_DT[i] + delta; if (collision_sensitivity_time > 0.0) { for (j=1;j<=3;j++) { epoch_rvariant[j] = collision_sensitivity_r_variant[i+6][j]; epoch_rprimevariant[j] = collision_sensitivity_rprime_variant[i+6][j]; } } if (VIflag2) { /* This is either a standard search, or a VI search in which the nominal trajectory is a VI; suppress collision detection for variant trajectories. */ update(epoch_rvariant, epoch_rprimevariant, start_time, rvariant, rprimevariant, time, false, collision_variant); } else { /* This is a VI search in which the nominal trajectory was a near-miss; allow collision detection of the variant trajectories. Set VIflag3 temporarily to prevent a recursive call back to probability_of_collision. */ VIflag3 = true; update(epoch_rvariant, epoch_rprimevariant, start_time, rvariant, rprimevariant, time, true, collision_variant); VIflag3 = false; } for (j=1;j<=3;j++) { phi[j][i+6] = (rvariant[j] - rnominal[j])/delta; phi[j+3][i+6] = (rprimevariant[j] - rprimenominal[j])/delta; collision_sensitivity_r_variant[i+6][j] = rvariant[j]; collision_sensitivity_rprime_variant[i+6][j] = rprimevariant[j]; } phi[i+6][i+6] = 1.0; A1_A2_DT[i] = A1_A2_DT[i] - delta; } } } /* Skip if a VI search has resulted in a variant collision */ if (probability != -999.9) { /* Set the collision_sensitivity_time to the end time of integration, so we can pick up there if there is another event. */ collision_sensitivity_time = time; /* Determine error covariance at time using epoch error covariance matrix big_p */ /* Note that big_p is an attribute of object (instance variable), hence available to all methods */ for (i=1;i<=8;i++) { for (j=1;j<=8;j++) { inter[i][j] = 0; for (k=1;k<=8;k++) inter[i][j] = inter[i][j] + big_p[i][k]*phi[j][k]; } } for (i=1;i<=6;i++) { for (j=1;j<=6;j++) { covariance[i][j] = 0; for (k=1;k<=8;k++) covariance[i][j] = covariance[i][j] + phi[i][k]*inter[k][j]; } } /* Define the "impact plane", i.e., the plane perpendicular to the planet's relative velocity vector which contains the planet's position at time. */ for (i=1;i<=3;i++) normal[i] = rprime_at_time[i] - planet_rprime[planet][i]; mag_normal = Math.sqrt(Math.pow(normal[1],2) + Math.pow(normal[2],2) + Math.pow(normal[3],2)); for (i=1;i<=3;i++) { normal[i] = normal[i]/mag_normal; } /* Find the vector d from the planet to the object, and the angle between d and the normal to the impact plane */ for (i=1;i<=3;i++) d[i] = r_at_time[i] - planet_r[planet][i]; mag_d = Math.sqrt(Math.pow(d[1],2) + Math.pow(d[2],2) + Math.pow(d[3],2)); for (i=1;i<=3;i++) d[i] = d[i]/mag_d; angle_cos = d[1]*normal[1] + d[2]*normal[2] + d[3]*normal[3]; angle_sin = Math.sqrt(Math.pow((d[2]*normal[3] - d[3]*normal[2]),2) + Math.pow((d[1]*normal[3] - d[3]*normal[1]),2) + Math.pow((d[1]*normal[2] - d[2]*normal[1]),2)); angle = Math.atan2(angle_sin,angle_cos); /* Project the object position at time onto the impact plane */ for (i=1;i<=3;i++) { r_project[i] = r_at_time[i] - mag_d*angle_cos*normal[i]; } /* Create a new "a-b-c" coordinate system. Basis vector a will be in the impact plane, with its direction defined by the vector from the projected object position to the planet. Basis vector c will be the normal to the impact plane. Basis vector b will be in the impact plane, defined so as to complete the right-handed system. */ for (i=1;i<=3;i++) { vector_a[i] = planet_r[planet][i] - r_project[i]; vector_c[i] = normal[i]; } mag_a = Math.sqrt(Math.pow(vector_a[1],2) + Math.pow(vector_a[2],2) + Math.pow(vector_a[3],2)); for (i=1;i<=3;i++) vector_a[i] = vector_a[i]/mag_a; vector_b[1] = -(vector_a[2]*vector_c[3] - vector_a[3]*vector_c[2]); vector_b[2] = (vector_a[1]*vector_c[3] - vector_a[3]*vector_c[1]); vector_b[3] = -(vector_a[1]*vector_c[2] - vector_a[2]*vector_c[1]); /* Create the transition matrix, and transform the (position) error covariance matrix at time from "i-j-k" to "a-b-c". Note that we will only transform the upper 3x3 of the covariance matrix; this is all we need for our purposes. */ for (i=1;i<=3;i++) { transition_matrix[1][i] = vector_a[i]; transition_matrix[2][i] = vector_b[i]; transition_matrix[3][i] = vector_c[i]; } for (i=1;i<=3;i++) { for (j=1;j<=3;j++) { inter[i][j] = 0; for (k=1;k<=3;k++) inter[i][j] = inter[i][j] + covariance[i][k]*transition_matrix[j][k]; } } for (i=1;i<=3;i++) { for (j=1;j<=3;j++) { abc_covariance[i][j] = 0; for (k=1;k<=3;k++) abc_covariance[i][j] = abc_covariance[i][j] + transition_matrix[i][k]*inter[k][j]; } } /* Use the transformed covariances to determine the probability that the planet and the object physically overlap on the impact plane. Note that this requires integrating the normal density function over the planet's disk. However, we are making a pretty huge simplification already by projecting onto the "impact plane", in order to predict the probability of collision before or after the time of nominal closest approach. Therefore, I feel justified in approximating this integral. The integral would be over all a and b that cover the planet's disk. I will approximate this numerically three different ways. First, by evaluating the ndf at the planet's center and multiplying by the area of the planet's disk. Second, by evaluating the ndf at four points along the a and b axes, and multiplying by area/4. Third, by evaluating the ndf at four points along the a-b diagonals, and multiplying by area/4. And finally, I will average these approximations to arrive at a final result. You have to stop somewhere, and I feel that this covers the disk fairly well. */ probability1 = Math.pow(planet_radius[planet],2)*Math.exp(-Math.pow(mag_a,2)/(2*Math.abs(abc_covariance[1][1])))/(2*Math.sqrt(Math.abs(abc_covariance[1][1]))*Math.sqrt(Math.abs(abc_covariance[2][2]))); probability2 = Math.pow(planet_radius[planet],2)*(2*Math.exp(-Math.pow(mag_a,2)/(2*Math.abs(abc_covariance[1][1])))*Math.exp(-Math.pow(planet_radius[planet],2)/(4*Math.abs(abc_covariance[2][2]))) + Math.exp(-Math.pow((mag_a - planet_radius[planet]/Math.sqrt(2.0)),2)/(2*Math.abs(abc_covariance[1][1]))) + Math.exp(-Math.pow((mag_a + planet_radius[planet]/Math.sqrt(2.0)),2)/(2*Math.abs(abc_covariance[1][1]))))/(8*Math.sqrt(Math.abs(abc_covariance[1][1]))*Math.sqrt(Math.abs(abc_covariance[2][2]))); probability3 = Math.pow(planet_radius[planet],2)*(2*Math.exp(-Math.pow((mag_a - planet_radius[planet]/Math.sqrt(2.0)),2)/(2*Math.abs(abc_covariance[1][1])))*Math.exp(-Math.pow((planet_radius[planet]/Math.sqrt(2.0)),2)/(2*Math.abs(abc_covariance[2][2]))) + 2*Math.exp(-Math.pow((mag_a + planet_radius[planet]/Math.sqrt(2.0)),2)/(2*Math.abs(abc_covariance[1][1])))*Math.exp(-Math.pow((planet_radius[planet]/Math.sqrt(2.0)),2)/(2*Math.abs(abc_covariance[2][2]))))/(8*Math.sqrt(Math.abs(abc_covariance[1][1]))*Math.sqrt(Math.abs(abc_covariance[2][2]))); probability = (probability1 + probability2 + probability3)/3.0; } return probability; } void detect_collision(double object_r_at_n[],double object_rprime_at_n[],double object_r_at_n1[],double object_rprime_at_n1[],double jultime_n,double jultime_n1,double collision[][]) { /* Given the position and velocity of an object at jultime_n and jultime_n1, determine whether the object collides with a planet, or has a near miss with a planet, between time steps n and n+1 in procedure update. Inputs: object_r_at_n[], object_rprime_at_n[] = the barycentric equatorial posvel of the object at time step n; object_r_at_n1[], object_rprime_at_n1[] = the barycentric equatorial posvel of the object at time step n+1; jultime_n, jultime_n1 = the Julian dates corresponding to time steps n and n+1 The output array "collision[][]" has two indices. The first index accounts for multiple encounters, the total being stored in "collision[0][0]". For the i'th encounter, entry "collision[i][1]" indicates whether a collision ( = 1) or a near miss ( = 2) occurs. In either case, "collision[i][0]" gives the minimum distance in A.U.s. "collision[i][2]" gives the planet (1-9, Moon=10, Sun=0), "collision[i][3]" gives the Julian date, "collision[i][4]" gives the nominal distance in A.U.s., "collision[i][5]" gives the estimated probability of collision, "collision[i][6]" gives the object's x-covariance at the time of the event, "collision[i][7]" gives the object's y-covariance at the time of the event, and "collision[i][8]" gives the object's z-covariance at the time of the event. */ double[] rho = new double[4]; double[][] planet_r_at_n = new double[11][4]; double[][] planet_r_at_n1 = new double[11][4]; double[][] planet_rprime_at_n = new double[11][4]; double[][] planet_rprime_at_n1 = new double[11][4]; double[] planet_i_r_at_n = new double[4]; double[] planet_i_rprime_at_n = new double[4]; double[] planet_i_r_at_time = new double[4]; double[] planet_i_rprime_at_time = new double[4]; double[] r_at_time = new double[4]; double[] rprime_at_time = new double[4]; double[] r_at_n = new double[4]; double[] rprime_at_n = new double[4]; double[] r_at_n1 = new double[4]; double[] rprime_at_n1 = new double[4]; double[][] covariance = new double[7][7]; double[][] collision_variant = new double[2000][9]; double distance_n1 = 0, upper_time = 0, lower_time = 0, time = 0, distance = 0, distance_n = 0, deriv_1_at_n = 0, deriv_1_at_n1 = 0, deriv_1_at_time = 0; int i = 0, j = 0, collision_counter = 0, counter = 0; collision_counter = (int)(collision[0][0]); /* Get the ephemerides of the planets at each time step, without correction for light time-of-flight */ planetary_ephemeris(jultime_n,rho,0); /* Initialize barycentric planetary posvel arrays at time step n */ for (i=0;i<=10;i++) { for (j=1;j<=3;j++) { planet_r_at_n[i][j] = planet_r[i][j]; planet_rprime_at_n[i][j] = planet_rprime[i][j]; } } /* Initialize barycentric object posvel arrays at time step n */ for (j=1;j<=3;j++) { r_at_n[j] = object_r_at_n[j]; rprime_at_n[j] = object_rprime_at_n[j]; } planetary_ephemeris(jultime_n1,rho,0); /* Initialize barycentric planetary posvel arrays at time step n+1 */ for (i=0;i<=10;i++) { for (j=1;j<=3;j++) { planet_r_at_n1[i][j] = planet_r[i][j]; planet_rprime_at_n1[i][j] = planet_rprime[i][j]; } } /* Initialize barycentric object posvel arrays at time step n+1 */ for (j=1;j<=3;j++) { r_at_n1[j] = object_r_at_n1[j]; rprime_at_n1[j] = object_rprime_at_n1[j]; } /* Check each planet (and the Sun) for collision */ for (i=0;i<=10;i++) { if ((!earth_only) || ((earth_only) && (i == 3))) { /* Initialize barycentric posvel arrays for planet i */ for (j=1;j<=3;j++) { planet_i_r_at_n[j] = planet_r_at_n[i][j]; planet_i_rprime_at_n[j] = planet_rprime_at_n[i][j]; } /* Check for collision at time step n+1. */ distance_n1 = Math.sqrt(Math.pow((object_r_at_n1[1] - planet_r_at_n1[i][1]),2) + Math.pow((object_r_at_n1[2] - planet_r_at_n1[i][2]),2) + Math.pow((object_r_at_n1[3] - planet_r_at_n1[i][3]),2)); if (distance_n1 <= planet_radius[i]) { /* Collision. Increment collision counter and use bisection to find it. */ collision_counter = collision_counter + 1; collision[0][0] = collision[0][0] + 1; upper_time = jultime_n1; lower_time = jultime_n; collision[collision_counter][0] = 0; collision[collision_counter][1] = 1; collision[collision_counter][2] = i; collision[collision_counter][3] = jultime_n1; collision[collision_counter][4] = distance_n1; while ((Math.abs(upper_time - lower_time) > event_localization) && (counter < 100)) { counter = counter + 1; time = lower_time + (upper_time - lower_time)/2; /* Set VIflag5 to prevent this bisection call to update from resetting the VI_trajectory_r and _prime vectors */ VIflag5 = true; update(r_at_n,rprime_at_n,jultime_n,r_at_time,rprime_at_time,time,false,collision_variant); VIflag5 = false; planetary_ephemeris(time,rho,0); for (j=1;j<=3;j++) { planet_i_r_at_time[j] = planet_r[i][j]; planet_i_rprime_at_time[j] = planet_rprime[i][j]; } distance = Math.sqrt(Math.pow((r_at_time[1] - planet_i_r_at_time[1]),2) + Math.pow((r_at_time[2] - planet_i_r_at_time[2]),2) + Math.pow((r_at_time[3] - planet_i_r_at_time[3]),2)); if (distance <= planet_radius[i]) { collision[collision_counter][3] = time; collision[collision_counter][4] = distance; upper_time = time; } else /* Collision has not yet occurred. Increase lower limit. */ lower_time = time; } System.out.println("COLLISION with planet " + i + " at " + time + " distance = " + distance); if (!VIflag3 && !VIflag4) { /* At this point, we know that this is either a standard search, or it is a VI search and the (current) nominal trajectory is a VI; temporarily set VIflag2 to prevent probability_of_collision from attempting collision detection of the variant trajectories. */ VIflag2 = true; collision[collision_counter][5] = probability_of_collision(i,r_at_time,rprime_at_time,time,covariance,collision_counter); VIflag2 = false; collision[collision_counter][6] = Math.sqrt(Math.abs(covariance[1][1])); collision[collision_counter][7] = Math.sqrt(Math.abs(covariance[2][2])); collision[collision_counter][8] = Math.sqrt(Math.abs(covariance[3][3])); if (!VIflag1) collision[collision_counter][0] = 3.0*covariance_second_eigenvalue(covariance); } else { /* If VIflag3, this is a VI search where nominal trajectory is a near-miss and variant trajectories are now being tested for collision; avoid a recursive call to probability_of_collision. */ collision[collision_counter][5] = 0; collision[collision_counter][6] = 0; collision[collision_counter][7] = 0; collision[collision_counter][8] = 0; } if (VIflag4) { /* If VIflag4, we want collision detection without probability_of_collision and its six variant trajectories. We also want to archive the state vector at time. */ for (j=1;j<=3;j++) { collision_nominal_state[collision_counter][j] = r_at_time[j]; collision_nominal_state[collision_counter][j+3] = rprime_at_time[j]; } collision[collision_counter][5] = 0; collision[collision_counter][6] = 0; collision[collision_counter][7] = 0; collision[collision_counter][8] = 0; } } else { /* Not in collision at time step n+1; evaluate intermediate points */ /* Calculate the first derivative of distance at time steps n and n+1; check for a minimum that is less than the planetary radius for collision, or less than near_miss_threshold for a near miss. */ distance_n = Math.sqrt(Math.pow((object_r_at_n[1] - planet_r_at_n[i][1]),2) + Math.pow((object_r_at_n[2] - planet_r_at_n[i][2]),2) + Math.pow((object_r_at_n[3] - planet_r_at_n[i][3]),2)); deriv_1_at_n = ((object_r_at_n[1] - planet_r_at_n[i][1])*(object_rprime_at_n[1] - planet_rprime_at_n[i][1]) + (object_r_at_n[2] - planet_r_at_n[i][2])*(object_rprime_at_n[2] - planet_rprime_at_n[i][2]) + (object_r_at_n[3] - planet_r_at_n[i][3])*(object_rprime_at_n[3] - planet_rprime_at_n[i][3]))/distance_n; deriv_1_at_n1 = ((object_r_at_n1[1] - planet_r_at_n1[i][1])*(object_rprime_at_n1[1] - planet_rprime_at_n1[i][1]) + (object_r_at_n1[2] - planet_r_at_n1[i][2])*(object_rprime_at_n1[2] - planet_rprime_at_n1[i][2]) + (object_r_at_n1[3] - planet_r_at_n1[i][3])*(object_rprime_at_n1[3] - planet_rprime_at_n1[i][3]))/distance_n1; if ((deriv_1_at_n < 0) && (deriv_1_at_n1 > 0) && (distance_n < (1.2*near_miss_threshold))) { /* A relevant minimum has occurred. Use bisection on the first derivative of distance to find it. */ upper_time = jultime_n1; lower_time = jultime_n; while ((Math.abs(upper_time - lower_time) > event_localization) && (counter < 100)) { counter = counter + 1; time = lower_time + (upper_time - lower_time)/2; /* Set VIflag5 to prevent this bisection call to update from resetting the VI_trajectory_r and _prime vectors */ VIflag5 = true; update(r_at_n,rprime_at_n,jultime_n,r_at_time,rprime_at_time,time,false,collision_variant); VIflag5 = false; planetary_ephemeris(time,rho,0); for (j=1;j<=3;j++) { planet_i_r_at_time[j] = planet_r[i][j]; planet_i_rprime_at_time[j] = planet_rprime[i][j]; } distance = Math.sqrt(Math.pow((r_at_time[1] - planet_i_r_at_time[1]),2) + Math.pow((r_at_time[2] - planet_i_r_at_time[2]),2) + Math.pow((r_at_time[3] - planet_i_r_at_time[3]),2)); deriv_1_at_time = ((r_at_time[1] - planet_i_r_at_time[1])*(rprime_at_time[1] - planet_i_rprime_at_time[1]) + (r_at_time[2] - planet_i_r_at_time[2])*(rprime_at_time[2] - planet_i_rprime_at_time[2]) + (r_at_time[3] - planet_i_r_at_time[3])*(rprime_at_time[3] - planet_i_rprime_at_time[3]))/distance; if (deriv_1_at_time < 0) /* Minimum has not yet occurred. Increase lower limit. */ lower_time = time; else /* Minimum has already occurred. Reduce upper limit. */ upper_time = time; } /* We've now located the minimum to great accuracy. Is it a collision, or a near miss? */ if (distance <= planet_radius[i]) { /* Collision. Increment collision counter and use bisection to find it. */ collision_counter = collision_counter + 1; collision[0][0] = collision[0][0] + 1; upper_time = time; lower_time = jultime_n; collision[collision_counter][0] = 0; collision[collision_counter][1] = 1; collision[collision_counter][2] = i; collision[collision_counter][3] = time; collision[collision_counter][4] = distance; while ((Math.abs(upper_time - lower_time) > event_localization) && (counter < 200)) { counter = counter + 1; time = lower_time + (upper_time - lower_time)/2; /* Set VIflag5 to prevent this bisection call to update from resetting the VI_trajectory_r and _prime vectors */ VIflag5 = true; update(r_at_n,rprime_at_n,jultime_n,r_at_time,rprime_at_time,time,false,collision_variant); VIflag5 = false; planetary_ephemeris(time,rho,0); for (j=1;j<=3;j++) { planet_i_r_at_time[j] = planet_r[i][j]; planet_i_rprime_at_time[j] = planet_rprime[i][j]; } distance = Math.sqrt(Math.pow((r_at_time[1] - planet_i_r_at_time[1]),2) + Math.pow((r_at_time[2] - planet_i_r_at_time[2]),2) + Math.pow((r_at_time[3] - planet_i_r_at_time[3]),2)); if (distance <= planet_radius[i]) { collision[collision_counter][3] = time; collision[collision_counter][4] = distance; upper_time = time; } else /* Collision has not yet occurred. Increase lower limit. */ lower_time = time; } System.out.println("COLLISION with planet " + i + " at " + time + " distance = " + distance); if (!VIflag3 && !VIflag4) { /* At this point, we know that this is either a standard search, or it is a VI search and the (current) nominal trajectory is a VI; temporarily set VIflag2 to prevent probability_of_collision from attempting collision detection of the variant trajectories. */ VIflag2 = true; collision[collision_counter][5] = probability_of_collision(i,r_at_time,rprime_at_time,time,covariance,collision_counter); VIflag2 = false; collision[collision_counter][6] = Math.sqrt(Math.abs(covariance[1][1])); collision[collision_counter][7] = Math.sqrt(Math.abs(covariance[2][2])); collision[collision_counter][8] = Math.sqrt(Math.abs(covariance[3][3])); if (!VIflag1) collision[collision_counter][0] = 3.0*covariance_second_eigenvalue(covariance); } else { /* If VIflag3, this is a VI search where nominal trajectory is a near-miss and variant trajectories are now being tested for collision; avoid a recursive call to probability_of_collision. */ collision[collision_counter][5] = 0; collision[collision_counter][6] = 0; collision[collision_counter][7] = 0; collision[collision_counter][8] = 0; } if (VIflag4) { /* If VIflag4, we want collision detection without probability_of_collision and its six variant trajectories. We also want to archive the state vector at time. */ for (j=1;j<=3;j++) { collision_nominal_state[collision_counter][j] = r_at_time[j]; collision_nominal_state[collision_counter][j+3] = rprime_at_time[j]; } collision[collision_counter][5] = 0; collision[collision_counter][6] = 0; collision[collision_counter][7] = 0; collision[collision_counter][8] = 0; } } else if (Math.abs(distance - planet_radius[i]) < near_miss_threshold) { /* Near miss. Increment collision counter and record data. */ collision_counter = collision_counter + 1; collision[0][0] = collision[0][0] + 1; collision[collision_counter][1] = 2; collision[collision_counter][2] = i; collision[collision_counter][3] = time; collision[collision_counter][4] = distance; if (!VIflag3 && !VIflag4) { /* We know that this is a near-miss. If it is a standard search, temporarily set VIflag2 to prevent probability_of_collision from attempting collision detection of the variant trajectories. My initial feeling was that, if it is a VI search, do not set VIflag2; probability_of_collision will attempt collision detection of the variant trajectories. Upon further consideration, however, it is unlikely that the variants created by probability_of_collision would lie near the LOV, or have comparably minimized residuals. Additionally, collision detection takes time. I'll therefore set VIflag2 in all cases. */ /* if (!VIflag1) VIflag2 = true; */ VIflag2 = true; collision[collision_counter][5] = probability_of_collision(i,r_at_time,rprime_at_time,time,covariance,collision_counter); VIflag2 = false; collision[collision_counter][6] = Math.sqrt(Math.abs(covariance[1][1])); collision[collision_counter][7] = Math.sqrt(Math.abs(covariance[2][2])); collision[collision_counter][8] = Math.sqrt(Math.abs(covariance[3][3])); if (!VIflag1) collision[collision_counter][0] = 3.0*covariance_second_eigenvalue(covariance); if (collision[collision_counter][0] < 0.0) collision[collision_counter][0] = 0.0; if (collision[collision_counter][5] == -999.9) { /* A variant trajectory has resulted in an impact. Reset the number of events and the probability so update will note the event and overwrite the Monte Carlo trajectory with the variant VI. */ collision[0][0] = 1; collision[1][1] = 1; collision[1][5] = -999.9; } } else { /* If VIflag3, this is a VI search where nominal trajectory is a near-miss and variant trajectories are now being tested for collision; avoid a recursive call to probability_of_collision. */ collision[collision_counter][5] = 0; collision[collision_counter][6] = 0; collision[collision_counter][7] = 0; collision[collision_counter][8] = 0; } if (VIflag4) { /* If VIflag4, we want collision detection without probability_of_collision and its six variant trajectories. We also want to archive the state vector at time. */ for (j=1;j<=3;j++) { collision_nominal_state[collision_counter][j] = r_at_time[j]; collision_nominal_state[collision_counter][j+3] = rprime_at_time[j]; } collision[collision_counter][5] = 0; collision[collision_counter][6] = 0; collision[collision_counter][7] = 0; collision[collision_counter][8] = 0; } } } } } } } public void get_mp_state_vector(double classical_elements[], double time1, double r[], double rprime[]) { /* Method to convert the classical heliocentric ecliptic elements of a minor planet at TDB time1 into a barycentric equatorial state vector at that same epoch. Inputs: classical elements[1] - a in AU [2] - e [3] - i in radians [4] - w in radians [5] - omega in radians [6] - T expressed as a TDB julian day [7] - mean anomaly in radians [8] - q in AU time1 - TDB epoch to which the elements refer, expressed as a julday Outputs: r, rprime - barycentric equatorial position/velocity of the body at epoch Tested and verified 10 August 1999. */ double[] sun_r = new double[4]; double[] sun_rprime = new double[4]; int i = 0; /* First, convert classical heliocentric ecliptic elements at TDB time1 to pos/vel */ convert_Keplerian_elements_to_state_vector(classical_elements,r,rprime,time1); /* Two-body routines function in modified time; convert velocity to AU/day */ for (i=1;i<=3;i++) rprime[i] = 0.01720209895*rprime[i]; /* Rotate asteroid heliocentric state vector from ecliptic to equatorial */ ecliptic_to_equatorial(r); ecliptic_to_equatorial(rprime); /* Get Sun's position at TDB time1 */ get_planet_posvel(time1,11,sun_r,sun_rprime); /* Add Sun's posvel to convert asteroid state vector from helio- to barycentric */ for (i=1;i<=3;i++) { r[i] = r[i] + sun_r[i]; rprime[i] = rprime[i] + sun_rprime[i]; } } void get_observed_radec(double state_r[], double state_rprime[], double time1, boolean geocentric, double latitude, double longitude, double altitude, double ut1_minus_tai, double observer_r[], double observer_rprime[], double time2, double radec[]) { /* Method that uses the barycentric equatorial state vector of a body at TDB time1 to calculate its observed radec from a site at TDT time2. Inputs: state_r and state_rprime in AU and AU/day time1 is a TDB julian day geocentric is a boolean indicating whether the radec should be geocentric or relative to latitude/longitude/altitude latitude and longitude are in radians, with longitude measured eastward altitude is in meters ut1_minus_tai is in seconds observer_r and observer_rprime are arrays containing the barycentric equatorial rectangular position and velocity of the observing site, in AU and AU/day time2 is a TDT julian day Output: radec is an array containing ra and dec, in radians Tested and verified 10 August 1999. */ /* This method conforms to the outline of "Explanatory Supplement to the Astronomical Almanac", 1992, section 3.31 (Apparent-Place Algorithm for Planets). Note: time1 is the TDB epoch of the asteroid's orbital elements time2 is the TDT time of observation on Earth */ double dist = 0, time3 = 0, ra = 0, dec = 0, lat = 0, longit = 0; double rho_dot = 0, r_mag = 0, r_dot = 0, oldtime3 = 10; double dist_u = 0, dist_e = 0, dist_q = 0, tau = 0; double[] EOP = new double[8]; double[] r = new double[4]; double[] rprime = new double[4]; double[] r2 = new double[4]; double[] r2prime = new double[4]; double[] r3 = new double[4]; double[] r3prime = new double[4]; double[] rho = new double[4]; double[] latlong = new double[3]; double[] vector_u = new double[4]; double[] vector_q = new double[4]; double[] vector_e = new double[4]; double[][] collision = new double[2000][9]; double[] sun_r = new double[4]; double[] sun_rprime = new double[4]; int i = 0, j = 0, counter = 0; /* Section 3.311 */ /* Convert TDT time2 to its TDB equivalent; if (UT1-TAI) is zero (unknown), tdt_to_tdb algorithm will estimate it */ time2 = tdt_to_tdb(time2,latitude,longitude,altitude,ut1_minus_tai); /* Section 3.312 */ if (geocentric == true) { /* Geocentric case - get Earth barycentric equatorial position and velocity at TDB time2 */ planetary_ephemeris(time2,rho,0); for (i=1;i<=3;i++) { observer_r[i] = planet_r[3][i]; observer_rprime[i] = planet_rprime[3][i]; } } else { /* Topocentric case */ /* The observer coordinates determined when the observation was entered may have been passed as an argument; these may have used precise values of the EOPs, so we don't want to recalculate. On the other hand, the observer coordinates may be zero; this method may be creating an ephemeris, so precise values of the EOPs may not be available. We'll test and act accordingly. */ if (observer_r[1] == 0.0) { /* Observer coordinates are unavailable; calculate them. */ /* Set EOPs to zero so that method get_observer_position will use estimates and defaults. However, allow for the use of ut1_minus_tai if it is provided. */ for (i=1;i<=7;i++) EOP[i] = 0.0; EOP[3] = ut1_minus_tai; /* get observer barycentric equatorial position and velocity at TDB time2 */ get_observer_position(latitude,longitude,altitude,time2,EOP,observer_r,observer_rprime); } /* If observer_r[1] is not equal to 0.0, observer coordinates are available; use them. */ } /* Create and normalize Earth heliocentric position vector at TDB time2 */ get_planet_posvel(time2,11,sun_r,sun_rprime); for (i=1;i<=3;i++) vector_e[i] = observer_r[i] - sun_r[i]; dist_e = Math.sqrt(Math.pow(vector_e[1],2) + Math.pow(vector_e[2],2) + Math.pow(vector_e[3],2)); for (i=1;i<=3;i++) vector_e[i] = vector_e[i]/dist_e; /* Section 3.313 */ /* Get asteroid barycentric equatorial state vector at time of observation TDB time2 */ update(state_r,state_rprime,time1,r3,r3prime,time2,false,collision); /* Sections 3.314 and 3.315 */ /* Iterative correction for light time-of-flight */ for (i=1;i<=3;i++) { /* r[] and rprime[] will represent state vector at time3, the time at which light left asteroid to be observed on Earth at time2 */ r[i] = r3[i]; rprime[i] = r3prime[i]; } /* Create and normalize the line-of-sight vector from observer at time2 to asteroid at time3 */ for (i=1;i<=3;i++) vector_u[i] = r[i] - observer_r[i]; dist_u = Math.sqrt(Math.pow(vector_u[1],2) + Math.pow(vector_u[2],2) + Math.pow(vector_u[3],2)); /* First approximation: time3 = time2 - dist_u/speed_of_light */ time3 = time2 - dist_u/speed_of_light; /* Get asteroid barycentric equatorial state vector at TDB time3 */ update(r3,r3prime,time2,r,rprime,time3,false,collision); while ((Math.abs(time3 - oldtime3) > 1e-6) && (counter < 100)) { oldtime3 = time3; counter = counter + 1; /* Create and normalize the line-of-sight vector from observer at time2 to asteroid at time3 */ for (i=1;i<=3;i++) vector_u[i] = r[i] - observer_r[i]; dist_u = Math.sqrt(Math.pow(vector_u[1],2) + Math.pow(vector_u[2],2) + Math.pow(vector_u[3],2)); for (i=1;i<=3;i++) vector_u[i] = vector_u[i]/dist_u; /* Create and normalize asteroid's heliocentric position vector at TDB time3 */ get_planet_posvel(time3,11,sun_r,sun_rprime); for (i=1;i<=3;i++) vector_q[i] = r[i] - sun_r[i]; dist_q = Math.sqrt(Math.pow(vector_q[1],2) + Math.pow(vector_q[2],2) + Math.pow(vector_q[3],2)); for (i=1;i<=3;i++) vector_q[i] = vector_q[i]/dist_q; /* Recalculate time3 based on updated estimate of light time-of-flight */ tau = dist_u/speed_of_light + (2*Math.pow(kappa,2)/Math.pow(speed_of_light,3))*Math.log((dist_e + dist_u + dist_q)/(dist_e - dist_u + dist_q)); time3 = time2 - tau; /* Find asteroid position at time3 */ update(r3,r3prime,time2,r,rprime,time3,false,collision); } /* Section 3.3110 */ /* Calculate ra/dec */ radec[2] = Math.atan2(vector_u[3],Math.sqrt(Math.pow(vector_u[1],2) + Math.pow(vector_u[2],2))); radec[1] = Math.atan2(vector_u[2],vector_u[1]); } public void update(double r1[],double r1prime[],double time1,double r2[],double r2prime[],double time2,boolean check_for_collision,double collision[][]) { /* Procedure to calculate the position and velocity of a body at time2 given the position and velocity of the body at time1. The Cowell method is used, with integration performed using Runge-Kutta Prince-Dormand 8(7). Local error per unit step is controlled by variable step size. Positions and velocities are given in A.U.s and A.U.s/day, while the times are Julian dates. The procedure also checks for collisions with planets, if the Boolean check_for_collision is set. The output array "collision[][]" has two indices. The first index accounts for multiple encounters from time1 to time2, the total being stored in "collision[0][0]". For the i'th encounter, entry "collision[i][1]" indicates whether a collision ( = 1) or a near miss ( = 2) occurs. In either case, "collision[i][2]" gives the planet (1-9, Moon=10, Sun=0), "collision[i][3]" gives the Julian date, "collision[i][4]" gives the nominal distance in A.U.s., "collision[i][5]" gives the estimated probability of collision, "collision[i][6]" gives the object's x-covariance at the time of the event, "collision[i][7]" gives the object's y-covariance at the time of the event, and "collision[i][8]" gives the object's z-covariance at the time of the event. "collision[i][0]" gives the semi-width of the 3-sigma uncertainty region in A.U.s. Reference for the Cowell method is Bate, Mueller, and White, "Fundamentals of Astrodynamics", pp 386-390. Reference for RK dopri 8(7) is P.J. Prince and J.R. Dormand, "High-order embedded Runge-Kutta formulae", J. Comp. Appl. Math., vol 7, 1981, p. 67-75. The error control is adapted from the class notes of Math 551 taught at The Pennsylvania State University during Fall 1989 by Prof. Douglas Arnold. */ double[] fatn = new double[4]; double[] r = new double[4]; double[] rprime = new double[4]; double[] r_at_n = new double[4]; double[] rprime_at_n = new double[4]; double[] fatnp1 = new double[4]; double[] rho = new double[4]; double toler = 0, h = 0, n1time = 0, ntime = 0, localerror = 0, ratio = 0, test_distance = 0, min_distance = 0; int i = 0, j = 0, k = 0, istep = 0, collision_counter = 0, counter = 0, counter_limit = 0, candidate = 0; /* Initialize output collision array */ for (i=0;i<=1999;i++) { for (j=0;j<=8;j++) { collision[i][j] = 0; } } if (time1 == time2) { /* Return the state vector unchanged */ for (k=1;k<=3;k++) { r2[k] = r1[k]; r2prime[k] = r1prime[k]; } } else { /* Set local error per unit step tolerance */ /* For use in the asteroid_ephemeris procedure, we will use a fixed tolerance designed to limit integration error over 200 years */ /* toler = integration_error/Math.abs(time2 - time1); */ toler = integration_error/73050.0; /* Determine the maximum number of attempted steps to be permitted */ /* counter_limit = (int)(Math.abs(time2 - time1)/0.01); */ counter_limit = 100; /* Initialize variables for first time step */ /* istep is the step being calculated */ istep = 1; /* Set initial step size */ h = time2 - time1; /* Limit initial h to +/- 20 days */ if (Math.abs(h) > 20.0) h = 20.0*h/Math.abs(h); /* n1time is the time at the (n+1)'st step, ntime is the time at the n'th step */ n1time = time1; ntime = time1; /* Initialize position, velocity, and acceleration vectors */ for (k=1;k<=3;k++) { r[k] = r1[k]; r_at_n[k] = r1[k]; rprime[k] = r1prime[k]; rprime_at_n[k] = r1prime[k]; /* During VI search, a great many Monte Carlo trajectories are created in a normal distribution about the nominal epoch state vector; differential correction then attempts to force a VI. In calculating the differential partial derivatives, probability_of_collision must create variant trajectories centered on the Monte Carlo trajectory, rather than the nominal epoch state vector. Thus, we want to save the Monte Carlo trajectory for use by probability_of_collision. */ if ((VIflag1) && (!VIflag2) && (!VIflag3) && (!VIflag5)) { VI_trajectory_r[k] = r1[k]; VI_trajectory_rprime[k] = r1prime[k]; } } get_acceleration(ntime,r_at_n,rprime_at_n,fatn); /* Begin integration loop */ integrationloop: while ((((n1time < time2) && (time2 > time1)) || ((n1time > time2) && (time2 < time1))) && (collision[collision_counter][1] != 1.0) && (counter < counter_limit)) { /* So long as we've not reached time2 or slammed into a planet, ... */ /* Increment time to the (n+1)'st step */ n1time = ntime + h; counter = counter + 1; /* if (Math.abs(h) < 0.01) System.out.println("step size = " + h); */ /* System.out.println("step size = " + h); */ /* Integrate to get r and rprime at the (n+1)'st time step */ /* Use Runge-Kutta Dormand-Prince 8(7) */ localerror = dopri8(ntime,h,fatn,rprime_at_n,r_at_n,rprime,r); /* If a variant trajectory collides with a planet, the error may go to infinite or NaN. If so, end the integration. */ if (!(Math.abs(localerror) < 1e+10)) counter = counter_limit; /* Compute ratio of projected new step size to present step size */ if (Math.abs(localerror/h) > 0.0) /* Avoid a division by zero for small step sizes */ ratio = .8*Math.pow(Math.abs(toler/(localerror/h)),(1.0/7.0)); else /* Error and step size are very small; allow step size to double */ ratio = 2; /* Test to see if step will be accepted; if already at minimum step size, ignore */ if ((Math.abs(localerror/h) > toler) && (Math.abs(h) > 0.001)) { /* Step has failed */ /* Determine a new step size */ if ((ratio < .2) && (istep > 1)) ratio = .2; /* Note that new step size is limited to one-fifth of the previous one, except on the first step, which allows faster determination of an initial step size */ n1time = ntime; h = ratio*h; /* Enforce minimum step size to prevent ultrasmall steps near collision */ if (Math.abs(h) < 0.001) h = 0.001*(h/Math.abs(h)); continue integrationloop; } /* Step has been accepted, or we're at the minimum step size. In either case, reset the counter. */ counter = 0; /* Check for collision, if so desired */ /* Note: This check for collision must occur after a step has been accepted; otherwise, the check will be made more than once for failed steps. */ if (check_for_collision == true) { /* TESTING PURPOSES ONLY!! */ /* System.out.println("check for collision jultime = " + ntime); */ /* TESTING PURPOSES ONLY!! */ detect_collision(r_at_n,rprime_at_n,r,rprime,ntime,n1time,collision); collision_counter = (int)(collision[0][0]); if (collision[1][5] == -999.9) { /* During a VI search, a variant trajectory has resulted in an impact. Overwrite the Monte Carlo trajectory with the variant VI. Since detect_collision has set collision[1][1] = 1.0, this method will end, and we will return to the VI search algorithm. */ for (k=1;k<=3;k++) { r1[k] = VI_trajectory_r[k]; r1prime[k] = VI_trajectory_rprime[k]; } } } if ((check_for_collision == false) || ((check_for_collision == true) && (collision[collision_counter][1] != 1.0))) { /* Determine a new step size */ if (ratio > 2) ratio = 2; /* Note that new step size is limited to double the previous one */ h = ratio*h; /* Enforce minimum step size to prevent ultrasmall steps near collision */ if (Math.abs(h) < 0.001) h = 0.001*(h/Math.abs(h)); /* Calculate second derivative of r at time n+1 */ get_acceleration(n1time,r,rprime,fatnp1); /* Increment step counter */ istep = istep + 1; /* Check to prevent overrunning endpoint of integration */ if (((((n1time + h) > time2) && (time2 > time1)) || (((n1time + h) < time2) && (time2 < time1))) && (n1time != time2)) /* Reduce step size to avoid overrunning endpoint of integration */ h = time2 - n1time; /* If a variant trajectory collides with a planet, the acceleration may go to infinite or NaN. If so, end the integration. */ if (!(Math.abs(fatnp1[1]) < 1e+10) || !(Math.abs(fatnp1[2]) < 1e+10) || !(Math.abs(fatnp1[3]) < 1e+10)) counter = counter_limit; /* Update integration variables */ ntime = n1time; for (k=1;k<=3;k++) { r_at_n[k] = r[k]; rprime_at_n[k] = rprime[k]; fatn[k] = fatnp1[k]; } } } if (!(counter < counter_limit)) { System.out.println("WARNING: update exceeds max steps"); /* Proceeding on the assumption that the excessively high local error was caused by an imminent collision, make a guess as to which planet was responsible, and populate the collision array accordingly */ planetary_ephemeris(n1time,rho,0); min_distance = 1; for (k=0;k<=10;k++) { test_distance = Math.sqrt(Math.pow((r[1] - planet_r[k][1]),2) + Math.pow((r[2] - planet_r[k][2]),2) + Math.pow((r[3] - planet_r[k][3]),2)); if (test_distance < min_distance) { min_distance = test_distance; candidate = k; } } if (min_distance < 0.01) { collision[0][0] = collision[0][0] + 1; collision[(int)collision[0][0]][0] = 0; collision[(int)collision[0][0]][1] = 1; collision[(int)collision[0][0]][2] = candidate; collision[(int)collision[0][0]][3] = n1time; collision[(int)collision[0][0]][4] = min_distance; collision[(int)collision[0][0]][5] = 0; collision[(int)collision[0][0]][6] = 0; collision[(int)collision[0][0]][7] = 0; collision[(int)collision[0][0]][8] = 0; } } /* Integration complete; return new state vector, and reset collision_sensitivity state vectors, if necessary */ for (k=1;k<=3;k++) { r2[k] = r[k]; r2prime[k] = rprime[k]; if ((check_for_collision == true) && (!VIflag3)) { /* Don't reset the variant trajectories if a variant is being checked for collision during VI search */ for (j=1;j<=8;j++) { collision_sensitivity_r_variant[j][k] = 0; collision_sensitivity_rprime_variant[j][k] = 0; } collision_sensitivity_time = 0; } } } } void initialize_dopri8() { /* Initialization of the integration coefficients for dopri8 */ /* Initialize arrays */ dopri8_a[2][1] = 1.0/18.0; dopri8_a[3][1] = 1.0/48.0; dopri8_a[3][2] = 1.0/16.0; dopri8_a[4][1] = 1.0/32.0; dopri8_a[4][2] = 0; dopri8_a[4][3] = 3.0/32.0; dopri8_a[5][1] = 5.0/16.0; dopri8_a[5][2] = 0; dopri8_a[5][3] = -75.0/64.0; dopri8_a[5][4] = 75.0/64.0; dopri8_a[6][1] = 3.0/80.0; dopri8_a[6][2] = 0; dopri8_a[6][3] = 0; dopri8_a[6][4] = 3.0/16.0; dopri8_a[6][5] = 3.0/20.0; dopri8_a[7][1] = 29443.841/614563.906; dopri8_a[7][2] = 0; dopri8_a[7][3] = 0; dopri8_a[7][4] = 77736.538/692538.347; dopri8_a[7][5] = -28693.883/1125000.000; dopri8_a[7][6] = 23124.283/1800000.000; dopri8_a[8][1] = 16016.141/946692.911; dopri8_a[8][2] = 0; dopri8_a[8][3] = 0; dopri8_a[8][4] = 61564.180/158732.637; dopri8_a[8][5] = 22789.713/633445.777; dopri8_a[8][6] = 545815.736/2771057.229; dopri8_a[8][7] = -180193.667/1043307.555; dopri8_a[9][1] = 39632.708/573591.083; dopri8_a[9][2] = 0; dopri8_a[9][3] = 0; dopri8_a[9][4] = -433636.366/683701.615; dopri8_a[9][5] = -421739.975/2616292.301; dopri8_a[9][6] = 100302.831/723423.059; dopri8_a[9][7] = 790204.164/839813.087; dopri8_a[9][8] = 800635.310/3783071.287; dopri8_a[10][1] = 246121.993/1340847.787; dopri8_a[10][2] = 0; dopri8_a[10][3] = 0; dopri8_a[10][4] = -37695042.795/15268766.246; dopri8_a[10][5] = -309121.744/1061227.803; dopri8_a[10][6] = -12992.083/490766.935; dopri8_a[10][7] = 6005943.493/2108947.869; dopri8_a[10][8] = 393006.217/1396673.457; dopri8_a[10][9] = 123872.331/1001029.789; dopri8_a[11][1] = -1028468.189/846180.014; dopri8_a[11][2] = 0; dopri8_a[11][3] = 0; dopri8_a[11][4] = 8478235.783/508512.852; dopri8_a[11][5] = 1311729.495/1432422.823; dopri8_a[11][6] = -10304129.995/1701304.382; dopri8_a[11][7] = -48777925.059/3047939.560; dopri8_a[11][8] = 15336726.248/1032824.649; dopri8_a[11][9] = -45442868.181/3398467.696; dopri8_a[11][10] = 3065993.473/597172.653; dopri8_a[12][1] = 185892.177/718116.043; dopri8_a[12][2] = 0; dopri8_a[12][3] = 0; dopri8_a[12][4] = -3185094.517/667107.341; dopri8_a[12][5] = -477755.414/1098053.517; dopri8_a[12][6] = -703635.378/230739.211; dopri8_a[12][7] = 5731566.787/1027545.527; dopri8_a[12][8] = 5232866.602/850066.563; dopri8_a[12][9] = -4093664.535/808688.257; dopri8_a[12][10] = 3962137.247/1805957.418; dopri8_a[12][11] = 65686.358/487910.083; dopri8_a[13][1] = 403863.854/491063.109; dopri8_a[13][2] = 0; dopri8_a[13][3] = 0; dopri8_a[13][4] = -5068492.393/434740.067; dopri8_a[13][5] = -411421.997/543043.805; dopri8_a[13][6] = 652783.627/914296.604; dopri8_a[13][7] = 11173962.825/925320.556; dopri8_a[13][8] = -13158990.841/6184727.034; dopri8_a[13][9] = 3936647.629/1978049.680; dopri8_a[13][10] = -160528.059/685178.525; dopri8_a[13][11] = 248638.103/1413531.060; dopri8_a[13][12] = 0; dopri8_c[1] = 0; dopri8_c[2] = 1.0/18.0; dopri8_c[3] = 1.0/12.0; dopri8_c[4] = 1.0/8.0; dopri8_c[5] = 5.0/16.0; dopri8_c[6] = 3.0/8.0; dopri8_c[7] = 59.0/400.0; dopri8_c[8] = 93.0/200.0; dopri8_c[9] = 5490023.248/9719169.821; dopri8_c[10] = 13.0/20.0; dopri8_c[11] = 1201146.811/1299019.798; dopri8_c[12] = 1; dopri8_c[13] = 1; dopri8_b[1] = 13451.932/455176.623; dopri8_b[2] = 0; dopri8_b[3] = 0; dopri8_b[4] = 0; dopri8_b[5] = 0; dopri8_b[6] = -808719.846/976000.145; dopri8_b[7] = 1757004.468/5645159.321; dopri8_b[8] = 656045.339/265891.186; dopri8_b[9] = -3867574.721/1518517.206; dopri8_b[10] = 465885.868/322736.535; dopri8_b[11] = 53011.238/667516.719; dopri8_b[12] = 2.0/45.0; dopri8_b[13] = 0; dopri8_bhat[1] = 14005.451/335480.064; dopri8_bhat[2] = 0; dopri8_bhat[3] = 0; dopri8_bhat[4] = 0; dopri8_bhat[5] = 0; dopri8_bhat[6] = -59238.493/1068277.825; dopri8_bhat[7] = 181606.767/758867.731; dopri8_bhat[8] = 561292.985/797845.732; dopri8_bhat[9] = -1041891.430/1371343.529; dopri8_bhat[10] = 760417.239/1151165.299; dopri8_bhat[11] = 118820.643/751138.087; dopri8_bhat[12] = -528747.749/2220607.170; dopri8_bhat[13] = 0.25; } double dopri8(double ntime,double h,double fatn[],double rprimen[],double rn[],double rprime[],double r[]) { /* Procedure to integrate the equations of motion (Cowell's method) from ntime to ntime+h using the Runge-Kutta-type Prince-Dormand 8(7) method (see procedure update). Also finds local error in the integration step, which is returned. Reference: P.J. Prince and J.R. Dormand, "High-order embedded Runge-Kutta formulae", J. Comp. Appl. Math., vol 7, 1981, p. 67-75. Also, E. Hairer, S.P. Norsett, G. Wanner, "Solving Ordinary Differential Equations I (Non-stiff problems), Springer-Verlag, 1987, pp. 132-133, 167-171, 193-195, and 260. */ double[][] dopri8_kprime = new double[14][4]; double[][] dopri8_k = new double[14][4]; double[] low_r = new double[4]; double[] rvariant = new double[4]; double[] rprimevariant = new double[4]; double[] kprime_temp = new double[4]; double localerror = 0; int i = 0, j = 0, n = 0; /* If necessary, initialize the integration coefficients */ if (dopri8_a[2][1] != 1.0/18.0) initialize_dopri8(); /* Begin the integration from ntime to ntime+h */ /* Stage 1 */ for (n=1;n<=3;n++) { /* Initialize the Runge-Kutta variables, and the pos/vel vectors. */ dopri8_kprime[1][n] = fatn[n]; dopri8_k[1][n] = rprimen[n]; /* r is the 8'th order position solution; low_r is the 7'th order solution */ r[n] = rn[n] + h*dopri8_bhat[1]*dopri8_k[1][n]; low_r[n] = rn[n] + h*dopri8_b[1]*dopri8_k[1][n]; rprime[n] = rprimen[n] + h*dopri8_bhat[1]*dopri8_kprime[1][n]; } /* Stages 2-13 */ for (i=2;i<=13;i++) { for (n=1;n<=3;n++) { /* rvariant and rprimevariant are the intermediate pos/vel used to evaluate the acceleration function at each stage */ rvariant[n] = rn[n]; rprimevariant[n] = rprimen[n]; dopri8_k[i][n] = rprimen[n]; for (j=1;j<=(i-1);j++) { /* Sum the terms in k, rvariant, and rprimevariant needed to evaluate the acceleration function and determine kprime at this stage */ rvariant[n] = rvariant[n] + h*(dopri8_a[i][j]*dopri8_k[j][n]); rprimevariant[n] = rprimevariant[n] + h*(dopri8_a[i][j]*dopri8_kprime[j][n]); dopri8_k[i][n] = dopri8_k[i][n] + h*(dopri8_a[i][j]*dopri8_kprime[j][n]); } } get_acceleration(ntime+dopri8_c[i]*h,rvariant,rprimevariant,kprime_temp); for (n=1;n<=3;n++) { dopri8_kprime[i][n] = kprime_temp[n]; r[n] = r[n] + h*dopri8_bhat[i]*dopri8_k[i][n]; low_r[n] = low_r[n] + h*dopri8_b[i]*dopri8_k[i][n]; rprime[n] = rprime[n] + h*dopri8_bhat[i]*dopri8_kprime[i][n]; } } /* Compute the local error in position */ localerror = 0; for (n=1;n<=3;n++) localerror = localerror + Math.pow((r[n] - low_r[n]),2); localerror = Math.sqrt(localerror); return localerror; } void get_acceleration(double jultime,double r[],double rprime[],double acceleration[]) { /* Procedure to calculate the accelerations on a body for use in the differential equations of motion in Cowell's method. First, we calculate the gravitational perturbations due to the planets, using the barycentric relativistic equation for point-mass acceleration (Ref: P. Kenneth Seidelmann, "Explanatory Supplement to the Astronomical Almanac", 1992, p. 281), and then the perturbing acceleration on a comet due to reactive water vaporization (Refs: Brian G. Marsden, "Comets and Non-Gravitational Forces II", The Astronomical Journal Vol 74 No. 5, pp. 720-734; B.G. Marsden, Z. Sekanina, D.K. Yeomans, "Comets and Non-Gravitational Forces V", The Astronomical Journal, Vol 78 No. 2, pp. 211-225; D.K. Yeomans and P.W. Chodas, "An Asymmetric Outgassing Model for Cometary Nongravitational Accelerations", The Astronomical Journal, Vol 98 No. 3, pp. 1083-1093). Note that jultime is a Julian date, r[] is the barycentric position vector in A.U.s of the body in question, rprime[] is the barycentric velocity in AU/day, and acceleration() is the acceleration vector on the body given in AU/day^3. Note also that planet_r(I,j) represent the j'th component of the position of the I'th solar system body, planet_rprime(I,j) represent the j'th component of the velocity of t