#include <cmath>
#include <string>
#include <map>
using namespace std;
class Frem{
public:
Frem ( map<string, double> data ) : t(0.0), m(data.at("m")), r(data.at("r")), x(0.0), y(0.0), z(data.at("h")),
vx(data.at("v")*cos(data.at("phi")*(M_PI/180))*sin(data.at("theta")*(M_PI/180))), vy(data.at("v") * sin(data.at("phi")*(M_PI/180))*sin(data.at("theta")*(M_PI/180))), vz(data.at("v") * cos(data.at("theta")*(M_PI/180))),
kappa(data.at("kappa")), g(data.at("g")),
wx(data.at("w") * cos(data.at("alpha")*(M_PI/180))), wy(data.at("w")*sin(data.at("alpha")*(M_PI/180))), wz(0.0){};
double get_t() const {return t;}
double get_x() const {return x;};
double get_y() const {return y;};
double get_z() const {return z;};
void make_step(double dt);
private:
double g;
double t;
double m, r;
double x, y, z;
double vx, vy, vz;
double wx, wy ,wz;
double kappa;
double ax(double vx, double vy, double vz, double wx, double wy, double wz) const;
double ay(double vx, double vy, double vz, double wx, double wy, double wz) const;
double az(double vx, double vy, double vz, double wx, double wy, double wz) const;
void Frem::make_step(double dt)
{
double kx1=vx*dt;
double ky1=vy*dt;
double kz1=vz*dt;
double kvx1=ax(vx, vy, vz, wx, wy, wz)*dt;
double kvy1=ay(vx, vy, vz, wx, wy, wz)*dt;
double kvz1=az(vx, vy, vz, wx, wy, wz)*dt;
double kx2=(vx+kvx1/2)*dt;
double ky2=(vy+kvy1/2)*dt;
double kz2=(vz+kvz1/2)*dt;
double kvx2=ax(vx+kvx1/2, vy+kvy1/2, vz+kvz1/2, wx, wy,wz)*dt;
double kvy2=ay(vx+kvx1/2, vy+kvy1/2, vz+kvz1/2, wx, wy,wz)*dt;
double kvz2=az(vx+kvx1/2, vy+kvy1/2, vz+kvz1/2, wx ,wy,wz)*dt;
double kx3=(vx+kvx2/2)*dt;
double ky3=(vy+kvy2/2)*dt;
double kz3=(vz+kvz2/2)*dt;
double kvx3=ax(vx+kvx2/2, vy+kvy2/2,vz+kvz2/2 ,wx ,wy,wz)*dt;
double kvy3=ay(vx+kvx2/2 ,vy+kvy2/2,vz+kvz2/2 ,wx ,wy,wz)*dt;
double kvz3=az(vx+kvx2/2 ,vy+kvy2/2,vz+kvz2/2 ,wx ,wy,wz)*dt;
double kx4=(vx+kvx3)*dt;
double ky4=(vy+kvy3)*dt;
double kz4=(vz+kvz3)*dt;
double kvx4=ax(vx+kvx3 ,vy+kvy3,vz+kvz3 ,wx ,wy,wz)*dt;
double kvy4=ay(vx+kvx3 ,vy+kvy3,vz+kvz3 ,wx ,wy,wz)*dt;
double kvz4=az(vx+kvx3 ,vy+kvy3,vz+kvz3 ,wx ,wy,wz)*dt;
t+=dt;
x+=(kx1 + 2*kx2 + 2*kx3 + kx4)/6;
y+=(ky1 + 2*ky2 + 2*ky3 + ky4)/6;
z+=(kz1 + 2*kz2 + 2*kz3 + kz4)/6;
vx+=(kvx1 + 2*kvx2 + 2*kvx3 + kvx4)/6;
vy+=(kvy1 + 2*kvy2 + 2*kvy3 + kvy4)/6;
vz+=(kvz1 + 2*kvz2 + 2*kvz3 + kvz4)/6;
}
double Frem::ax(double vx, double vy, double vz, double wx, double wy, double wz) const
{
return -kappa * 1.29 * M_PI * r * r /(2*m) * (vx - wx) * sqrt(pow(vx - wx, 2) + pow(vy - wy, 2) + pow(vz, 2));
}
double Frem::ay(double vx, double vy, double vz, double wx, double wy, double wz) const
{
return - kappa * 1.29 * M_PI * r * r /(2*m) * (vy - wy) * sqrt(pow(vx - wx, 2) + pow(vy - wy, 2) + pow(vz, 2));
}
double Frem::az(double vx, double vy, double vz, double wx, double wy,double wz) const
{
return -g -kappa * 1.29 * M_PI * r * r /(2*m) * (vz - wz) * sqrt(pow(vx - wx, 2) + pow(vy - wy, 2) + pow(vz, 2));
}
};