Accurate 3D Planet Simulation (no media required) by Ditch8th Aug 2006 14:31
|
|---|
|
Summary A fairly accurate 3D simulation of our solar system based on Keplerian elements. Will show planet positions as of todays date. Free to use and modify. Code is heavily commented. Description Please modify, upgrade, and improve. Please send me a copy of your updates, I'd love to see it..... Code ` This code was downloaded from The Game Creators
` It is reproduced here with full permission
` http://www.thegamecreators.com
`************************************************************************************************
`***************************[ Project: 3D Planets ]*****************************
`***************************[ 8/07/2006 ]*****************************
`***************************[ ]*****************************
`***************************[ W. Jean Floyd ]*****************************
`***************************[ ]*****************************
`***************************[ DarkBasic Pro V1.062 ]*****************************
`***************************[ ]*****************************
`***************************[ External Code Used: ]*****************************
`***************************[ Darthster's "Space Flight Code" ]*****************************
`***************************[ ]*****************************
`***************************[ Data Sources: ]*****************************
`***************************[ Keplerian Elements JPL/Caltech ]*****************************
`***************************[ http://ssd.jpl.nasa.gov ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
`* *
`* A Short List of Needed Improvements: *
`* *
`* 1. Add a correction factor for the rectangular display so circular orbits *
`* actually appear circular. For now, the display window is squared to 1024 x 1024. *
`* 2. Add all the bells and whistles. Actual planetary graphics, sounds, spinning planets, *
`* moons, and rings. This first release was purposely written without external media. *
`* 3. Better orbit tracks. These rectangles are a bit weak. *
`* 4. Add pull down menus and help screens. *
`* 5. Better mouse functionality. *
`* 6. More and improved data display options. *
`* *
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ INITIALIZATION STATEMENTS ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
sync on
sync rate 0
autocam off
Gosub checkdisplay
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ DIMENSIONS, VARIABLES, ]*****************************
`***************************[ AND ARRAYS ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
TYPE planetdata `An assigned data type for the orbital
`elements and other calculated values
planetname$ as string `0 Planet Name
a# as float `1 Semi-Major Axis
aCy# as float `2 Semi-Major Axis Change per Julian Day
e# as float `3 Eccentricity
eCy# as float `4 Eccentricity Change per Julian Day
i# as float `5 Inclination
iCy# as float `6 Inclination Change per Julian Day
l# as float `7 Mean Longitude
lCy# as float `8 Mean Longitude Change per Julian Day
w# as float `9 Long Perihelion
wCy# as float `10 Long Perihelion Change per Julian Day
o# as float `11 Long Ascending Node
oCy# as float `12 Long Ascending Node Change per Julian Day
dia# as float `13 Planet Diameter
M# as float `14 Mean Anomaly
bigE# as float `15 Eccentric Anomaly
V# as float `16 True Anomaly
R# as float `17 Heliocentric Radius
x# as float `18 X-axis Heliocentric Coordinate
y# as float `19 Y-Axis Heliocentric Coordinate
z# as float `20 Z-axis Heliocentric Coordinate
orbperiod# as float `21 Orbital Period (Earth Years)
endtype
global dim labels$(21) as string `An array to hold the 21 labels [0-20] for the
`orbital elements and calculated values
global dim planets(21,9) as planetdata `An array to hold the orbital elements and
`calculated values. Holds 21 data
`points [0-20] for the 9 planets [1-9]
`with field [0] left blank.
global dim KeyHold(256) `Dimensions the arrays used in the KeyHold()
global dim LastKeyHold(256) `and LastKeyHold() functions.
global D# `Julian Day calculated by the CalcJulian()
`function
global Cy# `Julian Century calculated by the CalcJulian()
`function [D#/36525.0]
global sizer=8 `Scale factor for planet placement on the
`computer screen
global month
global day
global year
global MidWidth = 512 `Horizontal center of a 1024 x 768 screen
global MidHeight = 320 `Vertical center of a 1024 x 768 screen
global theta# = 0 `Camera rotation values for pitch and yaw,
global phi# = 0 `needed for the Navigation SubRoutine.
global thruststep# = 0.005 `Camera thrust value.
global restartflag = 1 `Initial flag setting for spacebar keypress
global speed# = 1.0 `Sets the initial Julian Day increment to 1
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ CALL FUNCTIONS AND ]*****************************
`***************************[ SUBROUTINES ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
Gosub loadarrays `Load the data arrays with read statements
CalcJulian() `Call the CalcJulian() function to calculate
`the Julian day (D#) based on today's date as
`retrieved by the DBPro 'Get Date$()' command.
`Refer to the function for additional comments.
starter#=D# `The Julian day (D#) will be incremented in
`the main loop, so this step simply assigns
`today's Julian day to 'starter' to remember
`our starting point
Gosub SetupCameras `Routine to set up and postion the camera
Gosub MakeLights `Routine to make a couple of lights
Gosub MakeObjects `Routine to make orbital track objects
Gosub MakeBackDrop `Routine to make the backdrop
`************************************************************************************************
`************************************************************************************************
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ MAIN LOOP ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
`************************************************************************************************
`************************************************************************************************
do
Gosub CalculateOrbits `Calulates the orbital elements based on time
Gosub AnimatePlanets `Moves the planets based on the new elements
Gosub KeyPresses `Looks for keypresses and processes action
Gosub Navigation `Free flyer code using mouse and AWSD keys
Gosub PrintRoutine `Routine to print data to the screen
inc D#,speed# `Increments the Julian Day by your set speed
sync `Sync the screen
loop `Loop through it again
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ CHECK DISPLAY ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
checkdisplay:
if check display mode(1024,768,32) `Ensures screen can handle 1024 x 768
set display mode 1024,768,32
set window on `Window mode
set window size 1024,1024 `Square the screen (ratio in-work)
set text font "arial"
set text size 16
cls
text 1,0,current graphics card$() `Just for information, not really needed
text 1,20,"Setting Display Mode 1024 x 768, Full Screen"
text 1,40,"Setting Text Font to 'Arial' and Font Size to '16'"
sync:sync
wait 1500 `Pause to read the screen
else
if check display mode(640,480,16) `If screen can't handle 1024 x 768
set display mode 640,480,32 `it gives you a graceful exit, sort of.
cls
text 1,1,"Can't set required display mode of 1024 x 768.....press any key to quit."
sync:sync
wait key
else
end
endif
end
endif
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ LOAD ARRAYS ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
loadarrays:
for j=0 to 21 `Loads the 21 data labels for 'planetdata'
read labels$(j) `array using the 'read' command
next j
for k=1 to 9 `Loads Orbital Elements for each of the planets
read planets(0,k).planetname$ `as well as the planet diameter (dia#) and
read planets(1,k).a# `place the data into the first 14 fields in
read planets(2,k).aCy# `the 'planetdata' array and the last field,
read planets(3,k).e# `which is the orbital period. All other
read planets(4,k).eCy# `fields are calculated values and will be
read planets(5,k).i# `loaded during the calculation routine.
read planets(6,k).iCy#
read planets(7,k).l#
read planets(8,k).lCy#
read planets(9,k).w#
read planets(10,k).wCy#
read planets(11,k).o#
read planets(12,k).oCy#
read planets(13,k).dia#
read planets(21,k).orbperiod# `Not really used for anything at this point.
next k
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ SETUP CAMERAS ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
SetupCameras:
set camera range 0,1,5000 `Expands the range on the default camera
camx#=-236 `Initial placement of the camera
camy#=180
camz#=-334
phi#=21 `Initial rotation of the camera
theta#=36
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ Make and Place Lights ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
MakeLights:
set ambient light 25 `Tones down the ambient light
make light 1 `Add a few lights to brighten up the
make light 2 `orbit tracks
make light 3
position light 1,-700,0,0
position light 2,700,0,0
position light 3, 0,-700,0
set light range 1, 2000
set light range 2, 2000
set light range 3, 2000
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ MAKE OBJECTS SUBROUTINE ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
MakeObjects:
for k = 1 to 9
make object sphere k,(planets(13,k).dia#) `Make 9 planets with appropriate diameter
make object cylinder 500+k,100 `Nine cylinders for orbit tracks
rotate object 500+k,90,0,0 `Rotate tracks into the proper plane
set object cull 500+k,0 `Set cull to 0
next k
color object 501,rgb(255,0,0) `Each object is a"rough" orbit track, which
scale object 501,6,3,6 `means they are handcrafted and not calculated
position object 501,0,-.5,0 `Placement, size, and rotation is approximate
`especially for Pluto. They are used to give
color object 502,rgb(0,255,0) `you a visual indication of the orbital track.
scale object 502,11.99,3,11.44 `The planet locatio itself is more accurate.
position object 502,0,0,0
color object 503,rgb(0,0,255)
scale object 503,16.69,3,16.0
position object 503,0,0,0
color object 504,rgb(255,105,180)
scale object 504,23.85,3,24.65
position object 504,-1,.12,0
color object 505,rgb(255,255,0)
scale object 505,82.6,5,82.36
position object 505,-1.4,0.2,0
color object 506,rgb(0,255,255)
scale object 506,151,5,152
position object 506,-.07,-2.35,0
color object 507,rgb(255,0,0)
scale object 507,306.5959,6,308.59
position object 507,7.36,0,0
color object 508,rgb(0,255,0)
scale object 508,480.93,7,480.4
position object 508,0.3,1.54,0
color object 509,rgb(0,0,255)
scale object 509,627,8,619
position object 509,55.5,50.22,-23.159
rotate object 509,84.7,16.7,2.6
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ BackDrop Routine ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
MakeBackDrop: `Grid Routine, makes a grid design and
cls `wraps it around a sphere
ink rgb(37,37,74),rgb(0,0,0) `Change ink color to dark grey
for i = 64 to 704 step 64 `Draw vertical and horizontal lines
line 0,i,1024,i
next i
for i = 0 to 1024 step 64
line i,0,i,768,
next i
get image 1,0,0,1024,768 `Capture the screen to an image file
make object sphere 100,-3000,40,40 `Make a large sphere
texture object 100,1 `Texture it with the image
MakeLine3D(101,-1500,0,0,1500,0,0) `Draw three 3-D lines to represent the
MakeLine3D(102,0,-1500,0,0,1500,0) `x, y and z axis
MakeLine3D(103,0,0,-1500,0,0,1500)
ink rgb(255,255,255),rgb(0,0,0) `Change ink color back to white on black
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ CALCULATE ORBITS ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
CalculateOrbits: `Main Calculation section
Cy#=D#/36525.0 `Converts the Julian day to a Julian Century
for k=1 to 9 `Cycle through each planet and calculate
a#= planets(1,k).a#+ planets(2,k).aCy#*Cy#`orbital elements based on elapsed time
e#= planets(3,k).e#+ planets(4,k).eCy#*Cy#
i#= planets(5,k).i#+ planets(6,k).iCy#*Cy#/3600
l#= planets(7,k).l#+ planets(8,k).lCy#*Cy#/3600
w#= planets(9,k).w#+ planets(10,k).wCy#*Cy#/3600
o#= planets(11,k).o#+planets(12,k).oCy#*Cy#/3600
M#=wrapvalue(l#-w#) `Ensure this value is between 0 and 360
bigE#=M#+e#*sin(M#)*(1.0+e#*cos(M#)) `Estimate bigE#
Repeat `Interate until difference is small (<.1)
bigE2#=bigE#
bigE#=bigE2#-(bigE2#-e#*sin(bigE2#)-M#)/(1-e#*cos(bigE2#))
Until abs(bigE2#-bigE#) < .1
V#=wrapvalue(2*atan(sqrt((1+e#)/(1-e#))*tan(0.5*bigE#)))
R#=(a#*(1-e#^2))/(1+e#*cos(V#))
planets(14,k).M#=M# `Load calculated values into the array
planets(15,k).bigE#=bigE#
planets(16,k).V#=V#
planets(17,k).R#=R#
planets(18,k).x#=R#*(cos(o#)*cos(V#+w#-o#)-sin(o#)*sin(V#+w#-o#)*cos(i#))
planets(19,k).y#=R#*(sin(o#)*cos(V#+w#-o#)+cos(o#)*sin(V#+w#-o#)*cos(i#))
planets(20,k).z#=R#*(sin(V#+w#-o#)*sin(i#))
next k
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ Animate Objects ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
AnimatePlanets:
for k=1 to 9 `Cycle through each planet and place it in the
`calculated position. Uses a multiplier,
`sizer, to distribute it better on the screen
position object k,sizer*planets(18,k).x#,sizer*planets(19,k).y#,sizer*planets(20,k).z#
if object in screen(k) `Print the planet name to the screen
text object screen x(k)+8,object screen y(k)-6,planets(0,k).planetname$
endif
next k
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ Navigation Routine ]*****************************
`***************************[ Darthster's "Space Flight Code" ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
Navigation:
`Mouse Control for rotations x and y
theta#=curveangle(theta#+mousemovex(),theta#,8)
phi#=curveangle(phi#+mousemovey(),phi#,8)
if KeyHold(18) then thrust#=thrust#+.25 `E Key for High Speed Forward Thrusters
if KeyHold(46) then thrust#=thrust#-.25 `C Key for High Speed Reverse Thrusters
if KeyHold(17) then thrust#=thrust#+thruststep# `W
if KeyHold(30) then xaxisstrafe#=xaxisstrafe#-thruststep# `A
if KeyHold(31) then thrust#=thrust#-thruststep# `S
if KeyHold(32) then xaxisstrafe#=xaxisstrafe#+thruststep# `D
thrust#=thrust#*0.199 `Free flight code borrowed from
xaxisstrafe#=xaxisstrafe#*0.199 `Darthster's Space Flight Code (Thanks!)
yaxisstrafe#=yaxisstrafe#*0.199
xvel#=xvel#+(sin(theta#)*cos(phi#)*thrust#) `Velocity Calculations
yvel#=yvel#-(sin(phi#)*thrust#)
zvel#=zvel#+(cos(theta#)*cos(phi#)*thrust#)
xvel#=xvel#+(sin(theta#+90)*xaxisstrafe#)
zvel#=zvel#+(cos(theta#+90)*xaxisstrafe#)
xvel#=xvel#+(sin(theta#)*cos(phi#+90)*yaxisstrafe#)
yvel#=yvel#-(sin(phi#+90)*yaxisstrafe#)
zvel#=zvel#+(cos(theta#)*cos(phi#+90)*yaxisstrafe#)
xvel#=xvel#*0.99 `Decay Velocities and Sum
yvel#=yvel#*0.99
zvel#=zvel#*0.99
oldx#=camx#
oldy#=camy#
oldz#=camz#
camx#=camx#+xvel#
camy#=camy#+yvel#
camz#=camz#+zvel#
position camera 0,camx#,camy#,camz# `Reposition the camera
rotate camera 0,phi#,theta#,0 `Rotate the camera
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ Print Routine ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
PrintRoutine:
box 0,0,1024,25,rgb(0,0,100),rgb(0,0,0),rgb(0,0,100),rgb(0,0,0)
box 0,743,1024,768,rgb(0,0,0),rgb(0,0,100),rgb(0,0,0),rgb(0,0,100)
line 0,25,1024,25
line 0,743,1024,743
text 5,3,"Starting Date"+" "+str$(month)+"/"+str$(day)+"/"+str$(year)
text 630,3,"Days Since Starting Date: "+str$(int(D#-starter#))+" Years Since Starting Date: "+str$((int(D#-starter#)/365))
text 5,30,"FPS:" :text 135,30,str$(screen fps())
text 5,45,"Camera X,Y,Z:" :text 135,45,str$(camera position x(0))+", "+str$(camera position y(0))+", "+str$(camera position z(0))
text 5,60,"Camera Phi#:" :text 135,60,str$(phi#)
text 5,75,"Camera Theta#:" :text 135,75,str$(theta#)
text 5,90,"KeyState:" :text 135,90,str$(scancode())
text 5,105,"Speed(days per sec):" :text 135,105,str$(screen fps()*speed#)
text 5,748," 1,2,3,4,5 = Views * -key = Slower * +key = Faster * Spacebar = Stop and View for Today's Date/Resume * W,S,A,D,Mouse = Navigate"
Return
`************************************************************************************************
`***************************[ ]*****************************
`***************************[ Routine to Process Key Presses ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
KeyPresses:
UpdateKeyHold()
if KeyHit(2) `keypress=1
camx#=0
camy#=600
camz#=0
phi#=90
theta#=0
endif
if KeyHit(3) `keypress=2
camx#=-600
camy#=0
camz#=0
phi#=0
theta#=90
endif
if KeyHit(4) `keypress=3
camx#=0
camy#=-600
camz#=0
phi#=280
theta#=0
endif
if KeyHit(5) `keypress=4
camx#=-236
camy#=180
camz#=-334
phi#=21
theta#=36
endif
if KeyHit(6) `keypress=5
camx#=0
camy#=0
camz#=-91
phi#=0
theta#=0
endif
if KeyHold(74) `keypress=-
speed#=speed#-0.001
endif
if KeyHold(78) `keypress=+
speed#=speed#+0.001
endif
if KeyHit(57) `keypress=spacebar
restartflag=restartflag*-1 `Freezes the screen at today's date
if restartflag=-1
remember#=D#
D#=starter#
speed#=0
else `Press spacebar again and it resumes
D#=remember# `from your last point at speed#=1
speed#=1
endif
endif
Return
`************************************************************************************************
`***************************[ CalcJulian() ]*****************************
`***************************[ ]*****************************
`***************************[ Calculates Julian Based on GetDate ]*****************************
`***************************[ ]*****************************
`***************************[ D#= Days since J2000 ]*****************************
`***************************[ Cy#= Centuries since J2000 ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
function CalcJulian()
today$=get date$()
month=val(left$(today$,2))
day=val(mid$(today$,4))*10 + val(mid$(today$,5))
year=val(right$(today$,2))+2000
D#=int(367*year-(7*(year+((month+9)/12)))/4+(275*month)/9+day-730530)
Cy#=D#/36525.0
endfunction
`************************************************************************************************
`***************************[ UpdateKeyHold() ]*****************************
`***************************[ ]*****************************
`***************************[ Checks for a "held" key press ]*****************************
`************************************************************************************************
function UpdateKeyHold()
for t=0 to 256
LastKeyHold(t) = KeyHold(t)
next t
for t=0 to 256
KeyHold(t) = Keystate(t)
next t
endfunction
`************************************************************************************************
`***************************[ KeyHit(t) ]*****************************
`***************************[ ]*****************************
`***************************[ Checks for a "hit" key press ]*****************************
`************************************************************************************************
function KeyHit(t)
retval=0
if not LastKeyHold(t)
if KeyHold(t)
retval=1
endif
endif
endfunction retval
`************************************************************************************************
`***************************[ Makes a 3D Line Using Triangles ]*****************************
`***************************[ ]*****************************
`***************************[ ]*****************************
`************************************************************************************************
function MakeLine3D(number, x#, y#, z#, tox#, toy#, toz#)
make object triangle number, x#, y#, z#, x#, y#, z#, tox#, toy#, toz#
set object wireframe number,1
endfunction
`************************************************************************************************
`***************************[ Data ]*****************************
`***************************[ ]*****************************
`***************************[ Data Tags and Orbital Elements ]*****************************
`************************************************************************************************
dataset:
Data "Planet","SemiMajor Axis","SemiMajor Axis Cy","Eccentricity","Eccentricity Cy","Inclination","Inclination Cy","Mean Longitude","Mean Longitude Cy","Long Perihelion","Long Perihelion Cy","Long Ascending Node","Long Ascending Node Cy","Diameter","Mean Anomaly","Eccentric Anomaly","True Anomaly", "Heliocentric Radius","X","Y","Z","Orbital Period"
Data "Mercury",0.38709893,0.00000066,0.20563069,0.00002527,7.00487,-23.51,252.25084,538101628.3,77.45645,573.57,48.33167,-446.3,0.3825,.2408
Data "Venus",0.72333199,0.00000092,0.00677323,-0.00004938,3.39471,-2.86,181.97973,210664136.1,131.53298,-108.8,76.68069,-996.89,.9489,.6152
Data "Earth",1.00000011,-0.00000005,0.01671022,-0.00003804,0.00005,-46.94,100.46435,129597740.6,102.94719,1198.28,-11.26064,-18228.25,1.0,1.0
Data "Mars",1.52366231,-0.00007221,0.09341233,0.00011902,1.85061,-25.47,355.45332,68905103.78,336.04084,1560.78,49.57854,-1020.19,.5326,1.88
Data "Jupiter",5.20336301,0.00060737,0.04839266,-0.0001288,1.3053,-4.15,34.40438,10925078.35,14.75385,839.93,100.55615,1217.17,11.021,11.862
Data "Saturn",9.53707032,-0.0030153,0.0541506,-0.00036762,2.48446,6.11,49.94432,4401052.95,92.43194,-1948.89,113.71504,-1591.05,9.45,29.458
Data "Uranus",19.19126393,0.00152025,0.04716771,-0.0001915,0.76986,-2.09,313.23218,1542547.79,170.96424,1312.56,74.22988,-1681.4,4.00,84.01
Data "Neptune",30.06896348,-0.00125196,0.00858587,0.0000251,1.76917,-3.64,304.88003,786449.21,44.97135,-844.43,131.72169,-151.25,3.88,164.79
Data "Pluto",39.48168677,-0.00076912,0.24880766,0.00006465,17.14175,11.07,238.92881,522747.9,224.06676,-132.25,110.30347,-37.33,.1874,248.54
|