Classes | |
class | WindowParameters |
class | OrbitStream |
class | Analysis |
class | Lyapunov |
class | MeanMotionResonance |
class | PoincareSurfaceOfSection |
class | BodyConstants |
class | Body |
class | BodyWithParameter |
base element for intepolation More... | |
class | BodyWithEpoch |
class | JPLBody |
class | ConfigItem |
class | Config |
class | Vector |
class | VectorWithParameter |
class | Debug |
struct | gsl_d1_minimization_parameters |
class | FFTDataStructure |
class | FFTDataStream |
class | FFTPowerSpectrumBaseElement |
class | FFTPowerSpectrum |
class | Peak |
class | FFT |
FFT analysis. More... | |
class | File |
File base class. More... | |
class | ReadFile |
Read-only files class. More... | |
class | WriteFile |
Write-only files class. More... | |
class | ReadWriteFile |
Read and write files class. More... | |
class | Mercury5IntegrationFile |
Mercury 5 integration input files. More... | |
class | RadauModIntegrationFile |
Modified Radau input files. More... | |
class | SWIFTFile |
SWIFT integration file. More... | |
class | LocationFile |
Locations of the observatories. More... | |
class | MPCObsFile |
MPC observation file. More... | |
class | RWOFile |
AstDyS observation file, usually with a .rwo extension. More... | |
class | AsteroidDatabaseFile |
Read-only asteroid files class. More... | |
class | AstDySMatrixFile |
NEODyS and AstDyS .ctc and .ctm files. More... | |
class | NEODYSCAT |
NEODyS and AstDyS .cat file. More... | |
class | JPLDastcomNumFile |
class | JPLDastcomUnnumFile |
class | JPLDastcomCometFile |
class | AstorbFile |
Lowell asteroids database file. More... | |
class | MPCOrbFile |
MPC asteroids database file. More... | |
class | MPCCometFile |
MPC comets database file. More... | |
class | OrsaFile |
orsa default input-output file More... | |
class | OrsaConfigFile |
orsa configuration file More... | |
class | OrsaPaths |
class | TLEFile |
class | JPLFile |
JPL ephem file. More... | |
class | JPLCache |
class | Frame |
class | Integrator |
This is the interface for all the Integrator classes. More... | |
class | FixedTimestepIntegrator |
class | MultistepIntegrator |
class | VariableTimestepIntegrator |
class | ModifiedMidpoint |
Advances using a number of substeps (midpoints). For conservative and non-conservative Interaction. More... | |
class | Stoer |
Like the ModifiedMidpoint, but for conservative Interaction. More... | |
class | RungeKutta |
class | DissipativeRungeKutta |
class | Radau15 |
class | Leapfrog |
class | Interaction |
class | MappedTable |
class | Legendre |
class | Newton |
class | GravitationalTree |
class | Relativistic |
class | ArmonicOscillator |
class | GalacticPotentialAllen |
class | GalacticPotentialAllenPlusNewton |
class | JPLPlanetsNewton |
class | Observation |
class | Sky |
class | Orbit |
class | OrbitWithEpoch |
class | OptimizedOrbitPositions |
class | Asteroid |
class | AsteroidDatabase |
class | ObservationCandidate |
class | Location |
class | OrbitWithCovarianceMatrixGSL |
class | PreliminaryOrbit |
class | CloseApproach |
class | UnitBaseScale |
class | Units |
class | TimeStep |
class | Date |
class | UniverseTypeAwareTime |
class | UniverseTypeAwareTimeStep |
class | Angle |
class | Evolution |
This class collects all the Frames of an integration, sampled at a fixed sample_period. More... | |
class | Universe |
Typedefs | |
typedef orsa::gsl_d1_minimization_parameters | gsl_d1_minimization_parameters |
Enumerations | |
enum | ConfigEnum { JPL_EPHEM_FILE, JPL_DASTCOM_NUM, JPL_DASTCOM_UNNUM, JPL_DASTCOM_COMET, LOWELL_ASTORB, MPC_MPCORB, MPC_COMET, MPC_NEA, MPC_DAILY, MPC_DISTANT, MPC_PHA, MPC_UNUSUALS, ASTDYS_ALLNUM_CAT, ASTDYS_ALLNUM_CTC, ASTDYS_ALLNUM_CTM, ASTDYS_UFITOBS_CAT, ASTDYS_UFITOBS_CTC, ASTDYS_UFITOBS_CTM, NEODYS_CAT, NEODYS_CTC, OBSCODE, TLE_NASA, TLE_GEO, TLE_GPS, TLE_ISS, TLE_KEPELE, TLE_VISUAL, TLE_WEATHER, TEXTURE_SUN, TEXTURE_MERCURY, TEXTURE_VENUS, TEXTURE_EARTH, TEXTURE_MOON, TEXTURE_MARS, TEXTURE_JUPITER, TEXTURE_SATURN, TEXTURE_URANUS, TEXTURE_NEPTUNE, TEXTURE_PLUTO, NO_CONFIG_ENUM } |
enum | FFTSearch { HK = 0, D, PQ, NODE, ANOMALY, ANOMALY_PHASE, A_M, TESTING } |
signal types More... | |
enum | FFTSearchAmplitude { A, E, I, SIN_I, TAN_I_2, ONE } |
enum | FFTSearchPhase { OMEGA_NODE, OMEGA_PERICENTER, OMEGA_TILDE, MM, LAMBDA, ZERO } |
enum | FFTAlgorithm { algo_FFT, algo_FFTB, algo_MFT, algo_FMFT1, algo_FMFT2 } |
enum | FILE_STATUS { CLOSE, OPEN_R, OPEN_W } |
enum | M5COLS { C7, C10 } |
enum | OrsaFileDataType { OFDT_END_OF_FILE = 0, OFDT_UNIVERSE = 1, OFDT_EVOLUTION = 2, OFDT_FRAME = 3, OFDT_BODY = 4 } |
enum | JPL_planets { NONE = 0, MERCURY = 1, VENUS = 2, EARTH = 3, MARS = 4, JUPITER = 5, SATURN = 6, URANUS = 7, NEPTUNE = 8, PLUTO = 9, MOON = 10, SUN = 11, SOLAR_SYSTEM_BARYCENTER = 12, EARTH_MOON_BARYCENTER = 13, NUTATIONS = 14, LIBRATIONS = 15, EARTH_AND_MOON = 1000 } |
enum | IntegratorType { STOER = 1, BULIRSCHSTOER = 2, RUNGEKUTTA = 3, DISSIPATIVERUNGEKUTTA = 4, RA15 = 5, LEAPFROG = 6 } |
enum | InteractionType { NEWTON = 1, ARMONICOSCILLATOR = 2, GALACTIC_POTENTIAL_ALLEN = 3, GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON = 4, JPL_PLANETS_NEWTON = 5, GRAVITATIONALTREE = 6, NEWTON_MPI = 7, RELATIVISTIC = 8 } |
enum | CovarianceMatrixElements { Osculating, Equinoctal } |
enum | time_unit { YEAR = 1, DAY = 2, HOUR = 3, MINUTE = 4, SECOND = 5 } |
enum | length_unit { MPARSEC = 1, KPARSEC = 2, PARSEC = 3, LY = 4, AU = 5, EARTHMOON = 6, REARTH = 7, RMOON = 8, KM = 9, M = 10, CM = 11, LD = EARTHMOON, ER = REARTH, MR = RMOON } |
enum | mass_unit { MSUN = 1, MJUPITER = 2, MEARTH = 3, MMOON = 4, KG = 5, GRAM = 6 } |
enum | TimeScale { UTC = 1, UT = 2, TAI = 3, TDT = 4, GPS = 5, UT1 = UT, ET = TDT, TT = TDT } |
TimeScale enum, useful only when using a Real Universe. More information can be obtained here: http://www.hartrao.ac.za/nccsdoc/slalib/sun67.htx/node217.html. More... | |
enum | ReferenceSystem { EQUATORIAL = 1, ECLIPTIC = 2 } |
enum | UniverseType { Real = 1, Simulated = 2 } |
This enum is used to classify the Universes: a Real Universe is composed. More... | |
Functions | |
static double | local_mass (const JPL_planets planet) |
static double | local_J2 (const JPL_planets planet) |
static double | local_J3 (const JPL_planets planet) |
static double | local_J4 (const JPL_planets planet) |
static double | local_C22 (const JPL_planets planet) |
static double | local_C31 (const JPL_planets planet) |
static double | local_C32 (const JPL_planets planet) |
static double | local_C33 (const JPL_planets planet) |
static double | local_C41 (const JPL_planets planet) |
static double | local_C42 (const JPL_planets planet) |
static double | local_C43 (const JPL_planets planet) |
static double | local_C44 (const JPL_planets planet) |
static double | local_S31 (const JPL_planets planet) |
static double | local_S32 (const JPL_planets planet) |
static double | local_S33 (const JPL_planets planet) |
static double | local_S41 (const JPL_planets planet) |
static double | local_S42 (const JPL_planets planet) |
static double | local_S43 (const JPL_planets planet) |
static double | local_S44 (const JPL_planets planet) |
void | Interpolate (const vector< BodyWithParameter > &b_in, const double x, Body &b_out, Body &err_b_out) |
void | print (const Body &b) |
bool | operator== (const Body &b1, const Body &b2) |
bool | operator!= (const Body &b1, const Body &b2) |
void | Interpolate (const std::vector< BodyWithParameter > &b_in, const double x, Body &b_out, Body &err_b_out) |
double | KineticEnergy (const Body &) |
string | Label (const ConfigEnum e) |
void | Interpolate (const vector< VectorWithParameter > vx_in, const double x, Vector &v_out, Vector &err_v_out) |
Vector | operator * (const double f, const Vector &v) |
Vector | operator * (const Vector &v, const double f) |
Vector | operator/ (const Vector &v, const double f) |
Vector | operator+ (const Vector &u, const Vector &v) |
Vector | operator- (const Vector &u, const Vector &v) |
Vector | ExternalProduct (const Vector &u, const Vector &v) |
Vector | Cross (const Vector &u, const Vector &v) |
double | operator * (const Vector &u, const Vector &v) |
bool | operator== (const Vector &v1, const Vector &v2) |
bool | operator!= (const Vector &v1, const Vector &v2) |
void | Interpolate (const std::vector< VectorWithParameter > v_in, const double x, Vector &v_out, Vector &err_v_out) |
double | norm (const fftw_complex z) |
double | norm_sq (const fftw_complex z) |
fftw_complex | phi (double omega, fftw_complex in[], const int size) |
Discrete Fourier Transform. | |
fftw_complex | phi_Hanning (double omega, fftw_complex in[], const int size) |
Discrete Fourier Transform with Hanning windowing. | |
double | phi_amp (double omega, fftw_complex in[], const int size) |
amplitude for spectrum, without windowing | |
double | phi_Hanning_amp (double omega, fftw_complex in[], const int size) |
amplitude for spectrum, with Hanning windowing | |
double | phi_gsl (double x, void *params) |
double | phi_gsl_two (double x, void *params) |
double | phi_Hanning_gsl (double x, void *params) |
int | compare_binamp (const binamp *a, const binamp *b) |
sort binamp struct from the bigger to the smaller... | |
double | psd_max_again (const fftw_complex *transformed_signal, const int size) |
void | psd_max_again_many (const fftw_complex *transformed_signal, const int size, vector< double > &candidate, const unsigned int nfreq) |
double | psd_max (const fftw_complex *transformed_signal, const int size) |
void | apply_window (fftw_complex *signal_win, fftw_complex *signal, int size) |
void | amph (double *amp, double *phase, fftw_complex *phiR, fftw_complex *phiL, double freq, fftw_complex *in, int size) |
double | accurate_peak (double left, double center, double right, fftw_complex *in, int size) |
double | dQ (double y) |
int | SWIFTRawReadBinaryFile (FILE_TYPE file, const int version=2) |
static void | swap_4 (void *ptr) |
static void | swap_8 (void *ptr) |
static void | swap (void *ptr, unsigned int size) |
void | convert (OrsaFileDataType &ofdt, const unsigned int i) |
void | remove_leading_trailing_spaces (std::string &s) |
void | SetupSolarSystem (Frame &frame, const list< JPL_planets > &l, const UniverseTypeAwareTime &t) |
The Sun is automatically included. | |
string | JPL_planet_name (const JPL_planets p) |
double | radius (const JPL_planets p) |
void | convert (JPL_planets &jp, const unsigned int i) |
void | SetupSolarSystem (Frame &, const std::list< JPL_planets > &, const UniverseTypeAwareTime &) |
void | print (const Frame &f) |
static double | modified_mu (const vector< Body > &f, const unsigned int i) |
Vector | CenterOfMass (const vector< Body > &f) |
Vector | CenterOfMassVelocity (const vector< Body > &f) |
Vector | Barycenter (const vector< Body > &f) |
Vector | BarycenterVelocity (const vector< Body > &f) |
Vector | RelativisticBarycenter (const vector< Body > &f) |
Vector | RelativisticBarycenterVelocity (const vector< Body > &f) |
Vector | AngularMomentum (const vector< Body > &f, const Vector ¢er) |
Vector | BarycentricAngularMomentum (const vector< Body > &f) |
double | KineticEnergy (const vector< Body > &f) |
double | KineticEnergy (const std::vector< Body > &) |
Vector | CenterOfMass (const std::vector< Body > &) |
Vector | CenterOfMassVelocity (const std::vector< Body > &) |
Vector | Barycenter (const std::vector< Body > &) |
Vector | BarycenterVelocity (const std::vector< Body > &) |
Vector | RelativisticBarycenter (const std::vector< Body > &) |
Vector | RelativisticBarycenterVelocity (const std::vector< Body > &) |
Vector | AngularMomentum (const std::vector< Body > &, const Vector &) |
Vector | BarycentricAngularMomentum (const std::vector< Body > &) |
void | make_new_integrator (Integrator **i, const IntegratorType type) |
void | convert (IntegratorType &it, const unsigned int i) |
std::string | label (const IntegratorType it) |
std::string | label (const InteractionType it) |
void | make_new_interaction (Interaction **i, const InteractionType type) |
void | convert (InteractionType &it, const unsigned int i) |
double | delta_function (const unsigned int i, const unsigned int j) |
Vector | ComputeAcceleration (const list< Body >::const_iterator body_it, const list< TreeNode >::const_iterator node_domain_it, const bool compute_quadrupole=true) |
double | RMS_residuals (const vector< Observation > &obs, const OrbitWithEpoch &orbit) |
The RMS of the residuals in arcsec units. | |
double | residual (const Observation &, const OrbitWithEpoch &) |
The RMS of the residuals in arcsec units. | |
Sky | PropagatedSky_J2000 (const OrbitWithEpoch &orbit, const UniverseTypeAwareTime &final_time, const std::string &obscode, const bool integrate, const bool light_time_corrections) |
static Vector | unit_vector (const Angle &ra, const Angle &dec) |
double | RMS_residuals (const std::vector< Observation > &, const OrbitWithEpoch &) |
double | Weight (const Observation &obs1, const Observation &obs2, const double optimal_dt=FromUnits(1, DAY)) |
double | Weight (vector< Observation > &obs, const double optimal_dt=FromUnits(1, DAY)) |
void | BuildObservationTriplets (const vector< Observation > &obs, vector< ThreeObservations > &triplets, const double optimal_dt=FromUnits(1, DAY)) |
double | MOID_f (const gsl_vector *v, void *params) |
double | MOID (const Orbit &, const Orbit &, Vector &, Vector &) |
Minimal Orbit Intersection Distance between two orbits. | |
double | MOID2RB_f (const gsl_vector *v, void *params) |
double | MOID2RB (const Vector &, const Vector &, const Orbit &, const Orbit &, Vector &, Vector &) |
Minimal Orbit Intersection Distance between two orbits with different reference bodies. | |
double | E1 (void *p) |
double | M1 (void *p1, void *p2) |
void | S1 (const gsl_rng *r, void *p, double step_size) |
int | least_sq_f (const gsl_vector *v, void *params, gsl_vector *f) |
double | least_sq_diff_f (double x, void *params) |
int | least_sq_df (const gsl_vector *v, void *params, gsl_matrix *J) |
int | least_sq_fdf (const gsl_vector *v, void *params, gsl_vector *f, gsl_matrix *J) |
void | OrbitDifferentialCorrectionsLeastSquares (OrbitWithCovarianceMatrixGSL &orbit, const vector< Observation > &obs) |
int | gauss_v_f (const gsl_vector *v, void *params, gsl_vector *f) |
double | gauss_v_diff_f (double x, void *params) |
int | gauss_v_df (const gsl_vector *v, void *params, gsl_matrix *J) |
int | gauss_v_fdf (const gsl_vector *v, void *params, gsl_vector *f, gsl_matrix *J) |
void | gauss_v (const Vector &r, Vector &v, const vector< Observation > &obs) |
OrbitWithCovarianceMatrixGSL | Compute (const vector< Observation > &obs) |
double | poly_8 (double x, void *params) |
double | poly_8_df (double x, void *params) |
void | poly_8_fdf (double x, void *params, double *y, double *dy) |
void | poly_8_gsl_solve (poly_8_params ¶ms, vector< poly_8_solution > &solutions) |
void | Compute_Gauss (const vector< Observation > &obs_in, vector< PreliminaryOrbit > &preliminary_orbits) |
void | Compute_TestMethod (const vector< Observation > &obs_in, vector< PreliminaryOrbit > &preliminary_orbits) |
double | CloseApproaches_gsl_f (const gsl_vector *v, void *params) |
void | SearchCloseApproaches (const Evolution *evol, const unsigned int obj_index, const unsigned int index, std::vector< CloseApproach > &clapp, const double distance_threshold, const double time_accuracy) |
OrbitWithCovarianceMatrixGSL | Compute (const std::vector< Observation > &) |
General interface to the Orbit computation from a set of observations. | |
void | Compute_Gauss (const std::vector< Observation > &, std::vector< PreliminaryOrbit > &) |
void | Compute_TestMethod (const std::vector< Observation > &, std::vector< PreliminaryOrbit > &) |
void | OrbitDifferentialCorrectionsLeastSquares (OrbitWithCovarianceMatrixGSL &, const std::vector< Observation > &) |
double | secure_pow (double x, double y) |
double | secure_log (double x) |
double | secure_log10 (double x) |
double | secure_atan2 (double x, double y) |
double | secure_asin (double x) |
double | secure_acos (double x) |
double | secure_sqrt (double x) |
bool | operator== (const TAI_minus_UTC_element &x, const TAI_minus_UTC_element &y) |
bool | operator!= (const TAI_minus_UTC_element &x, const TAI_minus_UTC_element &y) |
bool | operator== (const ET_minus_UT_element &x, const ET_minus_UT_element &y) |
bool | operator!= (const ET_minus_UT_element &x, const ET_minus_UT_element &y) |
string | TimeScaleLabel (TimeScale ts) |
UniverseTypeAwareTimeStep | operator * (const double x, const UniverseTypeAwareTimeStep &ts) |
UniverseTypeAwareTimeStep | operator * (const UniverseTypeAwareTimeStep &ts, const double x) |
UniverseTypeAwareTimeStep | operator * (const int i, const UniverseTypeAwareTimeStep &ts) |
UniverseTypeAwareTimeStep | operator * (const UniverseTypeAwareTimeStep &ts, const int i) |
Angle | obleq (const Date &date) |
Angle | gmst (const Date &date) |
Angle | obleq_J2000 () |
void | EclipticToEquatorial (Vector &v, const Date &date) |
void | EquatorialToEcliptic (Vector &v, const Date &date) |
void | EclipticToEquatorial_J2000 (Vector &v) |
void | EquatorialToEcliptic_J2000 (Vector &v) |
static void | alpha_delta_meridian_moon (const Date &date, Angle &alpha_zero, Angle &delta_zero, Angle &W) |
void | alpha_delta_meridian (const JPL_planets p, const Date &date, Angle &alpha_zero, Angle &delta_zero, Angle &W) |
void | convert (time_unit &tu, const unsigned int i) |
void | convert (length_unit &lu, const unsigned int i) |
void | convert (mass_unit &mu, const unsigned int i) |
double | GetG () |
double | GetG_MKS () |
double | GetMSun () |
double | GetC () |
double | FromUnits (const double, const time_unit, const int=1) |
double | FromUnits (const double, const length_unit, const int=1) |
double | FromUnits (const double, const mass_unit, const int=1) |
std::string | TimeLabel () |
std::string | LengthLabel () |
std::string | MassLabel () |
void | convert (TimeScale &ts, const unsigned int i) |
double | sin (const Angle &alpha) |
double | cos (const Angle &alpha) |
double | tan (const Angle &alpha) |
void | sincos (const Angle &alpha, double &s, double &c) |
void | convert (ReferenceSystem &rs, const unsigned int i) |
Frame | StartFrame (const vector< BodyWithEpoch > &f, vector< JPL_planets > &jf, const Interaction *itr, const Integrator *itg, const UniverseTypeAwareTime &t) |
void | convert (UniverseType &ut, const unsigned int i) |
Frame | StartFrame (const std::vector< BodyWithEpoch > &, std::vector< JPL_planets > &, const Interaction *, const Integrator *, const UniverseTypeAwareTime &) |
A good frame to start an integration with. | |
Variables | |
const double | halfpi = PI/2 |
const double | pi = PI |
const double | twopi = PI+PI |
const double | pisq = PI*PI |
Config * | config |
The active configuration. | |
int | nast |
double | el [6] |
double | l_ts |
double | w_ts |
double | file_time |
OrsaPaths * | orsa_paths = 0 |
LocationFile * | location_file |
OrsaPaths * | orsa_paths |
JPLFile * | jpl_file |
JPLCache * | jpl_cache |
const gsl_siman_params_t | params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE, K, T_INITIAL, MU_T, T_MIN} |
const double | d_a = 0.1 |
const double | d_e = 0.01 |
const double | d_i = 0.1*(pi/180) |
const double | d_omega_node = 1.0*(pi/180) |
const double | d_omega_pericenter = 1.0*(pi/180) |
const double | d_M = 5.0*(pi/180) |
const double | c_a = secure_pow(1.0/d_a, 2) |
const double | c_e = secure_pow(1.0/d_e, 2) |
const double | c_i = secure_pow(1.0/d_i, 2) |
const double | c_omega_node = secure_pow(1.0/d_omega_node, 2) |
const double | c_omega_pericenter = secure_pow(1.0/d_omega_pericenter, 2) |
const double | c_M = secure_pow(1.0/d_M, 2) |
const double | G_MKS = 6.67259e-11 |
const double | MSUN_MKS = 1.9889e30 |
const double | MJUPITER_MKS = 1.8989e27 |
const double | MEARTH_MKS = 5.9742e24 |
const double | MMOON_MKS = 7.3483e22 |
const double | AU_MKS = 1.49597870660e11 |
const double | c_MKS = 299792458.0 |
const double | R_EARTH_MKS = 6378140.0 |
const double | R_MOON_MKS = 1737400.0 |
TimeScale | default_Date_timescale = TT |
const TAI_minus_UTC_element | TAI_minus_UTC_table_final_element = {0,0,0,0} |
const TAI_minus_UTC_element | TAI_minus_UTC_table [] |
const ET_minus_UT_element | ET_minus_UT_table_final_element = {0,0,0,0} |
const ET_minus_UT_element | ET_minus_UT_table [] |
Units * | units |
TimeScale | default_Date_timescale |
Units * | units = 0 |
Universe * | universe = 0 |
The active universe. | |
Config * | config = 0 |
The active configuration. | |
JPLFile * | jpl_file = 0 |
JPLCache * | jpl_cache = 0 |
LocationFile * | location_file = 0 |
Universe * | universe |
The active universe. |
enum ConfigEnum |
Definition at line 51 of file orsa_config.h.
00051 { 00052 JPL_EPHEM_FILE, 00053 JPL_DASTCOM_NUM, 00054 JPL_DASTCOM_UNNUM, 00055 JPL_DASTCOM_COMET, 00056 LOWELL_ASTORB, 00057 MPC_MPCORB, 00058 MPC_COMET, 00059 MPC_NEA, 00060 MPC_DAILY, 00061 MPC_DISTANT, 00062 MPC_PHA, 00063 MPC_UNUSUALS, 00064 ASTDYS_ALLNUM_CAT, 00065 ASTDYS_ALLNUM_CTC, 00066 ASTDYS_ALLNUM_CTM, 00067 ASTDYS_UFITOBS_CAT, 00068 ASTDYS_UFITOBS_CTC, 00069 ASTDYS_UFITOBS_CTM, 00070 NEODYS_CAT, 00071 NEODYS_CTC, 00072 OBSCODE, 00073 // TLE files 00074 TLE_NASA, 00075 TLE_GEO, 00076 TLE_GPS, 00077 TLE_ISS, 00078 TLE_KEPELE, 00079 TLE_VISUAL, 00080 TLE_WEATHER, 00081 // textures 00082 TEXTURE_SUN, 00083 TEXTURE_MERCURY, 00084 TEXTURE_VENUS, 00085 TEXTURE_EARTH, 00086 TEXTURE_MOON, 00087 TEXTURE_MARS, 00088 TEXTURE_JUPITER, 00089 TEXTURE_SATURN, 00090 TEXTURE_URANUS, 00091 TEXTURE_NEPTUNE, 00092 TEXTURE_PLUTO, 00093 // 00094 NO_CONFIG_ENUM 00095 };
Definition at line 35 of file orsa_orbit_gsl.h.
00035 { Osculating, // a,e,i,node,peri,M 00036 Equinoctal // a,k,h,q,p,lambda 00037 };
enum FFTAlgorithm |
Definition at line 101 of file orsa_fft.h.
00101 { 00102 algo_FFT, // Search only power spectrum peaks 00103 algo_FFTB, // Search only power spectrum peaks, on positive and negative frequencies 00104 algo_MFT, // Original method by Laskar 00105 algo_FMFT1, // MFT with linear freq. corrections (Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137) 00106 algo_FMFT2 // FMFT1 with addition non-linear corrections (Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137) 00107 };
enum FFTSearch |
enum FFTSearchAmplitude |
enum FFTSearchPhase |
Definition at line 99 of file orsa_fft.h.
00099 {OMEGA_NODE,OMEGA_PERICENTER,OMEGA_TILDE,MM,LAMBDA,ZERO};
enum FILE_STATUS |
enum IntegratorType |
Definition at line 36 of file orsa_integrator.h.
00036 { 00037 STOER=1, 00038 BULIRSCHSTOER=2, 00039 RUNGEKUTTA=3, 00040 DISSIPATIVERUNGEKUTTA=4, 00041 RA15=5, 00042 LEAPFROG=6 00043 };
enum InteractionType |
NEWTON | |
ARMONICOSCILLATOR | |
GALACTIC_POTENTIAL_ALLEN | |
GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON | |
JPL_PLANETS_NEWTON | |
GRAVITATIONALTREE | |
NEWTON_MPI | |
RELATIVISTIC |
Definition at line 46 of file orsa_interaction.h.
00046 { 00047 NEWTON=1, 00048 ARMONICOSCILLATOR=2, 00049 GALACTIC_POTENTIAL_ALLEN=3, 00050 GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON=4, 00051 JPL_PLANETS_NEWTON=5, 00052 GRAVITATIONALTREE=6, 00053 NEWTON_MPI=7, 00054 RELATIVISTIC=8 00055 };
enum JPL_planets |
NONE | |
MERCURY | |
VENUS | |
EARTH | |
MARS | |
JUPITER | |
SATURN | |
URANUS | |
NEPTUNE | |
PLUTO | |
MOON | |
SUN | |
SOLAR_SYSTEM_BARYCENTER | |
EARTH_MOON_BARYCENTER | |
NUTATIONS | |
LIBRATIONS | |
EARTH_AND_MOON |
Definition at line 30 of file orsa_file_jpl.h.
00030 { 00031 NONE=0, 00032 MERCURY=1, 00033 VENUS=2, 00034 EARTH=3, 00035 MARS=4, 00036 JUPITER=5, 00037 SATURN=6, 00038 URANUS=7, 00039 NEPTUNE=8, 00040 PLUTO=9, 00041 MOON=10, 00042 SUN=11, 00043 SOLAR_SYSTEM_BARYCENTER=12, 00044 EARTH_MOON_BARYCENTER=13, 00045 NUTATIONS=14, 00046 LIBRATIONS=15, 00047 EARTH_AND_MOON=1000 00048 };
enum length_unit |
enum M5COLS |
enum mass_unit |
enum OrsaFileDataType |
Definition at line 381 of file orsa_file.h.
00381 { 00382 OFDT_END_OF_FILE=0, 00383 OFDT_UNIVERSE=1, 00384 OFDT_EVOLUTION=2, 00385 OFDT_FRAME=3, 00386 OFDT_BODY=4 00387 };
enum ReferenceSystem |
enum time_unit |
enum TimeScale |
TimeScale enum, useful only when using a Real Universe. More information can be obtained here: http://www.hartrao.ac.za/nccsdoc/slalib/sun67.htx/node217.html.
Definition at line 251 of file orsa_units.h.
00251 { 00252 UTC=1, 00253 UT=2, 00254 TAI=3, 00255 TDT=4, 00256 GPS=5, 00257 // aliases 00258 UT1=UT, 00259 ET=TDT, 00260 TT=TDT 00261 };
enum UniverseType |
This enum is used to classify the Universes: a Real Universe is composed.
Definition at line 149 of file orsa_universe.h.
double orsa::accurate_peak | ( | double | left, | |
double | center, | |||
double | right, | |||
fftw_complex * | in, | |||
int | size | |||
) |
Definition at line 479 of file orsa_fft.cc.
References norm(), phi(), phi_gsl_two(), gsl_d1_minimization_parameters::pointer_points_sequence, and gsl_d1_minimization_parameters::size.
00479 { 00480 00481 // check values 00482 // if ( (center<0) || (center>0.5) ) return center; // ERROR... 00483 if ( (center<(-0.5)) || (center>0.5) ) { 00484 cerr << "warning!! Peak out of range!" << endl; 00485 return center; // ERROR... 00486 } 00487 00488 // int c_iter=0, l_iter=0, r_iter=0, max_iter = 100; 00489 int max_iter = 100; 00490 double center_amp, left_amp, right_amp; 00491 00492 left_amp = norm(phi(left, in,size)); 00493 right_amp = norm(phi(right, in,size)); 00494 center_amp = norm(phi(center,in,size)); 00495 00496 // TODO: check for gsl peak search conditions!!! 00497 00498 // gsl stuff 00499 gsl_d1_minimization_parameters par; 00500 gsl_function F; 00501 gsl_min_fminimizer *s; 00502 const gsl_min_fminimizer_type *T; 00503 // 00504 T = gsl_min_fminimizer_goldensection; 00505 s = gsl_min_fminimizer_alloc(T); 00506 // 00507 par.size = size; 00508 par.pointer_points_sequence = &in[0]; 00509 // 00510 F.function = &phi_gsl_two; 00511 F.params = ∥ 00512 // 00513 gsl_min_fminimizer_set(s,&F,center,left,right); 00514 00515 double epsrel = 1.0e-10; // for the d=1 minimization 00516 double epsabs = 0.0; // for the d=1 minimization 00517 00518 int iter = 0; 00519 int status; 00520 do { 00521 iter++; 00522 00523 status = gsl_min_fminimizer_iterate (s); 00524 // 00525 center = gsl_min_fminimizer_minimum (s); 00526 left = gsl_min_fminimizer_x_lower (s); 00527 right = gsl_min_fminimizer_x_upper (s); 00528 // delta = range_stop - range_start; 00529 // minimize! 00530 status = gsl_min_test_interval(left,right,epsrel,epsabs); 00531 00532 // post-minimization... 00533 /* printf("%5d [%f, %f] %.7f %.7e bin:%i\n", 00534 iter,range_start,range_stop,test_peak.frequency,range_stop-range_start,bin); 00535 */ 00536 } while ((status == GSL_CONTINUE) && (iter < max_iter)); 00537 00538 return (center); 00539 }
Here is the call graph for this function:
void alpha_delta_meridian | ( | const JPL_planets | p, | |
const Date & | date, | |||
Angle & | alpha_zero, | |||
Angle & | delta_zero, | |||
Angle & | W | |||
) |
Definition at line 2046 of file orsa_units.cc.
References alpha_delta_meridian_moon(), cos(), EARTH, Date::GetJulian(), JUPITER, MARS, MERCURY, MOON, NEPTUNE, PLUTO, SATURN, Angle::SetDPS(), Date::SetJ2000(), Angle::SetRad(), sin(), SUN, URANUS, and VENUS.
02047 { 02048 02049 Date tmp_date; 02050 tmp_date.SetJ2000(); 02051 const Date date_J2000(tmp_date); 02052 02053 // CHECK TIMESCALES!! 02054 const double d = (date.GetJulian()-date_J2000.GetJulian()); 02055 const double T = d / 36525.0; 02056 02057 switch (p) { 02058 case SUN: 02059 alpha_zero.SetDPS(286.13,0,0); 02060 delta_zero.SetDPS(63.87,0,0); 02061 W.SetDPS(84.10+14.1844000*d,0,0); 02062 break; 02063 case MERCURY: 02064 alpha_zero.SetDPS(281.01-0.033*T,0,0); 02065 delta_zero.SetDPS(61.45-0.005*T,0,0); 02066 W.SetDPS(329.548+6.1385025*d,0,0); 02067 break; 02068 case VENUS: 02069 alpha_zero.SetDPS(272.76,0,0); 02070 delta_zero.SetDPS(67.16,0,0); 02071 W.SetDPS(160.20-1.4813688*d,0,0); 02072 break; 02073 case EARTH: 02074 alpha_zero.SetDPS(0.0-0.641*T,0,0); 02075 delta_zero.SetDPS(90.0-0.557*T,0,0); 02076 W.SetDPS(190.147+360.9856235*d,0,0); 02077 break; 02078 case MOON: 02079 alpha_delta_meridian_moon(date,alpha_zero,delta_zero,W); 02080 break; 02081 case MARS: 02082 alpha_zero.SetDPS(317.68143-0.1061*T,0,0); 02083 delta_zero.SetDPS(52.88650-0.0609*T,0,0); 02084 W.SetDPS(176.630+350.89198226*d,0,0); 02085 break; 02086 case JUPITER: 02087 alpha_zero.SetDPS(268.05-0.009*T,0,0); 02088 delta_zero.SetDPS(64.49+0.003*T,0,0); 02089 W.SetDPS(284.95+870.5366420*d,0,0); 02090 break; 02091 case SATURN: 02092 alpha_zero.SetDPS(40.589-0.036*T,0,0); 02093 delta_zero.SetDPS(83.537-0.004*T,0,0); 02094 W.SetDPS(38.90+810.7939024*d,0,0); 02095 break; 02096 case URANUS: 02097 alpha_zero.SetDPS(257.311,0,0); 02098 delta_zero.SetDPS(-15.175,0,0); 02099 W.SetDPS(203.81-501.1600928*d,0,0); 02100 break; 02101 case NEPTUNE: 02102 { 02103 Angle N; 02104 N.SetDPS(357.85+52.316*T,0,0); 02105 // 02106 alpha_zero.SetDPS(299.36+0.70*sin(N),0,0); 02107 delta_zero.SetDPS(43.46-0.51*cos(N),0,0); 02108 W.SetDPS(253.18+536.3128492*d-0.48*sin(N),0,0); 02109 } 02110 break; 02111 case PLUTO: 02112 alpha_zero.SetDPS(313.02,0,0); 02113 delta_zero.SetDPS(9.09,0,0); 02114 W.SetDPS(236.77-56.3623195*d,0,0); 02115 break; 02116 // 02117 default: 02118 alpha_zero.SetRad(0.0); 02119 delta_zero.SetRad(0.0); 02120 W.SetRad(0.0); 02121 break; 02122 } 02123 }
Here is the call graph for this function:
static void orsa::alpha_delta_meridian_moon | ( | const Date & | date, | |
Angle & | alpha_zero, | |||
Angle & | delta_zero, | |||
Angle & | W | |||
) | [static] |
Definition at line 1982 of file orsa_units.cc.
References cos(), E1(), Date::GetJulian(), Angle::SetDPS(), Date::SetJ2000(), and sin().
Referenced by alpha_delta_meridian().
01983 { 01984 Date tmp_date; 01985 tmp_date.SetJ2000(); 01986 const Date date_J2000(tmp_date); 01987 01988 const double d = (date.GetJulian()-date_J2000.GetJulian()); 01989 const double T = d / 36525.0; 01990 01991 Angle E1, E2, E3, E4, E5, E6, E7, E8, E9, E10, E11, E12, E13; 01992 01993 E1.SetDPS( 125.045 - 0.0529921*d,0,0); 01994 E2.SetDPS( 250.089 - 0.1059842*d,0,0); 01995 E3.SetDPS( 260.008 + 13.0120009*d,0,0); 01996 E4.SetDPS( 176.625 + 13.3407154*d,0,0); 01997 E5.SetDPS( 357.529 + 0.9856003*d,0,0); 01998 E6.SetDPS( 311.589 + 26.4057084*d,0,0); 01999 E7.SetDPS( 134.963 + 13.0649930*d,0,0); 02000 E8.SetDPS( 276.617 + 0.3287146*d,0,0); 02001 E9.SetDPS( 34.226 + 1.7484877*d,0,0); 02002 E10.SetDPS( 15.134 - 0.1589763*d,0,0); 02003 E11.SetDPS(119.743 + 0.0036096*d,0,0); 02004 E12.SetDPS(239.961 + 0.1643573*d,0,0); 02005 E13.SetDPS( 25.053 + 12.9590088*d,0,0); 02006 02007 alpha_zero.SetDPS(269.9949 02008 + 0.0031*T 02009 - 3.8787*sin(E1) 02010 - 0.1204*sin(E2) 02011 + 0.0700*sin(E3) 02012 - 0.0172*sin(E4) 02013 + 0.0072*sin(E6) 02014 - 0.0052*sin(E10) 02015 + 0.0043*sin(E13),0,0); 02016 02017 delta_zero.SetDPS(66.5392 02018 + 0.0130*T 02019 + 1.5419*cos(E1) 02020 + 0.0239*cos(E2) 02021 - 0.0278*cos(E3) 02022 + 0.0068*cos(E4) 02023 - 0.0029*cos(E6) 02024 + 0.0009*cos(E7) 02025 + 0.0008*cos(E10) 02026 - 0.0009*cos(E13),0,0); 02027 02028 W.SetDPS(38.3213 02029 + 13.17635815*d 02030 - 1.4e-12*d*d 02031 + 3.5610*sin(E1) 02032 + 0.1208*sin(E2) 02033 - 0.0642*sin(E3) 02034 + 0.0158*sin(E4) 02035 + 0.0252*sin(E5) 02036 - 0.0066*sin(E6) 02037 - 0.0047*sin(E7) 02038 - 0.0046*sin(E8) 02039 + 0.0028*sin(E9) 02040 + 0.0052*sin(E10) 02041 + 0.0040*sin(E11) 02042 + 0.0019*sin(E12) 02043 - 0.0044*sin(E13),0,0); 02044 }
Here is the call graph for this function:
void orsa::amph | ( | double * | amp, | |
double * | phase, | |||
fftw_complex * | phiR, | |||
fftw_complex * | phiL, | |||
double | freq, | |||
fftw_complex * | in, | |||
int | size | |||
) |
Definition at line 467 of file orsa_fft.cc.
References norm(), phi(), and secure_atan2().
00467 { 00468 00469 *phiR = phi(freq, in,size); 00470 *phiL = phi(-freq,in,size); 00471 00472 *amp = norm(*phiR); 00473 *phase = secure_atan2(phiR->im,phiR->re); 00474 00475 // cerr << "amph() DUMP: amp = " << *amp << " phase = " << *phase << endl; 00476 00477 }
Here is the call graph for this function:
Vector orsa::AngularMomentum | ( | const std::vector< Body > & | , | |
const Vector & | ||||
) |
Vector orsa::AngularMomentum | ( | const vector< Body > & | f, | |
const Vector & | center | |||
) |
Definition at line 188 of file orsa_frame.cc.
References ExternalProduct().
Referenced by Frame::AngularMomentum(), and BarycentricAngularMomentum().
00188 { 00189 Vector vec_sum; 00190 unsigned int k; 00191 for (k=0;k<f.size();++k) { 00192 if (f[k].mass() > 0) { 00193 vec_sum += f[k].mass()*ExternalProduct(f[k].position()-center,f[k].velocity()); 00194 } 00195 } 00196 return (vec_sum); 00197 }
Here is the call graph for this function:
void orsa::apply_window | ( | fftw_complex * | signal_win, | |
fftw_complex * | signal, | |||
int | size | |||
) |
Definition at line 449 of file orsa_fft.cc.
00449 { 00450 00451 int j; 00452 double win_coeff,arg; 00453 00454 for(j=0;j<size;j++) { 00455 00456 arg = twopi*j/size; 00457 win_coeff = (1. - cos(arg)); 00458 00459 signal_win[j].re = signal[j].re * win_coeff; 00460 signal_win[j].im = signal[j].im * win_coeff; 00461 00462 } 00463 00464 }
Here is the call graph for this function:
Vector orsa::Barycenter | ( | const std::vector< Body > & | ) |
Vector orsa::Barycenter | ( | const vector< Body > & | f | ) |
Definition at line 148 of file orsa_frame.cc.
References CenterOfMass().
Referenced by Frame::Barycenter().
00148 { 00149 return (orsa::CenterOfMass(f)); 00150 }
Here is the call graph for this function:
Vector orsa::BarycenterVelocity | ( | const std::vector< Body > & | ) |
Vector orsa::BarycenterVelocity | ( | const vector< Body > & | f | ) |
Definition at line 152 of file orsa_frame.cc.
References CenterOfMassVelocity().
Referenced by Frame::BarycenterVelocity().
00152 { 00153 return (orsa::CenterOfMassVelocity(f)); 00154 }
Here is the call graph for this function:
Vector orsa::BarycentricAngularMomentum | ( | const std::vector< Body > & | ) |
Vector orsa::BarycentricAngularMomentum | ( | const vector< Body > & | f | ) |
Definition at line 203 of file orsa_frame.cc.
References AngularMomentum(), and Frame::Barycenter().
Referenced by Frame::BarycentricAngularMomentum().
00203 { 00204 return (orsa::AngularMomentum(f,Barycenter(f))); 00205 }
Here is the call graph for this function:
void orsa::BuildObservationTriplets | ( | const vector< Observation > & | obs, | |
vector< ThreeObservations > & | triplets, | |||
const double | optimal_dt = FromUnits(1,DAY) | |||
) |
Definition at line 122 of file orsa_orbit_gsl.cc.
References Observation::date, DAY, FromUnits(), and Date::GetJulian().
Referenced by Compute().
00122 { 00123 // sort(obs.begin(),obs.end()); 00124 triplets.clear(); 00125 unsigned int l,m,n; 00126 const unsigned int s = obs.size(); 00127 // build only triplets with observations delays smaller than and obs_delay_max 00128 const double obs_delay_max = FromUnits(180,DAY); 00129 // max ration between dalays 00130 const double obs_delay_ration_max = 10.0; 00131 // 00132 Observation ol,om,on; 00133 double dt12, dt23, dt_ratio; 00134 ThreeObservations triobs; triobs.resize(3); 00135 l=0; 00136 while (l<s) { 00137 ol=obs[l]; 00138 triobs[0] = ol; 00139 m=l+1; 00140 while (m<s) { 00141 om=obs[m]; 00142 triobs[1] = om; 00143 n=m+1; 00144 while (n<s) { 00145 on=obs[n]; 00146 triobs[2] = on; 00147 00148 triobs.ComputeWeight(optimal_dt); 00149 00150 dt12 = FromUnits(om.date.GetJulian()-ol.date.GetJulian(),DAY); 00151 dt23 = FromUnits(on.date.GetJulian()-om.date.GetJulian(),DAY); 00152 00153 if (dt12>dt23) 00154 dt_ratio = dt12/dt23; 00155 else 00156 dt_ratio = dt23/dt12; 00157 00158 if ((dt12<obs_delay_max) && (dt23<obs_delay_max) && (dt_ratio<obs_delay_ration_max)) triplets.push_back(triobs); 00159 00160 // debug 00161 std::cerr << l << "-" << m << "-" << n << std::endl; 00162 00163 ++n; 00164 } 00165 ++m; 00166 } 00167 ++l; 00168 } 00169 00170 sort(triplets.begin(),triplets.end()); 00171 }
Here is the call graph for this function:
Vector orsa::CenterOfMass | ( | const std::vector< Body > & | ) |
Vector orsa::CenterOfMass | ( | const vector< Body > & | f | ) |
Definition at line 122 of file orsa_frame.cc.
Referenced by Barycenter(), and Frame::CenterOfMass().
00122 { 00123 Vector sum_vec(0,0,0); 00124 double sum_mass = 0.0; 00125 for (unsigned int k=0; k<f.size(); ++k) { 00126 const double mass = f[k].mass(); 00127 if (mass > 0.0) { 00128 sum_vec += mass*f[k].position(); 00129 sum_mass += mass; 00130 } 00131 } 00132 return (sum_vec/sum_mass); 00133 }
Vector orsa::CenterOfMassVelocity | ( | const std::vector< Body > & | ) |
Vector orsa::CenterOfMassVelocity | ( | const vector< Body > & | f | ) |
Definition at line 135 of file orsa_frame.cc.
Referenced by BarycenterVelocity(), and Frame::CenterOfMassVelocity().
00135 { 00136 Vector sum_vec(0,0,0); 00137 double sum_mass = 0.0; 00138 for (unsigned int k=0; k<f.size(); ++k) { 00139 const double mass = f[k].mass(); 00140 if (mass > 0.0) { 00141 sum_vec += mass*f[k].velocity(); 00142 sum_mass += mass; 00143 } 00144 } 00145 return (sum_vec/sum_mass); 00146 }
double orsa::CloseApproaches_gsl_f | ( | const gsl_vector * | v, | |
void * | params | |||
) |
Definition at line 2530 of file orsa_orbit_gsl.cc.
Referenced by SearchCloseApproaches().
02530 { 02531 02532 CloseApproaches_gsl_parameters * p = (CloseApproaches_gsl_parameters *) params; 02533 02534 Frame f = p->f; 02535 02536 const UniverseTypeAwareTime gsl_time(gsl_vector_get(v,0)); 02537 02538 p->e->clear(); 02539 p->e->push_back(f); 02540 p->e->SetSamplePeriod(f - gsl_time); 02541 p->e->SetMaxUnsavedSubSteps(10000); 02542 p->e->Integrate(gsl_time,true); 02543 02544 f = (*(p->e))[p->e->size()-1]; 02545 02546 return (f[p->obj_index].position() - f[p->index].position()).Length(); 02547 }
int orsa::compare_binamp | ( | const binamp * | a, | |
const binamp * | b | |||
) |
sort binamp struct from the bigger to the smaller...
Definition at line 153 of file orsa_fft.cc.
Referenced by psd_max_again_many().
00153 { 00154 00155 if (( (*a).amp - (*b).amp) < 0) 00156 return 1; 00157 if (( (*a).amp - (*b).amp) > 0) 00158 return -1; 00159 00160 return 0; 00161 }
OrbitWithCovarianceMatrixGSL orsa::Compute | ( | const std::vector< Observation > & | ) |
General interface to the Orbit computation from a set of observations.
Referenced by OrbitWithEpoch::Compute().
OrbitWithCovarianceMatrixGSL orsa::Compute | ( | const vector< Observation > & | obs | ) |
Definition at line 1396 of file orsa_orbit_gsl.cc.
References BuildObservationTriplets(), Compute_Gauss(), Compute_TestMethod(), DAY, Sky::delta_arcsec(), FromUnits(), Universe::GetUniverseType(), OrbitDifferentialCorrectionsLeastSquares(), ORSA_ERROR, pi, OptimizedOrbitPositions::PropagatedSky_J2000(), Real, and universe.
01396 { 01397 01398 // vector<Observation> obs(obs_in); 01399 01400 if (obs.size() < 3) { 01401 ORSA_ERROR("at least three observations are needed."); 01402 OrbitWithCovarianceMatrixGSL o; 01403 return o; 01404 } 01405 01406 if (universe->GetUniverseType() != Real) { 01407 ORSA_ERROR("trying to compute an orbit using the wrong universe type: return."); 01408 OrbitWithCovarianceMatrixGSL o; 01409 return o; 01410 } 01411 01412 // test 01413 /* 01414 { 01415 OrbitWithEpoch orbit; 01416 01417 orbit.a = 1.75; 01418 orbit.e = 0.39705; 01419 orbit.i = 54.52*(pi/180); 01420 orbit.omega_node = 47.46 *(pi/180); 01421 orbit.omega_pericenter = 95.47 *(pi/180); 01422 orbit.M = 322.918 *(pi/180); 01423 // 01424 Date date; 01425 date.SetJulian(52950+2400000.5); 01426 orbit.epoch.SetDate(date); 01427 // 01428 OrbitDifferentialCorrectionsLeastSquares(orbit,obs); 01429 exit(0); 01430 } 01431 */ 01432 01433 // needed? 01434 // sort(obs.begin(),obs.end()); 01435 01436 // cerr << "Orbit::Compute() observations: " << obs.size() << endl; 01437 01438 // NOTE: check the reference system somewhere near here!! 01439 01440 /* 01441 if (0) { 01442 unsigned int k=0; 01443 while (k<obs.size()) { 01444 printf("obs. %i: JD = %f\n",k,obs[k].date.GetJulian()); 01445 ++k; 01446 } 01447 } 01448 */ 01449 01450 // vector<OrbitWithEpoch> preliminary_orbits; 01451 vector<PreliminaryOrbit> preliminary_orbits; 01452 01453 if (0) { 01454 vector<Observation> test_obs; 01455 // 01456 test_obs.push_back(obs[0]); 01457 test_obs.push_back(obs[obs.size()/2]); 01458 test_obs.push_back(obs[obs.size()-1]); 01459 // 01460 Compute_TestMethod(test_obs,preliminary_orbits); 01461 } 01462 01463 if (0) { 01464 Compute_TestMethod(obs,preliminary_orbits); 01465 } 01466 01467 // yet another test... 01468 if (0) { 01469 vector<Observation> test_obs; 01470 // 01471 test_obs.push_back(obs[0]); 01472 test_obs.push_back(obs[obs.size()/2]); 01473 test_obs.push_back(obs[obs.size()-1]); 01474 // 01475 Compute_Gauss(test_obs,preliminary_orbits); 01476 } 01477 01478 if (1) { 01479 vector<ThreeObservations> triplet; 01480 triplet.clear(); 01481 // BuildObservationTriplets(obs, triplet, FromUnits(2,HOUR)); 01482 // BuildObservationTriplets(obs, triplet, FromUnits(20,DAY)); 01483 BuildObservationTriplets(obs, triplet, FromUnits(270,DAY)); 01484 // 01485 cerr << "triplets: " << triplet.size() << endl; 01486 01487 unsigned int trial=0; 01488 while (trial<triplet.size()) { 01489 // cerr << "using triplet with weight = " << triplet[trial].Weight() << endl; 01490 // 01491 cerr << "delay 1->2: " << triplet[trial][1].date.GetJulian()-triplet[trial][0].date.GetJulian() << " days" << endl; 01492 cerr << "delay 2->3: " << triplet[trial][2].date.GetJulian()-triplet[trial][1].date.GetJulian() << " days" << endl; 01493 // 01494 01495 Compute_Gauss(triplet[trial],preliminary_orbits); 01496 // Compute_TestMethod(triplet[trial],preliminary_orbits); 01497 01498 ++trial; 01499 // 01500 // if (trial > 1) break; 01501 // if (trial > 1) break; 01502 // if (preliminary_orbits.size() > 7) break; 01503 if (trial > 500) break; 01504 01505 cerr << "----trial: " << trial << endl; 01506 } 01507 01508 } 01509 01510 // compute rms for each preliminary orbit 01511 { // don't remove THIS!! 01512 unsigned int r=0; 01513 while (r<preliminary_orbits.size()) { 01514 preliminary_orbits[r].ComputeRMS(obs); 01515 printf("%20.12f %20.12f %20.12f %20.12f\n",preliminary_orbits[r].a,preliminary_orbits[r].e,preliminary_orbits[r].i*(180/pi),preliminary_orbits[r].GetRMS()); 01516 ++r; 01517 } 01518 } 01519 // 01520 // sort preliminary orbits 01521 sort(preliminary_orbits.begin(),preliminary_orbits.end()); 01522 01523 { 01524 cerr << "total prelim. orbits: " << preliminary_orbits.size() << endl; 01525 unsigned int r=0; 01526 while (r<preliminary_orbits.size()) { 01527 printf("%20.12f %20.12f %20.12f %20.12f\n",preliminary_orbits[r].a,preliminary_orbits[r].e,preliminary_orbits[r].i*(180/pi),preliminary_orbits[r].GetRMS()); 01528 ++r; 01529 } 01530 } 01531 01532 // only the first 01533 OrbitDifferentialCorrectionsLeastSquares(preliminary_orbits[0],obs); 01534 01535 if (0) { 01536 unsigned int r=0; 01537 // double rms; 01538 while (r<preliminary_orbits.size()) { 01539 // rms = RMS_residuals(obs,preliminary_orbits[r]); 01540 // printf("%20.12f %20.12f %20.12f %20.12f\n",preliminary_orbits[r].a,preliminary_orbits[r].e,preliminary_orbits[r].i*(180/pi),rms); 01541 // if (rms < 1.0*3600) OrbitSimulatedAnnealing(preliminary_orbits[r],obs); 01542 // if (rms < 500.0) OrbitDifferentialCorrectionsMultimin(preliminary_orbits[r],obs); 01543 // if (rms < 200.0) OrbitDifferentialCorrectionsLeastSquares(preliminary_orbits[r],obs); 01544 01545 OrbitDifferentialCorrectionsLeastSquares(preliminary_orbits[r],obs); 01546 01547 // test... 01548 if (0) { 01549 01550 OptimizedOrbitPositions opt(preliminary_orbits[r]); 01551 01552 vector<Observation> good_obs; 01553 unsigned int k; 01554 for (k=0;k<obs.size();++k) { 01555 // if (PropagatedSky(preliminary_orbits[r],obs[k].date,obs[k].obscode).delta_arcsec(obs[k]) < 60.0) { // 5.0 01556 // if (PropagatedSky(preliminary_orbits[r],obs[k].date,obs[k].obscode).delta_arcsec(obs[k]) < 60.0) { // 5.0 01557 if (opt.PropagatedSky_J2000(obs[k].date,obs[k].obscode).delta_arcsec(obs[k]) < 300.0) { // 5.0, 60.0 01558 good_obs.push_back(obs[k]); 01559 } 01560 if (good_obs.size()>12) break; 01561 } 01562 01563 cerr << "good obs.: " << good_obs.size() << endl; 01564 01565 // improve orbit using 'good observations' 01566 if (good_obs.size() > 3) OrbitDifferentialCorrectionsLeastSquares(preliminary_orbits[r],good_obs); 01567 } 01568 01569 ++r; 01570 } 01571 01572 } 01573 01574 // WARNING: check on the preliminary_orbits dimension... 01575 // to be refined! 01576 return preliminary_orbits[0]; 01577 }
Here is the call graph for this function:
void orsa::Compute_Gauss | ( | const std::vector< Observation > & | , | |
std::vector< PreliminaryOrbit > & | ||||
) |
Referenced by Compute().
void orsa::Compute_Gauss | ( | const vector< Observation > & | obs_in, | |
vector< PreliminaryOrbit > & | preliminary_orbits | |||
) |
Definition at line 1695 of file orsa_orbit_gsl.cc.
References Orbit::a, A, OrbitWithEpoch::Compute(), Sky::Compute_J2000(), cos(), Cross(), DAY, Sky::dec(), Orbit::e, EARTH, ECLIPTIC, FromUnits(), GetC(), GetG(), JPLCache::GetJPLBody(), GetMSun(), Angle::GetRad(), Universe::GetReferenceSystem(), Orbit::i, jpl_cache, Vector::Length(), Vector::LengthSquared(), location_file, Vector::Normalized(), obleq_J2000(), LocationFile::ObsPos(), ORSA_ERROR, params, pi, poly_8_gsl_solve(), Body::position(), Sky::ra(), RMS_residuals(), Vector::rotate(), secure_pow(), sin(), SUN, universe, Vector::y, and Vector::z.
01695 { 01696 01697 cerr << "inside Compute_Gauss(...)" << endl; 01698 01699 // see A. E. Roy, "Orbital Motion", section 13.4 01700 01701 PreliminaryOrbit orbit; 01702 01703 vector<Observation> obs(obs_in); 01704 01705 if (obs.size() != 3) { 01706 ORSA_ERROR("Compute_Gauss(): three observations needed."); 01707 return; 01708 } 01709 01710 unsigned int j; 01711 01712 // TODO: the obs. are sorted in time? 01713 01714 sort(obs.begin(),obs.end()); 01715 01716 // debug 01717 /* { 01718 int l; 01719 for (l=0;l<3;++l) { 01720 printf("observation %i: julian=%14.5f ra=%g dec=%g\n",l,obs[l].date.GetJulian(),obs[l].ra.GetRad(),obs[l].dec.GetRad()); 01721 } 01722 } 01723 */ 01724 01725 double tau[3]; 01726 const double sqrt_GM = sqrt(GetG()*GetMSun()); 01727 // 01728 tau[2] = sqrt_GM*FromUnits(obs[1].date.GetJulian()-obs[0].date.GetJulian(),DAY); 01729 tau[0] = sqrt_GM*FromUnits(obs[2].date.GetJulian()-obs[1].date.GetJulian(),DAY); 01730 tau[1] = sqrt_GM*FromUnits(obs[2].date.GetJulian()-obs[0].date.GetJulian(),DAY); 01731 01732 /* 01733 Config conf; 01734 OrsaConfigFile ocf(&conf); 01735 ocf.Read(); 01736 ocf.Close(); 01737 */ 01738 01739 Vector r; 01740 Vector R[3]; // heliocentric earth position 01741 { 01742 /* 01743 Vector tmp_v; 01744 JPLFile jf(config->paths[JPL_EPHEM_FILE]->GetValue().c_str()); 01745 for (j=0;j<3;++j) { 01746 jf.GetEph(obs[j].date,EARTH,SUN,r,tmp_v); 01747 R[j] = r; 01748 } 01749 */ 01750 // 01751 for (j=0;j<3;++j) { 01752 R[j] = jpl_cache->GetJPLBody(EARTH,obs[j].date).position() - jpl_cache->GetJPLBody(SUN ,obs[j].date).position(); 01753 } 01754 } 01755 01756 // add the observer location to the earth position 01757 { 01758 /* 01759 LocationFile lf; 01760 lf.SetFileName(config->paths[OBSCODE]->GetValue().c_str()); 01761 lf.Read(); 01762 lf.Close(); 01763 */ 01764 for (j=0;j<3;++j) { 01765 R[j] += location_file->ObsPos(obs[j].obscode,obs[j].date); 01766 } 01767 } 01768 01769 Vector u_rho[3]; // object direction from the observer (unit vectors) 01770 for (j=0;j<3;++j) { 01771 u_rho[j].x = cos(obs[j].dec) * cos(obs[j].ra); 01772 u_rho[j].y = cos(obs[j].dec) * sin(obs[j].ra); 01773 u_rho[j].z = sin(obs[j].dec); 01774 01775 /* 01776 if (universe->GetReferenceSystem() == ECLIPTIC) { 01777 u_rho[j].rotate(0.0,-obleq(obs[j].date).GetRad(),0.0); 01778 } 01779 */ 01780 01781 /* 01782 if (universe->GetReferenceSystem() == ECLIPTIC) { 01783 EquatorialToEcliptic_J2000(u_rho[j]); 01784 } 01785 */ 01786 01787 if (universe->GetReferenceSystem() == ECLIPTIC) { 01788 // EquatorialToEcliptic_J2000(u_rho[j]); 01789 Angle obl = obleq_J2000(); 01790 u_rho[j].rotate(0.0,-obl.GetRad(),0.0); 01791 } 01792 01793 } 01794 01795 Vector cross = Cross(u_rho[0],u_rho[2]); 01796 // const Vector f = cross/cross.Length(); 01797 const Vector f = cross.Normalized(); 01798 01799 const double rho_1_f = u_rho[1]*f; 01800 01801 const double R_0_f = R[0]*f; 01802 const double R_1_f = R[1]*f; 01803 const double R_2_f = R[2]*f; 01804 01805 const double A = (tau[0]/tau[1]*R_0_f + tau[2]/tau[1]*R_2_f - R_1_f)/rho_1_f; 01806 // const double B = (tau[0]/tau[1]*(secure_pow(tau[1],2)-secure_pow(tau[0],2))*R_0_f + tau[2]/tau[1]*(secure_pow(tau[1],2)-secure_pow(tau[2],2))*R_2_f)/rho_1_f/6.0; 01807 const double B = (tau[0]/tau[1]*(tau[1]*tau[1]-tau[0]*tau[0])*R_0_f + tau[2]/tau[1]*(tau[1]*tau[1]-tau[2]*tau[2])*R_2_f)/rho_1_f/6.0; 01808 01809 const double Xl_Ym_Zn = R[1]*u_rho[1]; 01810 01811 poly_8_params params; 01812 params.coeff_6 = R[1].LengthSquared() + A*A + 2*A*Xl_Ym_Zn; 01813 params.coeff_3 = 2*A*B + 2*B*Xl_Ym_Zn; 01814 params.coeff_0 = B*B; 01815 vector<poly_8_solution> solutions; 01816 poly_8_gsl_solve(params,solutions); 01817 01818 if (solutions.size() > 0) { 01819 01820 Vector rho[3]; 01821 Vector r[3]; 01822 Vector v; // v[0] 01823 double c[3]; 01824 01825 double tmp_length; 01826 double tmp_value; 01827 unsigned int p; 01828 for (p=0;p<solutions.size();++p) { 01829 01830 // rho[1] = u_rho[1]*(A+(B/secure_pow(solutions[p].value,3))); 01831 // check 01832 // tmp_length = A + (B/secure_pow(solutions[p].value,3)); 01833 tmp_value = solutions[p].value; 01834 tmp_length = A + (B/(tmp_value*tmp_value*tmp_value)); 01835 // cerr << "tmp_length: " << tmp_length << endl; 01836 if (tmp_length <= 0.0) { 01837 // cerr << "out..." << endl; 01838 continue; 01839 } 01840 // 01841 rho[1] = u_rho[1]*tmp_length; 01842 01843 r[1] = R[1] + rho[1]; 01844 01845 // standard relation 01846 for (j=0;j<3;++j) { 01847 // c[j] = tau[j]/tau[1]*(1+(secure_pow(tau[1],2)-secure_pow(tau[j],2))/(6*secure_pow(r[1].Length(),3))); 01848 c[j] = tau[j]/tau[1]*(1+(tau[1]*tau[1]-tau[j]*tau[j])/(6*secure_pow(r[1].Length(),3))); 01849 // printf("OLD c[%i] = %g\n",j,c[j]); 01850 } 01851 01852 { 01853 01854 const Vector v_k = rho[1] - (c[0]*R[0] + c[2]*R[2] - R[1]); 01855 const double k = v_k.Length(); 01856 // const Vector u_k = v_k/v_k.Length(); 01857 const Vector u_k = v_k.Normalized(); 01858 01859 const double s02 = u_rho[0]*u_rho[2]; 01860 const double s0k = u_rho[0]*u_k; 01861 const double s2k = u_rho[2]*u_k; 01862 01863 // rho[0] = u_rho[0]*(k*(s0k-s02*s2k)/(1-secure_pow(s02,2)))/c[0]; 01864 // tmp_length = (k*(s0k-s02*s2k)/(1-secure_pow(s02,2)))/c[0]; 01865 tmp_length = (k*(s0k-s02*s2k)/(1-s02*s02))/c[0]; 01866 if (tmp_length <= 0.0) { 01867 // cerr << "out..." << endl; 01868 continue; 01869 } 01870 // 01871 rho[0] = u_rho[0]*tmp_length; 01872 01873 // rho[2] = u_rho[2]*(k*(s2k-s02*s0k)/(1-secure_pow(s02,2)))/c[2]; 01874 // tmp_length = (k*(s2k-s02*s0k)/(1-secure_pow(s02,2)))/c[2]; 01875 tmp_length = (k*(s2k-s02*s0k)/(1-s02*s02))/c[2]; 01876 if (tmp_length <= 0.0) { 01877 // cerr << "out..." << endl; 01878 continue; 01879 } 01880 // 01881 rho[2] = u_rho[2]*tmp_length; 01882 01883 r[0] = rho[0] + R[0]; 01884 r[2] = rho[2] + R[2]; 01885 01886 } 01887 01888 if (0) { 01889 // refinements 01890 int iter_gibbs=0; 01891 double old_c[3]; 01892 double delta_c = 1.0; 01893 while ((iter_gibbs<10) && (delta_c > 1.0e-6)) { 01894 01895 // Gibbs relation (see i.e. F. Zagar, "Astronomia Sferica e teorica", chap. XVIII) 01896 double Gibbs_B[3]; 01897 Gibbs_B[0] = (tau[1]*tau[2]-secure_pow(tau[0],2))/12; 01898 Gibbs_B[1] = (tau[0]*tau[2]-secure_pow(tau[1],2))/12; 01899 Gibbs_B[2] = (tau[0]*tau[1]-secure_pow(tau[2],2))/12; 01900 double Brm3[3]; 01901 for (j=0;j<3;++j) { 01902 Brm3[j] = Gibbs_B[j]/secure_pow(r[j].Length(),3); 01903 } 01904 delta_c = 0.0; 01905 for (j=0;j<3;++j) { 01906 old_c[j] = c[j]; 01907 c[j] = tau[j]/tau[1]*(1+Brm3[j])/(1-Brm3[1]); 01908 delta_c += secure_pow(c[j]-old_c[j],2); 01909 // printf("NEW c[%i] = %g\n",j,c[j]); 01910 } 01911 delta_c /= 3; 01912 delta_c = sqrt(delta_c); 01913 01914 { 01915 const Vector v_k = rho[1] - (c[0]*R[0] + c[2]*R[2] - R[1]); 01916 const double k = v_k.Length(); 01917 // const Vector u_k = v_k/Length(v_k); 01918 const Vector u_k = v_k.Normalized(); 01919 01920 const double s02 = u_rho[0]*u_rho[2]; 01921 const double s0k = u_rho[0]*u_k; 01922 const double s2k = u_rho[2]*u_k; 01923 01924 rho[0] = u_rho[0]*(k*(s0k-s02*s2k)/(1-secure_pow(s02,2)))/c[0]; 01925 rho[2] = u_rho[2]*(k*(s2k-s02*s0k)/(1-secure_pow(s02,2)))/c[2]; 01926 01927 r[0] = rho[0] + R[0]; 01928 r[2] = rho[2] + R[2]; 01929 } 01930 01931 ++iter_gibbs; 01932 } 01933 01934 } 01935 01936 // try a simpler rule 01937 // v = (r[1]-r[0])/(FromUnits(obs[0].date.GetJulian()-obs[1].date.GetJulian(),DAY)); 01938 /* v = ( (r[1]-r[0])/(FromUnits(obs[1].date.GetJulian()-obs[0].date.GetJulian(),DAY)) + 01939 (r[2]-r[1])/(FromUnits(obs[2].date.GetJulian()-obs[1].date.GetJulian(),DAY)) ) / 2.0; 01940 */ 01941 // v = velocity at the epoch of the first obs, obs[0] 01942 // Vector v = (r[1]-r[0])/(FromUnits(obs[0].date.GetJulian()-obs[1].date.GetJulian(),DAY)); 01943 Vector v = (r[1]-r[0])/(FromUnits(obs[1].date.GetJulian()-obs[0].date.GetJulian(),DAY)); 01944 01945 // light-time correction [to be checked!] 01946 r[0] += v*(r[0]-R[0]).Length()/GetC(); 01947 01948 // orbit.ref_body = Body("Sun",GetMSun(),Vector(0,0,0),Vector(0,0,0)); 01949 // orbit.Compute(r[1],v,GetG()*GetMSun()); 01950 orbit.Compute(r[0],v,GetG()*GetMSun(),obs[0].date); 01951 01952 // cerr << "HHH ref body mass : " << ref_body.mass() << endl; 01953 01954 // orbit.epoch = obs[1].date; 01955 // orbit.epoch = obs[0].date; 01956 01957 // quick test 01958 if (orbit.e < 1.0) { 01959 01960 // const Vector old_v = v; 01961 // test! 01962 // 01963 // test... maybe is needed... 01964 // gauss_v(r[0],v,obs); 01965 // 01966 // 01967 // update the orbit with the improved velocity 01968 // orbit.Compute(r[1],v,GetG()*GetMSun()); 01969 orbit.Compute(r[0],v,GetG()*GetMSun(),obs[0].date); 01970 01971 // debug 01972 cerr << " ---> orbit.a = " << orbit.a << endl; 01973 cerr << " ---> orbit.e = " << orbit.e << endl; 01974 cerr << " ---> orbit.i = " << orbit.i*(180/pi) << endl; 01975 01976 // Sky sky; 01977 // UniverseTypeAwareTime utat; utat.SetDate(obs[1].date); 01978 // sky.Compute(rho[1],utat); 01979 // cerr << obs[1].ra.GetRad() << " " << obs[1].dec.GetRad() << " " << sky.ra().GetRad() << " " << sky.dec().GetRad() << " " << RMS_residuals(obs) << " " << obs[1].date.GetJulian()-obs[0].date.GetJulian() << " " << obs[2].date.GetJulian()-obs[1].date.GetJulian() << endl; 01980 01981 // cerr << "better..." << endl; 01982 /* 01983 int j; 01984 Sky sky; 01985 UniverseTypeAwareTime utat; 01986 for (j=0;j<3;++j) { 01987 utat.SetDate(obs[j].date); 01988 // sky.Compute(rho[j],utat); 01989 sky.Compute(r[j]-R[j],utat); 01990 fprintf(stderr,"SKY ra = %30.20f dec = %30.20f\n",obs[j].ra.GetRad(),obs[j].dec.GetRad()); 01991 fprintf(stderr,"COMP ra = %30.20f dec = %30.20f\n",sky.ra().GetRad(),sky.dec().GetRad()); 01992 } 01993 */ 01994 01995 fprintf(stderr,"limited RMS: %g\n",RMS_residuals(obs,orbit)); 01996 01997 if (0) { 01998 UniverseTypeAwareTime utat; 01999 Sky sky; 02000 double hand_RMS=0; 02001 for (j=0;j<3;++j) { 02002 // sky.Compute(r[j]-R[j],utat); 02003 sky.Compute_J2000(r[j]-R[j]); 02004 hand_RMS += ( secure_pow((obs[j].ra.GetRad()-sky.ra().GetRad())*cos(obs[j].dec.GetRad()),2) + 02005 secure_pow(obs[j].dec.GetRad()-sky.dec().GetRad(),2) ); 02006 } 02007 hand_RMS /= 3.0; 02008 hand_RMS = sqrt(hand_RMS); 02009 hand_RMS *= (180/pi)*3600; 02010 fprintf(stderr,"RMS by hands: %g\n",hand_RMS); 02011 02012 // tests 02013 if (hand_RMS > 1.0) { 02014 cerr << "controlled exit..." << endl; 02015 exit(0); 02016 } 02017 02018 } 02019 02020 } else { 02021 // cerr << "orbit eccentricity > 1.0: a=" << orbit.a << " e=" << orbit.e << " i=" << orbit.i*180/pi << endl; 02022 } 02023 02024 /* 02025 if (orbit.e > 0.8) { 02026 cerr << "prelim. orbit: a = " << orbit.a << " e = " << orbit.e << " i = " << orbit.i*(180/pi) << endl; 02027 } else { 02028 cerr << "prelim. orbit: a = " << orbit.a << " e = " << orbit.e << " i = " << orbit.i*(180/pi) << " RMS residuals (arcsec): " << RMS_residuals(obs,orbit) << endl; 02029 } 02030 */ 02031 02032 if (orbit.e < 1.0) preliminary_orbits.push_back(orbit); 02033 } 02034 } 02035 02036 cerr << "out of Compute_Gauss(...)" << endl; 02037 }
Here is the call graph for this function:
void orsa::Compute_TestMethod | ( | const std::vector< Observation > & | , | |
std::vector< PreliminaryOrbit > & | ||||
) |
Referenced by Compute().
void orsa::Compute_TestMethod | ( | const vector< Observation > & | obs_in, | |
vector< PreliminaryOrbit > & | preliminary_orbits | |||
) |
Definition at line 2042 of file orsa_orbit_gsl.cc.
References Orbit::a, AU, OrbitWithEpoch::Compute(), config, cos(), DAY, Orbit::e, EARTH, ECLIPTIC, FromUnits(), gauss_v(), GetG(), GetMSun(), Universe::GetReferenceSystem(), Orbit::i, JPL_EPHEM_FILE, location_file, obleq_J2000(), LocationFile::ObsPos(), Config::paths, pi, Vector::rotate(), sin(), SUN, universe, Vector::y, and Vector::z.
02042 { 02043 02044 cerr << "inside Compute_TestMethod(...)" << endl; 02045 02046 PreliminaryOrbit orbit; 02047 02048 vector<Observation> obs(obs_in); 02049 02050 unsigned int j; 02051 02052 sort(obs.begin(),obs.end()); 02053 02054 /* 02055 Config conf; 02056 OrsaConfigFile ocf(&conf); 02057 ocf.Read(); 02058 ocf.Close(); 02059 */ 02060 02061 const unsigned int size = obs.size(); 02062 02063 Vector R[size]; // heliocentric earth position 02064 { 02065 Vector tmp_r, tmp_v; 02066 JPLFile jf(config->paths[JPL_EPHEM_FILE]->GetValue().c_str()); 02067 for (j=0;j<size;++j) { 02068 jf.GetEph(obs[j].date,EARTH,SUN,tmp_r,tmp_v); 02069 R[j] = tmp_r; 02070 } 02071 } 02072 02073 // add the observer location to the earth position 02074 { 02075 02076 /* 02077 LocationFile lf; 02078 lf.SetFileName(config->paths[OBSCODE]->GetValue().c_str()); 02079 lf.Read(); 02080 lf.Close(); 02081 */ 02082 for (j=0;j<size;++j) { 02083 R[j] += location_file->ObsPos(obs[j].obscode,obs[j].date); 02084 } 02085 } 02086 02087 Vector u_rho[size]; // object direction from the observer (unit vectors) 02088 for (j=0;j<size;++j) { 02089 u_rho[j].x = cos(obs[j].dec) * cos(obs[j].ra); 02090 u_rho[j].y = cos(obs[j].dec) * sin(obs[j].ra); 02091 u_rho[j].z = sin(obs[j].dec); 02092 02093 /* 02094 if (universe->GetReferenceSystem() == ECLIPTIC) { 02095 u_rho[j].rotate(0.0,-obleq(obs[j].date).GetRad(),0.0); 02096 } 02097 */ 02098 02099 /* 02100 if (universe->GetReferenceSystem() == ECLIPTIC) { 02101 EquatorialToEcliptic_J2000(u_rho[j]); 02102 } 02103 */ 02104 02105 if (universe->GetReferenceSystem() == ECLIPTIC) { 02106 // EquatorialToEcliptic_J2000(u_rho[j]); 02107 Angle obl = obleq_J2000(); 02108 u_rho[j].rotate(0.0,-obl.GetRad(),0.0); 02109 } 02110 02111 } 02112 02113 Vector rho[size]; 02114 Vector r[size]; 02115 Vector v; // v[0] 02116 02117 /* 02118 double tmp_rho = FromUnits(0.2,EARTHMOON); 02119 while (tmp_rho < FromUnits(50,AU)) { 02120 */ 02121 02122 /* 02123 double tmp_rho = FromUnits(0.5 ,AU); 02124 while (tmp_rho < FromUnits(0.55,AU)) { 02125 */ 02126 02127 double tmp_rho = FromUnits(0.535,AU); 02128 while (tmp_rho < FromUnits(0.536,AU)) { 02129 02130 cerr << "[******] tmp_rho: " << tmp_rho << endl; 02131 02132 for (j=0;j<size;++j) { 02133 rho[j] = u_rho[j]*tmp_rho; 02134 r[j] = R[j] + rho[j]; 02135 } 02136 02137 v = (r[1]-r[0])/(FromUnits(obs[0].date.GetJulian()-obs[1].date.GetJulian(),DAY)); 02138 02139 // least sq. v 02140 gauss_v(r[0],v,obs); 02141 02142 orbit.Compute(r[0],v,GetG()*GetMSun(),obs[0].date); 02143 // orbit.epoch = obs[0].date; 02144 02145 // debug 02146 cerr << " ---> orbit.a = " << orbit.a << endl; 02147 cerr << " ---> orbit.e = " << orbit.e << endl; 02148 cerr << " ---> orbit.i = " << orbit.i*(180/pi) << endl; 02149 02150 if (orbit.e < 1.0) preliminary_orbits.push_back(orbit); 02151 02152 // tmp_rho *= 1.3; 02153 tmp_rho += 0.005; 02154 } 02155 02156 }
Here is the call graph for this function:
Vector orsa::ComputeAcceleration | ( | const list< Body >::const_iterator | body_it, | |
const list< TreeNode >::const_iterator | node_domain_it, | |||
const bool | compute_quadrupole = true | |||
) |
Definition at line 96 of file orsa_interaction_tree.cc.
References Vector::IsZero(), Vector::LengthSquared(), secure_pow(), Vector::x, Vector::y, and Vector::z.
00096 { 00097 00098 Vector a; 00099 00100 if (node_domain_it->node_mass()==0.0) return a; 00101 00102 Vector d = node_domain_it->node_center_of_mass() - body_it->position(); 00103 00104 // monopole 00105 00106 const double l2 = d.LengthSquared(); 00107 00108 if (d.IsZero()) { 00109 cerr << "*** Warning: two objects in the same position! (" << l2 << ")" << endl; 00110 // continue; 00111 return a; 00112 } 00113 00114 a += d * secure_pow(l2,-1.5) * node_domain_it->node_mass(); 00115 00116 if (!compute_quadrupole) { 00117 return a; 00118 } 00119 00120 // quadrupole 00121 00122 double x[3]; 00123 00124 x[0] = d.x; 00125 x[1] = d.y; 00126 x[2] = d.z; 00127 00128 double coefficient = 0.0; 00129 unsigned int i,j; 00130 double c_node_quadrupole[3][3]; 00131 memcpy((void*)c_node_quadrupole, (const void*)node_domain_it->node_quadrupole(), 3*3*sizeof(double)); // works? 00132 for (i=0;i<3;++i) { 00133 for (j=0;j<3;++j) { 00134 coefficient += c_node_quadrupole[i][j]*x[i]*x[j]; 00135 } 00136 } 00137 00138 a += d * secure_pow(l2,-3.0) * coefficient; 00139 00140 return a; 00141 }
Here is the call graph for this function:
void orsa::convert | ( | UniverseType & | ut, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 154 of file orsa_universe.h.
References ORSA_ERROR, Real, and Simulated.
00154 { 00155 switch(i) { 00156 case 1: ut = Real; break; 00157 case 2: ut = Simulated; break; 00158 // 00159 default: 00160 ORSA_ERROR("conversion problem: i = %i",i); 00161 break; 00162 } 00163 }
void orsa::convert | ( | ReferenceSystem & | rs, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 591 of file orsa_units.h.
References ECLIPTIC, EQUATORIAL, and ORSA_ERROR.
00591 { 00592 switch(i) { 00593 case 1: rs = EQUATORIAL; break; 00594 case 2: rs = ECLIPTIC; break; 00595 // 00596 default: 00597 ORSA_ERROR("conversion problem: i = %i",i); 00598 break; 00599 } 00600 }
void orsa::convert | ( | TimeScale & | ts, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 263 of file orsa_units.h.
References GPS, ORSA_ERROR, TAI, TDT, UT, and UTC.
00263 { 00264 switch(i) { 00265 case 1: ts = UTC; break; 00266 case 2: ts = UT; break; 00267 case 3: ts = TAI; break; 00268 case 4: ts = TDT; break; 00269 case 5: ts = GPS; break; 00270 // 00271 default: 00272 ORSA_ERROR("conversion problem: i = %i",i); 00273 break; 00274 } 00275 }
void orsa::convert | ( | mass_unit & | mu, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 108 of file orsa_units.h.
References GRAM, KG, MEARTH, MJUPITER, MMOON, MSUN, and ORSA_ERROR.
00108 { 00109 switch(i) { 00110 case 1: mu = MSUN; break; 00111 case 2: mu = MJUPITER; break; 00112 case 3: mu = MEARTH; break; 00113 case 4: mu = MMOON; break; 00114 case 5: mu = KG; break; 00115 case 6: mu = GRAM; break; 00116 // 00117 default: 00118 ORSA_ERROR("conversion problem: i = %i",i); 00119 break; 00120 } 00121 }
void orsa::convert | ( | length_unit & | lu, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 79 of file orsa_units.h.
References AU, CM, EARTHMOON, KM, KPARSEC, LY, M, MPARSEC, ORSA_ERROR, PARSEC, REARTH, and RMOON.
00079 { 00080 switch(i) { 00081 case 1: lu = MPARSEC; break; 00082 case 2: lu = KPARSEC; break; 00083 case 3: lu = PARSEC; break; 00084 case 4: lu = LY; break; 00085 case 5: lu = AU; break; 00086 case 6: lu = EARTHMOON; break; 00087 case 7: lu = REARTH; break; 00088 case 8: lu = RMOON; break; 00089 case 9: lu = KM; break; 00090 case 10: lu = M; break; 00091 case 11: lu = CM; break; 00092 // 00093 default: 00094 ORSA_ERROR("conversion problem: i = %i",i); 00095 break; 00096 } 00097 }
void orsa::convert | ( | time_unit & | tu, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 47 of file orsa_units.h.
References DAY, HOUR, MINUTE, ORSA_ERROR, SECOND, and YEAR.
00047 { 00048 switch(i) { 00049 case 1: tu = YEAR; break; 00050 case 2: tu = DAY; break; 00051 case 3: tu = HOUR; break; 00052 case 4: tu = MINUTE; break; 00053 case 5: tu = SECOND; break; 00054 // 00055 default: 00056 ORSA_ERROR("conversion problem: i = %i",i); 00057 break; 00058 } 00059 }
void orsa::convert | ( | InteractionType & | it, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 57 of file orsa_interaction.h.
References ARMONICOSCILLATOR, GALACTIC_POTENTIAL_ALLEN, GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON, GRAVITATIONALTREE, JPL_PLANETS_NEWTON, NEWTON, NEWTON_MPI, ORSA_ERROR, and RELATIVISTIC.
00057 { 00058 switch(i) { 00059 case 1: it = NEWTON; break; 00060 case 2: it = ARMONICOSCILLATOR; break; 00061 case 3: it = GALACTIC_POTENTIAL_ALLEN; break; 00062 case 4: it = GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON; break; 00063 case 5: it = JPL_PLANETS_NEWTON; break; 00064 case 6: it = GRAVITATIONALTREE; break; 00065 case 7: it = NEWTON_MPI; break; 00066 case 8: it = RELATIVISTIC; break; 00067 // 00068 default: 00069 ORSA_ERROR("conversion problem: i = %i",i); 00070 break; 00071 } 00072 }
void orsa::convert | ( | IntegratorType & | it, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 45 of file orsa_integrator.h.
References BULIRSCHSTOER, DISSIPATIVERUNGEKUTTA, LEAPFROG, ORSA_ERROR, RA15, RUNGEKUTTA, and STOER.
00045 { 00046 switch(i) { 00047 case 1: it = STOER; break; 00048 case 2: it = BULIRSCHSTOER; break; 00049 case 3: it = RUNGEKUTTA; break; 00050 case 4: it = DISSIPATIVERUNGEKUTTA; break; 00051 case 5: it = RA15; break; 00052 case 6: it = LEAPFROG; break; 00053 // 00054 default: 00055 ORSA_ERROR("conversion problem: i = %i",i); 00056 break; 00057 } 00058 }
void orsa::convert | ( | JPL_planets & | jp, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 71 of file orsa_file_jpl.h.
References EARTH, EARTH_AND_MOON, EARTH_MOON_BARYCENTER, JUPITER, LIBRATIONS, MARS, MERCURY, MOON, NEPTUNE, NONE, NUTATIONS, ORSA_ERROR, PLUTO, SATURN, SOLAR_SYSTEM_BARYCENTER, SUN, URANUS, and VENUS.
00071 { 00072 switch(i) { 00073 case 0: jp = NONE; break; 00074 // 00075 case 1: jp = MERCURY; break; 00076 case 2: jp = VENUS; break; 00077 case 3: jp = EARTH; break; 00078 case 4: jp = MARS; break; 00079 case 5: jp = JUPITER; break; 00080 case 6: jp = SATURN; break; 00081 case 7: jp = URANUS; break; 00082 case 8: jp = NEPTUNE; break; 00083 case 9: jp = PLUTO; break; 00084 case 10: jp = MOON; break; 00085 case 11: jp = SUN; break; 00086 case 12: jp = SOLAR_SYSTEM_BARYCENTER; break; 00087 case 13: jp = EARTH_MOON_BARYCENTER; break; 00088 case 14: jp = NUTATIONS; break; 00089 case 15: jp = LIBRATIONS; break; 00090 // 00091 case 1000: jp = EARTH_AND_MOON; break; 00092 // 00093 default: 00094 ORSA_ERROR("conversion problem: i = %i",i); 00095 break; 00096 } 00097 }
void orsa::convert | ( | OrsaFileDataType & | ofdt, | |
const unsigned int | i | |||
) | [inline] |
Definition at line 389 of file orsa_file.h.
References OFDT_BODY, OFDT_END_OF_FILE, OFDT_EVOLUTION, OFDT_FRAME, OFDT_UNIVERSE, and ORSA_ERROR.
Referenced by OrsaFile::Read().
00389 { 00390 switch(i) { 00391 case 0: ofdt = OFDT_END_OF_FILE; break; 00392 case 1: ofdt = OFDT_UNIVERSE; break; 00393 case 2: ofdt = OFDT_EVOLUTION; break; 00394 case 3: ofdt = OFDT_FRAME; break; 00395 case 4: ofdt = OFDT_BODY; break; 00396 // 00397 default: 00398 ORSA_ERROR("conversion problem: i = %i",i); 00399 break; 00400 } 00401 }
double orsa::cos | ( | const Angle & | alpha | ) | [inline] |
Definition at line 567 of file orsa_units.h.
References Angle::GetRad().
Referenced by alpha_delta_meridian(), alpha_delta_meridian_moon(), apply_window(), Orbit::Compute(), Compute_Gauss(), Compute_TestMethod(), dQ(), Orbit::GetE(), phi(), phi_Hanning(), Orbit::RelativePosVel(), Vector::rotate(), and sincos().
00567 { 00568 return std::cos(alpha.GetRad()); 00569 }
Here is the call graph for this function:
Vector orsa::Cross | ( | const Vector & | u, | |
const Vector & | v | |||
) | [inline] |
Definition at line 182 of file orsa_coord.h.
References Vector::x, Vector::y, and Vector::z.
Referenced by Compute_Gauss().
00182 { 00183 return Vector (u.y*v.z-u.z*v.y, 00184 u.z*v.x-u.x*v.z, 00185 u.x*v.y-u.y*v.x); 00186 }
double orsa::delta_function | ( | const unsigned int | i, | |
const unsigned int | j | |||
) |
double orsa::dQ | ( | double | y | ) |
Definition at line 1319 of file orsa_fft.cc.
References cos(), pi, pisq, and sin().
01319 { 01320 01321 // y==0 01322 if (fabs(y) < 1.0e-12) return (0.0); 01323 01324 // y==PI 01325 if (fabs(y-pi) < 1.0e-12) return (-3.0/(4*pi)); 01326 01327 // y==-PI 01328 if (fabs(y+pi) < 1.0e-12) return ( 3.0/(4*pi)); 01329 01330 double ysq=y*y; 01331 return (pisq/(y*(pisq-ysq))*(cos(y)+(sin(y)/y)*((3*ysq-pisq)/(pisq-ysq)))); 01332 01333 }
Here is the call graph for this function:
double orsa::E1 | ( | void * | p | ) |
Definition at line 573 of file orsa_orbit_gsl.cc.
References RMS_residuals().
Referenced by alpha_delta_meridian_moon().
00573 { 00574 data *d = (data*) p; 00575 return (RMS_residuals(*d->obs,*d->orbit)); 00576 }
Here is the call graph for this function:
void EclipticToEquatorial | ( | Vector & | v, | |
const Date & | date | |||
) |
Definition at line 1965 of file orsa_units.cc.
References Angle::GetRad(), obleq(), and Vector::rotate().
01965 { 01966 v.rotate(0,obleq(date).GetRad(),0); 01967 }
Here is the call graph for this function:
void EclipticToEquatorial_J2000 | ( | Vector & | v | ) |
Definition at line 1973 of file orsa_units.cc.
References Angle::GetRad(), obleq_J2000(), and Vector::rotate().
Referenced by Sky::Compute_J2000(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), and AstorbFile::Read().
01973 { 01974 v.rotate(0,obleq_J2000().GetRad(),0); 01975 }
Here is the call graph for this function:
void EquatorialToEcliptic | ( | Vector & | v, | |
const Date & | date | |||
) |
Definition at line 1969 of file orsa_units.cc.
References Angle::GetRad(), obleq(), and Vector::rotate().
01969 { 01970 v.rotate(0,-obleq(date).GetRad(),0); 01971 }
Here is the call graph for this function:
void EquatorialToEcliptic_J2000 | ( | Vector & | v | ) |
Definition at line 1977 of file orsa_units.cc.
References Angle::GetRad(), obleq_J2000(), and Vector::rotate().
01977 { 01978 v.rotate(0,-obleq_J2000().GetRad(),0); 01979 }
Here is the call graph for this function:
Vector orsa::ExternalProduct | ( | const Vector & | u, | |
const Vector & | v | |||
) | [inline] |
Definition at line 176 of file orsa_coord.h.
References Vector::x, Vector::y, and Vector::z.
Referenced by AngularMomentum(), Orbit::Compute(), and SearchCloseApproaches().
00176 { 00177 return Vector (u.y*v.z-u.z*v.y, 00178 u.z*v.x-u.x*v.z, 00179 u.x*v.y-u.y*v.x); 00180 }
double FromUnits | ( | const | double, | |
const | mass_unit, | |||
const | int = 1 | |||
) |
Definition at line 137 of file orsa_universe.cc.
References Units::FromUnits(), and units.
Referenced by BuildObservationTriplets(), Compute(), Compute_Gauss(), Compute_TestMethod(), poly_8_gsl_solve(), and Weight().
00137 { return orsa::units->FromUnits(x,u,power); };
Here is the call graph for this function:
double FromUnits | ( | const | double, | |
const | length_unit, | |||
const | int = 1 | |||
) |
Definition at line 136 of file orsa_universe.cc.
References Units::FromUnits(), and units.
00136 { return orsa::units->FromUnits(x,u,power); };
Here is the call graph for this function:
double FromUnits | ( | const | double, | |
const | time_unit, | |||
const | int = 1 | |||
) |
Definition at line 135 of file orsa_universe.cc.
References Units::FromUnits(), and units.
Referenced by Evolution::Evolution(), GalacticPotentialAllen::GalacticPotentialAllen(), TimeStep::GetDouble(), JPLFile::GetEph(), JPLFile::GetMass(), Date::GetTime(), Evolution::Integrate(), TimeStep::operator *=(), OptimizedOrbitPositions::PropagatedOrbit(), PropagatedSky_J2000(), radius(), Mercury5IntegrationFile::Read(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), RadauModIntegrationFile::Read(), LocationFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), UniverseTypeAwareTime::SetTime(), Date::SetTime(), StartFrame(), SWIFTRawReadBinaryFile(), Date::Time(), and TimeStep::TimeStep().
00135 { return orsa::units->FromUnits(x,u,power); };
Here is the call graph for this function:
void orsa::gauss_v | ( | const Vector & | r, | |
Vector & | v, | |||
const vector< Observation > & | obs | |||
) |
Definition at line 1320 of file orsa_orbit_gsl.cc.
References gauss_v_df(), gauss_v_f(), gauss_v_fdf(), ORSA_ERROR, Vector::x, Vector::y, and Vector::z.
Referenced by Compute_TestMethod().
01320 { 01321 01322 if (obs.size() < 2) { 01323 ORSA_ERROR("gauss_v(...) needs at least two observations..."); 01324 return; 01325 } 01326 01327 // cerr << "..gauss_v..." << endl; 01328 01329 g_v_class g_v; 01330 01331 g_v.obs = obs; 01332 g_v.r = r; 01333 01334 gsl_multifit_fdfsolver *s = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder,obs.size(),3); 01335 01336 gsl_multifit_function_fdf gauss_v_function; 01337 01338 gauss_v_function.f = &gauss_v_f; 01339 gauss_v_function.df = &gauss_v_df; 01340 gauss_v_function.fdf = &gauss_v_fdf; 01341 gauss_v_function.n = obs.size(); 01342 gauss_v_function.p = 3; 01343 gauss_v_function.params = &g_v; 01344 01345 gsl_vector *x = gsl_vector_alloc(3); 01346 01347 cerr << "..inside, Length(r): " << (r).Length() << endl; 01348 01349 gsl_vector_set(x,0,v.x); 01350 gsl_vector_set(x,1,v.y); 01351 gsl_vector_set(x,2,v.z); 01352 01353 gsl_multifit_fdfsolver_set(s,&gauss_v_function,x); 01354 01355 size_t iter = 0; 01356 int status; 01357 do { 01358 01359 ++iter; 01360 01361 status = gsl_multifit_fdfsolver_iterate (s); 01362 01363 printf ("itaration status = %s\n", gsl_strerror (status)); 01364 01365 /* 01366 printf ("iter: %3u %15.8f %15.8f %15.8f |f(x)| = %g\n", 01367 iter, 01368 gsl_vector_get (s->x, 0), 01369 gsl_vector_get (s->x, 1), 01370 gsl_vector_get (s->x, 2), 01371 gsl_blas_dnrm2 (s->f)); 01372 */ 01373 01374 // status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1.0e-3); 01375 status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 5.0e-2); 01376 01377 printf ("convergence status = %s\n", gsl_strerror (status)); 01378 01379 } while ((status == GSL_CONTINUE) && (iter < 10)); 01380 01381 printf ("status = %s\n", gsl_strerror (status)); 01382 01383 // save results 01384 v.x = gsl_vector_get (s->x,0); 01385 v.y = gsl_vector_get (s->x,1); 01386 v.z = gsl_vector_get (s->x,2); 01387 01388 // free 01389 gsl_multifit_fdfsolver_free (s); 01390 gsl_vector_free (x); 01391 }
Here is the call graph for this function:
int orsa::gauss_v_df | ( | const gsl_vector * | v, | |
void * | params, | |||
gsl_matrix * | J | |||
) |
Definition at line 1268 of file orsa_orbit_gsl.cc.
References OrbitWithEpoch::Compute(), OrbitWithEpoch::epoch, gauss_v_diff_f(), GetG(), GetMSun(), Vector::x, Vector::y, and Vector::z.
Referenced by gauss_v(), and gauss_v_fdf().
01268 { 01269 01270 Vector velocity; 01271 // 01272 velocity.x = gsl_vector_get(v,0); 01273 velocity.y = gsl_vector_get(v,1); 01274 velocity.z = gsl_vector_get(v,2); 01275 01276 g_v_class *par = (g_v_class*) params; 01277 01278 vector<Observation> &obs(par->obs); 01279 01280 OrbitWithEpoch orbit; 01281 orbit.Compute(par->r,velocity,GetG()*GetMSun(),obs[0].date); 01282 // orbit.epoch = obs[1].date; 01283 // orbit.epoch = obs[0].date; 01284 01285 gauss_v_diff_par_class diff_par; 01286 01287 diff_par.r = par->r; 01288 diff_par.velocity = velocity; 01289 // diff_par.orbit_epoch = obs[1].date; 01290 diff_par.orbit_epoch = orbit.epoch; 01291 01292 gsl_function diff_f; 01293 double result, abserr; 01294 01295 diff_f.function = &gauss_v_diff_f; 01296 diff_f.params = &diff_par; 01297 01298 size_t i,j; 01299 for (i=0;i<obs.size();++i) { 01300 diff_par.obs = obs[i]; 01301 for (j=0;j<3;++j) { 01302 diff_par.var_index = j; 01303 // gsl_diff_central (&diff_f, gsl_vector_get(v,j), &result, &abserr); 01304 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), 1.0e-9, &result, &abserr); 01305 gsl_matrix_set (J, i, j, result); 01306 // cerr << "diff: " << result << " +/- " << abserr << endl; 01307 } 01308 } 01309 01310 return GSL_SUCCESS; 01311 }
Here is the call graph for this function:
double orsa::gauss_v_diff_f | ( | double | x, | |
void * | params | |||
) |
Definition at line 1245 of file orsa_orbit_gsl.cc.
References OrbitWithEpoch::Compute(), Sky::delta_arcsec(), GetG(), GetMSun(), PropagatedSky_J2000(), Vector::x, Vector::y, and Vector::z.
Referenced by gauss_v_df().
01245 { 01246 01247 // cerr << "gauss_v_diff_f()..." << endl; 01248 01249 gauss_v_diff_par_class *diff_par = (gauss_v_diff_par_class*) params; 01250 01251 // cerr << "gauss_v_diff_f()... var_index = " << diff_par->var_index << endl; 01252 01253 Vector velocity(diff_par->velocity); 01254 01255 switch(diff_par->var_index) { 01256 case 0: velocity.x = x; break; 01257 case 1: velocity.y = x; break; 01258 case 2: velocity.z = x; break; 01259 } 01260 01261 OrbitWithEpoch orbit; 01262 orbit.Compute(diff_par->r,velocity,GetG()*GetMSun(),diff_par->orbit_epoch); 01263 // orbit.epoch = diff_par->orbit_epoch; 01264 01265 return (PropagatedSky_J2000(orbit,diff_par->obs.date,diff_par->obs.obscode).delta_arcsec(diff_par->obs)); 01266 }
Here is the call graph for this function:
int orsa::gauss_v_f | ( | const gsl_vector * | v, | |
void * | params, | |||
gsl_vector * | f | |||
) |
Definition at line 1189 of file orsa_orbit_gsl.cc.
References Orbit::a, OrbitWithEpoch::Compute(), Sky::delta_arcsec(), Orbit::e, GetG(), GetMSun(), Orbit::i, ORSA_DEBUG, pi, OptimizedOrbitPositions::PropagatedSky_J2000(), Vector::x, Vector::y, and Vector::z.
Referenced by gauss_v(), and gauss_v_fdf().
01189 { 01190 01191 Vector velocity; 01192 // 01193 velocity.x = gsl_vector_get(v,0); 01194 velocity.y = gsl_vector_get(v,1); 01195 velocity.z = gsl_vector_get(v,2); 01196 01197 g_v_class *par = (g_v_class*) params; 01198 01199 vector<Observation> &obs(par->obs); 01200 01201 // cerr << "Length(r): " << (par->r).Length() << endl; 01202 // cerr << "Length(v): " << (velocity).Length() << endl; 01203 01204 OrbitWithEpoch orbit; 01205 orbit.Compute(par->r,velocity,GetG()*GetMSun(),obs[0].date); 01206 // orbit.epoch = obs[1].date; 01207 // orbit.epoch = obs[0].date; 01208 01209 { 01210 // debug 01211 /* cerr << "gauss_v_f() ---> orbit.a = " << orbit.a << endl; 01212 cerr << "gauss_v_f() ---> orbit.e = " << orbit.e << endl; 01213 cerr << "gauss_v_f() ---> orbit.i = " << orbit.i*(180/pi) << endl; 01214 */ 01215 ORSA_DEBUG( 01216 "gauss_v_f():\n" 01217 "orbit.a = %g\n" 01218 "orbit.e = %g\n" 01219 "orbit.i = %g\n", 01220 orbit.a,orbit.e,orbit.i*(180/pi)); 01221 } 01222 01223 OptimizedOrbitPositions opt(orbit); 01224 01225 size_t i; 01226 double arcsec; 01227 for (i=0;i<obs.size();++i) { 01228 arcsec = opt.PropagatedSky_J2000(obs[i].date,obs[i].obscode).delta_arcsec(obs[i]); 01229 gsl_vector_set(f,i,arcsec); 01230 // cerr << "f[" << i << "] = " << arcsec << endl; 01231 } 01232 01233 return GSL_SUCCESS; 01234 }
Here is the call graph for this function:
int orsa::gauss_v_fdf | ( | const gsl_vector * | v, | |
void * | params, | |||
gsl_vector * | f, | |||
gsl_matrix * | J | |||
) |
Definition at line 1313 of file orsa_orbit_gsl.cc.
References gauss_v_df(), and gauss_v_f().
Referenced by gauss_v().
01313 { 01314 gauss_v_f (v,params,f); 01315 gauss_v_df(v,params,J); 01316 return GSL_SUCCESS; 01317 }
Here is the call graph for this function:
double GetC | ( | ) |
Definition at line 127 of file orsa_universe.cc.
References units.
Referenced by Compute_Gauss(), modified_mu(), and PropagatedSky_J2000().
00127 { return (orsa::units->GetC()); }
double GetG | ( | ) |
Definition at line 124 of file orsa_universe.cc.
References units.
Referenced by Orbit::Compute(), Compute_Gauss(), Compute_TestMethod(), GalacticPotentialAllen::GalacticPotentialAllen(), gauss_v_df(), gauss_v_diff_f(), gauss_v_f(), JPLFile::GetMass(), GravitationalTree::GravitationalTree(), JPLPlanetsNewton::JPLPlanetsNewton(), least_sq_df(), least_sq_f(), OptimizedOrbitPositions::PropagatedOrbit(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), and AstorbFile::Read().
00124 { return (orsa::units->GetG()); }
double GetG_MKS | ( | ) |
Definition at line 125 of file orsa_universe.cc.
References units.
Referenced by JPLFile::GetMEarth_MKS(), JPLFile::GetMJupiter_MKS(), JPLFile::GetMMoon_MKS(), and JPLFile::GetMSun_MKS().
00125 { return (orsa::units->GetG_MKS()); }
double GetMSun | ( | ) |
Definition at line 126 of file orsa_universe.cc.
References units.
Referenced by Compute_Gauss(), Compute_TestMethod(), gauss_v_df(), gauss_v_diff_f(), gauss_v_f(), least_sq_df(), least_sq_f(), OptimizedOrbitPositions::PropagatedOrbit(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), and AstorbFile::Read().
00126 { return (orsa::units->GetMSun()); }
Angle gmst | ( | const Date & | date | ) |
Definition at line 1896 of file orsa_units.cc.
References Date::GetDayFraction(), Date::GetJulian(), Angle::SetHMS(), and UT.
01896 { 01897 // tested using xephem, very good agreement! 01898 // T in centuries from JD 2451545.0 UT1=UT 01899 Date d = date; 01900 // d.ConvertToTimeScale(UT); // UT or UTC ?? 01901 // const double T = (d.GetJulian() - 2451545.0)/36525.0; 01902 const double T = (d.GetJulian(UT) - 2451545.0)/36525.0; 01903 Angle a; 01904 a.SetHMS(d.GetDayFraction()*24+6,41,50.54841+((-0.0000062*T+0.093104)*T+8640184.812866)*T); 01905 return a; 01906 }
Here is the call graph for this function:
void orsa::Interpolate | ( | const std::vector< VectorWithParameter > | v_in, | |
const double | x, | |||
Vector & | v_out, | |||
Vector & | err_v_out | |||
) |
void orsa::Interpolate | ( | const vector< VectorWithParameter > | vx_in, | |
const double | x, | |||
Vector & | v_out, | |||
Vector & | err_v_out | |||
) |
Definition at line 102 of file orsa_coord.cc.
References Vector::Set().
00102 { 00103 00104 // cerr << "Interpolate() call..." << endl; 00105 00106 unsigned int i,j,m; 00107 00108 unsigned int i_closest; 00109 00110 double diff,denom; 00111 double ho,hp; 00112 double tmp_double; 00113 00114 unsigned int n_points = vx_in.size(); 00115 00116 /* cerr << " -----> n_points: " << n_points << endl; 00117 for (j=0;j<n_points;j++) 00118 std::cerr << "vx_in[" << j << "].par = " << vx_in[j].par << "\n"; 00119 */ 00120 00121 if (n_points < 2) { 00122 cerr << "too few points..." << endl; 00123 for (j=0;j<n_points;j++) { 00124 v_out = vx_in[0]; 00125 err_v_out.Set(0,0,0); 00126 } 00127 return; 00128 } 00129 00130 Vector * c = new Vector[n_points]; 00131 Vector * d = new Vector[n_points]; 00132 00133 Vector w; 00134 00135 diff = fabs(x-vx_in[0].par); 00136 i_closest = 0; 00137 00138 for (j=0;j<n_points;j++) { 00139 00140 tmp_double = fabs(x-vx_in[j].par); 00141 if (tmp_double < diff) { 00142 diff = tmp_double; 00143 i_closest = j; 00144 } 00145 00146 c[j] = vx_in[j]; 00147 d[j] = vx_in[j]; 00148 } 00149 00150 v_out = vx_in[i_closest]; 00151 err_v_out = vx_in[i_closest]; 00152 00153 // cerr << "local.. v: " << Length(v_out) << " i_closest: " << i_closest << endl; 00154 00155 i_closest--; 00156 00157 for (m=1;m<=(n_points-2);m++) { 00158 for (i=0;i<(n_points-m);i++) { 00159 ho = vx_in[i].par - x; 00160 hp = vx_in[i+m].par - x; 00161 00162 denom = ho - hp; 00163 00164 if (denom == 0.0) { 00165 cerr << "interpolate() --> Error: divide by zero\n"; 00166 cerr << "i: " << i << " m: " << m << endl; 00167 for (j=0;j<n_points;j++) 00168 cerr << "vx_in[" << j << "].par = " << vx_in[j].par << "\n"; 00169 exit(0); 00170 } 00171 00172 w = c[i+1] - d[i]; 00173 00174 d[i] = hp*w/denom; 00175 c[i] = ho*w/denom; 00176 00177 // cerr << " ++++++++++ d[" << i << "]: " << Length(d[i]) << endl; 00178 // cerr << " ++++++++++ c[" << i << "]: " << Length(c[i]) << endl; 00179 } 00180 00181 if ( (2*i_closest) < (n_points-m) ) { 00182 err_v_out = c[i_closest+1]; 00183 } else { 00184 err_v_out = d[i_closest]; 00185 i_closest--; 00186 } 00187 00188 v_out += err_v_out; 00189 } 00190 00191 delete [] c; 00192 delete [] d; 00193 }
Here is the call graph for this function:
void orsa::Interpolate | ( | const std::vector< BodyWithParameter > & | b_in, | |
const double | x, | |||
Body & | b_out, | |||
Body & | err_b_out | |||
) |
void orsa::Interpolate | ( | const vector< BodyWithParameter > & | b_in, | |
const double | x, | |||
Body & | b_out, | |||
Body & | err_b_out | |||
) |
Definition at line 442 of file orsa_body.cc.
References Body::position(), Body::SetPosition(), Body::SetVelocity(), and Body::velocity().
00442 { 00443 00444 unsigned int n_points = b_in.size(); 00445 00446 Vector p_interpolated, err_p_interpolated; 00447 Vector v_interpolated, err_v_interpolated; 00448 00449 unsigned int j; 00450 00451 vector<VectorWithParameter> pp(n_points), vv(n_points); 00452 00453 for(j=0;j<n_points;j++) { 00454 pp[j].Set(b_in[j].position()); 00455 vv[j].Set(b_in[j].velocity()); 00456 pp[j].par = vv[j].par = b_in[j].par; 00457 } 00458 00459 Interpolate(pp,x,p_interpolated,err_p_interpolated); 00460 Interpolate(vv,x,v_interpolated,err_v_interpolated); 00461 00462 b_out = b_in[0]; 00463 b_out.SetPosition(p_interpolated); 00464 b_out.SetVelocity(v_interpolated); 00465 00466 err_b_out = b_in[0]; 00467 err_b_out.SetPosition(err_p_interpolated); 00468 err_b_out.SetVelocity(err_v_interpolated); 00469 00470 }
Here is the call graph for this function:
std::string JPL_planet_name | ( | const JPL_planets | p | ) |
Definition at line 480 of file orsa_file_jpl.cc.
References EARTH, EARTH_MOON_BARYCENTER, JUPITER, MARS, MERCURY, MOON, NEPTUNE, PLUTO, SATURN, SUN, URANUS, and VENUS.
00480 { 00481 00482 string name; 00483 00484 switch(p) { 00485 case SUN: name = "Sun"; break; 00486 case MERCURY: name = "Mercury"; break; 00487 case VENUS: name = "Venus"; break; 00488 case EARTH: name = "Earth"; break; 00489 case MARS: name = "Mars"; break; 00490 case JUPITER: name = "Jupiter"; break; 00491 case SATURN: name = "Saturn"; break; 00492 case URANUS: name = "Uranus"; break; 00493 case NEPTUNE: name = "Neptune"; break; 00494 case PLUTO: name = "Pluto"; break; 00495 // 00496 case MOON: name = "Moon"; break; 00497 case EARTH_MOON_BARYCENTER: name = "Earth-Moon barycenter"; break; 00498 // 00499 default: break; 00500 } 00501 00502 return (name); 00503 }
double orsa::KineticEnergy | ( | const std::vector< Body > & | ) |
double orsa::KineticEnergy | ( | const vector< Body > & | f | ) |
Definition at line 207 of file orsa_frame.cc.
References KineticEnergy().
00207 { 00208 if (f.size()==0) return (0.0); 00209 double energy=0.0; 00210 unsigned int k=0; 00211 while (k<f.size()) { 00212 energy += f[k].KineticEnergy(); 00213 k++; 00214 } 00215 return (energy); 00216 }
Here is the call graph for this function:
double orsa::KineticEnergy | ( | const Body & | ) |
Referenced by KineticEnergy().
std::string label | ( | const InteractionType | it | ) |
Definition at line 37 of file orsa_interaction.cc.
References ARMONICOSCILLATOR, GALACTIC_POTENTIAL_ALLEN, GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON, GRAVITATIONALTREE, JPL_PLANETS_NEWTON, NEWTON, NEWTON_MPI, and RELATIVISTIC.
00037 { 00038 00039 std::string s(""); 00040 00041 switch (it) { 00042 case NEWTON: s = "Newton"; break; 00043 case NEWTON_MPI: s = "Newton (MPI)"; break; 00044 case ARMONICOSCILLATOR: s = "Armonic Oscillator"; break; 00045 case GALACTIC_POTENTIAL_ALLEN: s = "Galactic Potential (Allen)"; break; 00046 case GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON: s = "Galactic Potential (Allen) + Newton"; break; 00047 case JPL_PLANETS_NEWTON: s = "JPL planets + Newton"; break; 00048 case GRAVITATIONALTREE: s = "Gravitational TreeCode"; break; 00049 case RELATIVISTIC: s = "Newton + Relativistic effects"; break; 00050 } 00051 00052 return s; 00053 }
std::string orsa::label | ( | const IntegratorType | it | ) | [inline] |
Definition at line 60 of file orsa_integrator.h.
References BULIRSCHSTOER, DISSIPATIVERUNGEKUTTA, LEAPFROG, RA15, RUNGEKUTTA, and STOER.
Referenced by Label(), Mercury5IntegrationFile::Read(), and SWIFTFile::Read().
00060 { 00061 std::string s; 00062 switch (it) { 00063 case STOER: s="Stoer"; break; 00064 case BULIRSCHSTOER: s="Bulirsch-Stoer"; break; 00065 case RUNGEKUTTA: s="Runge-Kutta 4th order"; break; 00066 case DISSIPATIVERUNGEKUTTA: s="Dissipative Runge-Kutta 4th order"; break; 00067 case RA15: s="Everhart's RADAU 15th order"; break; 00068 case LEAPFROG: s="Leapfrog 2nd order"; break; 00069 } 00070 return s; 00071 }
std::string Label | ( | const ConfigEnum | e | ) |
Definition at line 32 of file orsa_config.cc.
References ASTDYS_ALLNUM_CAT, ASTDYS_ALLNUM_CTC, ASTDYS_ALLNUM_CTM, ASTDYS_UFITOBS_CAT, ASTDYS_UFITOBS_CTC, ASTDYS_UFITOBS_CTM, JPL_DASTCOM_COMET, JPL_DASTCOM_NUM, JPL_DASTCOM_UNNUM, JPL_EPHEM_FILE, label(), LOWELL_ASTORB, MPC_COMET, MPC_DAILY, MPC_DISTANT, MPC_MPCORB, MPC_NEA, MPC_PHA, MPC_UNUSUALS, NEODYS_CAT, NEODYS_CTC, NO_CONFIG_ENUM, OBSCODE, TEXTURE_EARTH, TEXTURE_JUPITER, TEXTURE_MARS, TEXTURE_MERCURY, TEXTURE_MOON, TEXTURE_NEPTUNE, TEXTURE_PLUTO, TEXTURE_SATURN, TEXTURE_SUN, TEXTURE_URANUS, TEXTURE_VENUS, TLE_GEO, TLE_GPS, TLE_ISS, TLE_KEPELE, TLE_NASA, TLE_VISUAL, and TLE_WEATHER.
00032 { 00033 string label; 00034 switch(e) { 00035 case JPL_EPHEM_FILE: label="JPL ephemeris"; break; 00036 case JPL_DASTCOM_NUM: label="JPL asteroids database (NUM)"; break; 00037 case JPL_DASTCOM_UNNUM: label="JPL asteroids database (UNNUM)"; break; 00038 case JPL_DASTCOM_COMET: label="JPL comets database"; break; 00039 case LOWELL_ASTORB: label="Lowell asteroids database"; break; 00040 case MPC_MPCORB: label="MPC asteroids database"; break; 00041 case MPC_COMET: label="MPC comets database"; break; 00042 case MPC_NEA: label="MPC asteroids database (NEA)"; break; 00043 case MPC_DAILY: label="MPC asteroids database (DAILY)"; break; 00044 case MPC_DISTANT: label="MPC asteroids database (DISTANT)"; break; 00045 case MPC_PHA: label="MPC asteroids database (PHA)"; break; 00046 case MPC_UNUSUALS: label="MPC asteroids database (UNUSUALS)"; break; 00047 case ASTDYS_ALLNUM_CAT: label="AstDyS asteroids database (CAT)"; break; 00048 case ASTDYS_ALLNUM_CTC: label="AstDyS asteroids database (CTC)"; break; 00049 case ASTDYS_ALLNUM_CTM: label="AstDyS asteroids database (CTM)"; break; 00050 case ASTDYS_UFITOBS_CAT: label="AstDyS unnumbered asteroids database (CAT)"; break; 00051 case ASTDYS_UFITOBS_CTC: label="AstDyS unnumbered asteroids database (CTC)"; break; 00052 case ASTDYS_UFITOBS_CTM: label="AstDyS unnumbered asteroids database (CTM)"; break; 00053 case NEODYS_CAT: label="NEODyS asteroids database (CAT)"; break; 00054 case NEODYS_CTC: label="NEODyS asteroids database (CTC)"; break; 00055 case OBSCODE: label="Observatory codes"; break; 00056 // TLE 00057 case TLE_NASA: label="TLE (NASA)"; break; 00058 case TLE_GEO: label="TLE (GEO)"; break; 00059 case TLE_GPS: label="TLE (GPS)"; break; 00060 case TLE_ISS: label="TLE (ISS)"; break; 00061 case TLE_KEPELE: label="TLE (KEPELE)"; break; 00062 case TLE_VISUAL: label="TLE (VISUAL)"; break; 00063 case TLE_WEATHER: label="TLE (WEATHER)"; break; 00064 // textures 00065 case TEXTURE_SUN : label="Sun's texture"; break; 00066 case TEXTURE_MERCURY: label="Mercury's texture"; break; 00067 case TEXTURE_VENUS: label="Venus's texture"; break; 00068 case TEXTURE_EARTH: label="Earth's texture"; break; 00069 case TEXTURE_MOON: label="Moon's texture"; break; 00070 case TEXTURE_MARS: label="Mars's texture"; break; 00071 case TEXTURE_JUPITER: label="Jupiter's texture"; break; 00072 case TEXTURE_SATURN: label="Saturn's texture"; break; 00073 case TEXTURE_URANUS: label="Uranus's texture"; break; 00074 case TEXTURE_NEPTUNE: label="Neptune's texture"; break; 00075 case TEXTURE_PLUTO: label="Pluto's texture"; break; 00076 // 00077 case NO_CONFIG_ENUM: label="This shuld not be used!"; break; 00078 } 00079 return label; 00080 }
Here is the call graph for this function:
int orsa::least_sq_df | ( | const gsl_vector * | v, | |
void * | params, | |||
gsl_matrix * | J | |||
) |
Definition at line 892 of file orsa_orbit_gsl.cc.
References Orbit::a, Orbit::e, OrbitWithEpoch::epoch, GetG(), GetMSun(), Orbit::i, least_sq_diff_f(), Orbit::M, Orbit::mu, Orbit::omega_node, and Orbit::omega_pericenter.
Referenced by least_sq_fdf(), and OrbitDifferentialCorrectionsLeastSquares().
00892 { 00893 00894 OrbitWithEpoch orbit; 00895 00896 orbit.a = gsl_vector_get(v,0); 00897 orbit.e = gsl_vector_get(v,1); 00898 orbit.i = gsl_vector_get(v,2); 00899 orbit.omega_node = gsl_vector_get(v,3); 00900 orbit.omega_pericenter = gsl_vector_get(v,4); 00901 orbit.M = gsl_vector_get(v,5); 00902 00903 orbit.mu = GetG()*GetMSun(); 00904 00905 par_class * par = (par_class*) params; 00906 00907 vector<Observation> * obs = par->obs; 00908 00909 orbit.epoch = par->orbit_epoch; 00910 00911 least_sq_diff_par_class diff_par; 00912 00913 diff_par.orbit = orbit; 00914 00915 gsl_function diff_f; 00916 double result, abserr; 00917 00918 diff_f.function = &least_sq_diff_f; 00919 diff_f.params = &diff_par; 00920 00921 size_t i,j; 00922 for (i=0;i<obs->size();++i) { 00923 diff_par.obs = (*obs)[i]; 00924 for (j=0;j<6;++j) { 00925 diff_par.var_index = j; 00926 00927 /* 00928 // gsl_diff_central(&diff_f, gsl_vector_get(v,j), &result, &abserr); 00929 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), 1.0e-6, &result, &abserr); 00930 gsl_matrix_set (J, i, j, result); 00931 // 00932 fprintf(stderr,"[-6] diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 00933 i,j,result,abserr,gsl_vector_get(v,j)); 00934 00935 00936 00937 // gsl_diff_central(&diff_f, gsl_vector_get(v,j), &result, &abserr); 00938 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), 1.0e-8, &result, &abserr); 00939 gsl_matrix_set (J, i, j, result); 00940 // 00941 fprintf(stderr,"[-8] diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 00942 i,j,result,abserr,gsl_vector_get(v,j)); 00943 00944 00945 00946 // gsl_diff_central(&diff_f, gsl_vector_get(v,j), &result, &abserr); 00947 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), 1.0e-9, &result, &abserr); 00948 gsl_matrix_set (J, i, j, result); 00949 // 00950 fprintf(stderr,"[-9] diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 00951 i,j,result,abserr,gsl_vector_get(v,j)); 00952 00953 // gsl_diff_central(&diff_f, gsl_vector_get(v,j), &result, &abserr); 00954 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), 1.0e-12, &result, &abserr); 00955 gsl_matrix_set (J, i, j, result); 00956 // 00957 fprintf(stderr,"[-12]diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 00958 i,j,result,abserr,gsl_vector_get(v,j)); 00959 */ 00960 00961 /* 00962 // gsl_diff_central(&diff_f, gsl_vector_get(v,j), &result, &abserr); 00963 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), 1.0e-13, &result, &abserr); 00964 gsl_matrix_set (J, i, j, result); 00965 // 00966 fprintf(stderr,"[-13]diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 00967 i,j,result,abserr,gsl_vector_get(v,j)); 00968 00969 00970 // gsl_diff_central(&diff_f, gsl_vector_get(v,j), &result, &abserr); 00971 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), 1.0e-14, &result, &abserr); 00972 gsl_matrix_set (J, i, j, result); 00973 // 00974 fprintf(stderr,"[-14]diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 00975 i,j,result,abserr,gsl_vector_get(v,j)); 00976 */ 00977 00978 /* 00979 // a smarter one? (with numeric limits) 00980 // gsl_diff_central(&diff_f, gsl_vector_get(v,j), &result, &abserr); 00981 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), gsl_vector_get(v,j)*100.0*std::numeric_limits<double>::epsilon(), &result, &abserr); 00982 gsl_matrix_set (J, i, j, result); 00983 // 00984 fprintf(stderr,"[lim]diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 00985 i,j,result,abserr,gsl_vector_get(v,j)); 00986 */ 00987 00988 /* 00989 if (par->use_covar) { 00990 00991 cerr << "cov: j=" << j << " cov: " << sqrt(gsl_matrix_get(par->covar,j,j)) << endl; 00992 00993 // again, but using the covar matrix to have an estimate of the step 00994 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), 0.001*sqrt(gsl_matrix_get(par->covar,j,j)), &result, &abserr); 00995 gsl_matrix_set (J, i, j, result); 00996 // 00997 fprintf(stderr,"[cov]diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 00998 i,j,result,abserr,gsl_vector_get(v,j)); 00999 } else { 01000 */ 01001 01002 gsl_deriv_central(&diff_f, gsl_vector_get(v,j), gsl_vector_get(v,j)*1.0e4*std::numeric_limits<double>::epsilon(), &result, &abserr); 01003 gsl_matrix_set (J, i, j, result); 01004 // 01005 fprintf(stderr,"[lim]diff[%03i][%i] = %20f +/- %20f at %20.12f\n", 01006 i,j,result,abserr,gsl_vector_get(v,j)); 01007 01008 // } 01009 01010 } 01011 } 01012 01013 return GSL_SUCCESS; 01014 }
Here is the call graph for this function:
double orsa::least_sq_diff_f | ( | double | x, | |
void * | params | |||
) |
Definition at line 870 of file orsa_orbit_gsl.cc.
References Orbit::a, Sky::delta_arcsec(), Orbit::e, Orbit::i, Orbit::M, Orbit::omega_node, Orbit::omega_pericenter, and PropagatedSky_J2000().
Referenced by least_sq_df().
00870 { 00871 00872 least_sq_diff_par_class * diff_par = (least_sq_diff_par_class *) params; 00873 00874 OrbitWithEpoch orbit(diff_par->orbit); 00875 00876 switch(diff_par->var_index) { 00877 case 0: orbit.a = x; break; 00878 case 1: orbit.e = x; break; 00879 case 2: orbit.i = x; break; 00880 case 3: orbit.omega_node = x; break; 00881 case 4: orbit.omega_pericenter = x; break; 00882 case 5: orbit.M = x; break; 00883 } 00884 00885 /* 00886 ORSA_ERROR("least_sq_diff_f(...): diff_par->var_index = %i x = %26.20f",diff_par->var_index,x); 00887 */ 00888 00889 return (PropagatedSky_J2000(orbit,diff_par->obs.date,diff_par->obs.obscode).delta_arcsec(diff_par->obs)); 00890 }
Here is the call graph for this function:
int orsa::least_sq_f | ( | const gsl_vector * | v, | |
void * | params, | |||
gsl_vector * | f | |||
) |
Definition at line 828 of file orsa_orbit_gsl.cc.
References Orbit::a, Sky::delta_arcsec(), Orbit::e, OrbitWithEpoch::epoch, GetG(), GetMSun(), Orbit::i, Orbit::M, Orbit::mu, Orbit::omega_node, Orbit::omega_pericenter, ORSA_DEBUG, ORSA_ERROR, pi, and OptimizedOrbitPositions::PropagatedSky_J2000().
Referenced by least_sq_fdf(), and OrbitDifferentialCorrectionsLeastSquares().
00828 { 00829 00830 OrbitWithEpoch orbit; 00831 00832 orbit.a = gsl_vector_get(v,0); 00833 orbit.e = gsl_vector_get(v,1); 00834 orbit.i = gsl_vector_get(v,2); 00835 orbit.omega_node = gsl_vector_get(v,3); 00836 orbit.omega_pericenter = gsl_vector_get(v,4); 00837 orbit.M = gsl_vector_get(v,5); 00838 00839 { 00840 ORSA_DEBUG( 00841 "least_sq_f():\n" 00842 "orbit.a = %g\n" 00843 "orbit.e = %g\n" 00844 "orbit.i = %g\n", 00845 orbit.a,orbit.e,orbit.i*(180/pi)); 00846 } 00847 00848 orbit.mu = GetG()*GetMSun(); 00849 00850 par_class * par = (par_class *) params; 00851 00852 vector<Observation> * obs = par->obs; 00853 00854 orbit.epoch = par->orbit_epoch; 00855 00856 OptimizedOrbitPositions opt(orbit); 00857 00858 size_t i; 00859 double arcsec; 00860 for (i=0;i<obs->size();++i) { 00861 arcsec = opt.PropagatedSky_J2000((*obs)[i].date,(*obs)[i].obscode).delta_arcsec((*obs)[i]); 00862 gsl_vector_set(f,i,arcsec); 00863 00864 ORSA_ERROR("f[%02i] = %g",i,arcsec); 00865 } 00866 00867 return GSL_SUCCESS; 00868 }
Here is the call graph for this function:
int orsa::least_sq_fdf | ( | const gsl_vector * | v, | |
void * | params, | |||
gsl_vector * | f, | |||
gsl_matrix * | J | |||
) |
Definition at line 1016 of file orsa_orbit_gsl.cc.
References least_sq_df(), and least_sq_f().
Referenced by OrbitDifferentialCorrectionsLeastSquares().
01016 { 01017 least_sq_f (v,params,f); 01018 least_sq_df(v,params,J); 01019 return GSL_SUCCESS; 01020 }
Here is the call graph for this function:
string LengthLabel | ( | ) |
Definition at line 140 of file orsa_universe.cc.
References Units::LengthLabel(), and units.
Referenced by OrsaFile::Read().
00140 { return orsa::units->LengthLabel(); };
Here is the call graph for this function:
static double orsa::local_C22 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 278 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00278 { 00279 if (planet == MOON) { 00280 return jpl_file->GetTag("C22M"); 00281 } else { 00282 return 0.0; 00283 } 00284 }
Here is the call graph for this function:
static double orsa::local_C31 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 286 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00286 { 00287 if (planet == MOON) { 00288 return jpl_file->GetTag("C31M"); 00289 } else { 00290 return 0.0; 00291 } 00292 }
Here is the call graph for this function:
static double orsa::local_C32 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 294 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00294 { 00295 if (planet == MOON) { 00296 return jpl_file->GetTag("C32M"); 00297 } else { 00298 return 0.0; 00299 } 00300 }
Here is the call graph for this function:
static double orsa::local_C33 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 302 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00302 { 00303 if (planet == MOON) { 00304 return jpl_file->GetTag("C33M"); 00305 } else { 00306 return 0.0; 00307 } 00308 }
Here is the call graph for this function:
static double orsa::local_C41 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 310 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00310 { 00311 if (planet == MOON) { 00312 return jpl_file->GetTag("C41M"); 00313 } else { 00314 return 0.0; 00315 } 00316 }
Here is the call graph for this function:
static double orsa::local_C42 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 318 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00318 { 00319 if (planet == MOON) { 00320 return jpl_file->GetTag("C42M"); 00321 } else { 00322 return 0.0; 00323 } 00324 }
Here is the call graph for this function:
static double orsa::local_C43 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 326 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00326 { 00327 if (planet == MOON) { 00328 return jpl_file->GetTag("C43M"); 00329 } else { 00330 return 0.0; 00331 } 00332 }
Here is the call graph for this function:
static double orsa::local_C44 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 334 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00334 { 00335 if (planet == MOON) { 00336 return jpl_file->GetTag("C44M"); 00337 } else { 00338 return 0.0; 00339 } 00340 }
Here is the call graph for this function:
static double orsa::local_J2 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 247 of file orsa_body.cc.
References EARTH, JPLFile::GetTag(), jpl_file, MOON, and SUN.
00247 { 00248 double J2_ = 0.0; 00249 switch (planet) { 00250 case SUN: J2_ = jpl_file->GetTag("J2SUN"); break; 00251 case EARTH: J2_ = jpl_file->GetTag("J2E"); break; 00252 case MOON: J2_ = jpl_file->GetTag("J2M"); break; 00253 default: /*********************************/ break; 00254 } 00255 return J2_; 00256 }
Here is the call graph for this function:
static double orsa::local_J3 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 258 of file orsa_body.cc.
References EARTH, JPLFile::GetTag(), jpl_file, and MOON.
00258 { 00259 double J3_ = 0.0; 00260 switch (planet) { 00261 case EARTH: J3_ = jpl_file->GetTag("J3E"); break; 00262 case MOON: J3_ = jpl_file->GetTag("J3M"); break; 00263 default: /*********************************/ break; 00264 } 00265 return J3_; 00266 }
Here is the call graph for this function:
static double orsa::local_J4 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 268 of file orsa_body.cc.
References EARTH, JPLFile::GetTag(), jpl_file, and MOON.
00268 { 00269 double J4_ = 0.0; 00270 switch (planet) { 00271 case EARTH: J4_ = jpl_file->GetTag("J4E"); break; 00272 case MOON: J4_ = jpl_file->GetTag("J4M"); break; 00273 default: /*********************************/ break; 00274 } 00275 return J4_; 00276 }
Here is the call graph for this function:
static double orsa::local_mass | ( | const JPL_planets | planet | ) | [static] |
Definition at line 243 of file orsa_body.cc.
References JPLFile::GetMass(), and jpl_file.
00243 { 00244 return jpl_file->GetMass(planet); 00245 }
Here is the call graph for this function:
static double orsa::local_S31 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 342 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00342 { 00343 if (planet == MOON) { 00344 return jpl_file->GetTag("S31M"); 00345 } else { 00346 return 0.0; 00347 } 00348 }
Here is the call graph for this function:
static double orsa::local_S32 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 350 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00350 { 00351 if (planet == MOON) { 00352 return jpl_file->GetTag("S32M"); 00353 } else { 00354 return 0.0; 00355 } 00356 }
Here is the call graph for this function:
static double orsa::local_S33 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 358 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00358 { 00359 if (planet == MOON) { 00360 return jpl_file->GetTag("S33M"); 00361 } else { 00362 return 0.0; 00363 } 00364 }
Here is the call graph for this function:
static double orsa::local_S41 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 366 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00366 { 00367 if (planet == MOON) { 00368 return jpl_file->GetTag("S41M"); 00369 } else { 00370 return 0.0; 00371 } 00372 }
Here is the call graph for this function:
static double orsa::local_S42 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 374 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00374 { 00375 if (planet == MOON) { 00376 return jpl_file->GetTag("S42M"); 00377 } else { 00378 return 0.0; 00379 } 00380 }
Here is the call graph for this function:
static double orsa::local_S43 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 382 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00382 { 00383 if (planet == MOON) { 00384 return jpl_file->GetTag("S43M"); 00385 } else { 00386 return 0.0; 00387 } 00388 }
Here is the call graph for this function:
static double orsa::local_S44 | ( | const JPL_planets | planet | ) | [static] |
Definition at line 390 of file orsa_body.cc.
References JPLFile::GetTag(), jpl_file, and MOON.
00390 { 00391 if (planet == MOON) { 00392 return jpl_file->GetTag("S44M"); 00393 } else { 00394 return 0.0; 00395 } 00396 }
Here is the call graph for this function:
double orsa::M1 | ( | void * | p1, | |
void * | p2 | |||
) |
Definition at line 593 of file orsa_orbit_gsl.cc.
References c_a, c_e, c_i, c_M, c_omega_node, c_omega_pericenter, and secure_pow().
Referenced by MOID_f().
00593 { 00594 data *d1 = (data*) p1; 00595 data *d2 = (data*) p2; 00596 00597 double dist = sqrt( c_a * secure_pow( d1->orbit->a - d2->orbit->a, 2) + 00598 c_e * secure_pow( d1->orbit->e - d2->orbit->e, 2) + 00599 c_i * secure_pow( d1->orbit->i - d2->orbit->i, 2) + 00600 c_omega_node * secure_pow( d1->orbit->omega_node - d2->orbit->omega_node, 2) + 00601 c_omega_pericenter * secure_pow( d1->orbit->omega_pericenter - d2->orbit->omega_pericenter,2) + 00602 c_M * secure_pow( d1->orbit->M - d2->orbit->M, 2) ) / sqrt(6.0); 00603 00604 std::cerr << "M1: " << dist << endl; 00605 00606 return dist; 00607 }
Here is the call graph for this function:
void make_new_integrator | ( | Integrator ** | i, | |
const IntegratorType | type | |||
) |
Definition at line 29 of file orsa_integrator.cc.
References DISSIPATIVERUNGEKUTTA, LEAPFROG, RA15, RUNGEKUTTA, and STOER.
Referenced by OrsaFile::Read(), and Evolution::SetIntegrator().
00029 { 00030 00031 delete (*i); 00032 (*i) = 0; 00033 00034 switch (type) { 00035 case STOER: (*i) = new Stoer; break; 00036 // case BULIRSCHSTOER: (*i) = new BulirschStoer; break; 00037 case RUNGEKUTTA: (*i) = new RungeKutta; break; 00038 case DISSIPATIVERUNGEKUTTA: (*i) = new DissipativeRungeKutta; break; 00039 case RA15: (*i) = new Radau15; break; 00040 case LEAPFROG: (*i) = new Leapfrog; break; 00041 } 00042 }
void make_new_interaction | ( | Interaction ** | i, | |
const InteractionType | type | |||
) |
Definition at line 55 of file orsa_interaction.cc.
References ARMONICOSCILLATOR, GALACTIC_POTENTIAL_ALLEN, GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON, GRAVITATIONALTREE, JPL_PLANETS_NEWTON, NEWTON, NEWTON_MPI, ORSA_WARNING, and RELATIVISTIC.
Referenced by OrsaFile::Read(), and Evolution::SetInteraction().
00055 { 00056 00057 delete (*i); 00058 (*i) = 0; 00059 00060 switch (type) { 00061 case NEWTON: (*i) = new Newton; break; 00062 case NEWTON_MPI: 00063 #ifdef HAVE_MPI 00064 (*i) = new Newton_MPI; 00065 #else 00066 ORSA_WARNING("read NEWTON_MPI interaction from application without MPI support."); 00067 #endif 00068 break; 00069 case ARMONICOSCILLATOR: (*i) = new ArmonicOscillator(1,1); break; 00070 case GALACTIC_POTENTIAL_ALLEN: (*i) = new GalacticPotentialAllen; break; 00071 case GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON: (*i) = new GalacticPotentialAllenPlusNewton; break; 00072 case JPL_PLANETS_NEWTON: /*************************************/ break; 00073 case GRAVITATIONALTREE: (*i) = new GravitationalTree; break; 00074 case RELATIVISTIC: (*i) = new Relativistic; break; 00075 } 00076 }
string MassLabel | ( | ) |
Definition at line 141 of file orsa_universe.cc.
References Units::MassLabel(), and units.
Referenced by OrsaFile::Read().
00141 { return orsa::units->MassLabel(); };
Here is the call graph for this function:
static double orsa::modified_mu | ( | const vector< Body > & | f, | |
const unsigned int | i | |||
) | [static] |
Definition at line 89 of file orsa_frame.cc.
References GetC().
Referenced by RelativisticBarycenter(), and RelativisticBarycenterVelocity().
00089 { 00090 if (f[i].has_zero_mass()) return 0.0; 00091 const double one_over_two_c = 1.0/(2.0*GetC()); 00092 // const double g = GetG(); 00093 // 00094 vector<double> mu(f.size()); 00095 for (unsigned int j=0; j<f.size(); ++j) { 00096 if (f[j].has_zero_mass()) { 00097 mu[j] = 0.0; 00098 } else { 00099 mu[j] = f[j].mu(); 00100 } 00101 } 00102 // 00103 double mod_mu = 0.0; 00104 double tmp_sum = 0.0; 00105 for (unsigned int j=0; j<f.size(); ++j) { 00106 if (j == i) continue; 00107 tmp_sum += mu[j]/(f[j].position()-f[i].position()).Length(); 00108 } 00109 // 00110 mod_mu += mu[i]*(1.0+one_over_two_c*(f[i].velocity().LengthSquared()-tmp_sum)); 00111 return mod_mu; 00112 }
Here is the call graph for this function:
double MOID | ( | const Orbit & | o1_in, | |
const Orbit & | o2_in, | |||
Vector & | r1, | |||
Vector & | r2 | |||
) |
Minimal Orbit Intersection Distance between two orbits.
Definition at line 223 of file orsa_orbit_gsl.cc.
References Orbit::M, MOID_f(), Orbit::omega_node, Orbit::omega_pericenter, pi, Orbit::RelativePosVel(), and twopi.
00223 { 00224 00225 Orbit o1(o1_in); 00226 Orbit o2(o2_in); 00227 00228 MOID_parameters par; 00229 par.o1 = &o1; 00230 par.o2 = &o2; 00231 00232 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex,2); 00233 00234 gsl_multimin_function moid_function; 00235 00236 moid_function.f = &MOID_f; 00237 moid_function.n = 2; 00238 moid_function.params = ∥ 00239 00240 gsl_vector *x = gsl_vector_alloc(2); 00241 gsl_vector *step_size = gsl_vector_alloc(2); 00242 00243 // gsl_vector_set(step_size,0,pi); 00244 // gsl_vector_set(step_size,1,pi); 00245 00246 MOID_solution best_solution; 00247 00248 int gsl_solutions_iter=0,max_gsl_solutions_iter=7; 00249 while (gsl_solutions_iter<max_gsl_solutions_iter) { 00250 00251 // gsl_vector_set(x,0,o1.M); 00252 // gsl_vector_set(x,1,o2.M); 00253 // gsl_vector_set(x,0,0.0); 00254 // gsl_vector_set(x,1,0.0); 00255 // gsl_vector_set(x,0,fmod(10*twopi -o1.omega_node-o1.omega_pericenter,twopi)); 00256 // gsl_vector_set(x,1,fmod(10*twopi+pi-o2.omega_node-o2.omega_pericenter,twopi)); 00257 gsl_vector_set(x,0,fmod(10*twopi+(twopi/(2*max_gsl_solutions_iter))*gsl_solutions_iter -o1.omega_node-o1.omega_pericenter,twopi)); 00258 gsl_vector_set(x,1,fmod(10*twopi+(twopi/(2*max_gsl_solutions_iter))*gsl_solutions_iter+pi-o2.omega_node-o2.omega_pericenter,twopi)); 00259 00260 gsl_vector_set(step_size,0,pi); 00261 gsl_vector_set(step_size,1,pi); 00262 00263 gsl_multimin_fminimizer_set (s, &moid_function, x, step_size); 00264 00265 size_t iter = 0; 00266 int status; 00267 do { 00268 00269 iter++; 00270 00271 // iter status 00272 status = gsl_multimin_fminimizer_iterate(s); 00273 00274 // convergence status 00275 status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s),1.0e-6); 00276 00277 // if (status == GSL_SUCCESS) printf ("Minimum found at:\n"); 00278 00279 /* printf ("%5d %.5f %.5f %10.5f\n", iter, 00280 (180/pi)*gsl_vector_get (s->x, 0), 00281 (180/pi)*gsl_vector_get (s->x, 1), 00282 gsl_multimin_fminimizer_minimum(s)); 00283 */ 00284 00285 // } while (status == GSL_CONTINUE && iter < 500); 00286 } while (status == GSL_CONTINUE && iter < 200); 00287 00288 // const double moid = gsl_multimin_fminimizer_minimum(s); 00289 00290 // cerr << "MOID iters: " << iter << endl; 00291 00292 if (status == GSL_SUCCESS) { 00293 best_solution.Insert(gsl_multimin_fminimizer_minimum(s),gsl_vector_get(s->x,0),gsl_vector_get(s->x,1)); 00294 } 00295 00296 ++gsl_solutions_iter; 00297 } 00298 00299 const double moid = best_solution.moid(); 00300 o1.M = best_solution.M1(); 00301 o2.M = best_solution.M2(); 00302 Vector v1, v2; 00303 // update r1 and r2 for output 00304 o1.RelativePosVel(r1,v1); 00305 o2.RelativePosVel(r2,v2); 00306 00307 /* 00308 printf ("%5d %.5f %.5f %.5f %.5f %10.5f\n", iter, 00309 (180/pi)*gsl_vector_get (s->x, 0), (180/pi)*o1.M, 00310 (180/pi)*gsl_vector_get (s->x, 1), (180/pi)*o2.M, 00311 gsl_multimin_fminimizer_minimum(s)); 00312 */ 00313 00314 // free 00315 gsl_multimin_fminimizer_free(s); 00316 gsl_vector_free(x); 00317 gsl_vector_free(step_size); 00318 00319 // check 00320 // WARNING! some vars can be modified by 00321 // this check area --> comment out when not debugging! 00322 if (0) { 00323 00324 Orbit o1(o1_in); 00325 Orbit o2(o2_in); 00326 00327 // random number generator 00328 // const int random_seed = 124323; 00329 struct timeval tv; 00330 struct timezone tz; 00331 gettimeofday(&tv,&tz); 00332 const int random_seed = tv.tv_usec; 00333 gsl_rng * rnd = gsl_rng_alloc(gsl_rng_gfsr4); 00334 gsl_rng_set(rnd,random_seed); 00335 int j,k; 00336 Vector r1,v1,r2,v2; 00337 double d, min_d=moid; 00338 // 00339 /* 00340 double angle; 00341 j=0; 00342 while (j < 500) { 00343 // coordinated test 00344 angle = twopi*gsl_rng_uniform(rnd); 00345 o1.M = fmod(10*twopi+angle-o1.omega_node-o1.omega_pericenter,twopi); 00346 o2.M = fmod(10*twopi+angle-o2.omega_node-o2.omega_pericenter,twopi); 00347 o1.RelativePosVel(r1,v1); 00348 o2.RelativePosVel(r2,v2); 00349 d = (r1-r2).Length(); 00350 if (d < min_d) { 00351 cerr << "(C) "; 00352 min_d = d; 00353 } 00354 ++j; 00355 } 00356 */ 00357 // 00358 // random test 00359 j=0; 00360 while (j < 50) { 00361 o1.M = twopi*gsl_rng_uniform(rnd); 00362 o1.RelativePosVel(r1,v1); 00363 k=0; 00364 while (k < 50) { 00365 // o1.M = twopi*gsl_rng_uniform(rnd); 00366 o2.M = twopi*gsl_rng_uniform(rnd); 00367 // o1.RelativePosVel(r1,v1); 00368 o2.RelativePosVel(r2,v2); 00369 d = (r1-r2).Length(); 00370 if (d < min_d) { 00371 // cerr << "WARNING: found another solution for the MOID: " << d << " < " << moid << " delta = " << d-moid << endl; 00372 // printf("WARNING: found another solution for the MOID: %.10f < %.10f delta: %.10f\n",d,moid,d-moid); 00373 // exit(0); 00374 std::cerr << "(R) "; 00375 min_d = d; 00376 } 00377 ++k; 00378 } 00379 ++j; 00380 } 00381 // 00382 if (min_d < moid) { 00383 // cerr << "WARNING: found another solution for the MOID: " << d << " < " << moid << " delta = " << d-moid << endl; 00384 printf("\n WARNING: found another solution for the MOID: %.10f < %.10f delta: %.10f\n",min_d,moid,min_d-moid); 00385 // exit(0); 00386 } 00387 // free random number generator 00388 gsl_rng_free(rnd); 00389 00390 } 00391 00392 return moid; 00393 }
Here is the call graph for this function:
double MOID2RB | ( | const Vector & | rb1, | |
const Vector & | rb2, | |||
const Orbit & | o1_in, | |||
const Orbit & | o2_in, | |||
Vector & | r1, | |||
Vector & | r2 | |||
) |
Minimal Orbit Intersection Distance between two orbits with different reference bodies.
Definition at line 447 of file orsa_orbit_gsl.cc.
References Orbit::M, MOID2RB_f(), Orbit::omega_node, Orbit::omega_pericenter, pi, Orbit::RelativePosVel(), and twopi.
00447 { 00448 00449 Orbit o1(o1_in); 00450 Orbit o2(o2_in); 00451 00452 MOID2RB_parameters par; 00453 par.o1 = &o1; 00454 par.o2 = &o2; 00455 par.rb1 = rb1; 00456 par.rb2 = rb2; 00457 00458 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex,2); 00459 00460 gsl_multimin_function moid2rb_function; 00461 00462 moid2rb_function.f = &MOID2RB_f; 00463 moid2rb_function.n = 2; 00464 moid2rb_function.params = ∥ 00465 00466 gsl_vector *x = gsl_vector_alloc(2); 00467 gsl_vector *step_size = gsl_vector_alloc(2); 00468 00469 // gsl_vector_set(step_size,0,pi); 00470 // gsl_vector_set(step_size,1,pi); 00471 00472 MOID2RB_solution best_solution; 00473 00474 int gsl_solutions_iter=0,max_gsl_solutions_iter=7; 00475 while (gsl_solutions_iter<max_gsl_solutions_iter) { 00476 00477 gsl_vector_set(x,0,fmod(10*twopi+(twopi/(2*max_gsl_solutions_iter))*gsl_solutions_iter -o1.omega_node-o1.omega_pericenter,twopi)); 00478 gsl_vector_set(x,1,fmod(10*twopi+(twopi/(2*max_gsl_solutions_iter))*gsl_solutions_iter+pi-o2.omega_node-o2.omega_pericenter,twopi)); 00479 00480 gsl_vector_set(step_size,0,pi); 00481 gsl_vector_set(step_size,1,pi); 00482 00483 gsl_multimin_fminimizer_set (s, &moid2rb_function, x, step_size); 00484 00485 size_t iter = 0; 00486 int status; 00487 do { 00488 00489 iter++; 00490 00491 // iter status 00492 status = gsl_multimin_fminimizer_iterate(s); 00493 00494 // convergence status 00495 status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s),1.0e-6); 00496 00497 // if (status == GSL_SUCCESS) printf ("Minimum found at:\n"); 00498 00499 /* printf ("%5d %.5f %.5f %10.5f\n", iter, 00500 (180/pi)*gsl_vector_get (s->x, 0), 00501 (180/pi)*gsl_vector_get (s->x, 1), 00502 gsl_multimin_fminimizer_minimum(s)); 00503 */ 00504 00505 // } while (status == GSL_CONTINUE && iter < 500); 00506 } while (status == GSL_CONTINUE && iter < 200); 00507 00508 // const double moid2rb = gsl_multimin_fminimizer_minimum(s); 00509 00510 // cerr << "MOID2RB iters: " << iter << endl; 00511 00512 if (status == GSL_SUCCESS) { 00513 best_solution.Insert(gsl_multimin_fminimizer_minimum(s),gsl_vector_get(s->x,0),gsl_vector_get(s->x,1)); 00514 } 00515 00516 ++gsl_solutions_iter; 00517 } 00518 00519 const double moid2rb = best_solution.moid2rb(); 00520 o1.M = best_solution.M1(); 00521 o2.M = best_solution.M2(); 00522 Vector v1, v2; 00523 // update r1 and r2 for output 00524 o1.RelativePosVel(r1,v1); 00525 o2.RelativePosVel(r2,v2); 00526 // not here... 00527 // r1 += rb1; 00528 // r2 += rb2; 00529 00530 /* 00531 printf ("%5d %.5f %.5f %.5f %.5f %10.5f\n", iter, 00532 (180/pi)*gsl_vector_get (s->x, 0), (180/pi)*o1.M, 00533 (180/pi)*gsl_vector_get (s->x, 1), (180/pi)*o2.M, 00534 gsl_multimin_fminimizer_minimum(s)); 00535 */ 00536 00537 // free 00538 gsl_multimin_fminimizer_free(s); 00539 gsl_vector_free(x); 00540 gsl_vector_free(step_size); 00541 00542 return moid2rb; 00543 }
Here is the call graph for this function:
double orsa::MOID2RB_f | ( | const gsl_vector * | v, | |
void * | params | |||
) |
Definition at line 404 of file orsa_orbit_gsl.cc.
Referenced by MOID2RB().
00404 { 00405 00406 double M1 = gsl_vector_get(v,0); 00407 double M2 = gsl_vector_get(v,1); 00408 00409 MOID2RB_parameters *p = (MOID2RB_parameters*) params; 00410 00411 p->o1->M = M1; 00412 p->o2->M = M2; 00413 00414 Vector r1,r2,v1,v2; 00415 00416 p->o1->RelativePosVel(r1,v1); 00417 p->o2->RelativePosVel(r2,v2); 00418 00419 return ((p->rb1+r1)-(p->rb2+r2)).Length(); 00420 }
double orsa::MOID_f | ( | const gsl_vector * | v, | |
void * | params | |||
) |
Definition at line 180 of file orsa_orbit_gsl.cc.
References M1().
Referenced by MOID().
00180 { 00181 00182 double M1 = gsl_vector_get(v,0); 00183 double M2 = gsl_vector_get(v,1); 00184 00185 MOID_parameters *p = (MOID_parameters*) params; 00186 00187 p->o1->M = M1; 00188 p->o2->M = M2; 00189 00190 Vector r1,r2,v1,v2; 00191 00192 p->o1->RelativePosVel(r1,v1); 00193 p->o2->RelativePosVel(r2,v2); 00194 00195 return (r1-r2).Length(); 00196 }
Here is the call graph for this function:
double orsa::norm | ( | const fftw_complex | z | ) |
Definition at line 38 of file orsa_fft.cc.
Referenced by accurate_peak(), amph(), phi_gsl_two(), psd_max(), psd_max_again(), and psd_max_again_many().
double orsa::norm_sq | ( | const fftw_complex | z | ) |
Definition at line 42 of file orsa_fft.cc.
Referenced by phi_amp(), phi_Hanning_amp(), and psd_max().
Angle obleq | ( | const Date & | date | ) |
Definition at line 1851 of file orsa_units.cc.
References Date::GetJulian(), Angle::SetDPS(), and UT.
Referenced by EclipticToEquatorial(), EquatorialToEcliptic(), and obleq_J2000().
01851 { 01852 // T in centuries from J2000 01853 Date d = date; 01854 // d.ConvertToTimeScale(UT); // UT or UTC ?? 01855 // const double T = (d.GetJulian() - 2451545.0)/36525.0; 01856 // DOUBLE-CHECK this "UT"!!! 01857 const double T = (d.GetJulian(UT) - 2451545.0)/36525.0; 01858 Angle a; 01859 // updated Feb 2004 01860 a.SetDPS(23,26,21.448+((0.001813*T-0.00059)*T-46.8150)*T); 01861 return a; 01862 }
Here is the call graph for this function:
Angle obleq_J2000 | ( | ) |
Definition at line 1909 of file orsa_units.cc.
References obleq(), and Date::SetJ2000().
Referenced by Compute_Gauss(), Compute_TestMethod(), EclipticToEquatorial_J2000(), EquatorialToEcliptic_J2000(), JPLFile::GetEph(), and TLEFile::Read().
01909 { 01910 Date d; d.SetJ2000(); 01911 return obleq(d); 01912 }
Here is the call graph for this function:
UniverseTypeAwareTimeStep operator * | ( | const UniverseTypeAwareTimeStep & | ts, | |
const int | i | |||
) |
UniverseTypeAwareTimeStep operator * | ( | const int | i, | |
const UniverseTypeAwareTimeStep & | ts | |||
) |
UniverseTypeAwareTimeStep operator * | ( | const UniverseTypeAwareTimeStep & | ts, | |
const double | x | |||
) |
UniverseTypeAwareTimeStep operator * | ( | const double | x, | |
const UniverseTypeAwareTimeStep & | ts | |||
) |
double orsa::operator * | ( | const Vector & | u, | |
const Vector & | v | |||
) | [inline] |
Vector orsa::operator * | ( | const Vector & | v, | |
const double | f | |||
) | [inline] |
Vector orsa::operator * | ( | const double | f, | |
const Vector & | v | |||
) | [inline] |
bool orsa::operator!= | ( | const ET_minus_UT_element & | x, | |
const ET_minus_UT_element & | y | |||
) | [inline] |
bool orsa::operator!= | ( | const TAI_minus_UTC_element & | x, | |
const TAI_minus_UTC_element & | y | |||
) | [inline] |
bool orsa::operator!= | ( | const Vector & | v1, | |
const Vector & | v2 | |||
) | [inline] |
bool orsa::operator!= | ( | const Body & | b1, | |
const Body & | b2 | |||
) | [inline] |
Vector orsa::operator+ | ( | const Vector & | u, | |
const Vector & | v | |||
) | [inline] |
Vector orsa::operator- | ( | const Vector & | u, | |
const Vector & | v | |||
) | [inline] |
Vector orsa::operator/ | ( | const Vector & | v, | |
const double | f | |||
) | [inline] |
bool orsa::operator== | ( | const ET_minus_UT_element & | x, | |
const ET_minus_UT_element & | y | |||
) | [inline] |
Definition at line 280 of file orsa_units.cc.
00280 { 00281 if (x.day != y.day) return false; 00282 if (x.month != y.month) return false; 00283 if (x.year != y.year) return false; 00284 if (x.ET_minus_UT != y.ET_minus_UT) return false; 00285 return true; 00286 }
bool orsa::operator== | ( | const TAI_minus_UTC_element & | x, | |
const TAI_minus_UTC_element & | y | |||
) | [inline] |
Definition at line 225 of file orsa_units.cc.
00225 { 00226 if (x.day != y.day) return false; 00227 if (x.month != y.month) return false; 00228 if (x.year != y.year) return false; 00229 if (x.TAI_minus_UTC != y.TAI_minus_UTC) return false; 00230 return true; 00231 }
bool orsa::operator== | ( | const Vector & | v1, | |
const Vector & | v2 | |||
) | [inline] |
Definition at line 195 of file orsa_coord.h.
References Vector::x, Vector::y, and Vector::z.
00195 { 00196 if (v1.x != v2.x) return false; 00197 if (v1.y != v2.y) return false; 00198 if (v1.z != v2.z) return false; 00199 return true; 00200 }
bool orsa::operator== | ( | const Body & | b1, | |
const Body & | b2 | |||
) | [inline] |
Definition at line 209 of file orsa_body.h.
References Body::BodyId(), Body::mass(), Body::name(), Body::position(), and Body::velocity().
00209 { 00210 if (b1.BodyId() != b2.BodyId()) return false; 00211 if (b1.name() != b2.name()) return false; 00212 if (b1.mass() != b2.mass()) return false; 00213 if (b1.position() != b2.position()) return false; 00214 if (b1.velocity() != b2.velocity()) return false; 00215 return true; 00216 }
Here is the call graph for this function:
void orsa::OrbitDifferentialCorrectionsLeastSquares | ( | OrbitWithCovarianceMatrixGSL & | , | |
const std::vector< Observation > & | ||||
) |
Referenced by Compute().
void orsa::OrbitDifferentialCorrectionsLeastSquares | ( | OrbitWithCovarianceMatrixGSL & | orbit, | |
const vector< Observation > & | obs | |||
) |
Definition at line 1023 of file orsa_orbit_gsl.cc.
References Orbit::a, Sky::dec(), Sky::delta_arcsec(), Orbit::e, OrbitWithEpoch::epoch, ERR, FIT, Angle::GetRad(), Orbit::i, least_sq_df(), least_sq_f(), least_sq_fdf(), Orbit::M, Orbit::omega_node, Orbit::omega_pericenter, Osculating, pi, OptimizedOrbitPositions::PropagatedSky_J2000(), Sky::ra(), secure_pow(), and OrbitWithCovarianceMatrixGSL::SetCovarianceMatrix().
01023 { 01024 01025 cerr << "inside odcls..." << endl; 01026 01027 if (obs.size() < 6) { 01028 cerr << "at least 6 observations are needed for OrbitDifferentialCorrectionsLeastSquares..." << endl; 01029 return; 01030 } 01031 01032 par_class par; 01033 01034 par.orbit_epoch = orbit.epoch; 01035 01036 par.obs = const_cast< vector<Observation>* > (&obs); 01037 01038 par.covar = gsl_matrix_alloc(6,6); 01039 par.use_covar = false; 01040 01041 /* 01042 { // check par.obs 01043 cerr << "par.obs->size()=" << par.obs->size() << endl; 01044 int y,m,d; 01045 for (unsigned int k=0;k<par.obs->size();++k) { 01046 (*par.obs)[k].date.GetGregor(y,m,d); 01047 cerr << "(*par.obs)[" << k << "] date: " << y << " " << m << " " << d << " t: " << (*par.obs)[k].date.GetJulian() << endl; 01048 } 01049 } 01050 */ 01051 01052 gsl_multifit_fdfsolver * s = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, obs.size(), 6); 01053 // gsl_multifit_fdfsolver * s = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmder, obs.size(), 6); 01054 01055 gsl_multifit_function_fdf least_sq_function; 01056 01057 least_sq_function.f = &least_sq_f; 01058 least_sq_function.df = &least_sq_df; 01059 least_sq_function.fdf = &least_sq_fdf; 01060 least_sq_function.n = obs.size(); 01061 least_sq_function.p = 6; 01062 least_sq_function.params = ∥ 01063 01064 gsl_vector * x = gsl_vector_alloc(6); 01065 01066 gsl_vector_set(x, 0, orbit.a); 01067 gsl_vector_set(x, 1, orbit.e); 01068 gsl_vector_set(x, 2, orbit.i); 01069 gsl_vector_set(x, 3, orbit.omega_node); 01070 gsl_vector_set(x, 4, orbit.omega_pericenter); 01071 gsl_vector_set(x, 5, orbit.M); 01072 01073 gsl_multifit_fdfsolver_set(s,&least_sq_function,x); 01074 01075 size_t iter = 0; 01076 int status; 01077 do { 01078 01079 ++iter; 01080 01081 status = gsl_multifit_fdfsolver_iterate(s); 01082 01083 printf("itaration status = %s\n",gsl_strerror(status)); 01084 01085 /* 01086 printf ("iter: %3u %15.8f %15.8f %15.8f |f(x)| = %g\n", 01087 iter, 01088 gsl_vector_get (s->x, 0), 01089 gsl_vector_get (s->x, 1), 01090 gsl_vector_get (s->x, 2), 01091 gsl_blas_dnrm2 (s->f)); 01092 */ 01093 01094 // status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-4); 01095 // status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-5); 01096 // status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-6); 01097 status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1.0e-8); 01098 // status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-8); 01099 // status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-9); 01100 01101 printf("convergence status = %s\n",gsl_strerror(status)); 01102 01103 { 01104 gsl_matrix * covar = gsl_matrix_alloc(6,6); 01105 gsl_multifit_covar(s->J,0.0,covar); 01106 printf ("more info:\n" 01107 "a = %12.6f +/- %12.6f\n" 01108 "e = %12.6f +/- %12.6f\n" 01109 "i = %12.6f +/- %12.6f\n" 01110 "\n", 01111 gsl_vector_get(s->x,0),sqrt(gsl_matrix_get(covar,0,0)), 01112 gsl_vector_get(s->x,1),sqrt(gsl_matrix_get(covar,1,1)), 01113 gsl_vector_get(s->x,2)*(180/pi),sqrt(gsl_matrix_get(covar,2,2))*(180/pi)); 01114 } 01115 01116 gsl_multifit_covar(s->J,0.0,par.covar); 01117 par.use_covar = true; 01118 01119 } while ((status == GSL_CONTINUE) && (iter < 1024)); 01120 01121 gsl_matrix * covar = gsl_matrix_alloc(6,6); 01122 gsl_multifit_covar(s->J,0.0,covar); 01123 // gsl_matrix_fprintf(stdout, covar,"%g"); 01124 01125 #define FIT(i) gsl_vector_get(s->x, i) 01126 #define ERR(i) sqrt(gsl_matrix_get(covar,i,i)) 01127 01128 printf("a = %12.6f +/- %12.6f\n", FIT(0), ERR(0)); 01129 printf("e = %12.6f +/- %12.6f\n", FIT(1), ERR(1)); 01130 printf("i = %12.6f +/- %12.6f\n", FIT(2)*(180/pi), ERR(2)*(180/pi)); 01131 printf("node = %12.6f +/- %12.6f\n", FIT(3)*(180/pi), ERR(3)*(180/pi)); 01132 printf("peri = %12.6f +/- %12.6f\n", FIT(4)*(180/pi), ERR(4)*(180/pi)); 01133 printf("M = %12.6f +/- %12.6f\n", FIT(5)*(180/pi), ERR(5)*(180/pi)); 01134 01135 printf ("status = %s\n", gsl_strerror (status)); 01136 01137 // update orbit! 01138 orbit.a = gsl_vector_get (s->x,0); 01139 orbit.e = gsl_vector_get (s->x,1); 01140 orbit.i = gsl_vector_get (s->x,2); 01141 orbit.omega_node = gsl_vector_get (s->x,3); 01142 orbit.omega_pericenter = gsl_vector_get (s->x,4); 01143 orbit.M = gsl_vector_get (s->x,5); 01144 01145 // set covariance matrix 01146 { 01147 double covm[6][6]; 01148 unsigned int i,j; 01149 for (i=0;i<6;++i) { 01150 for (j=0;j<6;++j) { 01151 covm[i][j] = gsl_matrix_get(covar,i,j); 01152 } 01153 } 01154 orbit.SetCovarianceMatrix((double**)covm,Osculating); 01155 } 01156 01157 // print errors and RMS. 01158 { 01159 OptimizedOrbitPositions opt(orbit); 01160 double arcsec, rms=0.0; 01161 for (unsigned int k=0;k<obs.size();++k) { 01162 arcsec = opt.PropagatedSky_J2000(obs[k].date,obs[k].obscode).delta_arcsec(obs[k]); 01163 rms += secure_pow(arcsec,2); 01164 fprintf(stderr,"delta_arcsec obs[%02i]: %g (%+.2f,%+.2f)\n",k,arcsec, 01165 (opt.PropagatedSky_J2000(obs[k].date,obs[k].obscode).ra().GetRad() -obs[k].ra.GetRad() )*(180/pi)*3600.0, 01166 (opt.PropagatedSky_J2000(obs[k].date,obs[k].obscode).dec().GetRad()-obs[k].dec.GetRad())*(180/pi)*3600.0); 01167 } 01168 rms = sqrt(rms/obs.size()); 01169 fprintf(stderr,"\n RMS: %g\n\n",rms); 01170 } 01171 01172 // free 01173 gsl_multifit_fdfsolver_free (s); 01174 gsl_vector_free (x); 01175 01176 cerr << "out of odcls..." << endl; 01177 }
Here is the call graph for this function:
fftw_complex orsa::phi | ( | double | omega, | |
fftw_complex | in[], | |||
const int | size | |||
) |
Discrete Fourier Transform.
Definition at line 49 of file orsa_fft.cc.
References cos(), sin(), and twopi.
Referenced by accurate_peak(), amph(), phi_amp(), and phi_gsl_two().
00049 { 00050 00051 fftw_complex result; 00052 result.re = 0; 00053 result.im = 0; 00054 00055 double arg,c,s; 00056 int k; 00057 00058 for(k=0;k<size;k++) { 00059 arg = twopi*k*omega; 00060 c = cos(arg); 00061 s = sin(arg); 00062 result.re += in[k].re * c + in[k].im * s; 00063 result.im -= in[k].re * s - in[k].im * c; 00064 } 00065 00066 double scale = (double)size; 00067 result.re /= scale; 00068 result.im /= scale; 00069 00070 // OUT HERE! 00071 return result; 00072 }
Here is the call graph for this function:
double orsa::phi_amp | ( | double | omega, | |
fftw_complex | in[], | |||
const int | size | |||
) |
double orsa::phi_gsl | ( | double | x, | |
void * | params | |||
) |
Definition at line 119 of file orsa_fft.cc.
References phi_amp(), gsl_d1_minimization_parameters::pointer_points_sequence, and gsl_d1_minimization_parameters::size.
00119 { 00120 00121 struct gsl_d1_minimization_parameters * p = (gsl_d1_minimization_parameters *) params; 00122 00123 // return (-phi_Hanning_amp(x, p->pointer_points_sequence, p->size)); 00124 return (-phi_amp(x, p->pointer_points_sequence, p->size)); 00125 00126 }
Here is the call graph for this function:
double orsa::phi_gsl_two | ( | double | x, | |
void * | params | |||
) |
Definition at line 128 of file orsa_fft.cc.
References norm(), phi(), gsl_d1_minimization_parameters::pointer_points_sequence, and gsl_d1_minimization_parameters::size.
Referenced by accurate_peak().
00128 { 00129 00130 struct gsl_d1_minimization_parameters * p = (gsl_d1_minimization_parameters *) params; 00131 00132 // return (-phi_Hanning_amp(x, p->pointer_points_sequence, p->size)); 00133 return (-norm(phi(x, p->pointer_points_sequence, p->size))); 00134 }
Here is the call graph for this function:
fftw_complex orsa::phi_Hanning | ( | double | omega, | |
fftw_complex | in[], | |||
const int | size | |||
) |
Discrete Fourier Transform with Hanning windowing.
Definition at line 77 of file orsa_fft.cc.
References cos(), sin(), and twopi.
Referenced by phi_Hanning_amp().
00077 { 00078 00079 fftw_complex result; 00080 result.re = 0; 00081 result.im = 0; 00082 00083 double arg,c,s; 00084 int k; 00085 double window_factor; 00086 00087 for(k=0;k<size;k++) { 00088 arg = twopi*k*omega; 00089 c = cos(arg); 00090 s = sin(arg); 00091 window_factor = (1-cos(k*twopi/size)); 00092 result.re += window_factor * (in[k].re * c + in[k].im * s); 00093 result.im -= window_factor * (in[k].re * s - in[k].im * c); 00094 } 00095 00096 double scale = (double)size; 00097 result.re /= scale; 00098 result.im /= scale; 00099 00100 // OUT HERE! 00101 return result; 00102 }
Here is the call graph for this function:
double orsa::phi_Hanning_amp | ( | double | omega, | |
fftw_complex | in[], | |||
const int | size | |||
) |
amplitude for spectrum, with Hanning windowing
Definition at line 113 of file orsa_fft.cc.
References norm_sq(), and phi_Hanning().
Referenced by phi_Hanning_gsl().
00113 { 00114 return sqrt( norm_sq(phi_Hanning( omega,&in[0],size)) + 00115 norm_sq(phi_Hanning(-omega,&in[0],size)) ); 00116 }
Here is the call graph for this function:
double orsa::phi_Hanning_gsl | ( | double | x, | |
void * | params | |||
) |
Definition at line 136 of file orsa_fft.cc.
References phi_Hanning_amp(), gsl_d1_minimization_parameters::pointer_points_sequence, and gsl_d1_minimization_parameters::size.
00136 { 00137 00138 struct gsl_d1_minimization_parameters * p = (gsl_d1_minimization_parameters *) params; 00139 00140 return (-phi_Hanning_amp(x, p->pointer_points_sequence, p->size)); 00141 // return (-phi_amp(x, p->pointer_points_sequence, p->size)); 00142 00143 }
Here is the call graph for this function:
double poly_8 | ( | double | x, | |
void * | params | |||
) |
Definition at line 1588 of file orsa_orbit_gsl.cc.
References secure_pow().
Referenced by poly_8_gsl_solve().
01588 { 01589 struct poly_8_params *p = (struct poly_8_params *) params; 01590 return (secure_pow(x,8) - p->coeff_6*secure_pow(x,6) - p->coeff_3*secure_pow(x,3) - p->coeff_0); 01591 }
Here is the call graph for this function:
double poly_8_df | ( | double | x, | |
void * | params | |||
) |
Definition at line 1593 of file orsa_orbit_gsl.cc.
References secure_pow().
Referenced by poly_8_gsl_solve().
01593 { 01594 struct poly_8_params *p = (struct poly_8_params *) params; 01595 return (8*secure_pow(x,7) - 6*p->coeff_6*secure_pow(x,5) - 3*p->coeff_3*secure_pow(x,2)); 01596 }
Here is the call graph for this function:
void poly_8_fdf | ( | double | x, | |
void * | params, | |||
double * | y, | |||
double * | dy | |||
) |
Definition at line 1598 of file orsa_orbit_gsl.cc.
References secure_pow().
Referenced by poly_8_gsl_solve().
01598 { 01599 struct poly_8_params *p = (struct poly_8_params *) params; 01600 *y = (secure_pow(x,8) - p->coeff_6*secure_pow(x,6) - p->coeff_3*secure_pow(x,3) - p->coeff_0); 01601 *dy = (8*secure_pow(x,7) - 6*p->coeff_6*secure_pow(x,5) - 3*p->coeff_3*secure_pow(x,2)); 01602 }
Here is the call graph for this function:
void orsa::poly_8_gsl_solve | ( | poly_8_params & | params, | |
vector< poly_8_solution > & | solutions | |||
) |
Definition at line 1609 of file orsa_orbit_gsl.cc.
References AU, FromUnits(), params, poly_8(), poly_8_df(), and poly_8_fdf().
Referenced by Compute_Gauss().
01609 { 01610 01611 const int max_iter = 100; 01612 const double x_start = FromUnits(0.0,AU); 01613 const double x_incr = FromUnits(0.2,AU); 01614 01615 const double nominal_relative_accuracy = 1.0e-5; 01616 01617 solutions.clear(); 01618 01619 poly_8_solution tmp_solution; 01620 01621 double x, x0; // this value is not used 01622 int gsl_status; 01623 // const gsl_root_fdfsolver_type *T; 01624 01625 gsl_root_fdfsolver *s; 01626 gsl_function_fdf FDF; 01627 01628 /* 01629 cerr << " poly_8_gsl_solve(): params.coeff_6 = " << params->coeff_6 << endl; 01630 cerr << " poly_8_gsl_solve(): params.coeff_3 = " << params->coeff_3 << endl; 01631 cerr << " poly_8_gsl_solve(): params.coeff_0 = " << params->coeff_0 << endl; 01632 */ 01633 01634 FDF.f = &poly_8; 01635 FDF.df = &poly_8_df; 01636 FDF.fdf = &poly_8_fdf; 01637 FDF.params = ¶ms; 01638 01639 // T = gsl_root_fdfsolver_steffenson; 01640 // s = gsl_root_fdfsolver_alloc (T); 01641 s = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_steffenson); 01642 01643 int iter = 0; 01644 while (iter<max_iter) { 01645 01646 ++iter; 01647 x = x_start+iter*x_incr; 01648 gsl_root_fdfsolver_set (s, &FDF, x); 01649 01650 int iter_gsl=0; 01651 const int max_iter_gsl = 100; 01652 do { 01653 ++iter_gsl; 01654 gsl_status = gsl_root_fdfsolver_iterate(s); 01655 x0 = x; 01656 x = gsl_root_fdfsolver_root(s); 01657 gsl_status = gsl_root_test_delta (x, x0, nominal_relative_accuracy, 0); 01658 // if (gsl_status == GSL_SUCCESS) printf ("Converged:\n"); 01659 // printf ("%5d %10.7f %10.7f\n",iter_gsl, x, x - x0); 01660 } while ((gsl_status == GSL_CONTINUE) && (iter_gsl < max_iter_gsl)); 01661 01662 if (gsl_status == GSL_SUCCESS) { 01663 tmp_solution.value = x; 01664 // gsl_root_fdfsolver_iterate(s); tmp_solution.error = fabs(x-gsl_root_fdfsolver_root(s)); // GSL doc: ...the error can be estimated more accurately by taking the difference between the current iterate and next iterate rather than the previous iterate. 01665 tmp_solution.error = nominal_relative_accuracy; 01666 01667 unsigned int k=0; 01668 bool duplicate=false; 01669 while (k<solutions.size()) { 01670 if (fabs(solutions[k].value-tmp_solution.value) < (solutions[k].error+tmp_solution.error)) { 01671 duplicate= true; 01672 break; 01673 } 01674 ++k; 01675 } 01676 01677 if (!duplicate) { 01678 solutions.push_back(tmp_solution); 01679 } 01680 } 01681 } 01682 01683 gsl_root_fdfsolver_free(s); 01684 01685 if (0) { 01686 cerr << "solutions found: " << solutions.size() << endl; 01687 unsigned int k=0; 01688 while (k<solutions.size()) { 01689 cerr << k << ": " << solutions[k].value << " accuracy: " << solutions[k].error << endl; 01690 ++k; 01691 } 01692 } 01693 }
Here is the call graph for this function:
void print | ( | const Frame & | f | ) |
Definition at line 59 of file orsa_frame.cc.
References print(), Frame::size(), and UniverseTypeAwareTime::Time().
00059 { 00060 cout << "Frame time: " << f.Time() << endl; 00061 cout << "Frame size: " << f.size() << endl; 00062 unsigned int l; 00063 for (l=0;l<f.size();l++) print(f[l]); 00064 }
Here is the call graph for this function:
void print | ( | const Body & | b | ) |
Definition at line 472 of file orsa_body.cc.
References Body::mass(), Body::name(), Body::position(), Body::velocity(), Vector::x, Vector::y, and Vector::z.
Referenced by print().
00472 { 00473 cout << "body name: " << b.name() << " mass: " << b.mass() << endl; 00474 /* if (b.parent) { 00475 cout << "parent body: " << b.parent->name() << endl; 00476 } 00477 */ 00478 00479 cout << "position: " << b.position().x << " " << b.position().y << " " << b.position().z << endl; 00480 cout << "velocity: " << b.velocity().x << " " << b.velocity().y << " " << b.velocity().z << endl; 00481 }
Here is the call graph for this function:
Sky PropagatedSky_J2000 | ( | const OrbitWithEpoch & | orbit, | |
const UniverseTypeAwareTime & | final_time, | |||
const std::string & | obscode, | |||
const bool | integrate = true , |
|||
const bool | light_time_corrections = true | |||
) |
Definition at line 662 of file orsa_orbit.cc.
References Integrator::accuracy, Sky::Compute_J2000(), DAY, EARTH, EARTH_AND_MOON, OrbitWithEpoch::epoch, FromUnits(), GetC(), UniverseTypeAwareTime::GetDate(), UniverseTypeAwareTime::GetTime(), Evolution::Integrate(), JUPITER, Vector::Length(), location_file, MARS, MERCURY, NEPTUNE, LocationFile::ObsPos(), Evolution::push_back(), OrbitWithEpoch::RelativePosVel(), SATURN, UniverseTypeAwareTime::SetDate(), Evolution::SetIntegrator(), Evolution::SetInteraction(), Evolution::SetMaxUnsavedSubSteps(), Body::SetPosition(), Evolution::SetSamplePeriod(), SetupSolarSystem(), Body::SetVelocity(), Frame::size(), Evolution::size(), SUN, Integrator::timestep, URANUS, and VENUS.
Referenced by gauss_v_diff_f(), and least_sq_diff_f().
00662 { 00663 00664 if (!integrate) { 00665 00666 list<JPL_planets> l; 00667 // 00668 l.push_back(SUN); 00669 l.push_back(EARTH); 00670 00671 Frame frame; 00672 00673 SetupSolarSystem(frame,l,final_time); 00674 00675 /* 00676 LocationFile lf; 00677 lf.SetFileName(config->paths[OBSCODE]->GetValue().c_str()); 00678 lf.Read(); 00679 lf.Close(); 00680 */ 00681 00682 Vector geo = frame[1].position(); 00683 // geo += lf.ObsPos(obscode,final_time.GetDate()); 00684 geo += location_file->ObsPos(obscode,final_time.GetDate()); 00685 00686 Vector r,v; 00687 orbit.RelativePosVel(r,v,final_time); 00688 // 00689 r += frame[0].position(); 00690 v += frame[0].velocity(); 00691 00692 const Vector relative_position = r - geo; 00693 00694 Sky sky; 00695 00696 if (light_time_corrections) { 00697 sky.Compute_J2000(relative_position-v*relative_position.Length()/GetC()); 00698 } else { 00699 sky.Compute_J2000(relative_position); 00700 } 00701 00702 return sky; 00703 } 00704 00705 Frame frame; 00706 00707 list<JPL_planets> l; 00708 // 00709 l.push_back(SUN); 00710 l.push_back(MERCURY); 00711 l.push_back(VENUS); 00712 l.push_back(EARTH_AND_MOON); 00713 l.push_back(MARS); 00714 l.push_back(JUPITER); 00715 l.push_back(SATURN); 00716 l.push_back(URANUS); 00717 l.push_back(NEPTUNE); 00718 00719 SetupSolarSystem(frame,l,orbit.epoch); 00720 00721 /* 00722 Config conf; 00723 OrsaConfigFile ocf(&conf); 00724 ocf.Read(); 00725 ocf.Close(); 00726 */ 00727 00728 /* 00729 LocationFile lf; 00730 lf.SetFileName(config->paths[OBSCODE]->GetValue().c_str()); 00731 lf.Read(); 00732 lf.Close(); 00733 */ 00734 00735 { 00736 Body b("object",0.0); 00737 // RelativePosVel(b.position,b.velocity); 00738 Vector r,v; 00739 orbit.RelativePosVel(r,v); 00740 // 00741 // cerr << "PropagatedSky ---> RelativePosVel ---> r: " << Length(r) << endl; 00742 // 00743 // cerr << "b.position: " << b.position.x << " " << b.position.y << " " << b.position.z << endl; 00744 // add Sun postion and velocity 00745 r += frame[0].position(); 00746 v += frame[0].velocity(); 00747 // cerr << "b.position: " << b.position.x << " " << b.position.y << " " << b.position.z << endl; 00748 b.SetPosition(r); 00749 b.SetVelocity(v); 00750 frame.push_back(b); 00751 00752 // debug 00753 // cerr << "object r: " << Length(r) << endl; 00754 } 00755 00756 const Frame orbit_epoch_frame = frame; 00757 00758 UniverseTypeAwareTime obs_utat; 00759 00760 Radau15 itg; 00761 // itg.accuracy = 1.0e-12; 00762 itg.accuracy = 1.0e-12; 00763 itg.timestep = FromUnits(1.0,DAY); 00764 // 00765 Newton newton; 00766 00767 // Frame last_frame; 00768 // Vector relative_position; 00769 // Vector v; 00770 00771 Evolution evol; 00772 00773 // evol.interaction = &newton; 00774 evol.SetInteraction(&newton); 00775 // evol.interaction = &itr; 00776 // evol.integrator = &itg; 00777 evol.SetIntegrator(&itg); 00778 evol.push_back(orbit_epoch_frame); 00779 evol.SetMaxUnsavedSubSteps(100000); 00780 obs_utat.SetDate(final_time.GetDate()); 00781 // evol.sample_period = orbit_epoch_frame.Time() - obs_utat.Time(); // important, otherwise the final frame is not saved 00782 // evol.sample_period = (orbit_epoch_frame.Time() - obs_utat.Time())*(1.0-1.0e-8); // important, otherwise the final frame is not saved 00783 // evol.sample_period = (orbit_epoch_frame.Time() - obs_utat.Time()); 00784 evol.SetSamplePeriod(orbit_epoch_frame - obs_utat); 00785 evol.Integrate(obs_utat.GetTime(),true); 00786 Frame last_frame = evol[evol.size()-1]; 00787 00788 // cerr << "sample_period: " << evol.sample_period << endl; 00789 // fprintf(stderr,"final time JD: %20.8f last frame JD: %20.8f\n",final_time.GetDate().GetJulian(),last_frame.GetDate().GetJulian()); 00790 00791 // cerr << "..is Sun? -----> " << frame[0].name() << endl; 00792 // cerr << "..is Earth? --> " << frame[3].name() << endl; 00793 // 00794 Vector geo = last_frame[3].position()-last_frame[0].position(); // earth center position 00795 // Vector geo = frame[3].position; // earth center position 00796 00797 // IMPORTANT!!! 00798 geo += location_file->ObsPos(obscode,final_time.GetDate()); 00799 00800 // cerr << "obs. pos. length: " << Length(lf.ObsPos(obscode,final_time.GetDate())) << " " << LengthLabel() << endl; 00801 // cerr << "obs. pos. length: " << FromUnits(Length(lf.ObsPos(obscode,final_time.GetDate())),REARTH,-1) << " " << units.label(REARTH) << endl; 00802 00803 // cerr << "Length(geo): " << Length(geo) << endl; 00804 00805 const Vector relative_position = last_frame[last_frame.size()-1].position() - last_frame[0].position() - geo; // object - (earth + observer) 00806 const Vector v = last_frame[last_frame.size()-1].velocity(); 00807 00808 // cerr << "relative_position: " << relative_position.x << " " << relative_position.y << " " << relative_position.z << " " << endl; 00809 00810 // cerr << "relative distance: " << Length(relative_position) << " " << LengthLabel() << endl; 00811 00812 Sky sky; 00813 // sky.Compute(relative_position,obs_utat); 00814 // sky.Compute(relative_position,last_frame); 00815 // sky.Compute_J2000(relative_position); 00816 // 00817 if (light_time_corrections) { 00818 sky.Compute_J2000(relative_position-v*relative_position.Length()/GetC()); 00819 } else { 00820 sky.Compute_J2000(relative_position); 00821 } 00822 00823 // fprintf(stderr,"propagated sky: ra = %16.12f dec = %16.12f orbit JD: %14.5f final time JD: %14.5f\n",sky.ra().GetRad(),sky.dec().GetRad(),orbit.epoch.GetDate().GetJulian(),final_time.GetDate().GetJulian()); 00824 00825 // debug test 00826 /* 00827 { 00828 00829 // random number generator 00830 // const int random_seed = 124323; 00831 struct timeval tv; 00832 struct timezone tz; 00833 gettimeofday(&tv,&tz); 00834 const int random_seed = tv.tv_usec; 00835 gsl_rng *rnd; 00836 rnd = gsl_rng_alloc(gsl_rng_gfsr4); 00837 gsl_rng_set(rnd,random_seed); 00838 00839 Vector geo,rand; 00840 int oc; 00841 Sky sky; 00842 // double h,m,s,d,p; 00843 for (oc=1;oc<400;++oc) { 00844 // geo = last_frame[3].position + lf.ObsPos(oc,final_time.GetDate()); 00845 rand.Set(gsl_rng_uniform(rnd)*2-1, 00846 gsl_rng_uniform(rnd)*2-1, 00847 gsl_rng_uniform(rnd)*2-1); 00848 rand /= Length(rand); 00849 rand *= FromUnits(1.0,REARTH); 00850 geo = last_frame[3].position+rand; 00851 relative_position = last_frame[last_frame.size()-1].position - geo; 00852 sky.Compute(relative_position,obs_utat); 00853 // sky.ra().GetHMS(h,m,s); 00854 // printf("R.A.: %2f %2f %10.4f ",h,m,s); 00855 // sky.dec().GetDPS(d,p,s); 00856 // printf("dec.: %2f %2f %10.4f\n",d,p,s); 00857 printf("%20.12f %20.12f\n",sky.ra().GetRad(),sky.dec().GetRad()); 00858 } 00859 00860 gsl_rng_free(rnd); 00861 } 00862 */ 00863 00864 return sky; 00865 }
Here is the call graph for this function:
double orsa::psd_max | ( | const fftw_complex * | transformed_signal, | |
const int | size | |||
) |
Definition at line 375 of file orsa_fft.cc.
References norm(), and norm_sq().
00375 { 00376 00377 vector<double> psd; 00378 00379 int k; 00380 00381 psd.resize(size/2+1); 00382 psd[0] = norm(transformed_signal[0])/size; 00383 for (k=1;k<(size+1)/2;k++) 00384 psd[k] = sqrt(norm_sq(transformed_signal[k])+norm_sq(transformed_signal[size-k]))/size; 00385 if ( (size%2) == 0) // size is even 00386 psd[size/2] = norm(transformed_signal[size/2])/size; // Nyquist freq. 00387 00388 /* cerr << " === DUMP PSD ===" << endl; 00389 for (k=0;psd.size();++k) { 00390 printf("PSD[%05i]=%g\n",k,psd[k]); 00391 } 00392 */ 00393 00394 double maxamp=0; 00395 int bin=0; 00396 unsigned int l; 00397 // 00398 double found=false; 00399 if ( (psd[0] > psd[1]) && 00400 (psd[0] > maxamp) ) { 00401 maxamp = psd[0]; 00402 bin = 0; 00403 found = true; 00404 } 00405 // 00406 for (l=1;l<(psd.size()-1);l++) { 00407 if ( ( (psd[l] > psd[l-1]) && (psd[l] > psd[l+1]) ) && 00408 (psd[l] > maxamp) ) { 00409 maxamp = psd[l]; 00410 bin = l; 00411 found = true; 00412 } 00413 } 00414 00415 /* for (l=1;l<psd.size();l++) { 00416 if (psd[l] > maxamp) { 00417 maxamp = psd[l]; 00418 bin = l; 00419 } 00420 } 00421 */ 00422 00423 // DEBUG: DUMP DATA TO FILE 00424 if (0) { 00425 int filenum=0; 00426 char filename[20]; 00427 sprintf(filename,"dump_psd_%02i.dat",filenum); 00428 cerr << "data dump on file " << filename << endl; 00429 FILE *fp = fopen(filename,"w"); 00430 if (fp!=0) { 00431 for (l=0;l<psd.size();l++) { 00432 fprintf(fp,"%i %g\n",l,psd[l]); 00433 } 00434 } 00435 fclose(fp); 00436 filenum++; 00437 } 00438 00439 // cerr << " psd_max: returning " << ((double)bin/(double)size) << " amp: " << maxamp << " bin: " << bin << " found: " << found << endl; 00440 00441 if (found) return ((double)bin/(double)size); 00442 else { 00443 cerr << "\n *\n **\n ***\n ****\n*************\n************** WARNING!!! peak don't found... returning (-1)!!\n*************\n ****\n ***\n **\n *\n" << endl; 00444 return (-1); 00445 } 00446 }
Here is the call graph for this function:
double orsa::psd_max_again | ( | const fftw_complex * | transformed_signal, | |
const int | size | |||
) |
Definition at line 164 of file orsa_fft.cc.
References norm().
00164 { 00165 00166 vector<double> psd_plus((size-1)/2), psd_minus((size-1)/2); 00167 double Nyquist=0; 00168 00169 int k; 00170 00171 const double DC = norm(transformed_signal[0])/size; 00172 for (k=0;k<(size-1)/2;k++) 00173 psd_plus[k] = norm(transformed_signal[k+1])/size; 00174 for (k=0;k<(size-1)/2;k++) 00175 psd_minus[k] = norm(transformed_signal[size-(k+1)])/size; 00176 if ( (size%2) == 0) // size is even 00177 Nyquist = norm(transformed_signal[size/2])/size; // Nyquist freq. 00178 00179 // cerr << "...some test data: " << DC << " " << psd_plus[0] << " " << psd_plus[1] << " " << psd_plus[2] << " " << psd_plus[3] << " " << psd_plus[4] << endl; 00180 00181 // DEBUG: DUMP DATA TO FILE 00182 if (0) { 00183 static int filenum=0; 00184 char filename[20]; 00185 sprintf(filename,"dump_psd_%02i.dat",filenum); 00186 cerr << "data dump on file " << filename << endl; 00187 FILE *fp = fopen(filename,"w"); 00188 if (fp!=0) { 00189 int l; 00190 for (l=psd_minus.size()-1;l>=0;l--) fprintf(fp,"%i %g\n",-(l+1),psd_minus[l]); 00191 fprintf(fp,"%i %g\n",0,DC); 00192 for (l=0;l<(int)psd_plus.size(); l++) fprintf(fp,"%i %g\n", (l+1),psd_plus[l]); 00193 if ( (size%2) == 0) fprintf(fp,"%i %g\n",size/2,Nyquist); 00194 } 00195 fclose(fp); 00196 filenum++; 00197 } 00198 00199 // psd.resize(size); 00200 // 00201 /* psd[0] = norm(transformed_signal[0])/size; 00202 for (k=1;k<(size-1)/2;k++) 00203 psd[k] = norm(transformed_signal[k])/size; 00204 for (k=(size/2)+1;k<size;k++) 00205 psd[k] = norm(transformed_signal[size-k])/size; 00206 if ( (size%2) == 0) // size is even 00207 psd[size/2] = norm(transformed_signal[size/2])/size; // Nyquist freq. 00208 */ 00209 00210 double maxpow=DC; 00211 int bin=0; 00212 00213 for (k=1;k<(size-3)/2;k++) { 00214 if ( (psd_plus[k] > psd_plus[k-1]) && 00215 (psd_plus[k] > psd_plus[k+1]) && 00216 (psd_plus[k] > maxpow) ) { 00217 maxpow = psd_plus[k]; 00218 bin = k+1; 00219 } 00220 } 00221 00222 if ( (psd_plus[0] > DC) && 00223 (psd_plus[0] > psd_plus[1]) && 00224 (psd_plus[0] > maxpow) ) { 00225 maxpow = psd_plus[0]; 00226 bin = 1; 00227 } 00228 00229 for (k=1;k<(size-3)/2;k++) { 00230 if ( (psd_minus[k] > psd_minus[k-1]) && 00231 (psd_minus[k] > psd_minus[k+1]) && 00232 (psd_minus[k] > maxpow) ) { 00233 maxpow = psd_minus[k]; 00234 bin = -(k+1); 00235 } 00236 } 00237 00238 if ( (psd_minus[0] > DC) && 00239 (psd_minus[0] > psd_minus[1]) && 00240 (psd_minus[0] > maxpow) ) { 00241 maxpow = psd_minus[0]; 00242 bin = -1; 00243 } 00244 00245 if ( (DC > psd_plus[0]) && 00246 (DC > psd_minus[0]) && 00247 (DC > maxpow) ) { 00248 maxpow = DC; 00249 bin = 0; 00250 } 00251 00252 // TODO: Nyquist... 00253 00254 // cerr << " psd_max_again: returning " << ((double)bin/(double)size) << " amp: " << maxpow << " bin: " << bin << endl; 00255 00256 // DEBUG: DUMP DATA TO FILE 00257 if (0) { 00258 int filenum=0; 00259 char filename[20]; 00260 cerr << "DC: " << DC << endl; 00261 cerr << "size: " << size << endl; 00262 sprintf(filename,"dump_psd_%02i.dat",filenum); 00263 cerr << "data dump on file " << filename << endl; 00264 FILE *fp = fopen(filename,"w"); 00265 if (fp!=0) { 00266 int l; 00267 for (l=psd_minus.size()-1;l>=0;l--) fprintf(fp,"%i %g\n",-(l+1),psd_minus[l]); 00268 fprintf(fp,"%i %g\n",0,DC); 00269 for (l=0;l<(int)psd_plus.size(); l++) fprintf(fp,"%i %g\n", (l+1),psd_plus[l]); 00270 if ( (size%2) == 0) fprintf(fp,"%i %g\n",size/2,Nyquist); 00271 } 00272 fclose(fp); 00273 filenum++; 00274 } 00275 00276 if (maxpow>0) return ((double)bin/(double)size); 00277 /* else { 00278 cerr << "\n *\n **\n ***\n ****\n*************\n************** WARNING!!! peak don't found... returning (-1)!!\n*************\n ****\n ***\n **\n *\n" << endl; 00279 return (-1); 00280 } 00281 */ 00282 00283 return ((double)bin/(double)size); 00284 }
Here is the call graph for this function:
void orsa::psd_max_again_many | ( | const fftw_complex * | transformed_signal, | |
const int | size, | |||
vector< double > & | candidate, | |||
const unsigned int | nfreq | |||
) |
Definition at line 286 of file orsa_fft.cc.
References compare_binamp(), and norm().
00286 { 00287 00288 vector<double> psd_plus((size-1)/2), psd_minus((size-1)/2); 00289 double Nyquist=0; 00290 00291 int k; 00292 00293 const double DC = norm(transformed_signal[0])/size; 00294 for (k=0;k<(size-1)/2;k++) 00295 psd_plus[k] = norm(transformed_signal[k+1])/size; 00296 for (k=0;k<(size-1)/2;k++) 00297 psd_minus[k] = norm(transformed_signal[size-(k+1)])/size; 00298 if ( (size%2) == 0) // size is even 00299 Nyquist = norm(transformed_signal[size/2])/size; // Nyquist freq. 00300 00301 vector<binamp> vec_binamp; 00302 binamp ba; 00303 00304 for (k=1;k<(size-3)/2;k++) { 00305 if ( (psd_plus[k] > psd_plus[k-1]) && 00306 (psd_plus[k] > psd_plus[k+1]) ) { 00307 // candidate.push_back((double)(k+1)/(double)size); 00308 ba.bin = k+1; 00309 ba.amp = psd_plus[k]; 00310 vec_binamp.push_back(ba); 00311 } 00312 } 00313 00314 if ( (psd_plus[0] > DC) && 00315 (psd_plus[0] > psd_plus[1]) ) { 00316 // candidate.push_back(1.0/(double)size); 00317 ba.bin = 1; 00318 ba.amp = psd_plus[0]; 00319 vec_binamp.push_back(ba); 00320 } 00321 00322 for (k=1;k<(size-3)/2;k++) { 00323 if ( (psd_minus[k] > psd_minus[k-1]) && 00324 (psd_minus[k] > psd_minus[k+1]) ) { 00325 // candidate.push_back((double)-(k+1)/(double)size); 00326 ba.bin = -(k+1); 00327 ba.amp = psd_minus[k]; 00328 vec_binamp.push_back(ba); 00329 } 00330 } 00331 00332 if ( (psd_minus[0] > DC) && 00333 (psd_minus[0] > psd_minus[1]) ) { 00334 // candidate.push_back(-1.0/(double)size); 00335 ba.bin = -1; 00336 ba.amp = psd_minus[0]; 00337 vec_binamp.push_back(ba); 00338 } 00339 00340 if ( (DC > psd_plus[0]) && 00341 (DC > psd_minus[0]) ) { 00342 // candidate.push_back(0.0); 00343 ba.bin = 0; 00344 ba.amp = DC; 00345 vec_binamp.push_back(ba); 00346 } 00347 00348 00349 unsigned int bsize = vec_binamp.size(); 00350 binamp *ptr_binamp; 00351 ptr_binamp = (binamp *) malloc(bsize*sizeof(binamp)); 00352 { 00353 unsigned int k; 00354 for (k=0;k<bsize;k++) { 00355 ptr_binamp[k].bin = vec_binamp[k].bin; 00356 ptr_binamp[k].amp = vec_binamp[k].amp; 00357 } 00358 } 00359 gsl_heapsort(ptr_binamp,bsize,sizeof(binamp),(gsl_comparison_fn_t)compare_binamp); 00360 00361 candidate.clear(); 00362 00363 { 00364 unsigned int k; 00365 for (k=0;k<nfreq;k++) { 00366 candidate.push_back((double)ptr_binamp[k].bin/(double)size); 00367 } 00368 } 00369 00370 free(ptr_binamp); 00371 }
Here is the call graph for this function:
double radius | ( | const JPL_planets | p | ) |
Definition at line 507 of file orsa_file_jpl.cc.
References EARTH, FromUnits(), JUPITER, KM, MARS, MERCURY, MOON, NEPTUNE, PLUTO, SATURN, SUN, URANUS, and VENUS.
00507 { 00508 double r=0; 00509 switch(p) { 00510 case SUN: r = FromUnits(695000. ,KM); break; 00511 case MERCURY: r = FromUnits( 2440. ,KM); break; 00512 case VENUS: r = FromUnits( 6051.84,KM); break; 00513 case EARTH: r = FromUnits( 6371.01,KM); break; 00514 case MARS: r = FromUnits( 3389.92,KM); break; 00515 case JUPITER: r = FromUnits( 69911. ,KM); break; 00516 case SATURN: r = FromUnits( 58232. ,KM); break; 00517 case URANUS: r = FromUnits( 25362. ,KM); break; 00518 case NEPTUNE: r = FromUnits( 24624. ,KM); break; 00519 case PLUTO: r = FromUnits( 1151. ,KM); break; 00520 // 00521 case MOON: r = FromUnits( 1738. ,KM); break; 00522 // 00523 default: break; 00524 } 00525 return (r); 00526 }
Here is the call graph for this function:
Vector orsa::RelativisticBarycenter | ( | const std::vector< Body > & | ) |
Vector orsa::RelativisticBarycenter | ( | const vector< Body > & | f | ) |
Definition at line 156 of file orsa_frame.cc.
References modified_mu().
Referenced by Frame::RelativisticBarycenter().
00156 { 00157 Vector sum_vec(0,0,0); 00158 double sum_mu = 0.0; 00159 for (unsigned int i=0; i<f.size(); ++i) { 00160 const double mod_mu = modified_mu(f,i); 00161 if (mod_mu > 0.0) { 00162 sum_vec += mod_mu * f[i].position(); 00163 sum_mu += mod_mu; 00164 } 00165 } 00166 // 00167 return (sum_vec/sum_mu); 00168 }
Here is the call graph for this function:
Vector orsa::RelativisticBarycenterVelocity | ( | const std::vector< Body > & | ) |
Vector orsa::RelativisticBarycenterVelocity | ( | const vector< Body > & | f | ) |
Definition at line 170 of file orsa_frame.cc.
References modified_mu().
Referenced by Frame::RelativisticBarycenterVelocity().
00170 { 00171 Vector sum_vec(0,0,0); 00172 double sum_mu = 0.0; 00173 for (unsigned int i=0; i<f.size(); ++i) { 00174 const double mod_mu = modified_mu(f,i); 00175 if (mod_mu > 0.0) { 00176 sum_vec += mod_mu * f[i].velocity(); 00177 sum_mu += mod_mu; 00178 } 00179 } 00180 // 00181 return (sum_vec/sum_mu); 00182 }
Here is the call graph for this function:
void orsa::remove_leading_trailing_spaces | ( | std::string & | s | ) | [inline] |
Definition at line 505 of file orsa_file.h.
Referenced by TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), OrsaConfigFile::Read(), LocationFile::Read(), RWOFile::Read(), MPCObsFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), and AstorbFile::Read().
00505 { 00506 00507 const int first = s.find_first_not_of(" "); 00508 s.erase(0,first); 00509 00510 const int last = s.find_last_not_of(" "); 00511 s.erase(last+1,s.size()); 00512 }
double residual | ( | const Observation & | obs, | |
const OrbitWithEpoch & | orbit | |||
) |
The RMS of the residuals in arcsec units.
Definition at line 652 of file orsa_orbit.cc.
References Observation::date, Sky::delta_arcsec(), Observation::obscode, and OptimizedOrbitPositions::PropagatedSky_J2000().
00652 { 00653 OptimizedOrbitPositions opt(orbit); 00654 Sky sky = opt.PropagatedSky_J2000(obs.date,obs.obscode); 00655 return fabs(sky.delta_arcsec(obs)); 00656 }
Here is the call graph for this function:
double orsa::RMS_residuals | ( | const std::vector< Observation > & | , | |
const OrbitWithEpoch & | ||||
) |
Referenced by Compute_Gauss(), and E1().
double orsa::RMS_residuals | ( | const vector< Observation > & | obs, | |
const OrbitWithEpoch & | orbit | |||
) |
The RMS of the residuals in arcsec units.
Definition at line 620 of file orsa_orbit.cc.
References Sky::delta_arcsec(), OptimizedOrbitPositions::PropagatedSky_J2000(), and secure_sqrt().
00620 { 00621 00622 double sum_delta_arcsec_squared=0.0; 00623 unsigned int k=0; 00624 Sky sky; 00625 OptimizedOrbitPositions opt(orbit); 00626 while (k<obs.size()) { 00627 // sky = PropagatedSky(orbit,obs[k].date,obs[k].obscode); 00628 sky = opt.PropagatedSky_J2000(obs[k].date,obs[k].obscode); 00629 // sum_delta_arcsec_squared += secure_pow(PropagatedSky(orbit,obs[k].date,obs[k].obscode).delta_arcsec(obs[k]),2); 00630 // sum_delta_arcsec_squared += secure_pow(sky.delta_arcsec(obs[k]),2); 00631 const double delta_arcsec = sky.delta_arcsec(obs[k]); 00632 sum_delta_arcsec_squared += delta_arcsec*delta_arcsec; 00633 // 00634 /* if (k>0) { 00635 fprintf(stderr,"T-%i-%i: %g days\n",k,k-1,obs[k].date.GetJulian()-obs[k-1].date.GetJulian()); 00636 } 00637 */ 00638 00639 // debug 00640 /* 00641 printf("%16.12f %16.12f\n",obs[k].ra.GetRad(),obs[k].dec.GetRad()); 00642 printf("%16.12f %16.12f\n",sky.ra().GetRad(),sky.dec().GetRad()); 00643 */ 00644 // 00645 00646 ++k; 00647 } 00648 00649 return (secure_sqrt(sum_delta_arcsec_squared/obs.size())); 00650 }
Here is the call graph for this function:
void orsa::S1 | ( | const gsl_rng * | r, | |
void * | p, | |||
double | step_size | |||
) |
Definition at line 609 of file orsa_orbit_gsl.cc.
References d_a, d_e, d_i, d_M, d_omega_node, d_omega_pericenter, pi, and twopi.
00609 { 00610 00611 // double u = gsl_rng_uniform(r); 00612 // new_x = u * 2 * step_size - step_size + old_x; 00613 00614 // memcpy(xp, &new_x, sizeof(new_x)); 00615 00616 data *d = (data*) p; 00617 00618 // cerr << "d_a = " << d_a << endl; 00619 00620 /* 00621 data old_d; 00622 Orbit o = *d->orbit; old_d.orbit = &o; 00623 vector<Observation> obs = *d->obs; old_d.obs = &obs; 00624 */ 00625 00626 /* 00627 double u, total = step_size; 00628 u = 0.40*total*(gsl_rng_uniform(r)-0.5); total -= u; new_d.orbit->a += u; 00629 u = 0.05*total*(gsl_rng_uniform(r)-0.5); total -= u; new_d.orbit->e += u; 00630 u = 0.50*total*(gsl_rng_uniform(r)-0.5); total -= u; new_d.orbit->i += u; 00631 u = 1.00*total*(gsl_rng_uniform(r)-0.5); total -= u; new_d.orbit->omega_node += u; 00632 u = 1.00*total*(gsl_rng_uniform(r)-0.5); total -= u; new_d.orbit->omega_pericenter += u; 00633 u = 1.00*total*(gsl_rng_uniform(r)-0.5); total -= u; new_d.orbit->M += u; 00634 */ 00635 00636 // cerr << "requested step size: " << step_size << endl; 00637 00638 d->orbit->a += step_size*(gsl_rng_uniform(r)-0.5)*d_a; 00639 d->orbit->e += step_size*(gsl_rng_uniform(r)-0.5)*d_e; 00640 d->orbit->i += step_size*(gsl_rng_uniform(r)-0.5)*d_i; 00641 d->orbit->omega_node += step_size*(gsl_rng_uniform(r)-0.5)*d_omega_node; 00642 d->orbit->omega_pericenter += step_size*(gsl_rng_uniform(r)-0.5)*d_omega_pericenter; 00643 d->orbit->M += step_size*(gsl_rng_uniform(r)-0.5)*d_M; 00644 00645 // checks 00646 while (d->orbit->a < 0) d->orbit->a += 0.1*step_size*gsl_rng_uniform(r)*d_a; 00647 while (d->orbit->e < 0) d->orbit->e += 0.1*step_size*gsl_rng_uniform(r)*d_e; 00648 // 00649 do { d->orbit->i = fmod(d->orbit->i+pi,pi); } while (d->orbit->i < 0); 00650 do { d->orbit->omega_node = fmod(d->orbit->omega_node+twopi,twopi); } while (d->orbit->omega_node < 0); 00651 do { d->orbit->omega_pericenter = fmod(d->orbit->omega_pericenter+twopi,twopi); } while (d->orbit->omega_pericenter < 0); 00652 do { d->orbit->M = fmod(d->orbit->M+twopi,twopi); } while (d->orbit->M < 0); 00653 00654 // cerr << "proposed step size: " << M1(p,(void*)(&old_d)) << endl; 00655 00656 // *(d->orbit) = *(new_d.orbit); 00657 // *(d->obs) = *(new_d.obs); 00658 00659 }
void SearchCloseApproaches | ( | const Evolution * | evol, | |
const unsigned int | obj_index, | |||
const unsigned int | index, | |||
std::vector< CloseApproach > & | clapp, | |||
const double | distance_threshold, | |||
const double | time_accuracy = FromUnits(1,SECOND) | |||
) |
Definition at line 2746 of file orsa_orbit_gsl.cc.
References CloseApproaches_gsl_f(), CloseApproach::distance, CloseApproach::epoch, ExternalProduct(), UniverseTypeAwareTimeStep::GetDouble(), Evolution::GetSamplePeriod(), Vector::Length(), CloseApproach::name, Vector::Normalize(), Vector::Normalized(), CloseApproach::relative_velocity, UniverseTypeAwareTime::SetTime(), Evolution::size(), CloseApproach::tp_xy, CloseApproach::tp_z, Vector::x, Vector::y, and Vector::z.
02746 { 02747 02748 if (index == obj_index) return; 02749 02750 // gsl init 02751 CloseApproaches_gsl_parameters par; 02752 par.e = new Evolution(*evol); 02753 // 02754 gsl_multimin_fminimizer * s = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex,1); 02755 // 02756 gsl_multimin_function clapp_function; 02757 clapp_function.f = &CloseApproaches_gsl_f; 02758 clapp_function.n = 1; 02759 clapp_function.params = ∥ 02760 // 02761 gsl_vector * x = gsl_vector_alloc(1); 02762 gsl_vector * step_size = gsl_vector_alloc(1); 02763 02764 // cerr << "clapp [3]" << endl; 02765 02766 vector<double> distance; 02767 distance.resize(evol->size()); 02768 02769 CloseApproach ca; 02770 02771 for (unsigned int j=0;j<evol->size();++j) { 02772 distance[j] = ((*evol)[j][obj_index].position() - (*evol)[j][index].position()).Length(); 02773 } 02774 02775 // skip first and last 02776 for (unsigned int j=1;j<(evol->size()-1);++j) { 02777 02778 if ( (distance[j] <= 1.5*distance_threshold) && 02779 (distance[j] <= distance[j-1]) && 02780 (distance[j] <= distance[j+1]) ) { 02781 02782 // gsl stuff 02783 par.f = (*evol)[j]; 02784 par.obj_index = obj_index; 02785 par.index = index; 02786 02787 gsl_vector_set(x,0,(*evol)[j].GetTime()); 02788 02789 gsl_vector_set(step_size,0,evol->GetSamplePeriod().GetDouble()); 02790 02791 gsl_multimin_fminimizer_set(s, &clapp_function, x, step_size); 02792 02793 // cerr << "clapp [4]" << endl; 02794 02795 size_t iter = 0; 02796 int status; 02797 do { 02798 02799 ++iter; 02800 02801 // iter status 02802 status = gsl_multimin_fminimizer_iterate(s); 02803 02804 // convergence status 02805 // status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s),FromUnits(1,SECOND)); 02806 // status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s),FromUnits(10.0,MINUTE)); 02807 // status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s),FromUnits(1,HOUR)); 02808 status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s),time_accuracy); 02809 02810 } while (status == GSL_CONTINUE && iter < 200); 02811 02812 // const double minimum = gsl_multimin_fminimizer_minimum(s); 02813 02814 if (status == GSL_SUCCESS) { 02815 if (gsl_multimin_fminimizer_minimum(s) < distance_threshold) { 02816 ca.name = (*evol)[j][index].name(); 02817 ca.epoch.SetTime(gsl_vector_get(s->x,0)); 02818 ca.distance = gsl_multimin_fminimizer_minimum(s); 02819 02820 { 02821 Frame f = par.f; 02822 // 02823 const UniverseTypeAwareTime gsl_time(gsl_vector_get(s->x,0)); 02824 // 02825 par.e->clear(); 02826 par.e->push_back(f); 02827 par.e->SetSamplePeriod(f - gsl_time); 02828 par.e->SetMaxUnsavedSubSteps(10000); 02829 par.e->Integrate(gsl_time,true); 02830 // 02831 f = (*(par.e))[par.e->size()-1]; 02832 // 02833 ca.relative_velocity = (f[obj_index].velocity() - f[index].velocity()).Length(); 02834 02835 // target plane 02836 { 02837 const Vector unit_v = (f[obj_index].velocity() - f[index].velocity()).Normalize(); 02838 const Vector d = f[obj_index].position() - f[index].position(); 02839 const Vector unit_d = d.Normalized(); 02840 const Vector unit_q = ExternalProduct(unit_v,unit_d).Normalize(); 02841 02842 // d' = \alpha d + \beta q 02843 // and d'.z = 0.0 02844 02845 if (unit_d.z != 0.0) { 02846 const double beta = sqrt(+ unit_q.x*unit_q.x 02847 + unit_q.y*unit_q.y 02848 + unit_q.z*unit_q.z/(unit_d.z*unit_d.z)*(unit_d.x*unit_d.x+unit_d.y*unit_d.y) 02849 - 2*(unit_q.z/unit_d.z)*(unit_q.x*unit_d.y+unit_q.y*unit_d.x)); 02850 const double alpha = - beta * unit_q.z / unit_d.z; 02851 02852 const Vector new_unit_d = (alpha*unit_d+beta*unit_q).Normalize(); 02853 const Vector new_unit_q = ExternalProduct(unit_v,new_unit_d).Normalize(); 02854 02855 // this solves the ambiguity in beta, that has 2 solutions: +/- sqrt(...) 02856 if (new_unit_q.z > 0.0) { 02857 ca.tp_xy = d*new_unit_d; 02858 ca.tp_z = d*new_unit_q; 02859 } else { 02860 ca.tp_xy = -d*new_unit_d; 02861 ca.tp_z = -d*new_unit_q; 02862 } 02863 02864 /* 02865 { 02866 // debug 02867 cerr << "new_unit_d: " << new_unit_d.x << " " << new_unit_d.y << " " << new_unit_d.z << endl; 02868 cerr << "new_unit_q: " << new_unit_q.x << " " << new_unit_q.y << " " << new_unit_q.z << endl; 02869 cerr << "d*unit_v: " << d*unit_v << endl; 02870 cerr << "sqrt((d*new_unit_d)*(d*new_unit_d)+(d*new_unit_q)*(d*new_unit_q)): " << sqrt((d*new_unit_d)*(d*new_unit_d)+(d*new_unit_q)*(d*new_unit_q)) << endl; 02871 } 02872 */ 02873 02874 } else { 02875 ca.tp_xy = d.Length(); 02876 ca.tp_z = 0.0; 02877 } 02878 02879 // ca.tp_xy = d*unit_xy; 02880 // ca.tp_z = d*unit_z; 02881 02882 } 02883 02884 } 02885 // 02886 clapp.push_back(ca); 02887 } 02888 02889 } 02890 02891 // cerr << "clapp [5]" << endl; 02892 02893 } 02894 02895 } 02896 02897 // } 02898 02899 // gsl free 02900 gsl_multimin_fminimizer_free(s); 02901 gsl_vector_free(x); 02902 gsl_vector_free(step_size); 02903 02904 delete par.e; 02905 }
Here is the call graph for this function:
double secure_acos | ( | double | x | ) |
Definition at line 99 of file orsa_secure_math.cc.
References ORSA_DOMAIN_ERROR.
Referenced by Orbit::Compute(), and Sky::Compute_J2000().
00099 { 00100 if ((x>1.0) || (x<-1.0)) { 00101 // domain error 00102 ORSA_DOMAIN_ERROR("secure_acos(%g) is undefined!",x); 00103 return 1.0; // better value? 00104 } else { 00105 return acos(x); 00106 } 00107 }
double secure_asin | ( | double | x | ) |
Definition at line 88 of file orsa_secure_math.cc.
References ORSA_DOMAIN_ERROR.
00088 { 00089 if ((x>1.0) || (x<-1.0)) { 00090 // domain error 00091 ORSA_DOMAIN_ERROR("secure_asin(%g) is undefined!",x); 00092 return 1.0; // better value? 00093 } else { 00094 return asin(x); 00095 } 00096 }
double secure_atan2 | ( | double | x, | |
double | y | |||
) |
Definition at line 73 of file orsa_secure_math.cc.
References ORSA_DOMAIN_ERROR.
Referenced by amph(), Orbit::Compute(), Sky::Compute_J2000(), and AstDySMatrixFile::Read().
00073 { 00074 if (x==0.0) { 00075 if (y==0.0) { 00076 // domain error 00077 ORSA_DOMAIN_ERROR("secure_atan2(%g,%g) is undefined!",x,y); 00078 return 1.0; // better value? 00079 } else { 00080 return atan2(x,y); 00081 } 00082 } else { 00083 return atan2(x,y); 00084 } 00085 }
double secure_log | ( | double | x | ) |
Definition at line 53 of file orsa_secure_math.cc.
References ORSA_DOMAIN_ERROR.
Referenced by Orbit::Compute(), and Orbit::RelativePosVel().
00053 { 00054 if (x>0) { 00055 return log(x); 00056 } else { 00057 ORSA_DOMAIN_ERROR("secure_log(%g) is undefined!",x); 00058 return 1.0; // better value? 00059 } 00060 }
double secure_log10 | ( | double | x | ) |
Definition at line 63 of file orsa_secure_math.cc.
References ORSA_DOMAIN_ERROR.
00063 { 00064 if (x>0) { 00065 return log10(x); 00066 } else { 00067 ORSA_DOMAIN_ERROR("secure_log10(%g) is undefined!",x); 00068 return 1.0; // better value? 00069 } 00070 }
double secure_pow | ( | double | x, | |
double | y | |||
) |
Definition at line 38 of file orsa_secure_math.cc.
References ORSA_DOMAIN_ERROR.
Referenced by Compute_Gauss(), ComputeAcceleration(), Orbit::GetE(), M1(), OrbitDifferentialCorrectionsLeastSquares(), poly_8(), poly_8_df(), poly_8_fdf(), GalacticPotentialAllen::PotentialEnergy(), and Units::Recompute().
00038 { 00039 00040 if (x<0.0) { 00041 if (rint(y)!=y) { 00042 ORSA_DOMAIN_ERROR("secure_pow(%g,%g) is undefined!",x,y); 00043 return 1.0; // better value? 00044 } else { 00045 return pow(x,y); 00046 } 00047 } else { 00048 return pow(x,y); 00049 } 00050 }
double secure_sqrt | ( | double | x | ) |
Definition at line 110 of file orsa_secure_math.cc.
References ORSA_DOMAIN_ERROR.
Referenced by Orbit::Compute(), Orbit::Period(), Orbit::RelativePosVel(), and RMS_residuals().
00110 { 00111 if (x<0) { 00112 // domain error 00113 ORSA_DOMAIN_ERROR("secure_sqrt(%g) is undefined!",x); 00114 return sqrt(fabs(x)); // better value? 00115 } else { 00116 return sqrt(x); 00117 } 00118 }
void orsa::SetupSolarSystem | ( | Frame & | , | |
const std::list< JPL_planets > & | , | |||
const UniverseTypeAwareTime & | ||||
) |
void orsa::SetupSolarSystem | ( | Frame & | frame, | |
const list< JPL_planets > & | l, | |||
const UniverseTypeAwareTime & | t | |||
) |
The Sun is automatically included.
Definition at line 413 of file orsa_file_jpl.cc.
References EARTH, EARTH_AND_MOON, JPLFile::EphemEnd(), JPLFile::EphemStart(), JPLCache::GetJPLBody(), UniverseTypeAwareTime::GetTime(), jpl_cache, jpl_file, MOON, ORSA_WARNING, UniverseTypeAwareTime::SetTime(), and SUN.
Referenced by JPLPlanetsNewton::PotentialEnergy(), OptimizedOrbitPositions::PropagatedOrbit(), and PropagatedSky_J2000().
00413 { 00414 00415 frame.clear(); 00416 frame.SetTime(t); 00417 00418 /* 00419 string name; 00420 double mass; 00421 Vector r,v; 00422 */ 00423 00424 // date checks 00425 if (t < jpl_file->EphemStart()) { 00426 ORSA_WARNING("epoch requested is before ephem file start time! (%g < %g)",t.GetTime(),jpl_file->EphemStart().GetTime()); 00427 return; 00428 } 00429 // 00430 if (t > jpl_file->EphemEnd()) { 00431 ORSA_WARNING("epoch requested is after ephem file end time! (%g > %g)",t.GetTime(),jpl_file->EphemStart().GetTime()); 00432 return; 00433 } 00434 00435 // auto include the Sun 00436 // name="Sun"; SetMRV(mass,r,v,jpl_file,t,SUN); frame.push_back(Body(name,mass,r,v)); 00437 frame.push_back(jpl_cache->GetJPLBody(SUN,t)); 00438 00439 JPL_planets pl; 00440 list<JPL_planets>::const_iterator it = l.begin(); 00441 while (it != l.end()) { 00442 pl = *it; 00443 if (pl==SUN) { 00444 ++it; 00445 continue; // sun already included 00446 } else if (pl==EARTH_AND_MOON) { 00447 // pl = EARTH; name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl); frame.push_back(Body(name,mass,r,v)); 00448 frame.push_back(jpl_cache->GetJPLBody(EARTH,t)); 00449 // pl = MOON; name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl); frame.push_back(Body(name,mass,r,v)); 00450 frame.push_back(jpl_cache->GetJPLBody(MOON,t)); 00451 } else { 00452 // name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl); frame.push_back(Body(name,mass,r,v)); 00453 frame.push_back(jpl_cache->GetJPLBody(pl,t)); 00454 } 00455 ++it; 00456 } 00457 }
Here is the call graph for this function:
double orsa::sin | ( | const Angle & | alpha | ) | [inline] |
Definition at line 563 of file orsa_units.h.
References Angle::GetRad().
Referenced by alpha_delta_meridian(), alpha_delta_meridian_moon(), Orbit::Compute(), Compute_Gauss(), Compute_TestMethod(), dQ(), Orbit::GetE(), phi(), phi_Hanning(), Units::Recompute(), Orbit::RelativePosVel(), Vector::rotate(), and sincos().
00563 { 00564 return std::sin(alpha.GetRad()); 00565 }
Here is the call graph for this function:
void orsa::sincos | ( | const Angle & | alpha, | |
double & | s, | |||
double & | c | |||
) | [inline] |
Definition at line 575 of file orsa_units.h.
References cos(), Angle::GetRad(), and sin().
Referenced by Orbit::GetE(), Orbit::RelativePosVel(), Vector::rotate(), and unit_vector().
00575 { 00576 #ifdef HAVE_SINCOS 00577 ::sincos(alpha.GetRad(),&s,&c); 00578 #else // HAVE_SINCOS 00579 s = std::sin(alpha.GetRad()); 00580 c = std::cos(alpha.GetRad()); 00581 #endif // ! HAVE_SINCOS 00582 }
Here is the call graph for this function:
Frame orsa::StartFrame | ( | const std::vector< BodyWithEpoch > & | , | |
std::vector< JPL_planets > & | , | |||
const Interaction * | , | |||
const Integrator * | , | |||
const UniverseTypeAwareTime & | ||||
) |
A good frame to start an integration with.
Frame orsa::StartFrame | ( | const vector< BodyWithEpoch > & | f, | |
vector< JPL_planets > & | jf, | |||
const Interaction * | itr, | |||
const Integrator * | itg, | |||
const UniverseTypeAwareTime & | t | |||
) |
Definition at line 584 of file orsa_universe.cc.
References Evolution::clear(), FromUnits(), UniverseTypeAwareTime::GetDate(), Universe::GetUniverseType(), Evolution::Integrate(), Evolution::name, ORSA_WARNING, Evolution::push_back(), Real, SECOND, UniverseTypeAwareTime::SetDate(), Evolution::SetIntegrator(), Evolution::SetIntegratorTimeStep(), Evolution::SetInteraction(), Evolution::SetSamplePeriod(), UniverseTypeAwareTime::SetTime(), Simulated, Evolution::size(), Frame::size(), UniverseTypeAwareTime::Time(), Integrator::timestep, and universe.
00584 { 00585 // Frame StartFrame(const vector<BodyWithEpoch> & f, vector<JPL_planets> & jf, InteractionType interaction_type, IntegratorType integrator_type, const UniverseTypeAwareTime & t) { 00586 00587 // const double original_timestep = itg->timestep; 00588 const UniverseTypeAwareTimeStep original_timestep = itg->timestep; 00589 00590 Frame frame_start; 00591 00592 switch (universe->GetUniverseType()) { 00593 case Real: { 00594 frame_start.SetDate(t.GetDate()); 00595 Frame running_frame; 00596 Evolution running_evolution; 00597 // running_evolution.interaction = itr->clone(); 00598 // running_evolution.integrator = itg->clone(); 00599 running_evolution.SetIntegrator(itg); 00600 running_evolution.SetInteraction(itr); 00601 00602 // put the JPL bodies into the frame_start 00603 /* 00604 if (jf.size()) { 00605 Date ddf = t.GetDate(); 00606 unsigned int k; 00607 for (k=0;k<jf.size();++k) { 00608 jf[k].SetEpoch(ddf); 00609 Body b(jf[k].name(), 00610 jf[k].mass(), 00611 jf[k].radius(), 00612 jf[k].position(), 00613 jf[k].velocity()); 00614 frame_start.push_back(b); 00615 } 00616 } 00617 */ 00618 // 00619 if (jf.size()) { 00620 unsigned int k; 00621 for (k=0;k<jf.size();++k) { 00622 frame_start.push_back(JPLBody(jf[k],t.GetDate())); 00623 } 00624 } 00625 00626 if (f.size()) { 00627 00628 Frame last_frame; 00629 unsigned int base_bodies, l; 00630 00631 unsigned int j=0; 00632 while (j < f.size()) { 00633 00634 running_frame.clear(); 00635 // running_frame.SetDate(f[j].epoch.GetDate()); 00636 running_frame.SetDate(f[j].GetEpoch().GetDate()); 00637 00638 running_frame.push_back(f[j]); 00639 ++j; 00640 00641 // while ((j<f.size()) && (f[j].epoch.GetDate() == running_frame.GetDate())) { 00642 while ((j<f.size()) && (f[j].Epoch().GetDate() == running_frame.GetDate())) { 00643 running_frame.push_back(f[j]); 00644 ++j; 00645 } 00646 00647 base_bodies = running_frame.size(); 00648 00649 // JPL 00650 /* 00651 Date ddt = running_frame.GetDate(); 00652 unsigned int k; 00653 for (k=0;k<jf.size();++k) { 00654 jf[k].SetEpoch(ddt); 00655 Body b(jf[k].name(), 00656 jf[k].mass(), 00657 jf[k].radius(), 00658 jf[k].position(), 00659 jf[k].velocity()); 00660 running_frame.push_back(b); 00661 } 00662 */ 00663 { 00664 unsigned int k; 00665 for (k=0;k<jf.size();++k) { 00666 running_frame.push_back(JPLBody(jf[k],running_frame.GetDate())); 00667 } 00668 } 00669 00670 running_evolution.clear(); 00671 running_evolution.push_back(running_frame); 00672 // running_evolution.sample_period = fabs(running_frame.GetTime() - frame_start.GetTime()); 00673 running_evolution.SetSamplePeriod((running_frame - frame_start).absolute()); 00674 running_evolution.SetIntegratorTimeStep(original_timestep); 00675 running_evolution.Integrate(t); 00676 00677 last_frame = running_evolution[running_evolution.size()-1]; 00678 if (fabs(last_frame.Time()-frame_start.Time()) > FromUnits(1.0e-3,SECOND)) { 00679 /* 00680 fprintf(stderr,"!!! last_frame.Time() != frame_start.Time() --->> %30.20f != %30.20f\n",last_frame.Time(),frame_start.Time()); 00681 cerr << "|dT| = " << FromUnits(fabs(last_frame.Time()-frame_start.Time()),SECOND,-1) << " seconds" << endl; 00682 print(last_frame); 00683 */ 00684 ORSA_WARNING( 00685 "last_frame.Time() != frame_start.Time() --->> %30.20f != %30.20f\n" 00686 "|dT| = %g seconds\n", 00687 last_frame.Time(),frame_start.Time(),FromUnits(fabs(last_frame.Time()-frame_start.Time()),SECOND,-1)); 00688 ORSA_WARNING(" objects in frame: "); 00689 unsigned int k; 00690 for (k=0;k<last_frame.size();++k) { 00691 ORSA_WARNING(" %s", last_frame[k].name().c_str()); 00692 } 00693 } 00694 00695 for (l=0;l<base_bodies;++l) { 00696 frame_start.push_back(last_frame[l]); 00697 } 00698 } 00699 } 00700 break; 00701 } 00702 case Simulated: { 00703 frame_start.SetTime(0.0); 00704 unsigned int k; 00705 for (k=0;k<f.size();++k) { 00706 frame_start.push_back(f[k]); 00707 } 00708 break; 00709 } 00710 } 00711 00712 return frame_start; 00713 }
Here is the call graph for this function:
static void orsa::swap | ( | void * | ptr, | |
unsigned int | size | |||
) | [static] |
Definition at line 1511 of file orsa_file.cc.
References ORSA_WARNING, swap_4(), and swap_8().
Referenced by OrsaFile::Read().
01511 { 01512 switch(size) { 01513 case 4: swap_4(ptr); break; 01514 case 8: swap_8(ptr); break; 01515 default: 01516 ORSA_WARNING("called read_swap with size = %i",size); 01517 break; 01518 } 01519 }
Here is the call graph for this function:
static void orsa::swap_4 | ( | void * | ptr | ) | [static] |
Definition at line 1495 of file orsa_file.cc.
References SWAP_MACRO.
Referenced by swap().
01495 { 01496 char *tptr = (char *)ptr, tchar; 01497 01498 SWAP_MACRO( tptr[0], tptr[3], tchar); 01499 SWAP_MACRO( tptr[1], tptr[2], tchar); 01500 }
static void orsa::swap_8 | ( | void * | ptr | ) | [static] |
Definition at line 1502 of file orsa_file.cc.
References SWAP_MACRO.
Referenced by swap().
01502 { 01503 char *tptr = (char *)ptr, tchar; 01504 01505 SWAP_MACRO( tptr[0], tptr[7], tchar); 01506 SWAP_MACRO( tptr[1], tptr[6], tchar); 01507 SWAP_MACRO( tptr[2], tptr[5], tchar); 01508 SWAP_MACRO( tptr[3], tptr[4], tchar); 01509 }
int orsa::SWIFTRawReadBinaryFile | ( | FILE_TYPE | file, | |
const int | version = 2 | |||
) |
Definition at line 1167 of file orsa_file.cc.
References el, file_time, FromUnits(), l_ts, nast, READ_FILE, w_ts, and YEAR.
Referenced by SWIFTFile::AsteroidsInFile(), and SWIFTFile::Read().
01167 { 01168 01169 // version == 1 --> not filetered data 01170 // version == 2 --> filtered data 01171 01172 int good = 0; 01173 double dummy; 01174 01175 if ( version == 1 ) { 01176 READ_FILE(&dummy, sizeof(int), 1,file); 01177 READ_FILE(&nast, sizeof(int), 1,file); 01178 READ_FILE(&el, sizeof(double), 6,file); 01179 READ_FILE(&file_time, sizeof(double), 1,file); 01180 good = READ_FILE(&dummy, sizeof(int), 1,file); 01181 } else if ( version == 2 ) { 01182 READ_FILE(&dummy, sizeof(int), 1,file); 01183 READ_FILE(&nast, sizeof(int), 1,file); 01184 READ_FILE(&el, sizeof(double), 6,file); 01185 READ_FILE(&l_ts, sizeof(double), 1,file); 01186 READ_FILE(&w_ts, sizeof(double), 1,file); 01187 READ_FILE(&file_time, sizeof(double), 1,file); 01188 good = READ_FILE(&dummy, sizeof(int), 1,file); 01189 } 01190 01191 // Units conversions 01192 file_time = FromUnits(file_time,YEAR); 01193 01194 return (good); 01195 }
Here is the call graph for this function:
double orsa::tan | ( | const Angle & | alpha | ) | [inline] |
Definition at line 571 of file orsa_units.h.
References Angle::GetRad().
Referenced by Orbit::Compute().
00571 { 00572 return std::tan(alpha.GetRad()); 00573 }
Here is the call graph for this function:
string TimeLabel | ( | ) |
Definition at line 139 of file orsa_universe.cc.
References Units::TimeLabel(), and units.
Referenced by OrsaFile::Read().
00139 { return orsa::units->TimeLabel(); };
Here is the call graph for this function:
std::string TimeScaleLabel | ( | TimeScale | ts | ) |
Definition at line 896 of file orsa_units.cc.
References ET, GPS, TAI, TDT, UT, UT1, and UTC.
00896 { 00897 if (ts==UTC) return "UTC"; 00898 if (ts==UT) return "UT"; 00899 if (ts==UT1) return "UT1"; 00900 if (ts==TAI) return "TAI"; 00901 if (ts==TDT) return "TDT"; 00902 if (ts==ET) return "ET"; 00903 if (ts==GPS) return "GPS"; 00904 return ""; 00905 }
static Vector orsa::unit_vector | ( | const Angle & | ra, | |
const Angle & | dec | |||
) | [static] |
Definition at line 1333 of file orsa_orbit.cc.
References Sky::dec(), Sky::ra(), and sincos().
Referenced by Sky::delta_arcsec().
01333 { 01334 double sin_ra, cos_ra; 01335 sincos(ra,sin_ra,cos_ra); 01336 // 01337 double sin_dec, cos_dec; 01338 sincos(dec,sin_dec,cos_dec); 01339 // 01340 return Vector(cos_dec*sin_ra, 01341 cos_dec*cos_ra, 01342 sin_dec); 01343 }
Here is the call graph for this function:
double orsa::Weight | ( | vector< Observation > & | obs, | |
const double | optimal_dt = FromUnits(1,DAY) | |||
) |
Definition at line 89 of file orsa_orbit_gsl.cc.
00089 { 00090 if (obs.size() < 2) return 0; // error... 00091 double w=0; 00092 sort(obs.begin(),obs.end()); 00093 unsigned int k=1; 00094 while (k < obs.size()) { 00095 w += Weight(obs[k-1],obs[k],optimal_dt); 00096 ++k; 00097 } 00098 return w; 00099 }
double orsa::Weight | ( | const Observation & | obs1, | |
const Observation & | obs2, | |||
const double | optimal_dt = FromUnits(1,DAY) | |||
) |
Definition at line 69 of file orsa_orbit_gsl.cc.
References Observation::date, DAY, FromUnits(), and Date::GetJulian().
00069 { 00070 double w; 00071 double dt = fabs(FromUnits(obs1.date.GetJulian()-obs2.date.GetJulian(),DAY)); 00072 /* 00073 if (dt<optimal_dt) { 00074 w = optimal_dt/dt; 00075 } else { 00076 w = 1+dt/optimal_dt; 00077 } 00078 */ 00079 // 00080 if (dt<optimal_dt) { 00081 w = dt/optimal_dt; 00082 } else { 00083 w = 1+optimal_dt/dt; 00084 } 00085 // 00086 return w; 00087 }
Here is the call graph for this function:
const double AU_MKS = 1.49597870660e11 |
Definition at line 50 of file orsa_units.cc.
const double c_MKS = 299792458.0 |
Definition at line 51 of file orsa_units.cc.
const double c_omega_node = secure_pow(1.0/d_omega_node, 2) |
const double c_omega_pericenter = secure_pow(1.0/d_omega_pericenter, 2) |
The active configuration.
Definition at line 45 of file orsa_universe.cc.
Referenced by Compute_TestMethod(), OrsaConfigFile::Read(), and OrsaConfigFile::Write().
The active configuration.
Definition at line 45 of file orsa_universe.cc.
Referenced by Compute_TestMethod(), OrsaConfigFile::Read(), and OrsaConfigFile::Write().
const double d_a = 0.1 |
const double d_e = 0.01 |
const double d_omega_node = 1.0*(pi/180) |
const double d_omega_pericenter = 1.0*(pi/180) |
Definition at line 218 of file orsa_units.cc.
Definition at line 218 of file orsa_units.cc.
double el[6] |
Definition at line 1165 of file orsa_file.cc.
Referenced by SWIFTFile::Read(), and SWIFTRawReadBinaryFile().
const ET_minus_UT_element ET_minus_UT_table[] |
const ET_minus_UT_element ET_minus_UT_table_final_element = {0,0,0,0} |
double file_time |
Definition at line 1165 of file orsa_file.cc.
Referenced by SWIFTFile::Read(), and SWIFTRawReadBinaryFile().
const double G_MKS = 6.67259e-11 |
const double halfpi = PI/2 |
Definition at line 47 of file orsa_common.h.
Definition at line 51 of file orsa_universe.cc.
Referenced by Compute_Gauss(), and SetupSolarSystem().
Definition at line 51 of file orsa_universe.cc.
Referenced by Compute_Gauss(), and SetupSolarSystem().
Definition at line 48 of file orsa_universe.cc.
Referenced by JPLBody::JPLBody(), local_C22(), local_C31(), local_C32(), local_C33(), local_C41(), local_C42(), local_C43(), local_C44(), local_J2(), local_J3(), local_J4(), local_mass(), local_S31(), local_S32(), local_S33(), local_S41(), local_S42(), local_S43(), local_S44(), JPLBody::SetEpoch(), and SetupSolarSystem().
Definition at line 48 of file orsa_universe.cc.
Referenced by JPLBody::JPLBody(), local_C22(), local_C31(), local_C32(), local_C33(), local_C41(), local_C42(), local_C43(), local_C44(), local_J2(), local_J3(), local_J4(), local_mass(), local_S31(), local_S32(), local_S33(), local_S41(), local_S42(), local_S43(), local_S44(), JPLBody::SetEpoch(), and SetupSolarSystem().
double l_ts |
Definition at line 1165 of file orsa_file.cc.
Referenced by SWIFTFile::Read(), and SWIFTRawReadBinaryFile().
Definition at line 54 of file orsa_universe.cc.
Referenced by Compute_Gauss(), Compute_TestMethod(), and PropagatedSky_J2000().
Definition at line 54 of file orsa_universe.cc.
Referenced by Compute_Gauss(), Compute_TestMethod(), and PropagatedSky_J2000().
const double MEARTH_MKS = 5.9742e24 |
Definition at line 48 of file orsa_units.cc.
const double MJUPITER_MKS = 1.8989e27 |
Definition at line 47 of file orsa_units.cc.
const double MMOON_MKS = 7.3483e22 |
Definition at line 49 of file orsa_units.cc.
const double MSUN_MKS = 1.9889e30 |
Definition at line 46 of file orsa_units.cc.
int nast |
Definition at line 1164 of file orsa_file.cc.
Referenced by SWIFTFile::AsteroidsInFile(), SWIFTFile::Read(), and SWIFTRawReadBinaryFile().
OrsaPaths* orsa_paths = 0 |
const gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE, K, T_INITIAL, MU_T, T_MIN} |
Definition at line 566 of file orsa_orbit_gsl.cc.
Referenced by Compute_Gauss(), and poly_8_gsl_solve().
const double pi = PI |
Definition at line 48 of file orsa_common.h.
Referenced by Compute(), Orbit::Compute(), Compute_Gauss(), Sky::Compute_J2000(), Compute_TestMethod(), Sky::delta_arcsec(), dQ(), gauss_v_f(), Angle::GetDPS(), Orbit::GetE(), Angle::GetHMS(), least_sq_f(), MOID(), MOID2RB(), OrbitDifferentialCorrectionsLeastSquares(), Mercury5IntegrationFile::Read(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), LocationFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), Units::Recompute(), S1(), Angle::SetDPS(), and Angle::SetHMS().
const double pisq = PI*PI |
Definition at line 50 of file orsa_common.h.
Referenced by dQ(), Orbit::Period(), and TLEFile::Read().
const double R_EARTH_MKS = 6378140.0 |
Definition at line 52 of file orsa_units.cc.
const double R_MOON_MKS = 1737400.0 |
Definition at line 53 of file orsa_units.cc.
const TAI_minus_UTC_element TAI_minus_UTC_table[] |
Initial value:
{ {1,1,1972,10}, {1,7,1972,11}, {1,1,1973,12}, {1,1,1974,13}, {1,1,1975,14}, {1,1,1976,15}, {1,1,1977,16}, {1,1,1978,17}, {1,1,1979,18}, {1,1,1980,19}, {1,7,1981,20}, {1,7,1982,21}, {1,7,1983,22}, {1,7,1985,23}, {1,1,1988,24}, {1,1,1990,25}, {1,1,1991,26}, {1,7,1992,27}, {1,7,1993,28}, {1,7,1994,29}, {1,1,1996,30}, {1,7,1997,31}, {1,1,1999,32}, {0,0,0,0} }
Definition at line 248 of file orsa_units.cc.
Referenced by Date::delta_seconds().
const TAI_minus_UTC_element TAI_minus_UTC_table_final_element = {0,0,0,0} |
const double twopi = PI+PI |
Definition at line 49 of file orsa_common.h.
Referenced by apply_window(), Orbit::Compute(), Sky::Compute_J2000(), Orbit::GetE(), MOID(), MOID2RB(), phi(), phi_Hanning(), OptimizedOrbitPositions::PropagatedOrbit(), JPLDastcomCometFile::Read(), AstDySMatrixFile::Read(), MPCCometFile::Read(), OrbitWithEpoch::RelativePosVel(), Orbit::RelativePosVel(), and S1().
Definition at line 39 of file orsa_universe.cc.
Referenced by FromUnits(), GetC(), GetG(), GetG_MKS(), GetMSun(), LengthLabel(), MassLabel(), TimeLabel(), and OrsaFile::Write().
Definition at line 39 of file orsa_universe.cc.
Referenced by FromUnits(), GetC(), GetG(), GetG_MKS(), GetMSun(), LengthLabel(), MassLabel(), TimeLabel(), and OrsaFile::Write().
The active universe.
Definition at line 42 of file orsa_universe.cc.
Referenced by UniverseTypeAwareTimeStep::absolute(), Compute(), Compute_Gauss(), Sky::Compute_J2000(), Compute_TestMethod(), Frame::ForceJPLEphemerisData(), UniverseTypeAwareTimeStep::GetDouble(), JPLFile::GetEph(), UniverseTypeAwareTime::GetTime(), Evolution::Integrate(), UniverseTypeAwareTimeStep::IsZero(), JPLPlanetsNewton::JPLPlanetsNewton(), UniverseTypeAwareTime::operator+(), UniverseTypeAwareTimeStep::operator+(), UniverseTypeAwareTime::operator+=(), UniverseTypeAwareTime::operator-(), UniverseTypeAwareTimeStep::operator-(), UniverseTypeAwareTime::operator-=(), UniverseTypeAwareTime::operator<(), UniverseTypeAwareTimeStep::operator<(), UniverseTypeAwareTime::operator<=(), UniverseTypeAwareTime::operator==(), UniverseTypeAwareTimeStep::operator==(), UniverseTypeAwareTime::operator>(), UniverseTypeAwareTimeStep::operator>(), UniverseTypeAwareTime::operator>=(), AstorbFile::Read(), MPCOrbFile::Read(), MPCCometFile::Read(), OrsaFile::Read(), AstDySMatrixFile::Read(), JPLDastcomNumFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomCometFile::Read(), NEODYSCAT::Read(), TLEFile::Read(), StartFrame(), OrsaFile::Write(), and Universe::~Universe().
The active universe.
Definition at line 42 of file orsa_universe.cc.
Referenced by UniverseTypeAwareTimeStep::absolute(), Compute(), Compute_Gauss(), Sky::Compute_J2000(), Compute_TestMethod(), Frame::ForceJPLEphemerisData(), UniverseTypeAwareTimeStep::GetDouble(), JPLFile::GetEph(), UniverseTypeAwareTime::GetTime(), Evolution::Integrate(), UniverseTypeAwareTimeStep::IsZero(), JPLPlanetsNewton::JPLPlanetsNewton(), UniverseTypeAwareTimeStep::operator+(), UniverseTypeAwareTime::operator+(), UniverseTypeAwareTime::operator+=(), UniverseTypeAwareTimeStep::operator-(), UniverseTypeAwareTime::operator-(), UniverseTypeAwareTime::operator-=(), UniverseTypeAwareTimeStep::operator<(), UniverseTypeAwareTime::operator<(), UniverseTypeAwareTime::operator<=(), UniverseTypeAwareTimeStep::operator==(), UniverseTypeAwareTime::operator==(), UniverseTypeAwareTimeStep::operator>(), UniverseTypeAwareTime::operator>(), UniverseTypeAwareTime::operator>=(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), OrsaFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), StartFrame(), OrsaFile::Write(), and Universe::~Universe().
double w_ts |