'++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
'Draws three PS files to make a planisphere. The first two are
'part of the base and the last is the moveable plate. Change
'value of phi for different latitudes.

'The PS code was contributed by Robert Lane and is provided here
'with his permission. He is working on a more elegant set of
'routines.

'grid.ps - RA grid with RA and DEC scale, ecliptic circle and a
'scale of days on the outside rim (365.25 to be exact). The
'ecliptic circle allows you to estimate the times of sunrise and
'set for a given day. Just put a straight edge between the pivot
'and the date - where the edge crosses the ecliptic marks the
'position of the Sun for that day. By turning the horizon plate
'until the Sun is just on the Western or Eastern horizon will
'set the time of sunset or sunrise on the hours scale on the
'horizon plate.

'stars.ps - Star map with 700 odd of the 907 Yale bright stars
'down to roughly Mag 4.5 that happen to have a declination above -55
'degrees. The data file yale.dat must be in the same directory.
'Put whatever you want in yale.dat, but each object has three
'numbers on a line - RA in decimal hours, DEC in decimal degrees
'and visual magnitude in decimals.

'horizon.ps - Moveable plate showing the horizon 'window', and a
'point representing the zenith. A scale of hours runs around the
'outside of this plate.

'To do: labels for scales, constellation or asterism lines.

DECLARE FUNCTION asin! (x!)
DECLARE FUNCTION atan2! (y!, x!)
DECLARE SUB Init (a$)
DECLARE SUB dLine (x1!, y1!, x2!, y2!)
DECLARE SUB dpoint (x!, y!)
DECLARE SUB dCircle (x!, y!, r!)
DECLARE SUB dArc (x!, y!, r!, ang1!, ang2!)
DECLARE SUB dRect (x1!, y1!, x2!, y2!)
DECLARE SUB endsec ()
DECLARE SUB dClose ()
CLS
PRINT "Planisphere program"
PRINT "-------------------"
PRINT

'constants
pi = 4 * ATN(1)
rads = pi / 180
degs = 180 / pi
phi = 52.5 * rads       'Alter this value for different latitudes
cotphi = 1 / TAN(phi)

dxfname$ = "grid"
CALL Init(dxfname$)

'plot a square framing the drawing to make the PSP import work ok
CALL dRect(0, 0, 300, 300)

'azimuthal projection - plot circles centred on celestial pole
'spaced equally by declination - distortion near horizon...
FOR i = 10 TO 150 STEP 10
CALL dCircle(150, 150, i)
NEXT i

'draw the RA lines
FOR i = 0 TO 23
x = 150 * COS(15 * rads * i)
y = 150 * SIN(15 * rads * i)
CALL dLine(150, 150, 150 + x, 150 + y)
NEXT i

'Now add a circle showing the days on the outer
'rim of the grid
CALL dCircle(150, 150, 147)
dayang = 360 / 365.25
FOR i = 0 TO 365.25
ang = i * dayang * rads
x1 = 147 * COS(ang)
y1 = 147 * SIN(ang)
x = 150 * COS(ang)
y = 150 * SIN(ang)
CALL dLine(150 + x1, 150 + y1, 150 + x, 150 + y)
NEXT i

'Now add the ecliptic circle plot
coseps = COS(rads * 23.439)
sineps = SIN(rads * 23.439)
'values at lambda = 0
x = 90
y = 0
'Each value of ecliptic longitude is treated as a point with zero
'ecliptic latitude and converted to the appropriate RA and DEC and
'then plotted in the usual way. The J2000 value for obliquity of
'the ecliptic is used....
FOR lambda = 1 TO 360
alpha = atan2(SIN(lambda * rads) * coseps, COS(lambda * rads))
delta = asin(sineps * SIN(lambda * rads)) * degs
'now find the polar coords and plot the line
x1 = (90 - delta) * COS(alpha)
y1 = (90 - delta) * SIN(alpha)
CALL dLine(150 + x, 150 + y, 150 + x1, 150 + y1)
'save current point for next line segment
x = x1
y = y1
NEXT lambda

'close the grid file
PRINT "Grid done"
CALL dClose

'print the 'stars' layer
dxfname$ = "stars"
CALL Init(dxfname$)
'plot a square framing the drawing to make the PSP import work ok
CALL dRect(0, 0, 300, 300)

'plot a circle same size as the base grid circumference
CALL dCircle(150, 150, 150)

'plot a cross in the centre of the stars layer
CALL dLine(148, 150, 152, 150)
CALL dLine(150, 148, 150, 152)

'plot 907 stars read from a file called yale.dat
a$ = "YALE"
OPEN a$ + ".dat" FOR INPUT AS #2
DO WHILE NOT EOF(2)
   INPUT #2, r, d, m
        polar = 90 - d
        IF polar < 145 THEN
                ra = r * rads * 15
                x = -polar * COS(ra) 'fixes flipped map problem
                y = polar * SIN(ra)
                'what function is best for magnitude??
                magn = .5 * SQR(5 - m)
                CALL dCircle(150 + x, 150 + y, magn)
   END IF
LOOP
CLOSE #2
'complete the stars file
PRINT "Stars done"
CALL dClose

'Now draw the horizon window file
dxfname$ = "horizon"
CALL Init(dxfname$)
'plot a square framing the drawing to make the PSP import work ok
CALL dRect(0, 0, 300, 300)
'plot a circle same size as the plate circumference, and one
'to act as inner ring for hours
CALL dCircle(150, 150, 147)
CALL dCircle(150, 150, 143)

'plot a cross in the centre of the plate to help with
'assembling the planisphere
CALL dLine(148, 150, 152, 150)
CALL dLine(150, 148, 150, 152)

'plot another small cross marking the zenith for the
'observer. The zenith will always be at 0 RA and a polar
'angle equal to 90 - latitude.
x = -(90 - phi * degs) * 1
y = (90 - phi * degs) * 0
CALL dLine(148 + x, 150 + y, 152 + x, 150 + y)
CALL dLine(150 + x, 148 + y, 150 + x, 152 + y)

'plot line segments that outline the horizon...
'set up the zeroth point
tand = cotphi * COS(rads * 0)
dec = ATN(tand)
polar = 90 - dec * degs
x = polar * COS(rads * 0)
y = polar * SIN(rads * 0)
'calculate the rest
FOR i = 1 TO 360
tand = cotphi * COS(rads * i)
dec = ATN(tand)
polar = 90 - dec * degs
x1 = polar * COS(rads * i)
y1 = polar * SIN(rads * i)
CALL dLine(150 + x, 150 + y, 150 + x1, 150 + y1)
x = x1
y = y1
NEXT i

'Now add a circle showing the hours from midnight
'This could be done in the loop above, but its easier
'to see what's going on having a separate one...
'This circle is a little smaller than the base circle
FOR i = 0 TO 23
ang = rads * 15 * i
x1 = 143 * COS(ang)
y1 = 143 * SIN(ang)
x = 147 * COS(ang)
y = 147 * SIN(ang)
CALL dLine(150 + x1, 150 + y1, 150 + x, 150 + y)
NEXT i

'finish the horizon plate
PRINT "Horizon window for latitude"; STR$(phi * degs); " done"
CALL dClose
END
'+++++++++++++++++++++++++++++++++++++++++++++++++++++++++

'***************************************************************************

FUNCTION asin (x)
'returns the arcsin of x in radians
c = SQR(1 - x * x)
asin = ATN(x / c)
END FUNCTION

FUNCTION atan2 (y, x)
'returns the inverse tan atan(y/x) in the
'correct quadrant
a = ATN(y / x)
pi = 4 * ATN(1)
IF x < o THEN a = a + pi
IF y < 0 AND x > 0 THEN a = a + pi + pi
atan2 = a
END FUNCTION

'Robert Lane's Post Script routines - used with permission
'**************************************************************************
SUB dArc (centerX, centerY, radius, startAngle, endAngle)
PRINT #1, "newpath 0.2 setlinewidth"
PRINT #1, centerX, centerY, radius, startAngle, endAngle
PRINT #1, "arc stroke"
PRINT #1,
END SUB

SUB dCircle (centerX, centerY, radius)
'Arc (centerX as double, centerY as double, radius, startAngle, endAngle)
'start and end angles are measured counterclockwise from positive x axis,
'in degrees
PRINT #1, "newpath 0.20 setlinewidth"
PRINT #1, centerX, centerY, radius, 0, 360
PRINT #1, "arc stroke"
PRINT #1,
END SUB

SUB dClose
PRINT #1, "showpage"
CLOSE 1
END SUB

SUB dLine (origX, origY, endX, endY)
PRINT #1,
PRINT #1, "newpath 0.2 setlinewidth"
PRINT #1, origX, origY, "moveto"
PRINT #1, endX, endY, "lineto"
PRINT #1, "stroke"
PRINT #1,
END SUB

SUB dRect (botX, botY, topX, topY)

PRINT #1, "newpath 0.2 setlinewidth"
PRINT #1, botX, botY, "moveto"
PRINT #1, botX, topY, "lineto"
PRINT #1, topX, topY, "lineto"
PRINT #1, topX, botY, "lineto"
PRINT #1, botX, botY, "lineto"
PRINT #1, "closepath stroke"

END SUB

SUB Init (thefile$)

OPEN thefile$ + ".ps" FOR OUTPUT AS #1
PRINT #1, "%!PS-Adobe-1.0"
PRINT #1, "%%BoundingBox: 0 0 300 300"
END SUB