#include "draw.h"
#include <math.h>
double iTime=0;
//Initialize 2 points givng them some values.
CDraw::CDraw(){
ZeroMemory(point,sizeof(gPOINT)*11);
//Sun
point[0].dMass=1.99e30;
//Earth
point[1].dPosX=-2.65e10;
point[1].dPosY=1.32754E11;
point[1].dPosZ=5.7558e10;
point[1].dVeloX=-2.98E+04;
point[1].dVeloY=-4.78E+03;
point[1].dVeloZ=-2.06E+03;
point[1].dMass=5.98e24;
bFirst=true;
DistanceFactor=5e8;
TimeFactor=50000;
}
//Deinitialize the points, giving them 0's
CDraw::~CDraw(){
ZeroMemory(point,sizeof(gPOINT)*11);
}
//Eulers method.
void CDraw::Euler(){
hDC = BeginPaint(FindWindow(0,schAppName), &Ps);
for(unsigned short i=0;i<2;i++){
point.dForceX=0;
point.dForceY=0;
point.dForceZ=0;
}
for(unsigned short i=0;i<2;i++){
for(unsigned short j=0;j<2;j++){
if(i!=j){
point.dForceX+=ForceX(point,point[j]);
point.dForceY+=ForceY(point,point[j]);
point.dForceZ+=ForceZ(point,point[j]);
}
}
}
for(unsigned short i=0;i<2;i++){
point.dVeloX=Velocity(point.dVeloX,Acceleration(point.dForceX,point.dMass),Time(1));
point.dVeloY=Velocity(point.dVeloY,Acceleration(point.dForceY,point.dMass),Time(1));
point.dVeloZ=Velocity(point.dVeloZ,Acceleration(point.dForceZ,point.dMass),Time(1));
}
for(unsigned short i=0;i<2;i++){
point.dPosX+=point.dVeloX*Time(1);
point.dPosY+=point.dVeloY*Time(1);
point.dPosZ+=point.dVeloY*Time(1);
}
for(unsigned short i=0;i<2;i++){
Point(point,COLOR(0,255,0));
}
iTime+=TimeFactor;
cout<<"Position x:"<<point[1].dPosX<<endl;
cout<<"Position y:"<<point[1].dPosY<<endl;
cout<<"Position z:"<<point[1].dPosZ<<endl;
cout<<"Years:"<<iTime/3.1536e7<<endl;
// system("pause");
EndPaint(FindWindow(0,schAppName), &Ps );
}
//Do the time scaling using Time factor
double CDraw::Time(long double dTimeMs){
return (dTimeMs*TimeFactor);
}
//calculate the new velocity. newVelocity=oldVelocioty +AccelerationTime
double CDraw::Velocity(double dLastVelocity,double dAcceleration,double dTime){
if(bFirst){
bFirst=false;
return dLastVelocity+dAcceleration*0.5*dTime;
}
return dLastVelocity+dAcceleration*dTime;
}
//Do acceleration from f=ma
double CDraw::Acceleration(double dForce,double dMass){
return dForce/dMass;
}
//Calculate force in the X direction
//Fx=forceGravity*changeinX/radius
double CDraw::ForceX(gPOINT point,gPOINT point2){
return (Gravity(point,point2)*(point2.dPosX-point.dPosX)/Pytheg(point,point2));
}
//Calculate force in the X direction
//Fy=forceGravity*changeiny/radius
double CDraw::ForceY(gPOINT point,gPOINT point2){
return (Gravity(point,point2)*(point2.dPosY-point.dPosY)/Pytheg(point,point2));
}
//Calculate force in the X direction
//Fy=forceGravity*changeiny/radius
double CDraw::ForceZ(gPOINT point,gPOINT point2){
return (Gravity(point,point2)*(point2.dPosZ-point.dPosZ)/Pytheg(point,point2));
}
//Gravity, Fg=Gmm/r^2
double CDraw::Gravity(gPOINT point,gPOINT point2){
return ((6.67e-11*point.dMass*point2.dMass)/pow(Pytheg(point,point2),2));
}
//r^2=a^2+b^2
//r=sqrt(a^2+b^2)
double CDraw::Pytheg(gPOINT point,gPOINT point2){
return sqrt(pow(point2.dPosX-point.dPosX,2)+pow(point2.dPosY-point.dPosY,2)+pow(point2.dPosZ-point.dPosZ,2));
}
Gravity in 3d
Author
Hello.
I added a 3d aspect to my gravity program and i am running into some problems....the rotation doesn;t happen. I am pretty sure i transported the stuff correctly and all.
The data i am sure is correct.
The step size is small enough, i can increase it 100 fold in 2d and still get accurate results, therefor this should still work.
The problem is probably somewhere in the code and i am missing some small thing.
Here is the code:
I know that is alot of code but it is very simple and alot of repetitaion.
Any ideas woudl be well appriciated, becuase i think i transformed all of it correcty into 3d.
Thx
Rotational accelerations only result from the application of torques on a body which can be thought of as the rotational equivalent to forces. Torque results from the application of a force which is not centered on a bodies center of mass. Since gravity applies forces across a body which 'even out' they can always be considered to act on the center of mass, hence never creating rotational motions. So in your example where you consider only point masses and gravitational forces it is quite correct that you never get rotation involved.
Author
My mistake i was talking about revolutions. I have sun + earth system in 3d, and there is some error in the formula that i missed. Sorry for the confusion with rotation.
I was talking about revolutions.
I was talking about revolutions.
point.dPosZ+=point.dVeloY*Time(1);
Try using a vector class to avoid this sort of error prone repetition.
Try using a vector class to avoid this sort of error prone repetition.
Author
thx alot, this is the stuff most of my time is wasted on :D
however there is another error that makes alot of trouble by offsetingthe results by a large amount....
any ideas?
Thx alot for your help.
[Edited by - rever on December 31, 2005 1:44:17 AM]
however there is another error that makes alot of trouble by offsetingthe results by a large amount....
any ideas?
#include <math.h>#include <string.h>#include "draw.h"double iTime=0;//Initialize 2 points givng them some values.CDraw::CDraw(){ ZeroMemory(point,sizeof(gPOINT)*11); //Sun point[0].dMass=1.99e30; //Mercury point[1].dPosX=-2.65e10; point[1].dPosY=1.32754E11; point[1].dPosZ=5.7558e10; point[1].dVeloX=-2.98E+04; point[1].dVeloY=-4.78E+03; point[1].dVeloZ=-2.06E+03; point[1].dMass=5.98e24; bFirst=true; DistanceFactor=5e8; TimeFactor=10000;}//Deinitialize the points, giving them 0'sCDraw::~CDraw(){ ZeroMemory(point,sizeof(gPOINT)*11);}//Eulers method. void CDraw::Euler(){ hDC = BeginPaint(FindWindow(0,schAppName), &Ps); for(unsigned short i=0;i<2;i++){ point.dForceX=0; point.dForceY=0; point.dForceZ=0; } for(unsigned short i=0;i<2;i++){ for(unsigned short j=0;j<2;j++){ if(i!=j){ point.dForceX+=ForceX(point,point[j]); point.dForceY+=ForceY(point,point[j]); point.dForceZ+=ForceZ(point,point[j]); } } } for(unsigned short i=0;i<2;i++){ point.dVeloX=Velocity(point.dVeloX,Acceleration(point.dForceX,point.dMass),Time(1)); point.dVeloY=Velocity(point.dVeloY,Acceleration(point.dForceY,point.dMass),Time(1)); point.dVeloZ=Velocity(point.dVeloZ,Acceleration(point.dForceZ,point.dMass),Time(1)); } for(unsigned short i=0;i<2;i++){ point.dPosX+=point.dVeloX*Time(1); point.dPosY+=point.dVeloY*Time(1); point.dPosZ+=point.dVeloZ*Time(1); } iTime+=TimeFactor; cout<<"Position x:"<<point[1].dPosX<<endl; cout<<"Position y:"<<point[1].dPosY<<endl; cout<<"Position z:"<<point[1].dPosZ<<endl; cout<<"Velocity x:"<<point[1].dVeloX<<endl; cout<<"Velocity y:"<<point[1].dVeloY<<endl; cout<<"Velocity z:"<<point[1].dVeloZ<<endl; cout<<"Years:"<<iTime/3.1536e7<<endl; bFirst=false; EndPaint(FindWindow(0,schAppName), &Ps );}//Do the time scaling using Time factordouble CDraw::Time(long double dTimeMs){ double Time=(dTimeMs*TimeFactor); return Time;}//calculate the new velocity. newVelocity=oldVelocioty +AccelerationTimedouble CDraw::Velocity(double dLastVelocity,double dAcceleration,double dTime){ if(bFirst){ double v=(dLastVelocity+dAcceleration*0.5*dTime); return v; } double v2=(dLastVelocity+dAcceleration*dTime); return v2;}//Do acceleration from f=madouble CDraw::Acceleration(double dForce,double dMass){ double a=dForce/dMass; return a;}//Calculate force in the X direction//Fx=forceGravity*changeinX/radiusdouble CDraw::ForceX(gPOINT point,gPOINT point2){ double x=(Gravity(point,point2)*(point2.dPosX-point.dPosX)/Pytheg(point,point2)); return x;} //Calculate force in the Y direction//Fy=forceGravity*changeiny/radiusdouble CDraw::ForceY(gPOINT point,gPOINT point2){ double y=(Gravity(point,point2)*(point2.dPosY-point.dPosY)/Pytheg(point,point2)); return y;}//Calculate force in the Z direction//Fy=forceGravity*changeinz/radiusdouble CDraw::ForceZ(gPOINT point,gPOINT point2){ double z=(Gravity(point,point2)*(point2.dPosZ-point.dPosZ)/Pytheg(point,point2)); return z;}//Gravity, Fg=Gmm/r^2double CDraw::Gravity(gPOINT point,gPOINT point2){ double g=((6.67e-11*point.dMass*point2.dMass)/pow(Pytheg(point,point2),2)); return g;}//r^2=a^2+b^2//r=sqrt(a^2+b^2)double CDraw::Pytheg(gPOINT point,gPOINT point2){ double p=sqrt(pow(point2.dPosX-point.dPosX,2)+pow(point2.dPosY-point.dPosY,2)+pow(point2.dPosZ-point.dPosZ,2)); return p;}Thx alot for your help.
[Edited by - rever on December 31, 2005 1:44:17 AM]
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
