/
Orbit.cpp
91 lines (74 loc) · 1.64 KB
/
Orbit.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include "Orbit.h"
#include "DNest4/code/DNest4.h"
#include <fstream>
#include <iostream>
#include <cmath>
using namespace std;
Orbit::Orbit()
{
}
void Orbit::load(const char* filename)
{
fstream fin(filename, ios::in);
if(!fin)
{
cerr<<"# ERROR: Couldn't open file "<<filename<<"."<<endl;
return;
}
// Empty vectors
vx.clear();
vy.clear();
double temp1, temp2;
while(fin>>temp1 && fin>>temp2)
{
vx.push_back(temp1);
vy.push_back(temp2);
}
cout<<"# Loaded "<<vx.size()<<" points from orbit file ";
cout<<filename<<"."<<endl;
fin.close();
}
vector<double> Orbit::evaluate(const std::vector<double>& arg_to_cos,
double viewing_angle) const
{
// Create radial velocities
double C = cos(viewing_angle);
double S = sin(viewing_angle);
vector<double> vr = vx;
double vrmax = -1E300;
double vrmin = 1E300;
int highest = 0;
for(size_t i=0; i<vr.size(); i++)
{
vr[i] = C*vx[i] + S*vy[i];
if(vr[i] > vrmax)
vrmax = vr[i];
if(vr[i] < vrmin)
vrmin = vr[i];
if(vr[i] > vr[highest])
highest = i;
}
for(size_t i=0; i<vr.size(); i++)
vr[i] = 2*(vr[i] - vrmin)/(vrmax - vrmin) - 1.;
vector<double> result = arg_to_cos;
int index;
double cc = 2.*M_PI*((double)(highest))/vr.size();
for(size_t i=0; i<result.size(); i++)
{
index = (int)(vr.size()*
DNest4::mod(arg_to_cos[i] + cc, 2.*M_PI)/(2*M_PI));
result[i] = vr[index];
}
return result;
}
void Orbit::test()
{
Orbit o;
o.load("Orbits/orbits0.710.dat");
vector<double> t;
for(double tt=-10; tt <= 10; tt += 0.01)
t.push_back(tt);
vector<double> y = o.evaluate(t, 1.);
for(size_t i=0; i<y.size(); i++)
cout<<t[i]<<' '<<y[i]<<endl;
}