Go Back   Science Forums
View Single Post
Old 10-30-2007   #6 (permalink)
CraigD's Avatar
CraigD
Creating


Location:
Silver Spring, MD, USA
 
CraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond reputeCraigD has a reputation beyond repute
 



Not Ranked  0 score     
Thumbs up Let's compare simulators

Quote:
Originally Posted by TheBigDog View Post
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
Reply With Quote
 
» Advertisement
» Current Poll
Who's the sexiest man alive? Johnny Depp or Robert Pattinson?
Johnny Depp - 27.27%
3 Votes
Robert Pattinson - 0%
0 Votes
Someone else (please specify) - 45.45%
5 Votes
I'm too macho to think a guy is sexy - 27.27%
3 Votes
Total Votes: 11
You may not vote on this poll.


All times are GMT -8. The time now is 07:45 PM.

Hypography?

Hypography [n.]: A combination of "hyperlink" and "bibliography" - ie, a list of links to electronic documents. Comparable to discography and bibliography, but not cartography.

We have been online since May 2000, and aim to be the best place to find and share science-related content of all kinds.

Share the love!

Please add more science to your life. Use our RSS feeds on your blog, your portal, or your favorite feedreader!


Powered by vBulletin® Version 3.8.3
Copyright ©2000 - 2009, Jelsoft Enterprises Ltd.
Copyright © 2000-2009 Hypography
Part of the Hypography - Science for Everyone Network