'd-Cephei.BAS Modified 27 July 2000 13 Jun 2009
'Robert S. Fritzius 305 Hillside Drive Starkville, MS 39759 U.S.A.
'fritzius@bellsouth.net see http://www.datasync.com/~rsf1/binaries.htm
'Method for calculating apparent radial velocity for C+V
'(non-Doppler) has been added in section 1300
GOSUB 50 ' Initialize variables
10 CLS
PRINT
PRINT " D-CEPHEI.BAS"
PRINT
PRINT " Menu"
PRINT
PRINT " (1) Review program information"
PRINT " (2) Select Run Parameters"
PRINT " ESC Quit"
15 Z$=INKEY$:IF Z$="" THEN 15
IF Z$="1" THEN 4000 'Review program information
IF Z$="2" THEN 20 'Select run options
IF Z$=ESC$ THEN END
GOTO 15
20 GOSUB 170 'Select Options
IF ESC = 1 THEN 10
GOSUB 270 'Run it
GOTO 20
50 'Initialize variables
DIM Xmax(3),Xmin(3)
DIM Ymax(3),Ymin(3)
'Ephemeris data Variable component
EV# = 0 'Eccentricity
VID = 0 'Iota (Tilt) Degrees
VOD = 270 'Omega Degrees (Ascending node to perihelion)
VODD = VOD 'Switchable value of VOD for two components
PVD = 5.4 'Sidereal Period in days
VTAU = 0 'Tau of perihelion (Not Used)
VWD = 180 'omega degrees Angle asc node to perihelion
VWR = VWd * P180 'omega radians Angle asc node to perihelion
K1 = 2.8 'Ellipse Generating constant
K2 = 1.35 'Ellipse eccentricity speed factor
MULT = .0065 'Orbit Radius Multiple of one A.U.
MULT1 = 0 'Orbit Radius used in calculations
NumComps = 1 'Number of Components 1 or 2
NumPoints = 3 'Number of emission points on surface of star
AU = 1.49E+11 'Astronomical Unit in meters
Pi = 3.141593 'Pi
TwoPi = 2 * Pi '2 Pi
P180 = Pi / 180 'Pi / 180
SecsDay = 86400! 'Seconds in a day
YrDays = 365.26 'Year in Days
YRS# = YrDays * SecsDay'Year in seconds
C = 300000000# 'Speed of Light M/Sec
LY# = C * YRS# 'Light Year in Meters
Lambda = 5.5E-07 'Wavelength (Meters) yellow-green
Xc# = 0 'X of Component 1
Yc# = 0 'Y of Component 1
Zc# = 0 'Z of Component 1
Dx# = 0
Dy# = 0
R# = 0 'Distance CG to Ellipse
RLY# = 0 'Distance Star-to-Observer in LY
Vx# = 0 'X Component variable speed of light
Vy# = 0 'Y Component variable speed of light
Vxy# = 0 'Orbital Velocity
'Positions of body in motion around binary C.O.G.
XS = 65 'X origin for sun on screen
YS = 65 'Y origin for sun on screen
VS = -20 / AU 'Vertical Scale for binary system
HS = 50 / AU 'Horizontal Scale for binary system
'Common Constants used on brightness and radial velocity graphs
Xorigin = 60 'X origin
Xmult = 3.1 'Horizontal scale ''
'Brightness versus Time display constants
YBo = 140 'Y Brightness origin
VB = -.7 'Vertical scale
'Radial Velocity versus Time display constants
YVo = 163 'Y Velocity origin
VRV = .8 * SQR(MULT) 'Vertical scale
'Photon positions in vincinity of observer constants
XO = 400 'X origin
YO = 75 'Y origin
VO = -40 / LY# 'Vertical Scale
HO = 100 / LY# 'Horizontal Scale
MODE = 1: MODE$ = "C+V"
DA = 1: DA$ = "dTau/dt"
NC = 2.25 'Number of Cycles (Orbits)
NB = 520 'Number of Brightness time bins
NS = 360 'Number of steps per orbit period
Radv$ = "######.### "
a$ = STRING$(79, " ")
ESC$ = CHR$(27)
Ax# = 0 'Acceleration of binary component in x direction
Temit = 0 'Time of Emisson at binary component
TOA# = 0 'Time Of Arrival at observer
Xold# = 0 'Old X, used in calc of Vx
Yold# = 0 'Old Y, used in calc of Vy
Zold# = 0 'Old Z, used in calc of Vz
Frac = 1! 'Fraction of orbit 2 size compared to orbit 1
Theta = 0 'Eclipse tracking angle
EC = 0 'Eclipse Flag 0 = off 1 = on
EC$ = "OFF" 'Eclipse Message
LBL = 20 'Lower Bin Limit (for display)
DIM NPR(NB)
TSF$ = "OFF" 'Tandy 1400HD Scale Factor string
TSF = 1.1 'Tandy 1400HD Scale Factor
COUNT=0 'Counter used to plot time tics
DIM TOA#(1000),DTOA#(1000)' Used to calculate DTDT
DELTAU=0 'd Tau
DTauDT=0 'd Tau / dt slope of app rad vel curve
RETURN
'********************************************************
170 '* Select Run Parameters *
'********************************************************
ESC = 0
200 CLS : SCREEN 2
OrbRad = MULT * AU
PVr = 5 * SecsDay
PVS = PVD * SecsDay 'Period of Variable in Seconds
Vref = TwoPi * MULT * AU / PVr
V = TwoPi * OrbRad / PVS
VB = -4 * Vref / V
L# = ((PVS / 2) * (C ^ 2) / (2 * V)) / LY#'This is a dimensionless ratio
IF DIST = 0 THEN DIST = .25 * L#
SRF = (29.6 / V) ^ 2 'Scale Reducing Factor Radial Velocity
YrFrac = PVD / YrDays
LOCATE 1, 1
PRINT "Binary Star Brightness versus Apparent Radial Velocity"
PRINT
PRINT "Options"
PRINT USING "(1) Orbit constant K1 ###.### "; K1
PRINT USING "(2) Orbit Radius in A.U.'s ##.#### "; MULT
PRINT USING "(3) Iota (Tilt) of orbit degrees ####.#### "; VID
PRINT USING "(4) Eccentricity of orbit #.#### "; EV#
PRINT USING "(5) Orbit period in days ###.### "; PVD
PRINT USING " Average orbit speed in Km/sec#####.### "; V * .001
PRINT USING "(6) Omega degrees (asc node) ####.# "; VOD
PRINT USING "(7) omega degrees (perihelion) ####.# "; VWD
PRINT USING " de Sitter Catch up distance ######.#### LY "; L#
PRINT USING "(8) Observer Distance ######.#### LY "; DIST
PRINT USING "(C) Modeling mode \ \ "; MODE$
IF MODE = 2 THEN PRINT a$: GOTO 205
PRINT USING "(D) Doppler Only or d Tau /dt \ \"; DA$
PRINT USING "(E) Eclipse function \ \ "; EC$
205 PRINT USING "(N) Number of components 1 or 2 # "; NumComps
IF NumComps = 1 THEN 208
PRINT USING "(P) Fraction radius Orbit 2 / Orbit 1 ##.## "; Frac
208 PRINT USING "(T) Tandy 1400HD Horiz Scale Factor \ \"; TSF$
PRINT "(Z) Restore original parameters"
LOCATE 5, 48: PRINT "At end of the run the"
LOCATE 6, 48: PRINT "display can be printed"
LOCATE 7, 48: PRINT "by pressing "
LOCATE 8, 48: PRINT "SHFT + PRINT SCREEN."
LOCATE
LOCATE 10, 48: PRINT "MS-DOS GRAPHICS driver"
LOCATE 11, 48: PRINT "must be pre-selected."
LOCATE 13, 48: PRINT "Otherwise press any key"
LOCATE 14, 48: PRINT "to start new run."
LOCATE 22, 1: PRINT STRING$(70, 32);
LOCATE 22, 1: PRINT "Which option (ENTER to execute) (ESC for MENU) ?";
ESC = 0
210 Z$ = INKEY$: IF Z$ = "" THEN 210
IF Z$ = ESC$ THEN ESC = 1: GOTO 250
IF Z$ = CHR$(13) THEN 250
IF Z$ = "1" THEN LOCATE 21, 1: INPUT "New K1 "; K1: GOTO 200
IF Z$ = "2" THEN LOCATE 21, 1: INPUT "New Orbit Radius a.u "; MULT: GOTO 200
IF Z$ = "3" THEN LOCATE 21, 1: INPUT "New Tilt degrees "; VID: GOTO 200
IF Z$ = "4" THEN LOCATE 21, 1: INPUT "New Eccentricity "; EV#: GOTO 200
IF Z$ = "5" THEN LOCATE 21, 1: INPUT "New Orbit period days "; PVD: GOTO 200
IF Z$ = "6" THEN LOCATE 21, 1: INPUT "New Omega degrees "; VOD: GOTO 200
IF Z$ = "7" THEN LOCATE 21, 1: INPUT "New omega Degrees "; VWd: GOTO 200
IF Z$ = "8" THEN LOCATE 21, 1: INPUT "New Observer Dist (LY) "; DIST: GOTO 200
Z = ASC(Z$)
IF Z > 96 AND Z < 123 THEN Z$ = CHR$(Z - 32)
IF Z$ = "C" AND MODE = 1 THEN MODE = 2: MODE$ = " C ": GOTO 200
IF Z$ = "C" AND MODE = 2 THEN MODE = 1: MODE$ = "C+V": GOTO 200
IF Z$ = "D" AND DA = 1 THEN DA = 2: DA$ = "DOPPLER": GOTO 200
IF Z$ = "D" AND DA = 2 THEN DA = 1: DA$ = "dTau/dt": GOTO 200
IF Z$="E" AND EC =0 THEN EC =1:EC$ ="ON": GOTO 200
IF Z$="E" AND EC =1 THEN EC =0:EC$ = "OFF": GOTO 200
IF Z$ = "N" AND NumComps = 1 THEN NumComps = 2: GOTO 200
IF Z$ = "N" AND NumComps = 2 THEN NumComps = 1: GOTO 200
IF Z$ = "P" THEN LOCATE 21, 1: INPUT "New Fraction "; Frac: GOTO 200
IF Z$ = "T" AND TSF = 1.1 THEN TSF = .7: TSF$ = " ON": GOTO 200
IF Z$ = "T" AND TSF = .7 THEN TSF = 1.1: TSF$ = "OFF": GOTO 200
IF Z$ = "Z" THEN RUN
GOTO 200
250 RETURN
270 CLS
PRINT USING "ORBIT RADIUS = #.### AU "; MULT;
PRINT USING "PERIOD = ##.## DAYS "; PVD;
PRINT USING "AVG ORBIT SPEED = ###.## KM/SEC"; V * .001
PRINT USING "DE SITTER DIST = ###.# LY "; L#;
PRINT USING "OBSERVER DIST = ###.##### LY "; DIST;
ON MODE GOTO 275, 280
275 PRINT USING "\ \"; DA$
GOTO 300
280 PRINT "DOPPLER C=CONST"
300 'LOCATE 3, 1: PRINT "Apparent radial velocity = Km/Sec"
PRINT USING "ORBIT ECCENTRICITY = #.###"; EV#;
LOCATE 3,30: PRINT USING "Perihelion = ### deg";VWd;
LOCATE 3, 55: PRINT "Space bar to pause/restart"
LOCATE 12, 9: PRINT "1";
LOCATE 9, 1: PRINT STRING$(30, 32);
IF TSF$ = " ON" THEN LOCATE 9, 1: PRINT " 4 2 ----> OBS";
IF TSF$ = "OFF" THEN LOCATE 9, 1: PRINT "4 2 ----> OBS";
LOCATE 5, 9: PRINT "3";
LOCATE 15, 3: PRINT "BRT";
LOCATE 19, 3: PRINT "RECEDING";
LOCATE 21, 3: PRINT "RAD VEL";
LOCATE 23, 3: PRINT "APPROACHING";
'ADJUST HORIZONTAL SCALE FOR SMALLER ELIPSES
'FOR I = 1 TO 15
' X = 16 + INT((I - 1) * 7.5 - .5)
' IF I > 0 AND I < 5 THEN J = I: GOTO 305
' IF I > 4 AND I < 9 THEN J = I - 4: GOTO 305
' IF I > 8 THEN J = I - 8
' IF X > 78 THEN 306
305 'LOCATE 24, X: PRINT STR$(J);
306 'NEXT I
FOR I=1 TO NumComps
Xmax(i)=-1e12:Xmin(i)=1e12
Ymax(i)=-1e12:Ymin(i)=1e12
NEXT I
RunNo = 1
'Zero "pulse arrival time" scaling bins
FOR I = 0 TO NB: NPR(I) = 0: NEXT I
'Ecliptic to Graph
310 ON RunNo GOTO 315,320
315 MULT1 = MULT
VODD = VOD
GOTO 330
320 MULT1=MULT*Frac
330 OrbRad = MULT1 * AU
VODD = VOD - 180
'Variable Component
VIR = VID * P180 'Iota (orbit tilt) radians
VOR = VODD * P180 'Omega radians (ascending node)
VWR = VWD * P180 'omega radians (asc node to perihelion)
VCI = COS(VIR)
VSI = SIN(VIR)
VCO = COS(VOR)
VSO = SIN(VOR)
VCW = COS(VWR)
VSW = SIN(VWR)
VL1 = VCO * VCW - VSO * VSW * VCI' + -
VM1 = VSO * VCW + VCO * VSW * VCI' + +
VN1 = VSW * VSI ' +
VL2 = -VCO * VSW - VSO * VCW * VCI' - -
VM2 = -VSO * VSW + VCO * VCW * VCI' - +
VN2 = VCW * VSI ' +
'Surface rotational speed components can be added to VX# and VY#
400 Xobs# = DIST * LY# 'X coord of observer in Km (from C.O.G.)
Yobs# = 0 'Y coord of observer in Km
Zobs# = 0 'Z coord of observer in KM
440 'Draw horizontal line for Radial Velocity info
Vrad# = 0
FOR bin% = LBL TO NB STEP 2
GOSUB 2400
NEXT bin%
450 'Set up initial conditions for ellipse generator
PVS = PVD * SecsDay 'Period of Component in Seconds
Vorb = TwoPi * OrbRad / PVS 'Average Orbit Speed in Meters/sec
'lprint using "Observer distance = ####.####### LY";DIST
'lprint " Temit DTavg Favg Vrel"
452 Dmul=(.7e-3)/(DIST)^.7 ' used as scale adjuster for distance
' in lines 1300+
on RunNo goto 454,456
454 'For First run
Ks=OrbRad*(1-EV#) 'Starting X comp (in plane of ellipse)
Nu = 0 'Starting Y comp
KsV = 0 'Starting Vx
NuV = (1 + K2 * EV#^1.1) * Vorb 'Starting Vy
goto 460
456 'For Second run
Ks=-OrbRad*(1-EV#) 'Starting X comp (in plane of ellipse)
Nu = 0 'Starting Y comp
KsV = 0 'Starting Vx
NuV = -(1 + K2 * EV#^1.1) * Vorb 'Starting Vy
460 STP = 1 / NS 'time between plotted positions
dt# = PVS * STP
RMAX# = 0
RMIN# = 2
D$ = "BM"'This is used to skip to plotted point without drawing a line
COUNT=0
FOR Temit = 0 TO NC STEP STP
COUNT=COUNT+1
for ii=1 to 9000:next ii 'speed reduction
'Calculate ellipse X's and Y's ( Ks's and Nu's)
R2# = (Ks ^ 2 + Nu ^ 2)'R Squared
R# = SQR(R2#)
'locate 22,1:print using "R/AU = ###.#### ";R#/AU
KdV = 1E+29 * (K1 / NS) * (MULT1 ^ 3 / R2#)
KsV = KsV - KdV * Ks / R#'Update X- component of Velocity
NuV = NuV - KdV * Nu / R#'Update Y- component of Velocity
Ks = Ks + KsV * dt#
Nu = Nu + NuV * dt#
465 Xc# = VL1 * Ks + VL2 * Nu'X for emitting component
Yc# = VM1 * Ks + VM2 * Nu
Zc# = VN1 * Ks + VN2 * Nu
IF Xc# = 0 THEN Theta = Pi / 2 * SGN(Yc#): GOTO 467
i=RunNo
'For orbit purity check
IF Xc#>Xmax(i) then Xmax(i)=Xc#
IF Xc#Ymax(i) then Ymax(i)=Yc#
IF Yc# 0 THEN Theta = Pi + Theta: GOTO 467'2nd Quadrant
IF Xc# < 0 AND Yc# < 0 THEN Theta = Pi + Theta: GOTO 467'3rd Quadrant
IF Xc# > 0 AND Yc# < 0 THEN Theta = 2 * Pi + Theta'4th Quadrant
467 'LOCATE 4,57:PRINT USING "Theta = ###.## radians";Theta;
GOSUB 2000 'Plot Variable Star positions
RLY# = (Xobs# - Xc#)
'This will give bogus Vx# on first time thru
IF Temit <= STP THEN 470
Vx# = (Xc# - Xold#) / dt#
Vy# = (Yc# - Yold#) / dt#
'Vz# = (Zc# - Zold#) / dt#
Vxy#=(Vx#^2+Vy#^2)^.5
470 Xold# = Xc#: 'X-Old's
Yold# = Yc#: 'Y-Old's
'Zold# = Zc#: 'Z-Old's
'Calculate Ax# 'acceleration of source w resp to obs
IF Temit <= STP THEN 480
'This will give a bogus result the first time thru
dVx# = (Vx# - OVx#)
480 OVx# = Vx# 'Old Vx
Ax# = 150000! * dVx# / dt# ' why 15e4?
ON MODE GOTO 910, 940
910 'Ballistic theory of light
ON DA GOTO 920, 930
920 'dt / D tau (Ritz's hypothesis)
'Vrad# = .5 * Ax# : 'don't use acceleration
'Vrad will be computed in 1300..
GOTO 1000
930 'Doppler (Sekerin's Ritzian with Doppler)
Vrad# = Vx#
GOTO 1000 '
940 'Wave theory of light
Vrad# = Vx#
GOTO 1000
1000 'Calculate Travel times
ON MODE GOTO 1100, 1200 'Mode 1 = C+V; Mode 2 = C
1100 'Ballistic speed of light C+V X,Y components
Sum# = C + Vx#
T# = (RLY# / Sum#) 'Travel time to Observer in secs
GOTO 1300
1200 'Constant speed of light C
T# = RLY# / C 'Travel time to Observer in secs
1300 'Calculate Time of Arrival (TOA)
TOA# = Temit * YrFrac + T# / YRS# - DIST 'TOA at observer
'Calculate d Tau/dt (Tau is modulated arrival time)
TOA#(COUNT)=TOA#*SecsDay*YrDays
DELTAU#=TOA#(COUNT)-TOA#(COUNT-1)
DTauDT#=ABS((DELTAU#)/dt#) 'd Tau/dt
'DTauDT#= (DELTAU#)/dt# 'd Tau/dt
IF MODE=2 THEN 1302
IF DA=2 THEN 1302
'Ballistic theory D Tau / dt radial velocity computation
Vrad#=((1/DTauDT# - 1) * c)*Dmul ' Vrad# for c+v
'See line 452 for Dmul
1302 'Z$=INKEY$:IF Z$="" THEN 1302
'if Temit>1.144 then lprint chr$(12):end
'Show arrival time modulation effects for second orbit
'******IF Temit > NC/2 THEN Vrad#=Vrad#/DTauDT
'Don't get rid of following. It's a trouble shooting feature
'Locate 16,1:print using "Te*YrFrac = #####.#### ";Temit*YrFrac;
'Locate 17,1:print using "T/yrs = #####.#### ";T#/YRS#;
'Locate 18,1:print using "DIST = #####.#### ";DIST;
'Locate 19,1:print using "Tarr = #####.#### ";TOA#
BIN% = INT(20 + (365 / PVD) * 80 * TOA#)
'BIN% is used for both Intensity and Apparent Radial velocity plots
IF EC=0 THEN 1370 ' Eclipse FCT is OFF ignore 1360
1360 IF Theta>.9*Pi AND Theta <1.1*Pi THEN 1400
1370 IF bin%<0 then bin%=NB+bin%:GOTO 1380
IF bin%>NB then bin%=bin%-NB
1380 NPR(bin%) = NPR(bin%) + 1
IF bin% < LBL OR bin% > NB THEN 1500
YB = NPR(bin%)
IF Temit <= STP THEN 1500
'Prog doesn't get right figures first iteration
1400 GOSUB 2300 'Plot intensity bin heights as they fill
GOSUB 2400 'Plot points on radial velocity curve
1450 'This shows the TOA of every twentieth arriving "sample"
'IF COUNT/10 <> INT(COUNT/10) THEN 1470
'Yg= INT(25 + COUNT/20)
'DRAW "BM" + STR$(Xg) + "," + STR$(Yg)
'DRAW "M" + STR$(Xg) + "," + STR$(Yg)
1470 Z$ = INKEY$: IF Z$ = "" THEN 1500
IF Z$ = CHR$(27) THEN Temit = NC: GOTO 1500
IF Z$ <> " " THEN 1500
1490 Z$ = INKEY$: IF Z$ = "" THEN 1490
IF Z$ = ESC$ THEN Temit = NC
1500 NEXT Temit
'Play it again Sam for Second star of binary
IF NumComps = 1 OR RunNo = 2 THEN 1900
'Flip second component to opposition of first
RunNo = 2
GOTO 310
1900 'Draw intensity envelope
bin% = LBL: YB = NPR(LBL): GOSUB 2500
FOR bin% = LBL + 1 TO NB - 1
YB = (NPR(bin% - 1) + 4 * NPR(bin%) + NPR(bin% + 1)) / 6
GOSUB 2600
NEXT bin%
LOCATE 4, 43: PRINT STRING$(30, " ");
'Locate 5,45:print using "Ecentricity = #.###";EV;
'Locate 6,45
'Print using "2* OrbRad = ##.###e8";2*OrbRad*1e-8
'For i=1 to NumComps
' Dx(i)=(Xmax(i)-Xmin(i))*1e-8
' Dy(i)=(Ymax(i)-Ymin(i))*1e-8
' DxDy(i)=Dx(i)/Dy(i)
' Locate 7+i,45
' print using "dx=##.###e8 dy=##.###e8 dx/dy=#.####";dx(i);dy(i);DxDy(i);
'Next i
'Pause for picture
1920 Z$ = INKEY$: IF Z$ = "" THEN 1920
LOCATE 23, 1: PRINT STRING$(78, 32);
LOCATE 23, 5: PRINT USING "Extinction Dist LY = #####.### "; DIST;
LOCATE 23, 37: INPUT "New Val (ENTER to change mode) "; Z
IF Z = 999 OR Z = 0 THEN RETURN
DIST = Z: GOTO 270
RETURN
2000 'Plot x,y positions of Component 1
'Xg = INT(HS * Xc# * TSF / MULT1 + XS + 5 * Temit)
Xg = INT(HS * Xc# * .8*TSF / MULT1 + XS + 4 * Temit)
'.8*TSF is to squeeze CRT display so screen dump comes out OK
Yg = INT(VS * Yc# * 1.1 / MULT1 + YS)
IF Xg < 0 OR Xg > 600 THEN RETURN
IF Yg < 0 OR Yg > 300 THEN RETURN
DRAW "BM" + STR$(Xg) + "," + STR$(Yg)
DRAW "M" + STR$(Xg) + "," + STR$(Yg)
RETURN
2300 'Draw dots in intensity bins
Xg = INT(Xmult * bin% + Xorigin)
Yg = INT(VB * YB + YBo)
IF Xg < 0 OR Xg > 680 THEN RETURN
IF Yg < 0 OR Yg > 300 THEN RETURN
DRAW "BM" + STR$(Xg) + "," + STR$(Yg)
DRAW "M" + STR$(Xg) + "," + STR$(Yg)
RETURN
2400 'Draw Radial Velocity info
Xg = INT(Xmult * bin% + Xorigin)
Yg = INT(VRV * Vrad# * SRF * 3000! + YVo)
'Yg = INT(-VRV * DTauDT * SRF * 3e7 + YVo)
IF Xg < 0 OR Xg > 680 THEN RETURN
IF Yg < 30 OR Yg > 300 THEN RETURN
DRAW "BM" + STR$(Xg) + "," + STR$(Yg)
DRAW "M" + STR$(Xg) + "," + STR$(Yg)
RETURN
2500 'Draw first point of intensity envelope
Xg = INT(Xmult * bin% + Xorigin)
Yg = INT(VB * YB + YBo)
IF Xg < 0 OR Xg > 600 THEN RETURN
IF Yg < 0 OR Yg > 300 THEN RETURN
DRAW "BM" + STR$(Xg) + "," + STR$(Yg)
DRAW "M" + STR$(Xg) + "," + STR$(Yg)
RETURN
2600 'Draw remaining points of intensity envelope
Xg = INT(Xmult * bin% + Xorigin)
Yg = INT(VB * YB + YBo)
'IF Xg<0 OR Xg>600 THEN RETURN
'IF Yg<0 OR Yg>300 THEN RETURN
DRAW "M" + STR$(Xg) + "," + STR$(Yg)
RETURN
4000 'Print program info
CLS
PRINT " D-CEPHEI"
PRINT " Robert S. Fritzius 305 Hillside Drive Starkville MS 39759 USA"
PRINT
PRINT " This program models V.I. Sekerin's (1988) hypothetical (C+V) modulation"
PRINT "effects on the observed brightness and apparent velocity of spectroscopic"
PRINT "binary stars. It uses the relation (C-V)/D = T + (C+V)/D to calculate the"
PRINT "distance required for slower light (component receeding from observer) to"
PRINT "be overtaken by the faster light (component approaching observer) where"
PRINT "T = one half orbital period and D is the distance to the observer.
PRINT " Orbital elements for one or two visual components of a binary star"
PRINT "are used to generate 2 1/4 cycles of orbital travel for each component."
PRINT "At one degree intervals each visual component 'emits' a light pulse toward"
PRINT "an observer who is located at the extinction distance from the binary."
PRINT "This extinction length is the distance at which essentially all of the"
PRINT "initial starlight has been absorbed and reradiated by the interstellar"
PRINT "medium, which is assumed to be stationary with respect to the observer.
PRINT " The program is written as though the light is unaffected by extinction"
PRINT "out to the extinction distance."
PRINT " Source-to-observer travel times and apparent velocity variations"
PRINT "are calculated for each pulse. These times are filed in brightness "
PRINT
PRINT "N = Next page M = Menu ESC = Quit"
4020 Z$=INKEY$:IF Z$="" THEN 4020
IF Z$=ESC$ THEN END
IF Z$="N" OR Z$="n" THEN 4030
IF Z$="M" OR Z$="m" THEN 10
4030 CLS
PRINT "time of arrival (TOA) bins. The brightness bins also serve as accumulators"
PRINT "so that their final summations show what an observer would 'see' in real time."
PRINT "In the C+V mode, bin filling is NOT a time-wise linear sequence. "
PRINT
PRINT "Pre-extinction light speed is selectable either as C or C+V, where V "
PRINT "is the visible source speed in the direction of the observer. Arrival"
PRINT "times and apparent velocities are calculated for each pulse. Arrival times"
PRINT "are 'filed' in brightness time of arrival (TOA) bins. These bins also serve"
PRINT "as accumulators so that their final summations show what an observer would "
PRINT "'see' in real time. In the C+V mode, bin filling is NOT a time-wise linear sequence. "
PRINT
PRINT "You may select Sekerin's Doppler (velocity) effects or Ritz's acceleration"
PRINT "effects to calculate observed apparent velocity variations. Uniquely pulsar-"
PRINT "like presentations occur where the extinction distance approximates or exceeds"
PRINT "0.5 of the 'overtake' distance."
PRINT
PRINT "The program's ellipse generator goes strange if you select an ellipticity"
PRINT "that exceeds 0.2.
PRINT
PRINT "P = Previous page Any other key to return to Menu"
4040 Z$ = INKEY$: IF Z$ = "" THEN 4040
IF Z$ = "P" OR Z$="p" THEN 4000
GOTO 10