REM
REM GE Copeland
REM Dept of Physics
REM Old Dominion University
REM Fall 2003
REM
CLS
PRINT " FORM-IT forming your own solar system."
PRINT "We consider closed orbits only."
PRINT "This program allows you to select a number of planetary objects in"
PRINT "a solar system, by inputing several elements of their orbits."
PRINT "Like:"
PRINT , "a = semimajor axis"
PRINT , "e = eccentricity "
PRINT , "w = the angle of the perihelion "
PRINT , "i = the angle of inclination "
PRINT , "v = the true anonmaly "
PRINT , "theta = position of the object at epoch"
PRINT
PRINT "The output will be a data file in the form of xandxdot.dat"
PRINT " This file format is used by many of the planetary simulation programs."
PRINT " Output will be the x,y,z and vx, vy and vz data for each"
PRINT "object and units will be AU and AU/day."
REM PRINT "Input the mass of the central star in Solar mass units": INPUT Ms
PRINT "Input the number of objects to generate.": INPUT np
DIM x(np), y(np), z(np), vx(np), vy(np), vz(np)
DIM a(np), e(np), w(np), In(np), v(np), Theta(np)

PRINT " Do you want to input data by hand or by a file?"
PRINT "Input HAND or FILE": INPUT Ans$
IF Ans$ = "HAND" OR Ans$ = "hand" OR Ans$ = "Hand" THEN
GOTO 10
ELSE
PRINT "File name is ELEMENTS.DAT"
PRINT
GOTO 11
END IF
10 REM
FOR i = 1 TO np
PRINT " For object #"; i;
PRINT "Input a(au)": INPUT a(i)
4 PRINT "Input eccentricity e": INPUT e(i)
IF e(i) >= 1! OR e(i) < 0 THEN
PRINT "closed orbits  between zero and <1."
GOTO 4
ELSE

END IF

PRINT "Input the angle of the perihelion - measured from Vernal Equinox"
INPUT w(i)
5 PRINT "Input the Inclination in degrees": INPUT In(i)
IF In(i) >= 90! OR In(i) < 0 THEN
PRINT " Inclination is between 0 and 90 degrees."
GOTO 5
ELSE
END IF

PRINT "Input the angle phi - position at epoch in degrees": INPUT Theta(i)
NEXT i

GOTO 9
11 REM input dat from the file
OPEN "ELEMENTS.DAT" FOR INPUT AS #3
FOR i = 1 TO np
INPUT #3, a(i), e(i), w(i), In(i), Theta(i)
NEXT i

CLOSE #3



9 REM

pi = 4! * ATN(1!)
REM convert to radians
conv = pi / 180!
REM - the main loop
REM calculations from WM Smart Spherical Astronomy
REM              chapter V pages around 126
REM
REM and AE Roy's Orbital Motion 2nd edition
REM Section 4.11 and page 102
REM The Orbit in space
REM

  OPEN "xxdot.dat" FOR OUTPUT AS #4

FOR i = 1 TO np
w(i) = conv * w(i)
In(i) = conv * In(i)
Theta(i) = conv * Theta(i)
NEXT i

REM PRINT "#   x   y    z   Vx  Vy   VZ"
FOR i = 1 TO np
b = a(i) * SQR(1 - e(i) ^ 2)
a1 = a(i)
v = Theta(i) - w(i)
ratio = (1 + e(i)) / (1 - e(i))
ratios = SQR(ratio)
EccAnn = 2! * ATN((TAN(v(i) * .5)) / ratios)
ce = COS(EccAnn)
se = SIN(EccAnn)
t = Theta(i)
om = w(i)
REM do l1 m1 n1  and l2 m2 and n2 of AE ROy page 102 ew 4.114 and 4.115
l1 = COS(t) * COS(om) - SIN(t) * SIN(om) * COS(In(i))
l2 = -COS(t) * SIN(om) - SIN(t) * COS(om) * COS(In(i))
REM
m1 = SIN(t) * COS(om) + COS(t) * SIN(om) * COS(In(i))
m2 = -SIN(t) * SIN(om) + COS(t) * COS(om) * COS(In(i))
REM
n1 = SIN(om) * SIN(In(i))
n2 = COS(om) * SIN(In(i))
REM

x(i) = a1 * l1 * ce + b * l2 * se - a1 * e(i) * l1
y(i) = a1 * m1 * ce + b * m2 * se - a1 * e(i) * m1
z(i) = a1 * n1 * ce + b * n2 * se - a1 * e(i) * n1

REM now do the dots
REM
R = a(i) * (1 - e(i) * ce)
per = SQR(a(i) ^ 3)
N = 2! * pi / (365.26 * per)
coef = N * a(i) / R
vx(i) = coef * (b * l2 * ce - a(i) * l1 * se)
vy(i) = coef * (b * m2 * ce - a(i) * m1 * se)
vz(i) = coef * (b * n2 * ce - a(i) * n1 * se)
vel = SQR(vx(i) ^ 2 + vy(i) ^ 2 + vz(i) ^ 2)

PRINT "Object #"; i;
PRINT "Speed ="; (vel * 149596000) / (24! * 3600!); " km/s"

REM PRINT "X,Y Z  "; x(i), y(i), z(i)
REM PRINT "V x y z "; vx(i), vy(i), vz(i)
PRINT #4, x(i), y(i), z(i), vx(i), vy(i), vz(i)

NEXT i
CLOSE #4
PRINT " Data written to XXDOT.DAT"
END