orsa Namespace Reference


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)
double KineticEnergy (const Body &)
string Label (const ConfigEnum e)
std::string Label (const ConfigEnum)
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)
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 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 &center)
Vector BarycentricAngularMomentum (const vector< Body > &f)
double KineticEnergy (const vector< Body > &f)
void make_new_integrator (Integrator **i, const IntegratorType type)
void convert (IntegratorType &it, const unsigned int i)
std::string label (const IntegratorType it)
void make_new_integrator (Integrator **, const IntegratorType)
std::string label (const InteractionType it)
void make_new_interaction (Interaction **i, const InteractionType type)
void convert (InteractionType &it, const unsigned int i)
void make_new_interaction (Interaction **, const InteractionType)
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 &obs, const OrbitWithEpoch &orbit)
 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 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 &o1_in, const Orbit &o2_in, Vector &r1, Vector &r2)
 Minimal Orbit Intersection Distance between two orbits.
double MOID2RB_f (const gsl_vector *v, void *params)
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.
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)
 General interface to the Orbit computation from a set of observations.
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 &params, 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)
void SearchCloseApproaches (const Evolution *, const unsigned int, const unsigned int, std::vector< CloseApproach > &, const double, const double=FromUnits(1, SECOND))
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)
std::string TimeLabel ()
std::string LengthLabel ()
std::string MassLabel ()
void convert (TimeScale &ts, const unsigned int i)
UniverseTypeAwareTimeStep operator * (const int, const UniverseTypeAwareTimeStep &)
UniverseTypeAwareTimeStep operator * (const UniverseTypeAwareTimeStep &, const int)
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)
void alpha_delta_meridian (const JPL_planets, const Date &, Angle &alpha_zero, Angle &delta_zero, Angle &W)
double FromUnits (const double x, const time_unit u, const int power)
double FromUnits (const double x, const length_unit u, const int power)
double FromUnits (const double x, const mass_unit u, const int power)
Frame StartFrame (const vector< BodyWithEpoch > &f, vector< JPL_planets > &jf, const Interaction *itr, const Integrator *itg, const UniverseTypeAwareTime &t)
 A good frame to start an integration with.
void convert (UniverseType &ut, const unsigned int i)

Variables

const double halfpi = PI/2
const double pi = PI
const double twopi = PI+PI
const double pisq = PI*PI
Configconfig = 0
 The active configuration.
int nast
double el [6]
double l_ts
double w_ts
double file_time
OrsaPathsorsa_paths = 0
LocationFilelocation_file = 0
JPLFilejpl_file = 0
JPLCachejpl_cache = 0
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 []
Unitsunits = 0
Universeuniverse = 0
 The active universe.


Typedef Documentation

typedef struct orsa::gsl_d1_minimization_parameters gsl_d1_minimization_parameters
 


Enumeration Type Documentation

enum ConfigEnum
 

Enumeration values:
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 

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   };

enum CovarianceMatrixElements
 

Enumeration values:
Osculating 
Equinoctal 

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
 

Enumeration values:
algo_FFT 
algo_FFTB 
algo_MFT 
algo_FMFT1 
algo_FMFT2 

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
 

signal types

Enumeration values:
HK 
D 
PQ 
NODE 
ANOMALY 
ANOMALY_PHASE 
A_M 
TESTING 

Definition at line 96 of file orsa_fft.h.

enum FFTSearchAmplitude
 

Enumeration values:
A 
E 
I 
SIN_I 
TAN_I_2 
ONE 

Definition at line 98 of file orsa_fft.h.

00098 {A,E,I,SIN_I,TAN_I_2,ONE};

enum FFTSearchPhase
 

Enumeration values:
OMEGA_NODE 
OMEGA_PERICENTER 
OMEGA_TILDE 
MM 
LAMBDA 
ZERO 

Definition at line 99 of file orsa_fft.h.

enum FILE_STATUS
 

Enumeration values:
CLOSE 
OPEN_R 
OPEN_W 

Definition at line 81 of file orsa_file.h.

00081 {CLOSE,OPEN_R,OPEN_W};

enum IntegratorType
 

Enumeration values:
STOER 
BULIRSCHSTOER 
RUNGEKUTTA 
DISSIPATIVERUNGEKUTTA 
RA15 
LEAPFROG 

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
 

Enumeration values:
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
 

Enumeration values:
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
 

Enumeration values:
MPARSEC 
KPARSEC 
PARSEC 
LY 
AU 
EARTHMOON 
REARTH 
RMOON 
KM 
M 
CM 
LD 
ER 
MR 

Definition at line 60 of file orsa_units.h.

00060                    {
00061     MPARSEC=1,
00062     KPARSEC=2,
00063     PARSEC=3,
00064     LY=4,
00065     AU=5,
00066     EARTHMOON=6,
00067     REARTH=7,
00068     RMOON=8,
00069     KM=9,
00070     M=10,
00071     CM=11,
00072     // aliases
00073     LD=EARTHMOON,
00074     ER=REARTH,
00075     MR=RMOON
00076   };

enum M5COLS
 

Enumeration values:
C7 
C10 

Definition at line 146 of file orsa_file.h.

00146 {C7,C10};

enum mass_unit
 

Enumeration values:
MSUN 
MJUPITER 
MEARTH 
MMOON 
KG 
GRAM 

Definition at line 98 of file orsa_units.h.

00098                  {
00099     MSUN=1,
00100     MJUPITER=2,
00101     MEARTH=3,
00102     MMOON=4,
00103     KG=5,
00104     GRAM=6
00105   };

enum OrsaFileDataType
 

Enumeration values:
OFDT_END_OF_FILE 
OFDT_UNIVERSE 
OFDT_EVOLUTION 
OFDT_FRAME 
OFDT_BODY 

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
 

Enumeration values:
EQUATORIAL 
ECLIPTIC 

Definition at line 585 of file orsa_units.h.

00585                        {
00586     EQUATORIAL=1,
00587     ECLIPTIC=2
00588   };

enum time_unit
 

Enumeration values:
YEAR 
DAY 
HOUR 
MINUTE 
SECOND 

Definition at line 38 of file orsa_units.h.

00038                  {
00039     YEAR=1,
00040     DAY=2,
00041     HOUR=3,
00042     MINUTE=4,
00043     SECOND=5
00044   };

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.

Enumeration values:
UTC 
UT 
TAI 
TDT 
GPS 
UT1 
ET 
TT 

Definition at line 250 of file orsa_units.h.

00250                  {
00251     UTC=1,
00252     UT=2,
00253     TAI=3,
00254     TDT=4,
00255     GPS=5,
00256     // aliases
00257     UT1=UT,
00258     ET=TDT,
00259     TT=TDT
00260   };

enum UniverseType
 

This enum is used to classify the Universes: a Real Universe is composed.

Enumeration values:
Real 
Simulated 

Definition at line 149 of file orsa_universe.h.

00149                     {
00150     Real=1,
00151     Simulated=2
00152   };


Function Documentation

double 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   = &par;
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,
const 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
 

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.

Referenced by Newton::Acceleration().

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 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 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 AngularMomentum const vector< Body > &  f,
const Vector &  center
 

Definition at line 188 of file orsa_frame.cc.

References ExternalProduct().

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 apply_window fftw_complex *  signal_win,
fftw_complex *  signal,
int  size
 

Definition at line 449 of file orsa_fft.cc.

References cos(), and twopi.

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 Barycenter const vector< Body > &  f  ) 
 

Definition at line 148 of file orsa_frame.cc.

Referenced by BarycentricAngularMomentum().

00148                                             {
00149     return (orsa::CenterOfMass(f));
00150   }

Vector BarycenterVelocity const vector< Body > &  f  ) 
 

Definition at line 152 of file orsa_frame.cc.

00152                                                     {
00153     return (orsa::CenterOfMassVelocity(f));
00154   }

Vector BarycentricAngularMomentum const vector< Body > &  f  ) 
 

Definition at line 203 of file orsa_frame.cc.

References Barycenter().

00203                                                             {
00204     return (orsa::AngularMomentum(f,Barycenter(f)));
00205   }

Here is the call graph for this function:

void 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 CenterOfMass const vector< Body > &  f  ) 
 

Definition at line 122 of file orsa_frame.cc.

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 CenterOfMassVelocity const vector< Body > &  f  ) 
 

Definition at line 135 of file orsa_frame.cc.

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 CloseApproaches_gsl_f const gsl_vector *  v,
void *  params
 

Definition at line 2530 of file orsa_orbit_gsl.cc.

References params.

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 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 Compute const vector< Observation > &  obs  ) 
 

General interface to the Orbit computation from a set of observations.

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.

Referenced by OrbitWithEpoch::Compute().

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 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::x, Vector::y, and Vector::z.

Referenced by Compute().

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 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(), Angle::GetRad(), Universe::GetReferenceSystem(), Orbit::i, JPL_EPHEM_FILE, location_file, obleq_J2000(), LocationFile::ObsPos(), Config::paths, pi, Vector::rotate(), sin(), SUN, universe, Vector::x, Vector::y, and Vector::z.

Referenced by Compute().

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 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.

Referenced by GravitationalTree::Acceleration().

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 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 convert ReferenceSystem rs,
const unsigned int  i
[inline]
 

Definition at line 590 of file orsa_units.h.

References ECLIPTIC, EQUATORIAL, and ORSA_ERROR.

00590                                                                   {
00591     switch(i) {
00592     case 1: rs = EQUATORIAL; break;
00593     case 2: rs = ECLIPTIC; break;
00594       //
00595     default:
00596       ORSA_ERROR("conversion problem: i = %i",i);
00597       break;
00598     }
00599   }

void convert TimeScale ts,
const unsigned int  i
[inline]
 

Definition at line 262 of file orsa_units.h.

References GPS, ORSA_ERROR, TAI, TDT, UT, and UTC.

00262                                                             {
00263     switch(i) {
00264     case 1: ts = UTC; break;
00265     case 2: ts = UT;  break;
00266     case 3: ts = TAI; break;
00267     case 4: ts = TDT; break;
00268     case 5: ts = GPS; break;
00269       //
00270     default:
00271       ORSA_ERROR("conversion problem: i = %i",i);
00272       break;       
00273     }
00274   }

void convert mass_unit mu,
const unsigned int  i
[inline]
 

Definition at line 107 of file orsa_units.h.

References GRAM, KG, MEARTH, MJUPITER, MMOON, MSUN, and ORSA_ERROR.

00107                                                             {
00108     switch(i) {
00109     case 1: mu = MSUN;     break;
00110     case 2: mu = MJUPITER; break;
00111     case 3: mu = MEARTH;   break;
00112     case 4: mu = MMOON;    break;
00113     case 5: mu = KG;       break;
00114     case 6: mu = GRAM;     break;    
00115       //
00116     default:
00117       ORSA_ERROR("conversion problem: i = %i",i);
00118       break;       
00119     }
00120   }

void convert length_unit lu,
const unsigned int  i
[inline]
 

Definition at line 78 of file orsa_units.h.

References AU, CM, EARTHMOON, KM, KPARSEC, LY, M, MPARSEC, ORSA_ERROR, PARSEC, REARTH, and RMOON.

00078                                                               {
00079     switch(i) {
00080     case 1:  lu = MPARSEC;   break;
00081     case 2:  lu = KPARSEC;   break;
00082     case 3:  lu = PARSEC;    break;
00083     case 4:  lu = LY;        break;
00084     case 5:  lu = AU;        break;
00085     case 6:  lu = EARTHMOON; break;
00086     case 7:  lu = REARTH;    break;
00087     case 8:  lu = RMOON;     break;
00088     case 9:  lu = KM;        break;
00089     case 10: lu = M;         break;
00090     case 11: lu = CM;        break;
00091       //
00092     default:
00093       ORSA_ERROR("conversion problem: i = %i",i);
00094       break;       
00095     }
00096   }

void convert time_unit tu,
const unsigned int  i
[inline]
 

Definition at line 46 of file orsa_units.h.

References DAY, HOUR, MINUTE, ORSA_ERROR, SECOND, and YEAR.

00046                                                             {
00047     switch(i) {
00048     case 1: tu = YEAR;   break;
00049     case 2: tu = DAY;    break;
00050     case 3: tu = HOUR;   break;
00051     case 4: tu = MINUTE; break;
00052     case 5: tu = SECOND; break;
00053       //
00054     default:
00055       ORSA_ERROR("conversion problem: i = %i",i);
00056       break;       
00057     }
00058   }

void 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 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 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 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 cos const Angle &  alpha  )  [inline]
 

Definition at line 566 of file orsa_units.h.

Referenced by alpha_delta_meridian(), alpha_delta_meridian_moon(), apply_window(), Orbit::Compute(), Compute_Gauss(), Compute_TestMethod(), dQ(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), Orbit::GetE(), LocationFile::ObsPos(), phi(), phi_Hanning(), Orbit::RelativePosVel(), Vector::rotate(), and sincos().

00566                                          {
00567     return std::cos(alpha.GetRad());
00568   }

Vector Cross const Vector &  u,
const Vector &  v
[inline]
 

Definition at line 182 of file orsa_coord.h.

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 delta_function const unsigned int  i,
const unsigned int  j
 

Definition at line 91 of file orsa_interaction_tree.cc.

00091                                                                     {
00092     if (i==j) return 1.0;
00093     return 0.0;
00094   }

double 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 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 obleq().

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 obleq_J2000().

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 obleq().

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 obleq_J2000().

Referenced by LocationFile::ObsPos().

01977                                              {
01978     v.rotate(0,-obleq_J2000().GetRad(),0);  
01979   }

Here is the call graph for this function:

Vector ExternalProduct const Vector &  u,
const Vector &  v
[inline]
 

Definition at line 176 of file orsa_coord.h.

Referenced by Newton::Acceleration(), 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  x,
const mass_unit  u,
const int  power
 

Definition at line 137 of file orsa_universe.cc.

00137 { return orsa::units->FromUnits(x,u,power); };

double FromUnits const double  x,
const length_unit  u,
const int  power
 

Definition at line 136 of file orsa_universe.cc.

00136 { return orsa::units->FromUnits(x,u,power); };

double FromUnits const double  x,
const time_unit  u,
const int  power
 

Definition at line 135 of file orsa_universe.cc.

00135 { return orsa::units->FromUnits(x,u,power); };

double FromUnits const   double,
const   time_unit,
const   int = 1
 

Referenced by BuildObservationTriplets(), Compute(), Compute_Gauss(), Compute_TestMethod(), Evolution::Evolution(), GalacticPotentialAllen::GalacticPotentialAllen(), TimeStep::GetDouble(), JPLFile::GetEph(), JPLFile::GetMass(), Date::GetTime(), Evolution::Integrate(), TimeStep::operator *=(), poly_8_gsl_solve(), 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(), TimeStep::TimeStep(), and Weight().

void 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(), and ORSA_ERROR.

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 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(), params, 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 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(), params, and PropagatedSky_J2000().

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 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, params, 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 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(), gauss_v_f(), and params.

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.

Referenced by Compute_Gauss(), modified_mu(), Newton::Newton(), OptimizedOrbitPositions::PropagatedSky_J2000(), PropagatedSky_J2000(), and Relativistic::Relativistic().

00127 { return (orsa::units->GetC()); }

double GetG  ) 
 

Definition at line 124 of file orsa_universe.cc.

Referenced by BodyConstants::BodyConstants(), 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(), AstorbFile::Read(), and Relativistic::Relativistic().

00124 { return (orsa::units->GetG()); }

double GetG_MKS  ) 
 

Definition at line 125 of file orsa_universe.cc.

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.

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.

Referenced by LocationFile::ObsPos().

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 Interpolate const vector< VectorWithParameter >  vx_in,
const double  x,
Vector &  v_out,
Vector &  err_v_out
 

Definition at line 101 of file orsa_coord.cc.

00101                                                                                                               {
00102     
00103     // cerr << "Interpolate() call..." << endl;
00104     
00105     unsigned int i,j,m;
00106     
00107     unsigned int i_closest;
00108     
00109     double      diff,denom;
00110     double      ho,hp;
00111     double      tmp_double;
00112 
00113     unsigned int n_points = vx_in.size();
00114     
00115     /* cerr << " -----> n_points: " << n_points << endl;
00116        for (j=0;j<n_points;j++)
00117        std::cerr << "vx_in[" << j << "].par = " <<  vx_in[j].par << "\n";
00118     */
00119     
00120     if (n_points < 2) {
00121       cerr << "too few points..." << endl;
00122       for (j=0;j<n_points;j++) {
00123         v_out = vx_in[0];
00124         err_v_out.Set(0,0,0);
00125       }
00126       return;
00127     }
00128     
00129     Vector * c = new Vector[n_points];
00130     Vector * d = new Vector[n_points];
00131     
00132     Vector w;
00133     
00134     diff = fabs(x-vx_in[0].par);
00135     i_closest = 0;
00136     
00137     for (j=0;j<n_points;j++) { 
00138       
00139       tmp_double = fabs(x-vx_in[j].par);
00140       if (tmp_double <  diff) {
00141         diff = tmp_double;
00142         i_closest = j;
00143       }
00144       
00145       c[j] = vx_in[j];
00146       d[j] = vx_in[j]; 
00147     }
00148     
00149     v_out     = vx_in[i_closest];
00150     err_v_out = vx_in[i_closest];
00151     
00152     // cerr << "local..  v: " << Length(v_out) << "  i_closest: " << i_closest << endl;
00153     
00154     i_closest--;
00155     
00156     for (m=1;m<=(n_points-2);m++) {
00157       for (i=0;i<(n_points-m);i++) {
00158         ho = vx_in[i].par   - x;
00159         hp = vx_in[i+m].par - x;
00160 
00161         denom = ho - hp;
00162 
00163         if (denom == 0.0) {
00164           cerr << "interpolate() --> Error: divide by zero\n";
00165           cerr << "i: " << i << "   m: " << m << endl;
00166           for (j=0;j<n_points;j++)
00167             cerr << "vx_in[" << j << "].par = " <<  vx_in[j].par << "\n";
00168           exit(0);
00169         }
00170         
00171         w = c[i+1] - d[i];
00172 
00173         d[i] = hp*w/denom;
00174         c[i] = ho*w/denom;
00175 
00176         // cerr << " ++++++++++  d[" << i << "]: " << Length(d[i]) << endl;
00177         // cerr << " ++++++++++  c[" << i << "]: " << Length(c[i]) << endl;
00178       }
00179       
00180       if ( (2*i_closest) < (n_points-m) ) {
00181         err_v_out = c[i_closest+1]; 
00182       } else {
00183         err_v_out = d[i_closest];
00184         i_closest--;
00185       }
00186       
00187       v_out += err_v_out;    
00188     }
00189     
00190     delete [] c;
00191     delete [] d;
00192   }

void 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.

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   }

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.

Referenced by JPLBody::JPLBody().

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 KineticEnergy const vector< Body > &  f  ) 
 

Definition at line 207 of file orsa_frame.cc.

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   }

double KineticEnergy const Body &   ) 
 

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 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  ) 
 

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 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, Orbit::omega_pericenter, and params.

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 least_sq_diff_f double  x,
void *  params
 

Definition at line 870 of file orsa_orbit_gsl.cc.

References Sky::delta_arcsec(), params, 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 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, params, 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 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(), least_sq_f(), and params.

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.

Referenced by OrsaFile::Read().

00140 { return orsa::units->LengthLabel(); };

static double local_C22 const JPL_planets  planet  )  [static]
 

Definition at line 278 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_C31 const JPL_planets  planet  )  [static]
 

Definition at line 286 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_C32 const JPL_planets  planet  )  [static]
 

Definition at line 294 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_C33 const JPL_planets  planet  )  [static]
 

Definition at line 302 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_C41 const JPL_planets  planet  )  [static]
 

Definition at line 310 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_C42 const JPL_planets  planet  )  [static]
 

Definition at line 318 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_C43 const JPL_planets  planet  )  [static]
 

Definition at line 326 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_C44 const JPL_planets  planet  )  [static]
 

Definition at line 334 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 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.

Referenced by JPLBody::JPLBody().

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 local_J3 const JPL_planets  planet  )  [static]
 

Definition at line 258 of file orsa_body.cc.

References EARTH, JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_J4 const JPL_planets  planet  )  [static]
 

Definition at line 268 of file orsa_body.cc.

References EARTH, JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_mass const JPL_planets  planet  )  [static]
 

Definition at line 243 of file orsa_body.cc.

References JPLFile::GetMass(), and jpl_file.

Referenced by JPLBody::JPLBody().

00243                                                      {
00244     return jpl_file->GetMass(planet);
00245   }

Here is the call graph for this function:

static double local_S31 const JPL_planets  planet  )  [static]
 

Definition at line 342 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_S32 const JPL_planets  planet  )  [static]
 

Definition at line 350 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_S33 const JPL_planets  planet  )  [static]
 

Definition at line 358 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_S41 const JPL_planets  planet  )  [static]
 

Definition at line 366 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_S42 const JPL_planets  planet  )  [static]
 

Definition at line 374 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_S43 const JPL_planets  planet  )  [static]
 

Definition at line 382 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 local_S44 const JPL_planets  planet  )  [static]
 

Definition at line 390 of file orsa_body.cc.

References JPLFile::GetTag(), jpl_file, and MOON.

Referenced by JPLBody::JPLBody().

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 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 MOID2RB_f(), and 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 **  ,
const   IntegratorType
 

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 **  ,
const   InteractionType
 

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.

Referenced by OrsaFile::Read().

00141 { return orsa::units->MassLabel();   };

static double modified_mu const vector< Body > &  f,
const unsigned int  i
[static]
 

Definition at line 89 of file orsa_frame.cc.

References GetC(), and Vector::Length().

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 = &par;
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 = &par;
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 MOID2RB_f const gsl_vector *  v,
void *  params
 

Definition at line 404 of file orsa_orbit_gsl.cc.

References M1(), and params.

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   }

Here is the call graph for this function:

double MOID_f const gsl_vector *  v,
void *  params
 

Definition at line 180 of file orsa_orbit_gsl.cc.

References M1(), and params.

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 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().

00038                                       {
00039     return sqrt( z.re * z.re + z.im * z.im );
00040   }

double 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().

00042                                          {
00043     return ( z.re * z.re + z.im * z.im );
00044   }

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 Newton::Acceleration(), 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 &  ,
const   int
 

UniverseTypeAwareTimeStep operator * const   int,
const UniverseTypeAwareTimeStep & 
 

UniverseTypeAwareTimeStep operator * const UniverseTypeAwareTimeStep &  ts,
const int  i
 

Definition at line 1572 of file orsa_units.cc.

01572                                                                                            {
01573     UniverseTypeAwareTimeStep _ts(ts);
01574     _ts *= i;
01575     return _ts; 
01576   }

UniverseTypeAwareTimeStep operator * const int  i,
const UniverseTypeAwareTimeStep &  ts
 

Definition at line 1566 of file orsa_units.cc.

01566                                                                                            {
01567     UniverseTypeAwareTimeStep _ts(ts);
01568     _ts *= i;
01569     return _ts; 
01570   }

UniverseTypeAwareTimeStep operator * const UniverseTypeAwareTimeStep &  ts,
const double  x
 

Definition at line 1560 of file orsa_units.cc.

01560                                                                                               {
01561     UniverseTypeAwareTimeStep _ts(ts);
01562     _ts *= x;
01563     return _ts; 
01564   }

UniverseTypeAwareTimeStep operator * const double  x,
const UniverseTypeAwareTimeStep &  ts
 

Definition at line 1554 of file orsa_units.cc.

01554                                                                                               {
01555     UniverseTypeAwareTimeStep _ts(ts);
01556     _ts *= x;
01557     return _ts;
01558   }

double operator * const Vector &  u,
const Vector &  v
[inline]
 

Definition at line 189 of file orsa_coord.h.

00189                                                               {
00190     return (u.x*v.x+
00191             u.y*v.y+
00192             u.z*v.z);
00193   }  

Vector operator * const Vector &  v,
const double  f
[inline]
 

Definition at line 156 of file orsa_coord.h.

00156                                                               {
00157     return Vector(v.x*f, v.y*f, v.z*f);
00158   }

Vector operator * const double  f,
const Vector &  v
[inline]
 

Definition at line 152 of file orsa_coord.h.

00152                                                               {
00153     return Vector(v.x*f, v.y*f, v.z*f);
00154   }

bool operator!= const ET_minus_UT_element &  x,
const ET_minus_UT_element &  y
[inline]
 

Definition at line 288 of file orsa_units.cc.

00288                                                                                        {
00289     return (!(x==y));
00290   }

bool operator!= const TAI_minus_UTC_element &  x,
const TAI_minus_UTC_element &  y
[inline]
 

Definition at line 233 of file orsa_units.cc.

00233                                                                                            {
00234     return (!(x==y));
00235   }

bool operator!= const Vector &  v1,
const Vector &  v2
[inline]
 

Definition at line 202 of file orsa_coord.h.

00202                                                                {
00203     return !(v1 == v2);
00204   }

bool operator!= const Body &  b1,
const Body &  b2
[inline]
 

Definition at line 218 of file orsa_body.h.

00218 { return !(b1 == b2); }

Vector operator+ const Vector &  u,
const Vector &  v
[inline]
 

Definition at line 164 of file orsa_coord.h.

00164                                                               {
00165     return Vector(u.x+v.x,
00166                   u.y+v.y,
00167                   u.z+v.z);
00168   }

Vector operator- const Vector &  u,
const Vector &  v
[inline]
 

Definition at line 170 of file orsa_coord.h.

00170                                                               {
00171     return Vector(u.x-v.x,
00172                   u.y-v.y,
00173                   u.z-v.z);
00174   }

Vector operator/ const Vector &  v,
const double  f
[inline]
 

Definition at line 160 of file orsa_coord.h.

00160                                                               {
00161     return Vector(v.x/f, v.y/f, v.z/f); 
00162   }

bool 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 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 operator== const Vector &  v1,
const Vector &  v2
[inline]
 

Definition at line 195 of file orsa_coord.h.

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 operator== const Body &  b1,
const Body &  b2
[inline]
 

Definition at line 209 of file orsa_body.h.

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   } 

void OrbitDifferentialCorrectionsLeastSquares OrbitWithCovarianceMatrixGSL &  orbit,
const vector< Observation > &  obs
 

Definition at line 1023 of file orsa_orbit_gsl.cc.

References Sky::dec(), Sky::delta_arcsec(), ERR, FIT, Angle::GetRad(), least_sq_df(), least_sq_f(), least_sq_fdf(), Osculating, pi, OptimizedOrbitPositions::PropagatedSky_J2000(), Sky::ra(), and secure_pow().

Referenced by Compute().

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 = &par;
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 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 phi_amp double  omega,
fftw_complex  in[],
const int  size
 

amplitude for spectrum, without windowing

Definition at line 106 of file orsa_fft.cc.

References norm_sq(), and phi().

Referenced by phi_gsl().

00106                                                                   {
00107     return sqrt( norm_sq(phi( omega,&in[0],size)) + 
00108                  norm_sq(phi(-omega,&in[0],size)) );
00109   }

Here is the call graph for this function:

double phi_gsl double  x,
void *  params
 

Definition at line 119 of file orsa_fft.cc.

References params, 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 phi_gsl_two double  x,
void *  params
 

Definition at line 128 of file orsa_fft.cc.

References norm(), params, 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 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 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 phi_Hanning_gsl double  x,
void *  params
 

Definition at line 136 of file orsa_fft.cc.

References params, 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 params, and 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 params, and 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 params, and 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 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, Orbit::e, 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 = &params;
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().

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.

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   }

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, Radau15::e, EARTH, EARTH_AND_MOON, FromUnits(), GetC(), UniverseTypeAwareTime::GetTime(), Evolution::Integrate(), JUPITER, Vector::Length(), location_file, MARS, MERCURY, NEPTUNE, LocationFile::ObsPos(), Evolution::push_back(), 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 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 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 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.

Referenced by Body::Body(), BodyConstants::BodyConstants(), BodyWithEpoch::BodyWithEpoch(), and JPLBody::JPLBody().

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 RelativisticBarycenter const vector< Body > &  f  ) 
 

Definition at line 156 of file orsa_frame.cc.

References modified_mu().

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 RelativisticBarycenterVelocity const vector< Body > &  f  ) 
 

Definition at line 170 of file orsa_frame.cc.

References modified_mu().

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 remove_leading_trailing_spaces std::string &  s  )  [inline]
 

Definition at line 505 of file orsa_file.h.

Referenced by JPLFile::GetTag(), JPLFile::JPLFile(), 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 Sky::delta_arcsec(), 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 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().

Referenced by Compute_Gauss(), PreliminaryOrbit::ComputeRMS(), and E1().

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 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 *  ,
const unsigned  int,
const unsigned  int,
std::vector< CloseApproach > &  ,
const   double,
const   double = FromUnits(1, SECOND)
 

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
 

Definition at line 2746 of file orsa_orbit_gsl.cc.

References CloseApproaches_gsl_f(), CloseApproach::distance, CloseApproach::epoch, ExternalProduct(), Vector::Length(), CloseApproach::name, Vector::Normalize(), Vector::Normalized(), CloseApproach::relative_velocity, UniverseTypeAwareTime::SetTime(), 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 = &par;
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(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), 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 GalacticPotentialAllen::Acceleration(), 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(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), 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 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, and SUN.

Referenced by JPLPlanetsNewton::Acceleration(), 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 sin const Angle &  alpha  )  [inline]
 

Definition at line 562 of file orsa_units.h.

Referenced by alpha_delta_meridian(), alpha_delta_meridian_moon(), Orbit::Compute(), Compute_Gauss(), Compute_TestMethod(), dQ(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), Orbit::GetE(), LocationFile::ObsPos(), phi(), phi_Hanning(), Units::Recompute(), Orbit::RelativePosVel(), Vector::rotate(), and sincos().

00562                                          {
00563     return std::sin(alpha.GetRad());
00564   }

void sincos const Angle &  alpha,
double &  s,
double &  c
[inline]
 

Definition at line 574 of file orsa_units.h.

References cos(), and sin().

Referenced by Newton::Acceleration(), Vector::rotate(), and unit_vector().

00574                                                                   {
00575 #ifdef HAVE_SINCOS
00576     ::sincos(alpha.GetRad(),&s,&c); 
00577 #else // HAVE_SINCOS
00578     s = std::sin(alpha.GetRad());
00579     c = std::cos(alpha.GetRad());
00580 #endif // ! HAVE_SINCOS
00581   }

Here is the call graph for this function:

Frame StartFrame const vector< BodyWithEpoch > &  f,
vector< JPL_planets > &  jf,
const Interaction *  itr,
const Integrator *  itg,
const UniverseTypeAwareTime &  t
 

A good frame to start an integration with.

Definition at line 584 of file orsa_universe.cc.

References Evolution::clear(), FromUnits(), UniverseTypeAwareTime::GetDate(), Universe::GetUniverseType(), Evolution::Integrate(), 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(), 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 swap void *  ptr,
unsigned int  size
[static]
 

Definition at line 1510 of file orsa_file.cc.

References ORSA_WARNING, swap_4(), and swap_8().

Referenced by OrsaFile::GoodFile(), and OrsaFile::Read().

01510                                                  {
01511     switch(size) {
01512     case 4: swap_4(ptr); break;
01513     case 8: swap_8(ptr); break;
01514     default:
01515       ORSA_WARNING("called read_swap with size = %i",size);
01516       break;
01517     }
01518   }

Here is the call graph for this function:

static void swap_4 void *  ptr  )  [static]
 

Definition at line 1494 of file orsa_file.cc.

References SWAP_MACRO.

Referenced by swap().

01494                                 {
01495     char *tptr = (char *)ptr, tchar;
01496     
01497     SWAP_MACRO( tptr[0], tptr[3], tchar);
01498     SWAP_MACRO( tptr[1], tptr[2], tchar);
01499   }

static void swap_8 void *  ptr  )  [static]
 

Definition at line 1501 of file orsa_file.cc.

References SWAP_MACRO.

Referenced by swap().

01501                                 {
01502     char *tptr = (char *)ptr, tchar;
01503     
01504     SWAP_MACRO( tptr[0], tptr[7], tchar);
01505     SWAP_MACRO( tptr[1], tptr[6], tchar);
01506     SWAP_MACRO( tptr[2], tptr[5], tchar);
01507     SWAP_MACRO( tptr[3], tptr[4], tchar);
01508   }

int SWIFTRawReadBinaryFile FILE_TYPE  file,
const int  version = 2
 

Definition at line 1166 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().

01166                                                                     {
01167     
01168     // version == 1 --> not filetered data
01169     // version == 2 --> filtered data
01170     
01171     int good = 0;
01172     double dummy;
01173     
01174     if ( version == 1 ) {  
01175       READ_FILE(&dummy,        sizeof(int),    1,file);
01176       READ_FILE(&nast,         sizeof(int),    1,file);
01177       READ_FILE(&el,           sizeof(double), 6,file);
01178       READ_FILE(&file_time,    sizeof(double), 1,file);
01179       good = READ_FILE(&dummy, sizeof(int),    1,file);
01180     } else if ( version == 2 ) {
01181       READ_FILE(&dummy,        sizeof(int),    1,file);
01182       READ_FILE(&nast,         sizeof(int),    1,file);
01183       READ_FILE(&el,           sizeof(double), 6,file);
01184       READ_FILE(&l_ts,         sizeof(double), 1,file);
01185       READ_FILE(&w_ts,         sizeof(double), 1,file);
01186       READ_FILE(&file_time,    sizeof(double), 1,file);
01187       good = READ_FILE(&dummy, sizeof(int),    1,file);
01188     }
01189     
01190     // Units conversions
01191     file_time = FromUnits(file_time,YEAR);
01192     
01193     return (good);
01194   }

Here is the call graph for this function:

double tan const Angle &  alpha  )  [inline]
 

Definition at line 570 of file orsa_units.h.

Referenced by Orbit::Compute(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), and OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation().

00570                                          {
00571     return std::tan(alpha.GetRad());
00572   }

string TimeLabel  ) 
 

Definition at line 139 of file orsa_universe.cc.

Referenced by OrsaFile::Read().

00139 { return orsa::units->TimeLabel();   };

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 unit_vector const Angle &  ra,
const Angle &  dec
[static]
 

Definition at line 1333 of file orsa_orbit.cc.

References 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 Weight vector< Observation > &  obs,
const double  optimal_dt = FromUnits(1,DAY)
 

Definition at line 89 of file orsa_orbit_gsl.cc.

References Weight().

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   }

Here is the call graph for this function:

double 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 DAY, and FromUnits().

Referenced by Weight().

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:


Variable Documentation

const double AU_MKS = 1.49597870660e11
 

Definition at line 50 of file orsa_units.cc.

const double c_a = secure_pow(1.0/d_a, 2)
 

Definition at line 586 of file orsa_orbit_gsl.cc.

Referenced by M1().

const double c_e = secure_pow(1.0/d_e, 2)
 

Definition at line 587 of file orsa_orbit_gsl.cc.

Referenced by M1().

const double c_i = secure_pow(1.0/d_i, 2)
 

Definition at line 588 of file orsa_orbit_gsl.cc.

Referenced by M1(), and Vector::rotate().

const double c_M = secure_pow(1.0/d_M, 2)
 

Definition at line 591 of file orsa_orbit_gsl.cc.

Referenced by M1().

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)
 

Definition at line 589 of file orsa_orbit_gsl.cc.

Referenced by M1().

const double c_omega_pericenter = secure_pow(1.0/d_omega_pericenter, 2)
 

Definition at line 590 of file orsa_orbit_gsl.cc.

Referenced by M1().

Config * config = 0
 

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
 

Definition at line 579 of file orsa_orbit_gsl.cc.

Referenced by S1().

const double d_e = 0.01
 

Definition at line 580 of file orsa_orbit_gsl.cc.

Referenced by S1().

const double d_i = 0.1*(pi/180)
 

Definition at line 581 of file orsa_orbit_gsl.cc.

Referenced by S1().

const double d_M = 5.0*(pi/180)
 

Definition at line 584 of file orsa_orbit_gsl.cc.

Referenced by S1().

const double d_omega_node = 1.0*(pi/180)
 

Definition at line 582 of file orsa_orbit_gsl.cc.

Referenced by S1().

const double d_omega_pericenter = 1.0*(pi/180)
 

Definition at line 583 of file orsa_orbit_gsl.cc.

Referenced by S1().

TimeScale default_Date_timescale = TT
 

Definition at line 218 of file orsa_units.cc.

double el[6]
 

Definition at line 1164 of file orsa_file.cc.

Referenced by SWIFTFile::Read(), and SWIFTRawReadBinaryFile().

const ET_minus_UT_element ET_minus_UT_table[]
 

Definition at line 298 of file orsa_units.cc.

Referenced by Date::delta_seconds().

const ET_minus_UT_element ET_minus_UT_table_final_element = {0,0,0,0}
 

Definition at line 297 of file orsa_units.cc.

Referenced by Date::delta_seconds().

double file_time
 

Definition at line 1164 of file orsa_file.cc.

Referenced by SWIFTFile::Read(), and SWIFTRawReadBinaryFile().

const double G_MKS = 6.67259e-11
 

Definition at line 45 of file orsa_units.cc.

Referenced by Units::GetG_MKS().

const double halfpi = PI/2
 

Definition at line 47 of file orsa_common.h.

Referenced by Newton::Acceleration().

JPLCache * jpl_cache = 0
 

Definition at line 51 of file orsa_universe.cc.

Referenced by Compute_Gauss(), and SetupSolarSystem().

JPLFile * jpl_file = 0
 

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 1164 of file orsa_file.cc.

Referenced by SWIFTFile::Read(), and SWIFTRawReadBinaryFile().

LocationFile * location_file = 0
 

Definition at line 54 of file orsa_universe.cc.

Referenced by Compute_Gauss(), Compute_TestMethod(), OptimizedOrbitPositions::PropagatedSky_J2000(), 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 1163 of file orsa_file.cc.

Referenced by SWIFTFile::AsteroidsInFile(), SWIFTFile::Read(), and SWIFTRawReadBinaryFile().

OrsaPaths * orsa_paths = 0
 

Definition at line 3232 of file orsa_file.cc.

Referenced by OrsaPaths::OrsaPaths().

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 CloseApproaches_gsl_f(), Compute_Gauss(), gauss_v_df(), gauss_v_diff_f(), gauss_v_f(), gauss_v_fdf(), least_sq_df(), least_sq_diff_f(), least_sq_f(), least_sq_fdf(), MOID2RB_f(), MOID_f(), phi_gsl(), phi_gsl_two(), phi_Hanning_gsl(), poly_8(), poly_8_df(), poly_8_fdf(), 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}
 

Definition at line 247 of file orsa_units.cc.

Referenced by Date::delta_seconds().

const double twopi = PI+PI
 

Definition at line 49 of file orsa_common.h.

Referenced by apply_window(), Orbit::Compute(), Sky::Compute_J2000(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), Orbit::GetE(), MOID(), MOID2RB(), phi(), phi_Hanning(), OptimizedOrbitPositions::PropagatedOrbit(), JPLDastcomCometFile::Read(), AstDySMatrixFile::Read(), MPCCometFile::Read(), OrbitWithEpoch::RelativePosVel(), Orbit::RelativePosVel(), and S1().

Units * units = 0
 

Definition at line 39 of file orsa_universe.cc.

Referenced by OrsaFile::Write().

Universe * universe = 0
 

The active universe.

Definition at line 42 of file orsa_universe.cc.

Referenced by UniverseTypeAwareTimeStep::absolute(), Newton::Acceleration(), Compute(), Compute_Gauss(), Sky::Compute_J2000(), Compute_TestMethod(), Frame::ForceJPLEphemerisData(), UniverseTypeAwareTimeStep::GetDouble(), JPLFile::GetEph(), UniverseTypeAwareTime::GetTime(), Evolution::Integrate(), UniverseTypeAwareTimeStep::IsZero(), JPLPlanetsNewton::JPLPlanetsNewton(), LocationFile::ObsPos(), 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
 

Definition at line 1164 of file orsa_file.cc.

Referenced by SWIFTRawReadBinaryFile().


Generated on Tue Jan 11 15:27:50 2005 for liborsa by  doxygen 1.4.0