class DifferentiationFunction
{
public:virtual void Differentiate(float x, float y, float z, float t,float& outX, float& outY, float& outZ) = NULL;virtual float GetExtendT() const = NULL;virtual float GetStartX() const{return 0.0f;}virtual float GetStartY() const{return 0.0f;}virtual float GetStartZ() const{return 0.0f;}
};

// 行星的导数曲线
class Planet : public DifferentiationFunction
{
public:Planet(){m_star_weight = 1.0f;m_planet_x = 5.0f;m_planet_y = 8.0f;m_planet_z = 1.0f;m_planet_weight = 0.1f;m_planet_speed_x = 4.0f;m_planet_speed_y = 0.0f;m_planet_speed_z = 0.0f;m_g = 100.0f;m_ek = 0.5f*m_planet_weight*(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z); // 1/2*m*v*vfloat r = sqrt(m_planet_x*m_planet_x + m_planet_y*m_planet_y + m_planet_z*m_planet_z);m_ep = -/*0.5f**/m_g*m_star_weight*m_planet_weight/r;m_e = m_ek + m_ep;}void Differentiate(float x, float y, float z, float t,float& outX, float& outY, float& outZ){t = t*10.0f;float sqd = x*x + y*y + z*z;float d = sqrt(sqd);float a = m_g*m_star_weight/sqd;float ax = -a*x/d;float ay = -a*y/d;float az = -a*z/d;outX = x + m_planet_speed_x*t + 0.5f*ax*t*t;outY = y + m_planet_speed_y*t + 0.5f*ay*t*t;outZ = z + m_planet_speed_z*t + 0.5f*az*t*t;m_planet_speed_x += ax*t;m_planet_speed_y += ay*t;m_planet_speed_z += az*t;float r = sqrt(outX*outX + outY*outY + outZ*outZ);m_ep = -/*0.5f**/m_g*m_star_weight*m_planet_weight/r;m_ek = m_e - m_ep;if (m_ek < 0.0f){m_ek = 0.0f;}float v = sqrt(2*m_ek/m_planet_weight);float w = sqrt(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z);m_planet_speed_x *= v/w;m_planet_speed_y *= v/w;m_planet_speed_z *= v/w;}float GetExtendT() const{return 20.0f;}float GetStartX() const{return m_planet_x;}float GetStartY() const{return m_planet_y;}float GetStartZ() const{return m_planet_z;}public:float m_star_weight;float m_planet_x;float m_planet_y;float m_planet_z;float m_planet_weight;float m_planet_speed_x;float m_planet_speed_y;float m_planet_speed_z;float m_g;      // 万有引力系数float m_e;float m_ek;     // 动能float m_ep;     // 引力势能
};

// 三体问题的导数曲线
class ThreeBody : public DifferentiationFunction
{
public:ThreeBody(){m_star1_x = -10.0f;m_star1_y = 0.0f;m_star1_z = 0.0f;m_star1_weight = 1.0f;m_star2_x = 10.0f;m_star2_y = 0.0f;m_star2_z = 0.0f;m_star2_weight = 1.0f;m_planet_x = 5.0f;m_planet_y = 5.0f;m_planet_z = 0.1f;m_planet_weight = 0.1f;m_planet_speed_x = 0.0f;m_planet_speed_y = 2.0f;m_planet_speed_z = 0.0f;m_g = 50.0f;m_ek = 0.5f*m_planet_weight*(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z); // 1/2*m*v*vfloat d1x = m_star1_x - m_planet_x;float d1y = m_star1_y - m_planet_y;float d1z = m_star1_z - m_planet_z;float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z;float d1 = sqrt(sqd1);m_ep1 = -m_g*m_star1_weight*m_planet_weight/d1;float d2x = m_star2_x - m_planet_x;float d2y = m_star2_y - m_planet_y;float d2z = m_star2_z - m_planet_z;float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z;float d2 = sqrt(sqd2);m_ep2 = -m_g*m_star2_weight*m_planet_weight/d2;m_e = m_ek + m_ep1 + m_ep2;}void Differentiate(float x, float y, float z, float t,float& outX, float& outY, float& outZ){t = t*20.0f;float d1x = m_star1_x - x;float d1y = m_star1_y - y;float d1z = m_star1_z - z;float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z;float d1 = sqrt(sqd1);float d2x = m_star2_x - x;float d2y = m_star2_y - y;float d2z = m_star2_z - z;float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z;float d2 = sqrt(sqd2);float a1 = m_g*m_star1_weight/sqd1;float a1x = a1*d1x/d1;float a1y = a1*d1y/d1;float a1z = a1*d1z/d1;float a2 = m_g*m_star2_weight/sqd2;float a2x = a2*d2x/d2;float a2y = a2*d2y/d2;float a2z = a2*d2z/d2;outX = x + m_planet_speed_x*t + 0.5f*(a1x + a2x)*t*t;outY = y + m_planet_speed_y*t + 0.5f*(a1y + a2y)*t*t;outZ = z + m_planet_speed_z*t + 0.5f*(a1z + a2z)*t*t;m_planet_speed_x += (a1x + a2x)*t;m_planet_speed_y += (a1y + a2y)*t;m_planet_speed_z += (a1z + a2z)*t;{float d1x = m_star1_x - outX;float d1y = m_star1_y - outY;float d1z = m_star1_z - outZ;float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z;float d1 = sqrt(sqd1);m_ep1 = -m_g*m_star1_weight*m_planet_weight/d1;float d2x = m_star2_x - outX;float d2y = m_star2_y - outY;float d2z = m_star2_z - outZ;float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z;float d2 = sqrt(sqd2);m_ep2 = -m_g*m_star2_weight*m_planet_weight/d2;m_ek = m_e - m_ep1 - m_ep2;if (m_ek < 0.0f){m_ek = 0.0f;}float v = sqrt(2*m_ek/m_planet_weight);float w = sqrt(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z);m_planet_speed_x *= v/w;m_planet_speed_y *= v/w;m_planet_speed_z *= v/w;}}float GetExtendT() const{return 20.0f;}float GetStartX() const{return m_planet_x;}float GetStartY() const{return m_planet_y;}float GetStartZ() const{return m_planet_z;}public:float m_star1_x;float m_star1_y;float m_star1_z;float m_star1_weight;float m_star2_x;float m_star2_y;float m_star2_z;float m_star2_weight;float m_planet_x;float m_planet_y;float m_planet_z;float m_planet_weight;float m_planet_speed_x;float m_planet_speed_y;float m_planet_speed_z;float m_g;      // 万有引力系数float m_e;float m_ek;     // 动能float m_ep1;    // 引力势能float m_ep2;    // 引力势能
};