|
Creating
Location: Silver Spring, MD, USA
|
Not Ranked
:
+0 / -0
0 score
Let's compare simulators
Quote:
Originally Posted by TheBigDog
I have written a gravity simulator. It works in 3D (x,y,z coordinates) and runs on a 1 second time interval.
|
 Excelent! I look forward to seeing it.
I finished my latest, GRAVSIM4, last weekend. It's main features/changes: - Like my previous GRAVSIM*s, it’s a command line interpreter.
- I redid its syntax completely (the old ones’ were kinda clunky).
- Everything in the syntax is an object that can be viewed or set (eg: Earth.M=5.9736e24 #sets the mass of the Earth to its usual value)
- The $ object is the simulator itself
- It’s major enhancement of the old ones is that the bodies (Earth, Moon, the ship, etc.) have an property (.X) that allow you to give each a program consisting of MUMPS code, allowing them to do stuff like change their mass and acceleration per a rocket motor & fuel program
- It’s not as practical, yet, as GRAVSIM2. I don’t plan to make it much more, rather I plan to write .X code for stuff like graphics output. The command line itself, I plan to leave with only primitive output features.
Here’s a run of the simulator showing a 100 kg rocket with 900 kg of fuel and specific impulse 450 s – basically a really big model rocket – go up then down in a little less than 104 s
Code:
XGRAVSIM4>Earth.M=5.9736e24 Moon.M=7.347673e22 .L=363104000,0,0 .V=0,1082,0 #The usual Earth and Moon
XGRAVSIM4>Ship.M=1000 .L=6378137,0,0 .V=0,463.8,0 #1000 kg Ship on the surface of the equator
XGRAVSIM4>Ship.X="n (B,I) s Mf=B(I,""X"",""Mf""),SI=B(I,""X"",""SI""),M=-Mf*B(""TI"")+B(I,0) s:M<B(I,""X"",""M0"") Mf=0 s A=Mf*SI/(M+B(I,0))*2,B(I,21)=B(I,""X"",21)*A,B(I,22)=B(I,""X"",22)*A,B(I,23)=B(I,""X"",23)*A,B(I,0)=-Mf*B(""TI"")+B(I,0)" #Rocket motor
XGRAVSIM4>.X("SI")=250 .X("Mf")=50 .X("M0")=100 #Rocket motor specific impulse, burn rate, and 100 kg minimum vehicle mass
XGRAVSIM4>.X(21)=1 .X(22)=0 .X(23)=0 #thrust vector
XGRAVSIM4>$.T=1 Ship.M .L .V .A #
#SHIP.M=950 .L=6378138.510032063895,463.7998218340913345,0 .V=3.02006412778904451,463.7996436681826685,0 .A=12.82051282051282051,0,0
XGRAVSIM4>$.T=10 Ship.M .L .V .A
#SHIP.M=500 .L=6378414.894590005906,4637.880633501694489,0 .V=75.20675035288839365,463.7643692911025346,0 .A=23.80952380952380952,0,0
XGRAVSIM4>$.T=20 Ship.M .L .V .A
#SHIP.M=100 .L=6380677.051336701228,9275.04882631241004,0 .V=377.1969438729633817,463.6575375527910744,0 .A=0,0,0
#Note: engine shutoff
XGRAVSIM4>$.T=30 Ship.M .L .V .A
#SHIP.M=100 .L=6383959.572000057278,13910.79377670996689,0 .V=279.3238639845522164,463.4797680659022456,0 .A=0,0,0
XGRAVSIM4>$.T=40 Ship.M .L .V .A
#SHIP.M=100 .L=6386263.816969892319,18544.40719033907758,0 .V=181.5368644980712998,463.2312502756681295,0 .A=0,0,0
XGRAVSIM4>$.T=50 Ship.M .L .V .A
#SHIP.M=100 .L=6387590.497549352828,23175.1820474100567,0 .V=83.80606490867480755,462.9120625488694604,0 .A=0,0,0
XGRAVSIM4>$.T=60 Ship.M .L .V .A
#SHIP.M=100 .L=6387940.02674845683,27802.41167964278307,0 .V=-13.89832117965128978,462.5221968732724911,0 .A=0,0,0
#Note: rocket coming down
XGRAVSIM4>$.T=100 Ship.M .L .V .A
#SHIP.M=100 .L=6379564.474404304418,46261.7085752850397,0 .V=-405.0469612267389591,460.2527731771944981,0 .A=0,0,0
XGRAVSIM4>$.T=101 Ship.M .L .V .A
#SHIP.M=100 .L=6379154.529490345703,46721.92565277778461,0 .V=-414.8428666906910533,460.1813818082953234,0 .A=0,0,0
XGRAVSIM4>$.T=102 Ship.M .L .V .A
#SHIP.M=100 .L=6378734.788041651078,47182.07097860561975,0 .V=-424.6400306985589834,460.1092698473749626,0 .A=0,0,0
XGRAVSIM4>$.T=103 Ship.M .L .V .A
#SHIP.M=100 .L=6378305.248784582742,47642.14383192611936,0 .V=-434.4384834381129305,460.0364367936242636,0 .A=0,0,0
XGRAVSIM4>$.T=104 Ship.M .L .V .A
#SHIP.M=100 .L=6377865.910415304252,48102.14349139155726,0 .V=-444.2382551188663537,459.9628821372515378,0 .A=0,0,0
#note: splat! at 104 seconds
XGRAVSIM4>$.ti=0 #how to exit the simulator
Here’s the current MUMPS code of GRAVSIM4:
Code:
N (XCU,XQMP,XGRAVSIM3D1,XGRAVSIM4) x XGRAVSIM3D1(-1) s LO="$",WF=0 f q:'B("TI") w:WF ! s WF=0,WLO="" r "XGRAVSIM4>",R,! x XCU,XGRAVSIM4((R="?":0,1:1)) ;XGRAVSIM4: YA Gravity simulator
s N="" f s N=(XGRAVSIM4(0,N)) q:N="" w (XGRAVSIM4(0,N),";"),! ;XGRAVSIM4(0)
Enter - Object.property to display a property ;XGRAVSIM4(0,1)
- Object.property=value to set a property to a value ;XGRAVSIM4(0,2)
- #followed by any text to be ignored (a comment) ;XGRAVSIM4(0,2.5)
- Any of the above, separated by one or more spaces ;XGRAVSIM4(0,2.7)
Note - If object is omitted, the previously referenced object is assumed ;XGRAVSIM4(0,3)
- Occurances of the text "XGRAVSIM4>" at the beginning are ignored. ;XGRAVSIM4(0,3.5)
- The special object $ has the following properties: ;XGRAVSIM4(0,4)
-- T - elapsed time of simulation. Set to a greater value to advance ;XGRAVSIM4(0,5)
simulation ;XGRAVSIM4(0,6)
-- TI - simulation interval. Default 1. The interval of simulation ;XGRAVSIM4(0,7)
time between calculations. Set to 0 to end simulation. ;XGRAVSIM4(0,8)
-- ALL - Display only. Displays all properties of all objects. ;XGRAVSIM4(0,8.5)
- Other object have the following properties: ;XGRAVSIM4(0,9)
-- M - mass of object ;XGRAVSIM4(0,10)
-- L (=n,n,n) location of object ;XGRAVSIM4(0,11)
-- V (=n,n,n) velocity of object ;XGRAVSIM4(0,12)
-- A (=n,n,n) acceleration of object ;XGRAVSIM4(0,13)
-- X (may be followed by (subscripts...), ="xecute code") code to be ;XGRAVSIM4(0,14)
xecuted each simulation interval. Meaningful symbols: ;XGRAVSIM4(0,15)
--- B(I,0) - mass of object -OK to set ;XGRAVSIM4(0,16)
--- B(I,21) ,22) ,23) - acceleration of object -OK to set ;XGRAVSIM4(0,17)
--- B(OI(object),1) ,2) ,3) - location of any object -read only, please ;XGRAVSIM4(0,18)
--- B(OI(object),11) ,12) ,13) - velocity of any object -read only ;XGRAVSIM4(0,19)
x XGRAVSIM4(1,-1) s (R(1),CC)=R,R(2)=" ",R(3)=0 x XQMP f CP=1:1:R s R(1)=CC,R(2)=" ",R(3)=CP x XQMP i R]"" s C=R q:C?1"#".e s R(1)=R,R(2)="=",R(3)=1 x XQMP s O=(R,"."),O=(O="":LO,1:O),P=(R,".",2),V=C,(V,1,(R)+1)="" x XGRAVSIM4(1,V]"",O'="$") s LO=O ;XGRAVSIM4(1)
f q:R'?1.an1">".e s (R,">")="",(R)="" ;XGRAVSIM4(1,-1)
s N="" f s N=(XGRAVSIM4(1,0,0,1,N)) q:N="" x XGRAVSIM4(1,0,0,1,N) q:N="" ;XGRAVSIM4(1,0,0): show system property
i P="T" s V=B("T") x XGRAVSIM4(2) s N="" ;XGRAVSIM4(1,0,0,1,1)
i P="TI" s V=B("TI") x XGRAVSIM4(2) s N="" ;XGRAVSIM4(1,0,0,1,2)
i P="ALL" x XGRAVSIM4(1,0,0,1,3,0) s O="" f s O=(OI(O)) q:O="" x XGRAVSIM4(1,0,0,1,3,1),XGRAVSIM4(1,0,0,1,3,1.1) ;XGRAVSIM4(1,0,0,1,3)
s O="$" f P="T","TI" x XGRAVSIM4(1,0,0) ;XGRAVSIM4(1,0,0,1,3,0)
f P="M","L","V","A" x XGRAVSIM4(1,0,1) ;XGRAVSIM4(1,0,0,1,3,1)
s (Q,Q0)=(B(OI(O),"X")),P="X" x:(@Q)#10 XGRAVSIM4(1,0,1) f s (P,Q)=(@Q) q:(@Q,2)'=Q0 s (P,1,(Q0))="X(" x XGRAVSIM4(1,0,1) ;XGRAVSIM4(1,0,0,1,3,1.1)
s V="undefined",WLO="" x XGRAVSIM4(2) ;XGRAVSIM4(1,0,0,1,999)
s N="" f s N=(XGRAVSIM4(1,0,1,1,N)) q:N="" x XGRAVSIM4(1,0,1,1,N) i N="" x XGRAVSIM4(2) q ;XGRAVSIM4(1,0,1): show object property
s I=(OI(O)) i 'I s P="",V="undefined",WLO="" s N="" ;XGRAVSIM4(1,0,1,1,-1)
i P="M" s V=B(I,0) s N="" ;XGRAVSIM4(1,0,1,1,1)
i P="L" s V=B(I,1)_","_B(I,2)_","_B(I,3) s N="" ;XGRAVSIM4(1,0,1,1,2)
i P="V" s V=B(I,11)_","_B(I,12)_","_B(I,13) s N="" ;XGRAVSIM4(1,0,1,1,3)
i P="A" s V=B(I,21)_","_B(I,22)_","_B(I,23) s N="" ;XGRAVSIM4(1,0,1,1,4)
i P?1"X".1"(".e k N s V=P,(V,1,2)="B(I,""X"""_(P="X":")",1:","),N((@V,"undefined"))="",N=(N),V=(N,3,(N)-1) k N s N="" ;XGRAVSIM4(1,0,1,1,5)
n N s V=Q,(V,1,(Q0))="X",N((@Q,"undefined"))="",N=(N),V=V_"="_(N,3,(N)-1) ;XGRAVSIM4(1,0,1,1,5,1)
s V="undefined",WLO="" x XGRAVSIM4(2) ;XGRAVSIM4(1,0,1,1,999)
s N="" f s N=(XGRAVSIM4(1,1,0,1,N)) q:N="" x XGRAVSIM4(1,1,0,1,N) q:N="" ;XGRAVSIM4(1,1,0): set system property
i P="T" s N="" f q:B("T")'<V x XGRAVSIM4(1,1,0,1,1,1),XGRAVSIM3D1(0) ;XGRAVSIM4(1,1,0,1,1)
f I=1:1:B("I9") x (B(I,"X")) ;XGRAVSIM4(1,1,0,1,1,1)
i P="TI" s B("TI")=V s N="" ;XGRAVSIM4(1,1,0,1,2)
s V="undefined",WLO="" x XGRAVSIM4(2) ;XGRAVSIM4(1,1,0,1,999)
s N="" f s N=(XGRAVSIM4(1,1,1,1,N)) q:N="" x XGRAVSIM4(1,1,1,1,N) q:N="" ;XGRAVSIM4(1,1,1): set object property
s I=(OI(O)) i P="M" s:'I I=B("I9")+1,OI(O)=I s B(I,0)=V x XGRAVSIM3D1(-1) s N="" ;XGRAVSIM4(1,1,1,1,1)
i I,P="L" s B(I,1)=+V,B(I,2)=+(V,",",2),B(I,3)=+(V,",",3) x XGRAVSIM3D1(-1) s N="" ;XGRAVSIM4(1,1,1,1,2)
i I,P="V" s B(I,11)=+V,B(I,12)=+(V,",",2),B(I,13)=+(V,",",3) x XGRAVSIM3D1(-1) s N="" ;XGRAVSIM4(1,1,1,1,3)
i I,P="A" s B(I,21)=+V,B(I,22)=+(V,",",2),B(I,23)=+(V,",",3) x XGRAVSIM3D1(-1) s N="" ;XGRAVSIM4(1,1,1,1,4)
i I,P?1"X".1"(".e s (P,1,2)="B(I,""X"""_(P="X":")",1:","),@(P_"="_V) s N="" ;XGRAVSIM4(1,1,1,1,5)
s:'I P="" s V="undefined",WLO="" x XGRAVSIM4(2) ;XGRAVSIM4(1,1,1,1,999)
w (WF:" ",1:"#"),(O=WLO:"",1:O),(P]"":".",1:""),P,"=",V s WF=1,WLO=O ;XGRAVSIM4(2): output O.P=V
i P="T" s N="" f q:B("T")'<V x XGRAVSIM4(1,1,0,1,1,1),XGRAVSIM3D1(0) ;XGRAVSIM4(1,1,0,1,1)
f I=1:1:B("I9") x (B(I,"X")) ;XGRAVSIM4(1,1,0,1,1,1)
i P="TI" s B("TI")=V s N="" ;XGRAVSIM4(1,1,0,1,2)
s N="" f s N=(XGRAVSIM4(1,1,1,1,N)) q:N="" x XGRAVSIM4(1,1,1,1,N) q:N="" ;XGRAVSIM4(1,1,1): set object property
s I=(OI(O)) i P="M" s:'I I=B("I9")+1,OI(O)=I s B(I,0)=V x XGRAVSIM3D1(-1) s N="" ;XGRAVSIM4(1,1,1,1,1)
i I,P="L" s B(I,1)=+V,B(I,2)=+(V,",",2),B(I,3)=+(V,",",3) x XGRAVSIM3D1(-1) s N="" ;XGRAVSIM4(1,1,1,1,2)
i I,P="V" s B(I,11)=+V,B(I,12)=+(V,",",2),B(I,13)=+(V,",",3) x XGRAVSIM3D1(-1) s N="" ;XGRAVSIM4(1,1,1,1,3)
i I,P="A" s B(I,21)=+V,B(I,22)=+(V,",",2),B(I,23)=+(V,",",3) x XGRAVSIM3D1(-1) s N="" ;XGRAVSIM4(1,1,1,1,4)
i I,P?1"X".1"(".e s (P,1,2)="B(I,""X"""_(P="X":")",1:","),@(P_"="_V) s N="" ;XGRAVSIM4(1,1,1,1,5)
s:'I P="" s V="undefined",WLO="" x XGRAVSIM4(2) ;XGRAVSIM4(1,1,1,1,999)
w (WF:" ",1:"#"),(O=WLO:"",1:O),(P]"":".",1:""),P,"=",V s WF=1,WLO=O ;XGRAVSIM4(2): output O.P=V
n (R) s R=(R,"abcdefghijklmnopqrstuvwxyz","ABCDEFGHIJKLMNOPQRSTUVWXYZ") ;XCU: fold to uppercase
n (XQMP,R) s Q="""",D=R(2),(M,R)=D_R(1) x XQMP(1) s R=((R(3)):(R,((M,D,1,R(3)))+2,((M,D,1,R(3)+1))),1:(M,D)-1) ;XQMP: quote sensitive R=(R(1),R(2),R(3)) or R=(R(1),R(2)) if 'R(3)
s D1=((D)+1) f I=2:2:(R(1),Q) s (M,Q,I)=(("",((M,Q,I)))," ",D1) ;XQMP(1)
;XGRAVSIM3D1: 3-d gravity simulator core
N (XGRAVSIM3D1,B) s B("G")=(B("G"),6.6742e-11),B("TI")=(B("TI"),1),B("T")=(B("T"),0),B("S")=(B("S"),1),B("I9")=0 x XGRAVSIM3D1(-1,1) m B=B2 ;XGRAVSIM3D1(-1): Assure B valid for call to XGRAVSIM3D1(0)
f J=1:1 s I=(B("")) q:'I s B("I9")=J m B2(J)=B(I) k B(I) f K=0:1:3,11:1:13,21:1:23 s B2(J,K)=(B2(J,K),0) ;XGRAVSIM3D1(-1,1)
N (XGRAVSIM3D1,B) s I9=B("I9"),G=B("G"),TI=B("TI"),S=B("S") X XGRAVSIM3D1(2),XGRAVSIM3D1(3) M B=B1 s B("T")=TI*S+B("T") ;XGRAVSIM3D1(0): Increment B,T by 1*TI
F I2=1:1:I1-1,I1+1:1:I9 S DY=B(I2,1)-Y,DX=B(I2,2)-X,DZ=B(I2,3)-Z,D2=DY*DY+(DX*DX)+(DZ*DZ),D=D2**.5,A=G*B(I2,0)/D2*TI*S,VY=DY/D*A+VY,VX=DX/D*A+VX,VZ=DZ/D*A+VZ ;XGRAVSIM3D1(1): change body I1 velocity
F I1=1:1:I9 S VY=B(I1,11),VX=B(I1,12),VZ=B(I1,13),Y=VY*TI/2+B(I1,1),X=VX*TI/2+B(I1,2),Z=VZ*TI/2+B(I1,3) X XGRAVSIM3D1(1) S VY=B(I1,21)*TI*S+VY,VX=B(I1,22)*TI*S+VX,VZ=B(I1,23)*TI*S+VZ,B1(I1,11)=VY,B1(I1,12)=VX,B1(I1,13)=VZ ;XGRAVSIM3D1(2): change all velocities
F I=1:1:I9 S B(I,1)=B(I,11)+B1(I,11)/2*TI*S+B(I,1),B(I,2)=B(I,12)+B1(I,12)/2*TI*S+B(I,2),B(I,3)=B(I,13)+B1(I,13)/2*TI*S+B(I,3) ;XGRAVSIM3D1(3): change all positions
I’m curious to see how TBD’s VB.net-based simulator and my MUMPS-based simulator compare.
TBD, can you run yours for the following “standard”, a single geostationary orbit? Here are the parameters (all units kg, m, & s, calculated from this wikipedia article): - Gravitational constant = 6.6742e-11
- Orbital radius (from center of Earth) = 42164140.10012396182
- Orbital speed = 3074.661175977908351
- Sample position and velocity at 43082 and 86164 s (half and full orbit)
Here’re the results from GRAVSIM4
Code:
GRAVSIM4>Earth.M=5.9736e24 #usual Earth
GRAVSIM4>Ship.M=10 .L=42164140.10012396182,0,0 .V=0,3074.661175977908351,0 #Ship in geostationary orbit
GRAVSIM4>$.t=43082 ship.l .v
#SHIP.L=-42145157.25275492353,-59527.65293408068424,0 .V=4.343756684837521444,-3076.039916833619272,0
#note: 42145199.29245684292 3076.042983798489842
$.t=86164 ship.l .v
#SHIP.L=42163972.13018178319,119001.7055088261738,0 .V=-8.679701336015950322,3074.648927402258311,0
#note: 42164140.06238785073 3074.661178730290591
Note that the orbit seems stable – its speed changed by only .00000275238224 m/s – but isn’t perfectly circular, dipping from about 42164140 to 42145199, and goes more than 1 full orbit in 86164 s.
----------------
Moderator: Computers and Technology; Medical Science; Science Projects and Homework; Philosophy of Science; Physics and Mathematics; Environmental Studies 
Last edited by CraigD; 11-05-2007 at 07:06 AM..
Reason: Updated code to most current version
|