'++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 'Draws four DXF files to make a planisphere. The first three are 'part of the base and the last is the moveable plate. Change 'value of phi for different latitudes. 'grid.dxf - 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.dxf - 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. 'conline.dxf - Generates constellation lines imported from the file 'conlines.dat extracted from John Walker's excellent freeware program 'Home Planet. I'm throwing away the constellation name information 'at present.... 'If you do not have a graphics program that allows you to super- 'impose the three files above, then just modify the program below 'to print all these layers to the same dxf file. 'horizon.dxf - Moveable plate showing the horizon 'window', and a 'point representing the zenith. A scale of hours runs around the 'outside of this plate. This dxf file MUST be separate to the base 'and must be printed on transparent material. 'To do: labels for scales. 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 PRINT "Stars done" CALL dclose 'print the 'conlines' layer dxfname$ = "conline" 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 conlines layer CALL dline(148, 150, 152, 150) CALL dline(150, 148, 150, 152) 'attempt at constellation lines using the CSV file conlines.csv 'from Home Planet by John Walker a$ = "CONLINES" OPEN a$ + ".DAT" FOR INPUT AS #2 DO WHILE NOT EOF(2) INPUT #2, conname$, r1, d1, r2, d2 'very crude clipping criterion - if either end 'of the current line is below the bottom of the 'star chart limit, don't draw it... IF d1 > -5500 AND d2 > -5500 THEN r1 = r1 / 1000 * 15 * rads d1 = d1 / 100 r2 = r2 / 1000 * 15 * rads d2 = d2 / 100 p1 = 90 - d1 p2 = 90 - d2 x1 = -p1 * COS(r1) y1 = p1 * SIN(r1) x2 = -p2 * COS(r2) y2 = p2 * SIN(r2) CALL dline(150 + x1, 150 + y1, 150 + x2, 150 + y2) END IF LOOP CLOSE #2 'complete the stars file PRINT "conlines 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 SUB darc (x, y, r, ang1, ang2) 'adds arc of radius r and centre (x,y), starting at 'ang1 in degrees and ending at ang2 degrees to 'current list of entities PRINT #1, 0 PRINT #1, "ARC" PRINT #1, 10 PRINT #1, x PRINT #1, 20 PRINT #1, y PRINT #1, 30 PRINT #1, 0 PRINT #1, 40 PRINT #1, r PRINT #1, 50 PRINT #1, ang1 PRINT #1, 51 PRINT #1, ang2 END SUB SUB dcircle (x, y, r) 'adds a circle of radius r and centre (x,y) to current 'list of entities PRINT #1, 0 PRINT #1, "CIRCLE" PRINT #1, 10 PRINT #1, x PRINT #1, 20 PRINT #1, y PRINT #1, 30 PRINT #1, 0 PRINT #1, 40 PRINT #1, r END SUB SUB dclose 'closes the entities list (ENDSEC) and prints the end of file 'code (EOF), then closes the qbasic disc file. PRINT #1, 0 PRINT #1, "ENDSEC" PRINT #1, 0 PRINT #1, "EOF" CLOSE 1 END SUB SUB dline (x1, y1, x2, y2) ' adds a line from (x1,y1) to (x2,y2) to list of entities PRINT #1, 0 PRINT #1, "LINE" PRINT #1, 8 PRINT #1, "0" PRINT #1, 10 PRINT #1, x1 PRINT #1, 20 PRINT #1, y1 PRINT #1, 30 PRINT #1, 0! PRINT #1, 11 PRINT #1, x2 PRINT #1, 21 PRINT #1, y2 PRINT #1, 31 PRINT #1, 0! END SUB SUB dpoint (x, y) 'adds a point to entities at absolute coordinates (x,y) PRINT #1, 0 PRINT #1, "POINT" PRINT #1, 10 PRINT #1, x PRINT #1, 20 PRINT #1, y PRINT #1, 30 PRINT #1, 0 END SUB SUB drect (x1, y1, x2, y2) 'adds a rectangle of bottom left corner x1, y1 and 'top right corner x2, y2 to current list of entities CALL dline(x1, y1, x1, y2) CALL dline(x1, y2, x2, y2) CALL dline(x2, y2, x2, y1) CALL dline(x2, y1, x1, y1) END SUB SUB init (a$) 'opens the file and sets up the first few lines of the 'dxf file OPEN "o", 1, a$ + ".dxf" PRINT #1, 0 PRINT #1, "SECTION" PRINT #1, 2 PRINT #1, "ENTITIES" END SUB