Re: sundial program
- From: david.williams@xxxxxxxxxx (David Williams)
- Date: Wed, 25 Feb 2009 21:54:27 -0600
-> On Feb 22, 8:20 pm, robert.wo...@xxxxxxxxxxxxxx (Robert Wolfe) wrote:
-> > Richard, do you have the source code available to this somewhere to downlo
-> > I'd like to run this through PowerBASIC (both DOS and Windows version) to
-> > what happens :)
-> Surely you need the QBASIC source, as originally listed here by
-> David? I can't imagine my machine-translated BBC BASIC source being a
-> lot of use to you! Or have I misunderstood?
-> Richard.
-> http://www.rtrussell.co.uk/
If he wants my QBasic source, I'll append it below. However, you did
make some valuable changes, such as making the printing routine work
with a modern printer. Maybe he would like to include those in his
versions.
dow
--------------------------------------------------------
' SUNDIAL.BAS. David Williams. 2002
' This version 2009 Feb 06
DECLARE SUB Lcopy11 ()
DECLARE SUB KeyPause ()
DECLARE SUB Choice (N%, C%, L$)
DECLARE FUNCTION YesNo$ ()
DECLARE FUNCTION ET (D)
CLS
' initialize array used in graphing routine
DIM ML(1 TO 13)
FOR X = 1 TO 13
READ ML(X)
NEXT
DATA 33,29,33,31,33,31,33,33,31,33,31,33,0
' Month lengths adjusted for numbers of pixels on screen
DATA 7, 97, 343, -15, 20: ' plotting parameters for E. of T. graph
READ ML%, VT%, VB%, LL%, UL%
CLS
PRINT "This program designs a horizontal sundial, for use at any"
PRINT "location. Three diagrams are produced. The first is the"
PRINT "horizontal plate, showing the hour lines. This should be"
PRINT "mounted so the north-south line shown on the design is aligned"
PRINT "true north-south, pointing away from the nearer pole (so the"
PRINT "line should point southward in the northern hemisphere)."
PRINT "The second diagram is a template for the sundial's gnomon."
PRINT "This should be mounted in a north-south vertical plane, in a"
PRINT "line with the north-south line on the horizontal plate. The"
PRINT "left end of the base of the gnomon should be at the point"
PRINT "where all the hour-lines intersect. The sloping edge of the"
PRINT "gnomon will be pointing at the nearer celestial pole, i.e. at"
PRINT "the Pole Star in the northern hemisphere."
PRINT "The third diagram is produced only if you have selected that"
PRINT "the dial show clock time. It is a graph that shows the number"
PRINT "of minutes that must be subtracted from the dial reading to"
PRINT "make it agree with a clock."
PRINT "When the program pauses at any point (such as right now),"
PRINT "press any key to continue."
CALL KeyPause
CLS
PRINT "You can select whether to have the diagrams appear only on the"
PRINT "screen, or if they should also be printed by an Epson-type"
PRINT "printer, and/or saved to disk for later printing by some other"
PRINT "utility. If you select to save them to disk, you will be asked"
PRINT "for a filename. Suppose you select FNAME. The three files on"
PRINT "disk will be named FNAME.DG1, FNAME.DG2 and FNAME.DG3. The"
PRINT "files are BSAVEd from screen memory. The first seven bytes of"
PRINT "each file contain addresses, etc.. The rest is a string of"
PRINT "bytes that are the bytes in memory representing the SCREEN 11"
PRINT "image."
CALL KeyPause
CLS
PRINT "Should the dial show <c>lock or local <s>olar time? ";
DO
K$ = LCASE$(INKEY$)
LOOP UNTIL K$ = "c" OR K$ = "s"
PRINT K$
PI = 4 * ATN(1)
PI2 = PI / 2
DR = PI / 180
RH = PI / 12
PRINT "Use negative numbers for directions opposite from those shown."
INPUT "Latitude (deg. N)"; LT
IF K$ = "c" THEN
INPUT "Longitude (deg. E)"; LG
INPUT "Time Zone Offset (hours ahead of GMT)"; TZ
A0 = TZ * RH - DR * LG
ELSE
A0 = 0
END IF
ST = SIN(DR * LT)
C% = 0
PRINT "Should the dial be printed on Epson-type printer";
IF YesNo$ = "Y" THEN
C% = 1
BEEP
PRINT "Ensure printer is ready, then press any key."
CALL KeyPause
END IF
PRINT "Should screens be saved to disk, for later printing";
IF YesNo$ = "Y" THEN
DO
INPUT "Filename (8 chars max, no suffix)"; F$
IF LEN(F$) > 0 AND LEN(F$) < 9 AND INSTR(F$, ".") = 0 THEN
EXIT DO
ELSE
BEEP
PRINT "Illegal entry!"
END IF
LOOP
C% = C% OR 2
END IF
' draw horizontal sundial plate
SCREEN 11
CLS
IF K$ = "c" THEN
NH% = 12 + TZ - LG / 15
ELSE
NH% = 12
END IF
FOR X% = NH% - 7 TO NH% + 7
A = RH * (X% - 12) - A0
Q% = SGN(A * LT)
D = ABS(A) - PI2
IF ABS(D) < .0000001 THEN
AZ = PI2 * Q%
ELSE
AZ = ATN(ST * TAN(A))
IF D > 0 THEN AZ = AZ + PI * Q%
END IF
' AZ is direction of hour-line in radians clockwise from pole
SZ = SIN(AZ) / 2
CA = AZ + ATN(SZ / SQR(1 - SZ * SZ))
' CA is direction from centre of dial to end of hour line
CS = SIN(CA)
CC = COS(CA)
LOCATE 15.2 - 13.5 * CC, 40.5 + 27.5 * CS
PRINT MID$(STR$(X% + 12 * (X% > 12)), 2);
LINE (320, 340)-(320 + 200 * CS, 240 - 200 * CC)
NEXT
CIRCLE (320, 240), 200
CIRCLE (320, 240), 239
LINE (320, 340)-(320, 479)
LOCATE 26, 39
PRINT "N - S";
LOCATE 27, 38
PRINT "L I N E";
CALL Choice(1, C%, F$)
'draw gnomon template
'gnomon size makes shadow of tip just reach
'circumference of dial at noon on summer solstice.
CLS
AT = 23.5 * DR' axial tilt
R = 150 / COS(AT)' radius of circumscribed circle (locus of tip)
A = 2 * ABS(LT) * DR - AT' tip angle above horizontal wrt centre
PX% = R * COS(A)
PY% = 300 - R * (SIN(AT) + SIN(A))
SELECT CASE PX%
CASE IS > 150
ZX% = 150
CASE IS < -75
ZX% = -75
CASE ELSE
ZX% = PX%
END SELECT
XS% = 357 - ZX% \ 2
OX% = XS% - 150
PX% = PX% + XS%
ZX% = ZX% + XS%
PRINT TAB(27); "GNOMON TEMPLATE";
LINE (OX%, 300)-(PX%, PY%)
LINE (OX%, 300)-(ZX%, 300)
LINE (ZX%, 300)-(PX%, PY%)
L$ = "South"
R$ = "North"
IF LT < 0 THEN SWAP L$, R$
LOCATE 20, 21
PRINT L$; " <- Base -> "; R$
PRINT TAB(15); "Gnomon should be in North-South vertical plane, with "
PRINT TAB(15); L$; " end of base at intersection of hour lines.";
CALL Choice(2, C%, F$)
IF K$ = "c" THEN ' do equation of time if clock time selected
' draw equation of time graph
CLS
PRINT TAB(30); "EQUATION OF TIME"
LOCATE 3, 15
PRINT "Graph shows difference in minutes between clock"
PRINT TAB(15); "and sundial time. Positive difference means"
PRINT TAB(15); "sundial is ahead of clock, and vice versa."
LOCATE 16, 67
PRINT "-="
FOR T = LL% TO UL% STEP 5
LINE (137, 247 - 6.4 * T)-(530, 247 - 6.4 * T)
LOCATE 16 - T / 2.5, 14
PRINT RIGHT$(" " + STR$(T), 3);
IF T = 0 THEN PRINT " ="
NEXT
LOCATE ML%, 20
PRINT "Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec"
' Do graphing
T = 0
D = 0
MC = 0
FOR M = 0 TO 384
M1 = 147 + M
IF M = T THEN
' Draw vertical line
LINE (M1, VT%)-(M1, VB%)
MC = MC + 1
T = T + ML(MC)
ELSE
' Plot point(s)
YY = 247 - 6.4 * ET(D) ' function used
PSET (M1, YY)
D = D + 1
END IF
IF ML(MC) = 33 AND M = T - 18 THEN M = M + 1
NEXT
CALL Choice(3, C%, F$)
END IF
SCREEN 0
END
SUB Choice (N%, C%, F$)
IF C% AND 2 THEN
F1$ = F$ + ".DG" + RIGHT$(STR$(N%), 1)
DEF SEG = &HA000
BSAVE F1$, 0, 80& * 480
DEF SEG
END IF
IF C% AND 1 THEN
CALL Lcopy11
ELSE
CALL KeyPause
END IF
END SUB
FUNCTION ET (D) STATIC
' Calculates equation of time in minutes on day number D of year.
' (D = 0 on January 1.)
' STATIC means variables are preserved between calls of function
IF E2 = 0 THEN ' first call, initialize constants
PI = 4 * ATN(1)
W = 2 * PI / 365 ' earth's mean orbital angular speed in radians/day
DR = 180 / PI ' degree/radian factor
C = -23.45 / DR ' reverse angle of earth's axial tilt in radians
CT = COS(C) ' cosine of reverse tilt
E2 = 2 * .0167 ' twice earth's orbital eccentricity
SP = 12 * W ' 12 days from December solstice to perihelion
END IF
A = W * (D + 10) ' Solstice 10 days before Jan 1
B = A + E2 * SIN(A - SP)
C = (A - ATN(TAN(B) / CT)) / PI
ET = 720 * (C - INT(C + .5))
' in 720 minutes, earth rotates PI radians relative to sun
END FUNCTION
SUB KeyPause
DO WHILE INKEY$ <> ""
LOOP
DO WHILE INKEY$ = ""
LOOP
END SUB
SUB Lcopy11
' Copies monochrome SCREEN 11 image to printer, with printed image
' aligned with long axis across the page. Note: To fit on 8-inch
' wide paper, 32-pixel margins on both sides are not printed.
DEFINT A-Y
DEFLNG Z
DIM V(0 TO 7), A(0 TO 7, 4 TO 75)' col's 0-3 and 76-79 not printed
V(7) = 1 ' reVerse-order bits
FOR X = 7 TO 1 STEP -1
V(X - 1) = V(X) + V(X)
NEXT
F = FREEFILE
OPEN "PRN" FOR BINARY AS F
DEF SEG = &HA000
OUT &H3CE, 4
OUT &H3CF, 0
N$ = STRING$(8, 0) ' 8 nulls for empty square
E$ = CHR$(27) ' ESC character
L$ = E$ + "C" + CHR$(0) + CHR$(11) ' page length = 11 inches
PUT F, , L$
L$ = E$ + "A" + CHR$(8) ' set printer to 8/72 lpi
PUT F, , L$
Z = 4
FOR S = 0 TO 59 ' 60 stripes, each 8 rows deep
N = -1 ' null-stripe flag
FOR R = 0 TO 7 ' read 8 rows into array
FOR C = 4 TO 75
A(R, C) = PEEK(Z)
IF N THEN IF A(R, C) THEN N = 0
Z = Z + 1
NEXT
Z = Z + 8
NEXT
IF NOT N THEN ' not null stripe
L$ = E$ + "*" + CHR$(5) + MKI$(576) ' set plotter mode
PUT F, , L$
FOR C = 4 TO 75
Q = 0
FOR R = 0 TO 7
IF A(R, C) THEN Q = 1: EXIT FOR
NEXT
IF Q THEN ' not null square
FOR Y = 0 TO 7 ' analyse 8 bits per byte
B = 0
FOR R = 0 TO 7 ' 8 bytes to be analysed
IF A(R, C) AND V(Y) THEN B = B OR V(R)
NEXT R
L$ = CHR$(B) ' send byte to printer
PUT F, , L$
NEXT Y
ELSE ' null square
PUT F, , N$
END IF
NEXT C
END IF
L$ = CHR$(13) + CHR$(10) ' terminate line (stripe)
PUT F, , L$
NEXT S
L$ = CHR$(12) ' form feed
PUT F, , L$
L$ = E$ + "@" ' reset printer
PUT F, , L$
CLOSE F
DEF SEG
END SUB
DEFSNG A-Z
FUNCTION YesNo$
PRINT "? (y/N) ";
T$ = " YyNn" + CHR$(13)
DO
LOOP UNTIL INKEY$ = ""
DO
K$ = INKEY$
LOOP UNTIL INSTR(T$, K$) > 1
IF K$ = CHR$(13) THEN
K$ = "N"
ELSE
K$ = UCASE$(K$)
END IF
PRINT K$
YesNo$ = K$
END FUNCTION
------------------------------------------------------
.
- Follow-Ups:
- Re: sundial program
- From: Richard Russell
- Re: sundial program
- Prev by Date: Re: sundial program
- Next by Date: Re: sundial program
- Previous by thread: Re: sundial program
- Next by thread: Re: sundial program
- Index(es):
Relevant Pages
|
