'Basic program to see and understand
    'the movement of the sun and planets
    'in terms of x and y forces

    'John Erich Ebner 11-7-09          roseandjohn2001@yahoo.com

    nomainwin
    WindowWidth =800
    WindowHeight =800'780

    open "solarwobble110709" for graphics as #1
    #1,"font times new roman 18"
    '#1,"font arial 14 italic"
    #1,"home ; down ; posxy x y "
    #1,"fill cyan"  'window color
    #1,"backcolor red"  'fill color
    '#1,"backcolor green"  'fill color
    #1,"color blue "  'pencolor
    #1,"size 1"  'pen size

    #1,"rule notor"  'pen mode
    '#1,"rule notand"
    '#1,"rule xor"
    '#1,"rule "; _R2_NOTXORPEN

    s=0  'used in save
    size=1  'this is the size multiplier
    xo=380:yo=350  'start location of sun
    xoo=xo:yoo=yo  'copy of original xo an yo
     sun=0  'sun = 0 is a dot, sun = 1 is a disk
     gosub [var]  'load orbital parameters

     'various things to plot
     'size=1:t=p(1)/2:gosub [var]:gosub [orrary] ':gosub [save]
     'size=1:t=p(1)/2:gosub [wobble]  ':gosub [save]
     'size=1:t=p(1)/7:gosub [wobble]  ':gosub [save]
     size=1:t=p(1)/14:gosub [wobble]  ':gosub [save]
     'size=1:t=p(3)/12:gosub [wobble]  ':gosub [save]
     'size=15:t=4000001:gosub [wobble]  ':gosub [save]
     'size=20:t=p(3)/12:sun=0:gosub [wobble] ':gosub [save]
      'size=20:t=p(3)/12:sun=1:gosub [wobble]  :gosub [save]
     'size=4:t=p(5)/50:sun=0:gosub [wobble]  ':gosub [save]

     #1, "flush"
     wait
    close #1'
    end

[wobble]
xs=xo:ys=yo
r=38*size

'initial values of velocity vx and vy are necessary
'for the sun to loop in somewhat tight orbits and stay on the page
'rather than drift off, on a loopy diagonal path = a helical path,
'which it is inclined to do with the planets in tow
'around or across the galaxy

vx=10.4 'starting values of vx and vy in m/s
vy=9.4  'these two result in 10.75 m/s or 23.8 mph

if (t>.5E6 and t<1E6) then vx=10.1:vy=9.6
if (t>1E6 and t<2E6) then vx=10.2:vy=9.5
if (t>2E6 and t<4E6) then vx=10.5:vy=9.4
if (t>4E6) then vx=10.4:vy=9.4

scale=90E6/size  'ratio of meters/pixels

'one orbit of Pluto equals p(9)/p(1) = 1027.8008 orbits of Mercury

     'p(1)/2=3.8E6 seconds, mercurys orbit of 88.023 days/2
     'p(1)/7=1.08E6 seconds, mercurys orbit of 88.023 days/7
     'p(1)/14=543E3 seconds, mercurys orbit of 88.023 days/14
     'p(3)/12=2.6E6  monthly
     'p(3)/52=607E3   weekly


't=time increment, planet orbital period/something
't=(p(3)/12) is monthly
for zz= 1 to p(9)/t 'main loop, Pluto orbit seconds/t
   gosub [calc]
   gosub [plot]
   'for cc=1 to 20000:next cc' slow down loop
   'nn=nn+1:if nn=100 then nn=0:for rr=0 to 1000000:next rr:#1,"cls"
next zz

'gosub [cross]  'make a cross at the center of the average locations of the sun

return


[calc]
fx=0:fy=0
for p=1 to 9  'calc force and planetary angle
   fx=fx+f(p)*cos(a(p))  'sum x forces
   fy=fy+f(p)*sin(a(p))  'sum y forces
   a(p)=a(p)+6.283185/(p(p)/t)  'increment planetary angles
   if a(p)>6.263185 then a(p)=a(p)-6.283185
next p

'acc=force/mass   1.307E22/1.989E30=6.5711E-9_m/s^2
'velocity=vo+acc*time   vo+6.5711E-9*3802609=2.4987E-2_m/s
'dist=velocity*time   2.4987E-2*3802609=95017.6_m
'dist=(vo+acc*time)*time

'solar radius 696E6_m
'solar mass 1.989E30_kg

ax=fx/1.989E30 'force/mass = solar acceleration
ay=fy/1.989E30

xd=(vx+ax*t)*t 'x distance uses the old vx and the new ax
yd=(vy+ay*t)*t

vx=(vx+ax*t) 'save the velocity vx for the next loop
vy=(vy+ay*t)

xi=xd/scale  'x increment pixels
yi=yd/scale  'y increment pixels

xo=xo+xi:yo=yo+yi  'move the sun pixels

xs=xs+xo:ys=ys+yo  'sum x & y coordinates for average

return  'comment out "return" to print data

xnn=(xd^2+yd^2)^.5 'actual distance meters
xii=(xi^2+yi^2)^.5    'distance of increment pixels

#1,"place "; 400;" ";20
#1,"\";  xn;"   ";yn;"   ";xnn  'actual distance

#1," place "; 400;" ";50
#1,"\";  xi;"   ";yi;"   ";xii  'drawing distance

#1,"place "; 400;" ";80
#1,"\";  scale;"   ";696E6/scale  'the scale in meters per pixel and sun radius

#1,"place "; 400;" ";110
#1,"\";  vx;"  ";vy;"       "  'solar velocity
return


[plot]
#1,"place ";xo;"  ";yo' sun center
#1,"color yellow"
'nn=nn+1:if nn=12 then nn=0:#1,"backcolor blue":#1,"circlefilled ";3 else #1,"backcolor red":#1,"circlefilled ";3
if sun=0 then #1,"circlefilled "; 3
if sun=1 then #1,"circlefilled "; 696E6/scale 'sun as a disk to scale with the suns movement not with the planets

for p=1 to 9  'p is planet number
  xp=p*r*cos(a(p))  'p is planet number
  yp=p*r*sin(a(p))   'r is arbitrary radius used for scaling
  #1,"place ";xo+xp;"  ";yo+yp  'location of planet dot
 select case p  'get planet color
   case 1
   #1,"color darkred":gosub [dot]
   case 2
   #1,"color yellow":gosub [dot]
   case 3
   #1,"color darkblue":gosub [dot]
   case 4
   #1,"color blue":gosub [dot]
   case 5
   #1,"color lightgray":gosub [dot]
   case 6
   #1,"color darkgreen":gosub [dot]
   case 7
   #1,"color blue":gosub [dot]
   case 8
   #1,"color darkblue":gosub [dot]
   case 9
   #1,"color darkpink":gosub [dot]
 end select

 'nn=f(p)/f(9):#1,"circlefilled ";log(nn)+2;'size proportional to force

 '#1,"circle 1";
 '#1," goto ";xp;"  ";yp
next p
return

[dot]
 #1,"circlefilled ";3;  'make a planet dot
 'for cc=1 to 9000:next cc' slow down loop
 '#1,"circlefilled ";25;  'make a planet dot
return
 #1,"color blue"
 #1,"goto ";xo;"  ";yo 'draw a line to the sun
 xc=398:yc=240
 #1,"place ";xc;"  ";yc
 #1,"circlefilled ";5;'si
 #1,"place ";0;"  ";yc
 #1,"goto ";800;"  ";yc
 #1,"place ";xc;"  ";0
 #1,"goto ";xc;"  ";800
return

[cross] 'draw a cross at the average x,y position
xc=xs/zz  'average x location of the center of the sun
yc=ys/zz
#1,"place ";xc;"  ";yc
#1,"circlefilled ";5;'si
#1,"place ";0;"  ";yc
#1,"goto ";800;"  ";yc
#1,"place ";xc;"  ";0
#1,"goto ";xc;"  ";800
#1,"place ";xc;"  ";yc
#1,"circle ";696E6/scale
return  'comment out return to print data
#1,"place "; 20;" ";20
#1,"\";  xoo;"     ";yoo  'start

#1,"place "; 20;" ";50
pz=((xc-xoo)^2+(yc-yoo)^2)^.5
#1, "\";  xc;"   ";yc;"  ";pz

#1,"place "; 20;" ";80
#1,"\"; pz*scale
#1,"place "; 20;" ";110
#1,"\"; vxx;"  ";vyy
#1,"place "; 20;" ";140
#1,"\"; "Plutos orbit ";r(9)/scale;"  pixels or ";r(9)/scale/72/12;"  feet"
#1,"place "; 20;" ";170
#1,"\";  xc;"     ";yc  'start

'#1,"place ";xc;"  ";yc
'#1,"circle ";174
return

[orrary]
xo=350:yo=350
rad=35
    #1,"size 3"  'pen size
    #1,"backcolor yellow" ' fill color
    #1,"fill white"  'window color
    #1,"place ";xo;" ";yo;
    #1,"circlefilled 15"
for p=1 to 9
    #1,"color cyan"
    #1,"backcolor yellow" ' fill color
    #1,"circle ";p*rad
    #1,"color green"
    #1,"goto ";xo+p*rad*cos(a(p));" ";yo
    #1,"color blue"
    #1,"goto ";xo+p*rad*cos(a(p));" ";yo+p*rad*sin(a(p))
    #1,"color blue":#1," circlefilled 10"
    #1,"color red"
    #1,"goto ";xo;" ";yo;
    xs=xs+p*rad*cos(a(p))
    ys=ys+p*rad*sin(a(p))
next p
#1,"size 5"
#1,"place ";xo;" ";yo;
#1,"goto ";xs;" ";ys;
#1,"size 1"
return


[var]
#1,"cls"
#1,"fill cyan" 'window color

#1,"place ";xo;"  ";yo
'#1,"circlefilled ";10;   'start dot

a(1)= 5.74054'angm planetary angle to start  on 1-1944
a(2)=  4.94436  'angv
a(3)=  2.97685'ange
a(4)= 2.08145'angm
a(5)= 2.51852'angj
a(6)=  1.50738'angs
a(7)= 1.18698 'angu
a(8)= 3.19058 'angn
a(9)= 2.22995'angp

'radians    5.74054    4.94436    2.97685    2.08145    2.51852    1.50738    1.18698    3.19058    2.22995
'meters    59173816240    108858295040    148645462240    244486702240    801841967520    1351206274000    2893791908480    4527751173440    5658824592960

r(1)=59173816240'm distance to sun on 1-1944 1:00UT DST
r(2)=108858295040'v
r(3)=148645462240'e
r(4)=244486702240'm
r(5)=801841967520'j distance to sun
r(6)=1351206274000's
r(7)=2893791908480'u
r(8)=4527751173440'n
r(9)=5658824592960   'p

m(0)=1.989E30'mass sun
m(1)=3.30E23 'mass mercury
m(2)=4.87E24 'venus
m(3)=6.05E24 'earth + moon
m(4)=6.42E23 'mars
m(5)=1.90E27 'jupiter
m(6)=5.69E26 'saturn
m(7)=8.68E25 'uranus
m(8)=1.02E26 'neptune
m(9)=1.29E22 'pluto


f(1)=1.307E22'fm force is calculated, this is a check
f(2)=5.522E22'fv
f(3)=3.547E22'fe
f(4)=1.641E21'fm
f(5)=4.17E23'fj
f(6)=3.71E22'fs
f(7)=1.39E21'fu
f(8)=6.69E20'fn
f(9)=4.898E16'fn

p(1)=07605219'period seconds 88.023 days
p(2)=19407509
p(3)=31556925'earth year in seconds
p(4)=59358576
p(5)=374328244
p(6)=929540783
p(7)=2652990685
p(8)=5200896809
p(9)=7816650322'pluto 247.7 years
'mercury makes 1027.8008 orbits per pluto
' in 247.7 years
'2*mercury=2055.6 points
'7*mercury=7194.6 points
'14*mercury=14389.2 points
G=6.674E-11
'#1," place 20 20"
for p=1 to 9 'calculate the gravitational force for all the planets
f(p)=G*m(0)*m(p)/r(p)^2
'#1, "\";  n;"      ";f(n);
next p
return



[save]
    #1,"getbmp drawing ";0;" ";(0);" ";(800);" ";(800)'x y width height
   s$= str$(s):s$="000"+s$:s$=right$(s$,3)
   filename$ = "orbdoc" +s$+ ".bmp"
   'filedialog "Save as...", ".bmp", filename$
   'if filename$ = "" then wait
   bmpsave "drawing", filename$
   unloadbmp "drawing"
   s=s+1
return

a(1)= 5.18342'angm planetary angle to start  on 1-1948
a(2)=  0.04569 'angv
a(3)=  1.86983'ange
a(4)= 2.27295'angm
a(5)= 4.38746'angj
a(6)=  2.41035'angs
a(7)= 1.47316 'angu
a(8)= 3.33581'angn
a(9)= 2.33274'angp


r(1)=65743395520'm distance to sun on 1-1948 1:00UT DST
r(2)=108658623920'v
r(3)=147107005760'e
r(4)=246855139520'm
r(5)=797282279200'j distance to sun
r(6)=1371592191200's
r(7)=2854047362320'u
r(8)=4531245410560'n
r(9)=5510734176640   'p
'return
'radians    5.18342    0.04569    1.86983    2.27295    4.38746    2.41035    1.47316    3.33581    2.33274
'meters    65743395520    108658623920    147107005760    246855139520    797282279200    1371592191200    2854047362320    4531245410560    5510734176640

r(1)=5.79E10'm distance to sun on 8-19-09 1:00UT DST
r(2)=1.082E11'v
r(3)=1.496E11'e
r(4)=2.279E11'm
r(5)=7.78E11'j distance to sun
r(6)=1.43E12's
r(7)=2.87E12'u
r(8)=4.49E12'n
r(9)=5.913E12'p

a(1)= 4.3689'angm planetary angle to start  on 8-19-09
a(2)=  .991769'angv
a(3)=  .5767008'ange
a(4)= .8523509'angm
a(5)= 5.628'angj
a(6)= 3.036'angs
a(7)= 6.184'angu
a(8)= 5.674'angn
a(9)= 4.752657'angp