Program: Light.bas
Implementation Select: Implementation
Operation Select: Program: LIGHT.BAS
DECLARE SUB CALC1 ()
DECLARE SUB CALC2 (angle#, y#, display#)
DEFDBL A-Z
' LIGHT.BAS
' Revision 1.0 1 Dec 1995
' Added speed of gravity Revision 2.0 18 Jan 2003
' Configuration 1 is old configuration : Observer variable
' Configuration 2 Observer Rim-of-Sun and Star are on one line
DIM SHARED Rsun, c, dispscale, display, yscale, yoff, countdt
DIM SHARED scalefact, cg, vys0, rmin, earthdist, m
DIM SHARED count, sax, say, dt, xOB, yOB, voff, lnoff
DIM SHARED xsold, ysold ' Old display cord sun.
DIM SHARED xsave(10), ysave(10), isave
SCREEN 12
CONST ESC = 27, ENTER = 13
CONST pi = 3.141592653589#
CONST arcsec = 180# * 3600 / pi
c = 299792458 ' Speed of light m/sec
c = c / 1000 ' Speed of light in km/sec
earthdist = 149598023# ' Earth Sun distance
planet = 1 ' 1 = Sun 2=Jupiter
plname$ = "Sun" ' Sun
display = 0
lnoff = 4: ' ln offset
yoff = lnoff * 16 ' y offset of
dt = .001# ' delta t
countdt = INT(.1 / dt): ' delta t in counts
IF countdt < 1 THEN countdt = 1
IF countdt > 4000 THEN countdt = 4000
cg = 0 * c ' speed of gravity
scalefact = 1 ' scale factor (1 = normal)
vxs0 = 0: vys0 = 0 ' Initial speed Sun km/sec
conf = 1 ' Initial configuration
CLS
LOCATE 1, 25: PRINT " Light ray around the Sun"
COLOR 14: LOCATE 1, 60: PRINT "Configuration"; conf
GOTO c3
c0:
COLOR 14
ln = 2 ' line number = 2
LOCATE ln, 1: PRINT "Enter any key to continue; ESC to end";
DO: a$ = INKEY$: LOOP UNTIL a$ <> ""
IF a$ = CHR$(ESC) THEN END
c3:
ln = 2 ' line number = 2
CLS
COLOR 15
LOCATE 1, 25: PRINT " Light ray around the "; plname$
COLOR 14: LOCATE 1, 60: PRINT "Configuration"; conf
COLOR 15
LOCATE ln + 1, 1: PRINT "0 or Enter = Start Simulation "
COLOR 14
LOCATE ln + 2, 1: PRINT "1 = Configuration type "; TAB(50); "Configuration = "; conf
LOCATE ln + 3, 1: PRINT "2 = Scale factor "; TAB(50); "Scaletype = "; scalefact
LOCATE ln + 4, 1: PRINT "3 = Speed of gravity factor "; TAB(50); "cg = "; cg; "* c"
LOCATE ln + 5, 1: PRINT "4 = Initial Speed of "; plname$; TAB(50); "vys0 = "; vys0
LOCATE ln + 6, 1: PRINT "5 = Sun or Planet "; TAB(50); "planet = "; planet
COLOR 15
LOCATE ln + 7, 1: PRINT ">5 = End "
LOCATE ln, 1: PRINT SPACE$(70)
LOCATE ln, 1: INPUT "Enter parameter # "; typex ' y multiplication factor
SELECT CASE typex
CASE IS = 0
GOTO c5 ' Start Simulation
CASE IS = 1
LOCATE ln + 2, 5: INPUT "Enter configuration 1 or 2 "; conf1
IF conf1 = 1 OR conf1 = 2 THEN conf = conf1
GOTO c3
CASE IS = 2
LOCATE ln + 3, 5: INPUT "Scale factor >0 "; scalefactx
IF scalefactx > 0 THEN scalefact = scalefactx
GOTO c3
CASE IS = 3
LOCATE ln + 4, 5: INPUT "Speed of gravity factor "; cgx
IF cgx >= 0 THEN cg = cgx
GOTO c3
CASE IS = 4
LOCATE ln + 5, 5: PRINT "Initial Speed of "; plname$; : INPUT vys0
GOTO c3
CASE IS = 5
LOCATE ln + 6, 5: INPUT "Sun or Planet (1 or 2) "; planetx
IF planetx = 1 OR planetx = 2 THEN planet = planetx
IF planet = 1 THEN plname$ = "Sun "
IF planet = 2 THEN plname$ = "Jupiter"
GOTO c3
CASE ELSE
END
END SELECT
c5:
ln = 12
SELECT CASE planet
CASE 1
Rsun = 696000# ' Radius Sun
m = 6.672D-20 * 1.9891D+30 ' m = G times mass of Sun
dt = .001# ' delta t
countdt = INT(.1 / dt): ' delta t in counts
yscale = 100000 ' y scale
dangle = .000001# * pi ' configuration 2
CASE 2
Rsun = 143800# / 2 ' Radius Jupiter
m = 6.672D-20 * 1.8988D+27 ' m = G times mass of Jupiter
dt = .001# * Rsun / 696000 ' delta t
' countdt = INT(.1 / dt): ' delta t in counts
yscale = 10000000 ' y scale
dangle = .00000002# * pi ' configuration 2
END SELECT
CLS :
COLOR 15
LOCATE 1, 25: PRINT " Light ray around the "; plname$
COLOR 14: LOCATE 1, 60: PRINT "Configuration"; conf
yscale = yscale * scalefact
x = 0: y = 0
COLOR 7:
IF scalefact >= 1 THEN
CIRCLE (300 + x, 260 + yoff - y), 45 / scalefact' Display sun
ELSE
CIRCLE (300 + x, 260 + yoff - y), 45 ' Display sun
END IF
LOCATE 17 + lnoff, 40: PRINT RTRIM$(plname$); " (0,0)";
IF vys0 = 0 THEN PRINT " "
IF vys0 > 0 THEN PRINT " /"; CHR$(92) ' Print arrow head up
IF vys0 < 0 THEN PRINT " "; CHR$(92); "/" ' Print arrow head down
SELECT CASE conf
CASE 1
dispscale = 12000 ' common display scale
IF planet = 2 THEN dispscale = 1200
CALC1
GOTO c0
CASE 2
dispscale = 14000 ' common display scale
IF planet = 2 THEN dispscale = 1400
voff = 0
angle1 = pi
CALC2 angle1, y1, 0
angle3 = pi
c9:
a$ = INKEY$:
IF UCASE$(a$) = "E" THEN GOTO c0
IF a$ = CHR$(ESC) THEN GOTO c0
' Initial min and max angle calculation
angle3 = angle3 - dangle
CALC2 angle3, y3, 0
IF y3 > yOB THEN GOTO c10 ' Okay
angle1 = angle3: y1 = y3
GOTO c9
c10:
' Refinement
angle2 = (angle1 + angle3) / 2
CALC2 angle2, y2, 0
COLOR 14
' PRINT xOB, yOB
IF y1 < yOB AND y2 < yOB AND y1 <> y2 THEN
angle1 = angle2: y1 = y2
ELSE
IF y3 > yOB AND y2 > yOB AND y3 <> y2 THEN
angle3 = angle2: y3 = y2
ELSE
angle2 = (angle1 + angle3) / 2
CALC2 angle2, y2, 0
' COLOR 14
' LOCATE 2, 1: INPUT " Input any character"; a$
GOTO c0
END IF
END IF
angle2 = (angle1 + angle3) / 2
a$ = INKEY$:
IF UCASE$(a$) = "E" THEN GOTO c0
IF a$ = CHR$(ESC) THEN GOTO c0
GOTO c10
CASE ELSE
END
END SELECT
SUB CALC1
count = 0: sax = 0: say = 0:
vx0 = -c: vy0 = 0 ' Speed of light km/sec
x0 = 5 * Rsun: ' Initial x
y0 = scalefact * Rsun ' Initial y
' vxs0 = 0: vys0 = -1000 ' Initial speed Sun km/sec
r0 = SQR(x0 * x0 + y0 * y0)
IF cg > 0 THEN
ts = r0 / (cg * c) ' Initial time
ELSE
ts = 0
END IF
xs0 = 0: ys0 = -(r0 / c) * vys0 ' Initial position Sun
rmin = r0
x = x0: y = y0
vx = vx0: vy = vy0
xs = xs0: ys = ys0
vxs = vxs0: vys = vys0 ' Initialize speed of sun
COLOR 15: PSET (300 + xs0, 260 + yoff - ys0 / dispscale) ' Display sun
xsold = xs0: ysold = ys0 ' old cord. sun
end1 = 0
calc11:
count = count + 1
dy = y0 - y: dispy = y0 - yscale * dy
IF cg > 0 THEN
FOR i1 = 1 TO 2
xsv = xs - vxs * ts: ysv = ys - vys * ts ' virtual position sun
dx = x - xsv: dy = y - ysv
r2 = dx * dx + dy * dy:
r = SQR(r2) ' distance light Sun
ts = r / (cg * c)
NEXT i1
ELSE
dx = x - xs: dy = y - ys
r2 = dx * dx + dy * dy: r = SQR(r2) ' distance light Sun
END IF
a = m / r2
ax = -a * dx / r: ay = -a * dy / r
sax = sax + ax: say = say + ay ' ax and ay accu.
vx = vx0 + sax * dt: vy = vy0 + say * dt
v = SQR(vx * vx + vy * vy)
dx = vx * dt: dy = vy * dt
x = x + dx: y = y + dy ' position photon
dxs = vxs * dt: dys = vys * dt
xs = xs + dxs: ys = ys + dys ' position sun
IF ABS(x) > x0 AND count > 100 THEN end1 = 1
IF count MOD countdt = 0 THEN
COLOR 14
PSET (300 + x / dispscale, 260 + yoff - dispy / dispscale / scalefact)
END IF
IF count MOD countdt * 10 = 0 OR end1 = 1 THEN
COLOR 7
LOCATE 3, 1: PRINT USING "a ##.###########"; a;
PRINT " "; TAB(40); "v "; v; " ":
PRINT "x "; x; " "; TAB(40); "y "; y; " "
PRINT "r "; SQR(r2); " ";
al = dy / dx: al = ATN(al): alarcsec = al * arcsec: PRINT TAB(40); alarcsec; "arc sec "
COLOR 0: PSET (300 + xsold, 260 + yoff - ysold / dispscale) ' Display sun
COLOR 15: PSET (300 + xs, 260 + yoff - ys / dispscale) ' Display sun
xsold = xs: ysold = ys ' old cord. sun
END IF
COLOR 15
IF r < rmin THEN rmin = r: cond = 1
IF r > rmin AND cond = 1 THEN
x1 = x: y1 = y: cond = 2: count1 = count' Start of down hill
COLOR 0: PSET (300 + xsold, 260 + yoff - ysold / dispscale) ' Display sun
COLOR 15: PSET (300 + xs, 260 + yoff - ys / dispscale) ' Display sun
xsold = xs: ysold = ys ' old cord. sun
END IF
IF ABS(x) < x0 OR count < 100 THEN GOTO calc11
COLOR 0: PSET (300 + xsold, 260 + yoff - ysold / dispscale) ' Display sun
COLOR 15: PSET (300 + xs, 260 + yoff - ys / dispscale) ' Display sun
xsold = xs: ysold = ys ' old cord. sun
x2 = x: y2 = y: count2 = count ' Final position
PRINT "x0 "; x0; TAB(25); "y0 "; y0; TAB(50); "rmin ="; rmin
PRINT "x1 "; x1; TAB(25); "y1 "; y1; TAB(60); ' coord at min distance
PRINT USING "delta t = #.######"; dt
PRINT "x2 "; x2; TAB(25); "y2 "; y2; TAB(60); "cg "; cg ' coord at final distance
PRINT "dx "; dx; TAB(25);
PRINT USING "dy ##.##########"; dy;
PRINT TAB(60); "v sun "; vys0 ' angle
dist = earthdist * dy / dx
PRINT "Earth dist x"; earthdist; TAB(25); "dist y"; dist
' final angle
al = dy / dx: al = ATN(al): alarcsec = al * arcsec:
hoek = al * 180 / pi: hoek1 = INT(hoek)
hoek = 60 * (hoek - hoek1): hoek2 = INT(hoek)
hoek3 = 60 * (hoek - hoek2)
PRINT USING "alpha #.###########"; al; : PRINT " radian";
PRINT USING "#.##############"; TAB(30); alarcsec + ABS(alarcsecnew); : PRINT " arc sec "
dispy0 = y0
dispy1 = y0 - yscale * (y0 - y1)
dispy2 = y0 - yscale * (y0 - y2)
COLOR 7
' horizontal line x2,y0 - x1,y0
LINE (300 + x2 / dispscale, 260 + yoff - y0 / dispscale / scalefact)-(300 + x1 / dispscale, 260 + yoff - y0 / dispscale / scalefact)
' vertical line x2,y0 - x2,y2
LINE (300 + x2 / dispscale, 260 + yoff - y0 / dispscale / scalefact)-(300 + x2 / dispscale, 258 + yoff - dispy2 / dispscale / scalefact)
' helling x2,y2 - x1,y0
LINE (300 + x2 / dispscale, 258 + yoff - dispy2 / dispscale / scalefact)-(300 + x1 / dispscale, 260 + yoff - y0 / dispscale / scalefact)
COLOR 15
LOCATE 12 + lnoff, 9: PRINT "Earth dist x"
IF 4 / scalefact > 1 THEN
LOCATE 14 + lnoff, 26 - scalefact * 4: PRINT "alpha"
IF scalefact > .35 THEN
LOCATE 13 + 4 / scalefact + lnoff, 3: PRINT "dist y"
ELSE
LOCATE 28, 3: PRINT "dist y"
END IF
END IF
LOCATE 14 + lnoff, 72: PRINT "x0,y0"
LOCATE 15 + lnoff, 34: PRINT "x1,y1"
IF scalefact > .7 THEN LOCATE 14.5 + 8 / scalefact + lnoff, 1: PRINT "x2,y2"; ' start = 22
LOCATE 28, 30: PRINT "dist x = distance Earth and Sun"
LOCATE 29, 30: PRINT "y1, y2 and dist y are drawn on scale by"; INT(yscale);
' COLOR 14
' LOCATE 2, 1: INPUT " Input any character"; a$
END SUB
SUB CALC2 (angle, y, display)
count = 0: sax = 0: say = 0:
vx0 = c * COS(angle): vy0 = c * SIN(angle)
' PRINT vx0; vy0: END
x0 = 5 * Rsun: ' Initial x
y0 = scalefact * Rsun ' Initial y
xOB = -x0
yOB = y0
r0 = SQR(x0 * x0 + y0 * y0)
IF cg > 0 THEN
ts = r0 / (cg * c) ' Initial time
ELSE
ts = 0
END IF
xs0 = 0: ys0 = -(r0 / c) * vys0 + voff ' Initial position Sun
rmin = r0
x = x0: y = y0
vx = vx0: vy = vy0
xs = xs0: ys = ys0
vxs = vxs0: vys = vys0 ' Initialize speed of sun
COLOR 0: PSET (300 + xsold, 260 + yoff - ysold / dispscale) ' Display sun
COLOR 15: PSET (300 + xs0, 260 + yoff - ys0 / dispscale) ' Display sun
xsold = xs0: ysold = ys0 ' old cord. sun
dy = y0 - yOB: dispy = y0 - yscale * dy ' Observer
PSET (300 + xOB / dispscale, 260 + yoff - dispy / dispscale / scalefact)
anew = 0
end1 = 0
calc21:
count = count + 1
' PRINT x, y, r, r0; count
' IF count = 4 THEN END
dy = y0 - y: dispy = y0 - yscale * dy
IF cg > 0 THEN
FOR i1 = 1 TO 2
xsv = xs - vxs * ts: ysv = ys - vys * ts ' virtual position sun
dx = x - xsv: dy = y - ysv
r2 = dx * dx + dy * dy:
r = SQR(r2) ' distance light Sun
ts = r / (cg * c)
NEXT i1
ELSE
dx = x - xs: dy = y - ys
r2 = dx * dx + dy * dy: r = SQR(r2) ' distance light Sun
END IF
a = m / r2
ax = -a * dx / r: ay = -a * dy / r
sax = sax + ax: say = say + ay ' ax and ay accu.
vx = vx0 + sax * dt: vy = vy0 + say * dt
v = SQR(vx * vx + vy * vy)
dx = vx * dt: dy = vy * dt
x = x + dx: y = y + dy ' position photon
dxs = vxs * dt: dys = vys * dt
xs = xs + dxs: ys = ys + dys ' position sun
al = dy / dx: al = ATN(al): alarcsec = al * arcsec:
IF anew = 0 THEN anew = al: alarcsecnew = alarcsec
IF ABS(x) > x0 AND count > 100 THEN end1 = 1
IF count MOD countdt * 10 = 0 THEN
COLOR 0: PSET (300 + xsave(isave), 260 + ysave(isave)) ' Display sun
COLOR 14
PSET (300 + x / dispscale, 260 + yoff - dispy / dispscale / scalefact)
xsave(isave) = x / dispscale: ysave(isave) = yoff - dispy / dispscale / scalefact
isave = isave + 1
IF isave > 10 THEN isave = 0
END IF
IF count MOD countdt * 20 = 0 OR end1 = 1 THEN
COLOR 7
LOCATE 3, 1: PRINT USING "a ##.###########"; a;
PRINT " "; TAB(40); "v "; v; " ":
PRINT "x "; x; " "; TAB(40); "y "; y; " "
PRINT "r "; SQR(r2); " ";
PRINT TAB(40); ABS(alarcsec) + ABS(alarcsecnew); "arc sec "
COLOR 0: PSET (300 + xsold, 260 + yoff - ysold / dispscale) ' Display sun
COLOR 15: PSET (300 + xs, 260 + yoff - ys / dispscale) ' Display sun
xsold = xs: ysold = ys ' old cord. sun
END IF
COLOR 15
' PSET (300 + xs, 260 + yoff - ys / dispscale) ' Display sun
IF r < rmin THEN rmin = r: cond = 1
IF r > rmin AND cond = 1 THEN
x1 = x: y1 = y: cond = 2: count1 = count' Start of down hill
COLOR 0: PSET (300 + xsold, 260 + yoff - ysold / dispscale) ' Display sun
COLOR 15: PSET (300 + xs, 260 + yoff - ys / dispscale) ' Display sun
xsold = xs: ysold = ys ' old cord. sun
END IF
IF ABS(x) < x0 OR count < 100 THEN GOTO calc21
COLOR 14
COLOR 0: PSET (300 + xsold, 260 + yoff - ysold / dispscale) ' Display sun
COLOR 15: PSET (300 + xs, 260 + yoff - ys / dispscale) ' Display sun
xsold = xs: ysold = ys ' old cord. sun
x2 = x: y2 = y: count2 = count ' Final position
PRINT "x0 "; x0; TAB(25); "y0 "; y0; TAB(50); "rmin ="; rmin
PRINT "x1 "; x1; TAB(25); "y1 "; y1; TAB(60); ' coord at min distance
PRINT USING "delta t = #.######"; dt
PRINT "x2 "; x2; TAB(25); "y2 "; y2; TAB(60); "cg "; cg ' coord at final distance
PRINT "dx "; dx; TAB(25);
PRINT USING "dy ##.##########"; dy;
PRINT TAB(60); "v sun "; vys0 ' angle
' final angle
al = dy / dx: al = ATN(al): alarcsec = al * arcsec:
PRINT USING "alpha #.###########"; al + ABS(anew); : PRINT " radian";
PRINT USING "#.##############"; TAB(30); alarcsec + ABS(alarcsecnew); : PRINT " arc sec "
COLOR 7
LOCATE 14 + lnoff, 72: PRINT "x0,y0"
LOCATE 13 + lnoff, 35: PRINT "x1,y1"
LOCATE 14 + lnoff, 1: PRINT "x2,y2";
' LOCATE 17 + lnoff, 40: PRINT "Sun (0,0)"
LOCATE 29, 28: PRINT "y1, y2 and dist y are drawn on scale by"; INT(yscale);
END SUB
Back to my home page Contents of This Document