1 8a3b2ceb 2004-04-24 devnull #include <u.h>
2 8a3b2ceb 2004-04-24 devnull #include <libc.h>
3 8a3b2ceb 2004-04-24 devnull #include <bio.h>
4 8a3b2ceb 2004-04-24 devnull #include <draw.h>
5 8a3b2ceb 2004-04-24 devnull #include <event.h>
6 8a3b2ceb 2004-04-24 devnull #include <ctype.h>
7 8a3b2ceb 2004-04-24 devnull #include "map.h"
8 8a3b2ceb 2004-04-24 devnull #undef RAD
9 8a3b2ceb 2004-04-24 devnull #undef TWOPI
10 8a3b2ceb 2004-04-24 devnull #include "sky.h"
12 8a3b2ceb 2004-04-24 devnull /* convert to milliarcsec */
13 347aa646 2009-04-30 rsc static int32 c = MILLIARCSEC; /* 1 degree */
14 347aa646 2009-04-30 rsc static int32 m5 = 1250*60*60; /* 5 minutes of ra */
16 8a3b2ceb 2004-04-24 devnull DAngle ramin;
17 8a3b2ceb 2004-04-24 devnull DAngle ramax;
18 8a3b2ceb 2004-04-24 devnull DAngle decmin;
19 8a3b2ceb 2004-04-24 devnull DAngle decmax;
20 8a3b2ceb 2004-04-24 devnull int folded;
22 8a3b2ceb 2004-04-24 devnull Image *grey;
23 8a3b2ceb 2004-04-24 devnull Image *alphagrey;
24 8a3b2ceb 2004-04-24 devnull Image *green;
25 8a3b2ceb 2004-04-24 devnull Image *lightblue;
26 8a3b2ceb 2004-04-24 devnull Image *lightgrey;
27 8a3b2ceb 2004-04-24 devnull Image *ocstipple;
28 8a3b2ceb 2004-04-24 devnull Image *suncolor;
29 8a3b2ceb 2004-04-24 devnull Image *mooncolor;
30 8a3b2ceb 2004-04-24 devnull Image *shadowcolor;
31 8a3b2ceb 2004-04-24 devnull Image *mercurycolor;
32 8a3b2ceb 2004-04-24 devnull Image *venuscolor;
33 8a3b2ceb 2004-04-24 devnull Image *marscolor;
34 8a3b2ceb 2004-04-24 devnull Image *jupitercolor;
35 8a3b2ceb 2004-04-24 devnull Image *saturncolor;
36 8a3b2ceb 2004-04-24 devnull Image *uranuscolor;
37 8a3b2ceb 2004-04-24 devnull Image *neptunecolor;
38 8a3b2ceb 2004-04-24 devnull Image *plutocolor;
39 8a3b2ceb 2004-04-24 devnull Image *cometcolor;
41 8a3b2ceb 2004-04-24 devnull Planetrec *planet;
43 347aa646 2009-04-30 rsc int32 mapx0, mapy0;
44 347aa646 2009-04-30 rsc int32 mapra, mapdec;
45 8a3b2ceb 2004-04-24 devnull double mylat, mylon, mysid;
46 8a3b2ceb 2004-04-24 devnull double mapscale;
47 8a3b2ceb 2004-04-24 devnull double maps;
48 8a3b2ceb 2004-04-24 devnull int (*projection)(struct place*, double*, double*);
50 4590ad02 2006-03-19 devnull char *fontname = "/lib/font/bit/luc/unicode.6.font";
52 8a3b2ceb 2004-04-24 devnull /* types Coord and Loc correspond to types in map(3) thus:
53 8a3b2ceb 2004-04-24 devnull Coord == struct coord;
54 8a3b2ceb 2004-04-24 devnull Loc == struct place, except longitudes are measured
55 8a3b2ceb 2004-04-24 devnull positive east for Loc and west for struct place
58 8a3b2ceb 2004-04-24 devnull typedef struct Xyz Xyz;
59 8a3b2ceb 2004-04-24 devnull typedef struct coord Coord;
60 8a3b2ceb 2004-04-24 devnull typedef struct Loc Loc;
62 8a3b2ceb 2004-04-24 devnull struct Xyz{
63 8a3b2ceb 2004-04-24 devnull double x,y,z;
66 8a3b2ceb 2004-04-24 devnull struct Loc{
67 8a3b2ceb 2004-04-24 devnull Coord nlat, elon;
70 8a3b2ceb 2004-04-24 devnull Xyz north = { 0, 0, 1 };
73 8a3b2ceb 2004-04-24 devnull plotopen(void)
75 8a3b2ceb 2004-04-24 devnull if(display != nil)
76 8a3b2ceb 2004-04-24 devnull return 1;
77 8a3b2ceb 2004-04-24 devnull if(initdraw(drawerror, nil, nil) < 0){
78 8a3b2ceb 2004-04-24 devnull fprint(2, "initdisplay failed: %r\n");
79 8a3b2ceb 2004-04-24 devnull return -1;
81 8a3b2ceb 2004-04-24 devnull grey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x777777FF);
82 8a3b2ceb 2004-04-24 devnull lightgrey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0xAAAAAAFF);
83 8a3b2ceb 2004-04-24 devnull alphagrey = allocimage(display, Rect(0, 0, 2, 2), RGBA32, 1, 0x77777777);
84 8a3b2ceb 2004-04-24 devnull green = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x00AA00FF);
85 8a3b2ceb 2004-04-24 devnull lightblue = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x009EEEFF);
86 8a3b2ceb 2004-04-24 devnull ocstipple = allocimage(display, Rect(0, 0, 2, 2), CMAP8, 1, 0xAAAAAAFF);
87 8a3b2ceb 2004-04-24 devnull draw(ocstipple, Rect(0, 0, 1, 1), display->black, nil, ZP);
88 8a3b2ceb 2004-04-24 devnull draw(ocstipple, Rect(1, 1, 2, 2), display->black, nil, ZP);
90 8a3b2ceb 2004-04-24 devnull suncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFF77FF);
91 8a3b2ceb 2004-04-24 devnull mooncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAAAFF);
92 8a3b2ceb 2004-04-24 devnull shadowcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x00000055);
93 8a3b2ceb 2004-04-24 devnull mercurycolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFAAAAFF);
94 8a3b2ceb 2004-04-24 devnull venuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
95 8a3b2ceb 2004-04-24 devnull marscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFF5555FF);
96 8a3b2ceb 2004-04-24 devnull jupitercolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFFAAFF);
97 8a3b2ceb 2004-04-24 devnull saturncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAFFAAFF);
98 8a3b2ceb 2004-04-24 devnull uranuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77DDDDFF);
99 8a3b2ceb 2004-04-24 devnull neptunecolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77FF77FF);
100 8a3b2ceb 2004-04-24 devnull plutocolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x7777FFFF);
101 8a3b2ceb 2004-04-24 devnull cometcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
102 8a3b2ceb 2004-04-24 devnull font = openfont(display, fontname);
103 8a3b2ceb 2004-04-24 devnull if(font == nil)
104 8a3b2ceb 2004-04-24 devnull fprint(2, "warning: no font %s: %r\n", fontname);
105 8a3b2ceb 2004-04-24 devnull return 1;
109 8a3b2ceb 2004-04-24 devnull * Function heavens() for setting up observer-eye-view
110 8a3b2ceb 2004-04-24 devnull * sky charts, plus subsidiary functions.
111 8a3b2ceb 2004-04-24 devnull * Written by Doug McIlroy.
114 8a3b2ceb 2004-04-24 devnull /* heavens(nlato, wlono, nlatc, wlonc) sets up a map(3)-style
115 8a3b2ceb 2004-04-24 devnull coordinate change (by calling orient()) and returns a
116 8a3b2ceb 2004-04-24 devnull projection function heavens for calculating a star map
117 8a3b2ceb 2004-04-24 devnull centered on nlatc,wlonc and rotated so it will appear
118 8a3b2ceb 2004-04-24 devnull upright as seen by a ground observer at nlato,wlono.
119 8a3b2ceb 2004-04-24 devnull all coordinates (north latitude and west longitude)
120 8a3b2ceb 2004-04-24 devnull are in degrees.
124 8a3b2ceb 2004-04-24 devnull mkcoord(double degrees)
126 8a3b2ceb 2004-04-24 devnull Coord c;
128 8a3b2ceb 2004-04-24 devnull c.l = degrees*PI/180;
129 8a3b2ceb 2004-04-24 devnull c.s = sin(c.l);
130 8a3b2ceb 2004-04-24 devnull c.c = cos(c.l);
131 8a3b2ceb 2004-04-24 devnull return c;
135 8a3b2ceb 2004-04-24 devnull ptov(struct place p)
139 8a3b2ceb 2004-04-24 devnull v.z = p.nlat.s;
140 8a3b2ceb 2004-04-24 devnull v.x = p.nlat.c * p.wlon.c;
141 8a3b2ceb 2004-04-24 devnull v.y = -p.nlat.c * p.wlon.s;
142 8a3b2ceb 2004-04-24 devnull return v;
146 8a3b2ceb 2004-04-24 devnull dot(Xyz a, Xyz b)
148 8a3b2ceb 2004-04-24 devnull return a.x*b.x + a.y*b.y + a.z*b.z;
152 8a3b2ceb 2004-04-24 devnull cross(Xyz a, Xyz b)
156 8a3b2ceb 2004-04-24 devnull v.x = a.y*b.z - a.z*b.y;
157 8a3b2ceb 2004-04-24 devnull v.y = a.z*b.x - a.x*b.z;
158 8a3b2ceb 2004-04-24 devnull v.z = a.x*b.y - a.y*b.x;
159 8a3b2ceb 2004-04-24 devnull return v;
163 8a3b2ceb 2004-04-24 devnull len(Xyz a)
165 8a3b2ceb 2004-04-24 devnull return sqrt(dot(a, a));
168 8a3b2ceb 2004-04-24 devnull /* An azimuth vector from a to b is a tangent to
169 8a3b2ceb 2004-04-24 devnull the sphere at a pointing along the (shorter)
170 8a3b2ceb 2004-04-24 devnull great-circle direction to b. It lies in the
171 8a3b2ceb 2004-04-24 devnull plane of the vectors a and b, is perpendicular
172 8a3b2ceb 2004-04-24 devnull to a, and has a positive compoent along b,
173 8a3b2ceb 2004-04-24 devnull provided a and b span a 2-space. Otherwise
174 8a3b2ceb 2004-04-24 devnull it is null. (aXb)Xa, where X denotes cross
175 8a3b2ceb 2004-04-24 devnull product, is such a vector. */
178 8a3b2ceb 2004-04-24 devnull azvec(Xyz a, Xyz b)
180 8a3b2ceb 2004-04-24 devnull return cross(cross(a,b), a);
183 8a3b2ceb 2004-04-24 devnull /* Find the azimuth of point b as seen
184 8a3b2ceb 2004-04-24 devnull from point a, given that a is distinct
185 8a3b2ceb 2004-04-24 devnull from either pole and from b */
188 8a3b2ceb 2004-04-24 devnull azimuth(Xyz a, Xyz b)
190 8a3b2ceb 2004-04-24 devnull Xyz aton = azvec(a,north);
191 8a3b2ceb 2004-04-24 devnull Xyz atob = azvec(a,b);
192 8a3b2ceb 2004-04-24 devnull double lenaton = len(aton);
193 8a3b2ceb 2004-04-24 devnull double lenatob = len(atob);
194 8a3b2ceb 2004-04-24 devnull double az = acos(dot(aton,atob)/(lenaton*lenatob));
196 8a3b2ceb 2004-04-24 devnull if(dot(b,cross(north,a)) < 0)
197 8a3b2ceb 2004-04-24 devnull az = -az;
198 8a3b2ceb 2004-04-24 devnull return az;
201 8a3b2ceb 2004-04-24 devnull /* Find the rotation (3rd argument of orient() in the
202 8a3b2ceb 2004-04-24 devnull map projection package) for the map described
203 8a3b2ceb 2004-04-24 devnull below. The return is radians; it must be converted
204 8a3b2ceb 2004-04-24 devnull to degrees for orient().
208 8a3b2ceb 2004-04-24 devnull rot(struct place center, struct place zenith)
210 8a3b2ceb 2004-04-24 devnull Xyz cen = ptov(center);
211 8a3b2ceb 2004-04-24 devnull Xyz zen = ptov(zenith);
213 8a3b2ceb 2004-04-24 devnull if(cen.z > 1-FUZZ) /* center at N pole */
214 8a3b2ceb 2004-04-24 devnull return PI + zenith.wlon.l;
215 8a3b2ceb 2004-04-24 devnull if(cen.z < FUZZ-1) /* at S pole */
216 8a3b2ceb 2004-04-24 devnull return -zenith.wlon.l;
217 8a3b2ceb 2004-04-24 devnull if(fabs(dot(cen,zen)) > 1-FUZZ) /* at zenith */
218 8a3b2ceb 2004-04-24 devnull return 0;
219 8a3b2ceb 2004-04-24 devnull return azimuth(cen, zen);
223 8a3b2ceb 2004-04-24 devnull heavens(double zlatdeg, double zlondeg, double clatdeg, double clondeg))(struct place*, double*, double*)
225 8a3b2ceb 2004-04-24 devnull struct place center;
226 8a3b2ceb 2004-04-24 devnull struct place zenith;
228 8a3b2ceb 2004-04-24 devnull center.nlat = mkcoord(clatdeg);
229 8a3b2ceb 2004-04-24 devnull center.wlon = mkcoord(clondeg);
230 8a3b2ceb 2004-04-24 devnull zenith.nlat = mkcoord(zlatdeg);
231 8a3b2ceb 2004-04-24 devnull zenith.wlon = mkcoord(zlondeg);
232 8a3b2ceb 2004-04-24 devnull projection = stereographic();
233 8a3b2ceb 2004-04-24 devnull orient(clatdeg, clondeg, rot(center, zenith)*180/PI);
234 8a3b2ceb 2004-04-24 devnull return projection;
238 347aa646 2009-04-30 rsc maptoxy(int32 ra, int32 dec, double *x, double *y)
240 8a3b2ceb 2004-04-24 devnull double lat, lon;
241 8a3b2ceb 2004-04-24 devnull struct place pl;
243 8a3b2ceb 2004-04-24 devnull lat = angle(dec);
244 8a3b2ceb 2004-04-24 devnull lon = angle(ra);
245 8a3b2ceb 2004-04-24 devnull pl.nlat.l = lat;
246 8a3b2ceb 2004-04-24 devnull pl.nlat.s = sin(lat);
247 8a3b2ceb 2004-04-24 devnull pl.nlat.c = cos(lat);
248 8a3b2ceb 2004-04-24 devnull pl.wlon.l = lon;
249 8a3b2ceb 2004-04-24 devnull pl.wlon.s = sin(lon);
250 8a3b2ceb 2004-04-24 devnull pl.wlon.c = cos(lon);
251 8a3b2ceb 2004-04-24 devnull normalize(&pl);
252 8a3b2ceb 2004-04-24 devnull return projection(&pl, x, y);
255 8a3b2ceb 2004-04-24 devnull /* end of 'heavens' section */
258 347aa646 2009-04-30 rsc setmap(int32 ramin, int32 ramax, int32 decmin, int32 decmax, Rectangle r, int zenithup)
261 8a3b2ceb 2004-04-24 devnull double minx, maxx, miny, maxy;
263 8a3b2ceb 2004-04-24 devnull c = 1000*60*60;
264 8a3b2ceb 2004-04-24 devnull mapra = ramax/2+ramin/2;
265 8a3b2ceb 2004-04-24 devnull mapdec = decmax/2+decmin/2;
267 8a3b2ceb 2004-04-24 devnull if(zenithup)
268 8a3b2ceb 2004-04-24 devnull projection = heavens(mylat, mysid, (double)mapdec/c, (double)mapra/c);
270 8a3b2ceb 2004-04-24 devnull orient((double)mapdec/c, (double)mapra/c, 0.);
271 8a3b2ceb 2004-04-24 devnull projection = stereographic();
273 8a3b2ceb 2004-04-24 devnull mapx0 = (r.max.x+r.min.x)/2;
274 8a3b2ceb 2004-04-24 devnull mapy0 = (r.max.y+r.min.y)/2;
275 8a3b2ceb 2004-04-24 devnull maps = ((double)Dy(r))/(double)(decmax-decmin);
276 8a3b2ceb 2004-04-24 devnull if(maptoxy(ramin, decmin, &minx, &miny) < 0)
277 8a3b2ceb 2004-04-24 devnull return -1;
278 8a3b2ceb 2004-04-24 devnull if(maptoxy(ramax, decmax, &maxx, &maxy) < 0)
279 8a3b2ceb 2004-04-24 devnull return -1;
281 8a3b2ceb 2004-04-24 devnull * It's tricky to get the scale right. This noble attempt scales
282 8a3b2ceb 2004-04-24 devnull * to fit the lengths of the sides of the rectangle.
284 8a3b2ceb 2004-04-24 devnull mapscale = Dx(r)/(minx-maxx);
285 8a3b2ceb 2004-04-24 devnull if(mapscale > Dy(r)/(maxy-miny))
286 8a3b2ceb 2004-04-24 devnull mapscale = Dy(r)/(maxy-miny);
288 8a3b2ceb 2004-04-24 devnull * near the pole It's not a rectangle, though, so this scales
289 8a3b2ceb 2004-04-24 devnull * to fit the corners of the trapezoid (a triangle at the pole).
291 8a3b2ceb 2004-04-24 devnull mapscale *= (cos(angle(mapdec))+1.)/2;
292 8a3b2ceb 2004-04-24 devnull if(maxy < miny){
293 8a3b2ceb 2004-04-24 devnull /* upside down, between zenith and pole */
294 8a3b2ceb 2004-04-24 devnull fprint(2, "reverse plot\n");
295 8a3b2ceb 2004-04-24 devnull mapscale = -mapscale;
297 8a3b2ceb 2004-04-24 devnull return 1;
301 347aa646 2009-04-30 rsc map(int32 ra, int32 dec)
303 8a3b2ceb 2004-04-24 devnull double x, y;
304 8a3b2ceb 2004-04-24 devnull Point p;
306 8a3b2ceb 2004-04-24 devnull if(maptoxy(ra, dec, &x, &y) > 0){
307 8a3b2ceb 2004-04-24 devnull p.x = mapscale*x + mapx0;
308 8a3b2ceb 2004-04-24 devnull p.y = mapscale*-y + mapy0;
310 8a3b2ceb 2004-04-24 devnull p.x = -100;
311 8a3b2ceb 2004-04-24 devnull p.y = -100;
313 8a3b2ceb 2004-04-24 devnull return p;
317 8a3b2ceb 2004-04-24 devnull dsize(int mag) /* mag is 10*magnitude; return disc size */
319 8a3b2ceb 2004-04-24 devnull double d;
321 8a3b2ceb 2004-04-24 devnull mag += 25; /* make mags all positive; sirius is -1.6m */
322 8a3b2ceb 2004-04-24 devnull d = (130-mag)/10;
323 8a3b2ceb 2004-04-24 devnull /* if plate scale is huge, adjust */
324 8a3b2ceb 2004-04-24 devnull if(maps < 100.0/MILLIARCSEC)
325 8a3b2ceb 2004-04-24 devnull d *= .71;
326 8a3b2ceb 2004-04-24 devnull if(maps < 50.0/MILLIARCSEC)
327 8a3b2ceb 2004-04-24 devnull d *= .71;
328 8a3b2ceb 2004-04-24 devnull return d;
332 8a3b2ceb 2004-04-24 devnull drawname(Image *scr, Image *col, char *s, int ra, int dec)
334 8a3b2ceb 2004-04-24 devnull Point p;
336 8a3b2ceb 2004-04-24 devnull if(font == nil)
338 8a3b2ceb 2004-04-24 devnull p = addpt(map(ra, dec), Pt(4, -1)); /* font has huge ascent */
339 8a3b2ceb 2004-04-24 devnull string(scr, p, col, ZP, font, s);
343 8a3b2ceb 2004-04-24 devnull npixels(DAngle diam)
345 8a3b2ceb 2004-04-24 devnull Point p0, p1;
347 8a3b2ceb 2004-04-24 devnull diam = DEG(angle(diam)*MILLIARCSEC); /* diam in milliarcsec */
348 8a3b2ceb 2004-04-24 devnull diam /= 2; /* radius */
349 8a3b2ceb 2004-04-24 devnull /* convert to plate scale */
350 8a3b2ceb 2004-04-24 devnull /* BUG: need +100 because we sometimes crash in library if we plot the exact center */
351 8a3b2ceb 2004-04-24 devnull p0 = map(mapra+100, mapdec);
352 8a3b2ceb 2004-04-24 devnull p1 = map(mapra+100+diam, mapdec);
353 8a3b2ceb 2004-04-24 devnull return hypot(p0.x-p1.x, p0.y-p1.y);
357 8a3b2ceb 2004-04-24 devnull drawdisc(Image *scr, DAngle semidiam, int ring, Image *color, Point pt, char *s)
361 8a3b2ceb 2004-04-24 devnull d = npixels(semidiam*2);
362 8a3b2ceb 2004-04-24 devnull if(d < 5)
364 8a3b2ceb 2004-04-24 devnull fillellipse(scr, pt, d, d, color, ZP);
365 8a3b2ceb 2004-04-24 devnull if(ring == 1)
366 8a3b2ceb 2004-04-24 devnull ellipse(scr, pt, 5*d/2, d/2, 0, color, ZP);
367 8a3b2ceb 2004-04-24 devnull if(ring == 2)
368 8a3b2ceb 2004-04-24 devnull ellipse(scr, pt, d, d, 0, grey, ZP);
370 8a3b2ceb 2004-04-24 devnull d = stringwidth(font, s);
371 8a3b2ceb 2004-04-24 devnull pt.x -= d/2;
372 8a3b2ceb 2004-04-24 devnull pt.y -= font->height/2 + 1;
373 8a3b2ceb 2004-04-24 devnull string(scr, pt, display->black, ZP, font, s);
378 8a3b2ceb 2004-04-24 devnull drawplanet(Image *scr, Planetrec *p, Point pt)
380 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "sun") == 0){
381 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, suncolor, pt, nil);
384 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "moon") == 0){
385 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, mooncolor, pt, nil);
388 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "shadow") == 0){
389 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 2, shadowcolor, pt, nil);
392 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "mercury") == 0){
393 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, mercurycolor, pt, "m");
396 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "venus") == 0){
397 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, venuscolor, pt, "v");
400 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "mars") == 0){
401 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, marscolor, pt, "M");
404 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "jupiter") == 0){
405 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, jupitercolor, pt, "J");
408 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "saturn") == 0){
409 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 1, saturncolor, pt, "S");
413 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "uranus") == 0){
414 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, uranuscolor, pt, "U");
418 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "neptune") == 0){
419 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, neptunecolor, pt, "N");
423 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "pluto") == 0){
424 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, plutocolor, pt, "P");
428 8a3b2ceb 2004-04-24 devnull if(strcmp(p->name, "comet") == 0){
429 8a3b2ceb 2004-04-24 devnull drawdisc(scr, p->semidiam, 0, cometcolor, pt, "C");
433 8a3b2ceb 2004-04-24 devnull pt.x -= stringwidth(font, p->name)/2;
434 8a3b2ceb 2004-04-24 devnull pt.y -= font->height/2;
435 8a3b2ceb 2004-04-24 devnull string(scr, pt, grey, ZP, font, p->name);
439 8a3b2ceb 2004-04-24 devnull tolast(char *name)
441 8a3b2ceb 2004-04-24 devnull int i, nlast;
442 8a3b2ceb 2004-04-24 devnull Record *r, rr;
444 8a3b2ceb 2004-04-24 devnull /* stop early to simplify inner loop adjustment */
445 8a3b2ceb 2004-04-24 devnull nlast = 0;
446 8a3b2ceb 2004-04-24 devnull for(i=0,r=rec; i<nrec-nlast; i++,r++)
447 8a3b2ceb 2004-04-24 devnull if(r->type == Planet)
448 e9ad8994 2005-05-01 devnull if(name==nil || strcmp(r->u.planet.name, name)==0){
449 8a3b2ceb 2004-04-24 devnull rr = *r;
450 8a3b2ceb 2004-04-24 devnull memmove(rec+i, rec+i+1, (nrec-i-1)*sizeof(Record));
451 8a3b2ceb 2004-04-24 devnull rec[nrec-1] = rr;
454 8a3b2ceb 2004-04-24 devnull nlast++;
459 347aa646 2009-04-30 rsc bbox(int32 extrara, int32 extradec, int quantize)
462 8a3b2ceb 2004-04-24 devnull Record *r;
463 8a3b2ceb 2004-04-24 devnull int ra, dec;
464 8a3b2ceb 2004-04-24 devnull int rah, ram, d1, d2;
465 8a3b2ceb 2004-04-24 devnull double r0;
467 8a3b2ceb 2004-04-24 devnull ramin = 0x7FFFFFFF;
468 8a3b2ceb 2004-04-24 devnull ramax = -0x7FFFFFFF;
469 8a3b2ceb 2004-04-24 devnull decmin = 0x7FFFFFFF;
470 8a3b2ceb 2004-04-24 devnull decmax = -0x7FFFFFFF;
472 8a3b2ceb 2004-04-24 devnull for(i=0,r=rec; i<nrec; i++,r++){
473 8a3b2ceb 2004-04-24 devnull if(r->type == Patch){
474 8a3b2ceb 2004-04-24 devnull radec(r->index, &rah, &ram, &dec);
475 8a3b2ceb 2004-04-24 devnull ra = 15*rah+ram/4;
476 8a3b2ceb 2004-04-24 devnull r0 = c/cos(dec*PI/180);
477 8a3b2ceb 2004-04-24 devnull ra *= c;
478 8a3b2ceb 2004-04-24 devnull dec *= c;
479 8a3b2ceb 2004-04-24 devnull if(dec == 0)
480 8a3b2ceb 2004-04-24 devnull d1 = c, d2 = c;
481 8a3b2ceb 2004-04-24 devnull else if(dec < 0)
482 8a3b2ceb 2004-04-24 devnull d1 = c, d2 = 0;
484 8a3b2ceb 2004-04-24 devnull d1 = 0, d2 = c;
485 8a3b2ceb 2004-04-24 devnull }else if(r->type==SAO || r->type==NGC || r->type==Planet || r->type==Abell){
486 e9ad8994 2005-05-01 devnull ra = r->u.ngc.ra;
487 e9ad8994 2005-05-01 devnull dec = r->u.ngc.dec;
488 8a3b2ceb 2004-04-24 devnull d1 = 0, d2 = 0, r0 = 0;
490 8a3b2ceb 2004-04-24 devnull continue;
491 8a3b2ceb 2004-04-24 devnull if(dec+d2+extradec > decmax)
492 8a3b2ceb 2004-04-24 devnull decmax = dec+d2+extradec;
493 8a3b2ceb 2004-04-24 devnull if(dec-d1-extradec < decmin)
494 8a3b2ceb 2004-04-24 devnull decmin = dec-d1-extradec;
495 8a3b2ceb 2004-04-24 devnull if(folded){
496 8a3b2ceb 2004-04-24 devnull ra -= 180*c;
497 8a3b2ceb 2004-04-24 devnull if(ra < 0)
498 8a3b2ceb 2004-04-24 devnull ra += 360*c;
500 8a3b2ceb 2004-04-24 devnull if(ra+r0+extrara > ramax)
501 8a3b2ceb 2004-04-24 devnull ramax = ra+r0+extrara;
502 8a3b2ceb 2004-04-24 devnull if(ra-extrara < ramin)
503 8a3b2ceb 2004-04-24 devnull ramin = ra-extrara;
506 8a3b2ceb 2004-04-24 devnull if(decmax > 90*c)
507 8a3b2ceb 2004-04-24 devnull decmax = 90*c;
508 8a3b2ceb 2004-04-24 devnull if(decmin < -90*c)
509 8a3b2ceb 2004-04-24 devnull decmin = -90*c;
510 8a3b2ceb 2004-04-24 devnull if(ramin < 0)
511 8a3b2ceb 2004-04-24 devnull ramin += 360*c;
512 8a3b2ceb 2004-04-24 devnull if(ramax >= 360*c)
513 8a3b2ceb 2004-04-24 devnull ramax -= 360*c;
515 8a3b2ceb 2004-04-24 devnull if(quantize){
516 8a3b2ceb 2004-04-24 devnull /* quantize to degree boundaries */
517 8a3b2ceb 2004-04-24 devnull ramin -= ramin%m5;
518 8a3b2ceb 2004-04-24 devnull if(ramax%m5 != 0)
519 8a3b2ceb 2004-04-24 devnull ramax += m5-(ramax%m5);
520 8a3b2ceb 2004-04-24 devnull if(decmin > 0)
521 8a3b2ceb 2004-04-24 devnull decmin -= decmin%c;
523 8a3b2ceb 2004-04-24 devnull decmin -= c - (-decmin)%c;
524 8a3b2ceb 2004-04-24 devnull if(decmax > 0){
525 8a3b2ceb 2004-04-24 devnull if(decmax%c != 0)
526 8a3b2ceb 2004-04-24 devnull decmax += c-(decmax%c);
527 8a3b2ceb 2004-04-24 devnull }else if(decmax < 0){
528 8a3b2ceb 2004-04-24 devnull if(decmax%c != 0)
529 8a3b2ceb 2004-04-24 devnull decmax += ((-decmax)%c);
533 8a3b2ceb 2004-04-24 devnull if(folded){
534 8a3b2ceb 2004-04-24 devnull if(ramax-ramin > 270*c){
535 8a3b2ceb 2004-04-24 devnull fprint(2, "ra range too wide %ld°\n", (ramax-ramin)/c);
536 8a3b2ceb 2004-04-24 devnull return -1;
538 8a3b2ceb 2004-04-24 devnull }else if(ramax-ramin > 270*c){
539 8a3b2ceb 2004-04-24 devnull folded = 1;
540 8a3b2ceb 2004-04-24 devnull return bbox(0, 0, quantize);
543 8a3b2ceb 2004-04-24 devnull return 1;
547 8a3b2ceb 2004-04-24 devnull inbbox(DAngle ra, DAngle dec)
549 8a3b2ceb 2004-04-24 devnull int min;
551 8a3b2ceb 2004-04-24 devnull if(dec < decmin)
552 8a3b2ceb 2004-04-24 devnull return 0;
553 8a3b2ceb 2004-04-24 devnull if(dec > decmax)
554 8a3b2ceb 2004-04-24 devnull return 0;
555 8a3b2ceb 2004-04-24 devnull min = ramin;
556 8a3b2ceb 2004-04-24 devnull if(ramin > ramax){ /* straddling 0h0m */
557 8a3b2ceb 2004-04-24 devnull min -= 360*c;
558 8a3b2ceb 2004-04-24 devnull if(ra > 180*c)
559 8a3b2ceb 2004-04-24 devnull ra -= 360*c;
561 8a3b2ceb 2004-04-24 devnull if(ra < min)
562 8a3b2ceb 2004-04-24 devnull return 0;
563 8a3b2ceb 2004-04-24 devnull if(ra > ramax)
564 8a3b2ceb 2004-04-24 devnull return 0;
565 8a3b2ceb 2004-04-24 devnull return 1;
569 347aa646 2009-04-30 rsc gridra(int32 mapdec)
571 8a3b2ceb 2004-04-24 devnull mapdec = abs(mapdec)+c/2;
572 8a3b2ceb 2004-04-24 devnull mapdec /= c;
573 8a3b2ceb 2004-04-24 devnull if(mapdec < 60)
574 8a3b2ceb 2004-04-24 devnull return m5;
575 8a3b2ceb 2004-04-24 devnull if(mapdec < 80)
576 8a3b2ceb 2004-04-24 devnull return 2*m5;
577 8a3b2ceb 2004-04-24 devnull if(mapdec < 85)
578 8a3b2ceb 2004-04-24 devnull return 4*m5;
579 8a3b2ceb 2004-04-24 devnull return 8*m5;
582 8a3b2ceb 2004-04-24 devnull #define GREY (nogrey? display->white : grey)
585 8a3b2ceb 2004-04-24 devnull plot(char *flags)
587 8a3b2ceb 2004-04-24 devnull int i, j, k;
588 8a3b2ceb 2004-04-24 devnull char *t;
590 8a3b2ceb 2004-04-24 devnull int ra, dec;
592 8a3b2ceb 2004-04-24 devnull Point p, pts[10];
593 8a3b2ceb 2004-04-24 devnull Record *r;
594 8a3b2ceb 2004-04-24 devnull Rectangle rect, r1;
595 8a3b2ceb 2004-04-24 devnull int dx, dy, nogrid, textlevel, nogrey, zenithup;
596 8a3b2ceb 2004-04-24 devnull Image *scr;
597 8a3b2ceb 2004-04-24 devnull char *name, buf[32];
598 8a3b2ceb 2004-04-24 devnull double v;
600 8a3b2ceb 2004-04-24 devnull if(plotopen() < 0)
602 8a3b2ceb 2004-04-24 devnull nogrid = 0;
603 8a3b2ceb 2004-04-24 devnull nogrey = 0;
604 8a3b2ceb 2004-04-24 devnull textlevel = 1;
605 8a3b2ceb 2004-04-24 devnull dx = 512;
606 8a3b2ceb 2004-04-24 devnull dy = 512;
607 8a3b2ceb 2004-04-24 devnull zenithup = 0;
608 8a3b2ceb 2004-04-24 devnull for(;;){
609 8a3b2ceb 2004-04-24 devnull if(t = alpha(flags, "nogrid")){
610 8a3b2ceb 2004-04-24 devnull nogrid = 1;
611 8a3b2ceb 2004-04-24 devnull flags = t;
612 8a3b2ceb 2004-04-24 devnull continue;
614 8a3b2ceb 2004-04-24 devnull if((t = alpha(flags, "zenith")) || (t = alpha(flags, "zenithup")) ){
615 8a3b2ceb 2004-04-24 devnull zenithup = 1;
616 8a3b2ceb 2004-04-24 devnull flags = t;
617 8a3b2ceb 2004-04-24 devnull continue;
619 8a3b2ceb 2004-04-24 devnull if((t = alpha(flags, "notext")) || (t = alpha(flags, "nolabel")) ){
620 8a3b2ceb 2004-04-24 devnull textlevel = 0;
621 8a3b2ceb 2004-04-24 devnull flags = t;
622 8a3b2ceb 2004-04-24 devnull continue;
624 8a3b2ceb 2004-04-24 devnull if((t = alpha(flags, "alltext")) || (t = alpha(flags, "alllabel")) ){
625 8a3b2ceb 2004-04-24 devnull textlevel = 2;
626 8a3b2ceb 2004-04-24 devnull flags = t;
627 8a3b2ceb 2004-04-24 devnull continue;
629 8a3b2ceb 2004-04-24 devnull if(t = alpha(flags, "dx")){
630 8a3b2ceb 2004-04-24 devnull dx = strtol(t, &t, 0);
631 8a3b2ceb 2004-04-24 devnull if(dx < 100){
632 8a3b2ceb 2004-04-24 devnull fprint(2, "dx %d too small (min 100) in plot\n", dx);
635 8a3b2ceb 2004-04-24 devnull flags = skipbl(t);
636 8a3b2ceb 2004-04-24 devnull continue;
638 8a3b2ceb 2004-04-24 devnull if(t = alpha(flags, "dy")){
639 8a3b2ceb 2004-04-24 devnull dy = strtol(t, &t, 0);
640 8a3b2ceb 2004-04-24 devnull if(dy < 100){
641 8a3b2ceb 2004-04-24 devnull fprint(2, "dy %d too small (min 100) in plot\n", dy);
644 8a3b2ceb 2004-04-24 devnull flags = skipbl(t);
645 8a3b2ceb 2004-04-24 devnull continue;
647 8a3b2ceb 2004-04-24 devnull if((t = alpha(flags, "nogrey")) || (t = alpha(flags, "nogray"))){
648 8a3b2ceb 2004-04-24 devnull nogrey = 1;
649 8a3b2ceb 2004-04-24 devnull flags = skipbl(t);
650 8a3b2ceb 2004-04-24 devnull continue;
652 8a3b2ceb 2004-04-24 devnull if(*flags){
653 8a3b2ceb 2004-04-24 devnull fprint(2, "syntax error in plot\n");
658 8a3b2ceb 2004-04-24 devnull flatten();
659 8a3b2ceb 2004-04-24 devnull folded = 0;
661 8a3b2ceb 2004-04-24 devnull if(bbox(0, 0, 1) < 0)
663 8a3b2ceb 2004-04-24 devnull if(ramax-ramin<100 || decmax-decmin<100){
664 8a3b2ceb 2004-04-24 devnull fprint(2, "plot too small\n");
668 8a3b2ceb 2004-04-24 devnull scr = allocimage(display, Rect(0, 0, dx, dy), RGB24, 0, DBlack);
669 8a3b2ceb 2004-04-24 devnull if(scr == nil){
670 8a3b2ceb 2004-04-24 devnull fprint(2, "can't allocate image: %r\n");
673 8a3b2ceb 2004-04-24 devnull rect = scr->r;
674 8a3b2ceb 2004-04-24 devnull rect.min.x += 16;
675 8a3b2ceb 2004-04-24 devnull rect = insetrect(rect, 40);
676 8a3b2ceb 2004-04-24 devnull if(setmap(ramin, ramax, decmin, decmax, rect, zenithup) < 0){
677 8a3b2ceb 2004-04-24 devnull fprint(2, "can't set up map coordinates\n");
680 8a3b2ceb 2004-04-24 devnull if(!nogrid){
681 8a3b2ceb 2004-04-24 devnull for(x=ramin; ; ){
682 8a3b2ceb 2004-04-24 devnull for(j=0; j<nelem(pts); j++){
683 8a3b2ceb 2004-04-24 devnull /* use double to avoid overflow */
684 8a3b2ceb 2004-04-24 devnull v = (double)j / (double)(nelem(pts)-1);
685 8a3b2ceb 2004-04-24 devnull v = decmin + v*(decmax-decmin);
686 8a3b2ceb 2004-04-24 devnull pts[j] = map(x, v);
688 8a3b2ceb 2004-04-24 devnull bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
690 8a3b2ceb 2004-04-24 devnull if(folded){
691 8a3b2ceb 2004-04-24 devnull ra -= 180*c;
692 8a3b2ceb 2004-04-24 devnull if(ra < 0)
693 8a3b2ceb 2004-04-24 devnull ra += 360*c;
695 8a3b2ceb 2004-04-24 devnull p = pts[0];
696 8a3b2ceb 2004-04-24 devnull p.x -= stringwidth(font, hm5(angle(ra)))/2;
697 8a3b2ceb 2004-04-24 devnull string(scr, p, GREY, ZP, font, hm5(angle(ra)));
698 8a3b2ceb 2004-04-24 devnull p = pts[nelem(pts)-1];
699 8a3b2ceb 2004-04-24 devnull p.x -= stringwidth(font, hm5(angle(ra)))/2;
700 8a3b2ceb 2004-04-24 devnull p.y -= font->height;
701 8a3b2ceb 2004-04-24 devnull string(scr, p, GREY, ZP, font, hm5(angle(ra)));
702 8a3b2ceb 2004-04-24 devnull if(x == ramax)
704 8a3b2ceb 2004-04-24 devnull x += gridra(mapdec);
705 8a3b2ceb 2004-04-24 devnull if(x > ramax)
706 8a3b2ceb 2004-04-24 devnull x = ramax;
708 8a3b2ceb 2004-04-24 devnull for(y=decmin; y<=decmax; y+=c){
709 8a3b2ceb 2004-04-24 devnull for(j=0; j<nelem(pts); j++){
710 8a3b2ceb 2004-04-24 devnull /* use double to avoid overflow */
711 8a3b2ceb 2004-04-24 devnull v = (double)j / (double)(nelem(pts)-1);
712 8a3b2ceb 2004-04-24 devnull v = ramin + v*(ramax-ramin);
713 8a3b2ceb 2004-04-24 devnull pts[j] = map(v, y);
715 8a3b2ceb 2004-04-24 devnull bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
716 8a3b2ceb 2004-04-24 devnull p = pts[0];
717 8a3b2ceb 2004-04-24 devnull p.x += 3;
718 8a3b2ceb 2004-04-24 devnull p.y -= font->height/2;
719 8a3b2ceb 2004-04-24 devnull string(scr, p, GREY, ZP, font, deg(angle(y)));
720 8a3b2ceb 2004-04-24 devnull p = pts[nelem(pts)-1];
721 8a3b2ceb 2004-04-24 devnull p.x -= 3+stringwidth(font, deg(angle(y)));
722 8a3b2ceb 2004-04-24 devnull p.y -= font->height/2;
723 8a3b2ceb 2004-04-24 devnull string(scr, p, GREY, ZP, font, deg(angle(y)));
726 8a3b2ceb 2004-04-24 devnull /* reorder to get planets in front of stars */
727 8a3b2ceb 2004-04-24 devnull tolast(nil);
728 8a3b2ceb 2004-04-24 devnull tolast("moon"); /* moon is in front of everything... */
729 8a3b2ceb 2004-04-24 devnull tolast("shadow"); /* ... except the shadow */
731 8a3b2ceb 2004-04-24 devnull for(i=0,r=rec; i<nrec; i++,r++){
732 e9ad8994 2005-05-01 devnull dec = r->u.ngc.dec;
733 e9ad8994 2005-05-01 devnull ra = r->u.ngc.ra;
734 8a3b2ceb 2004-04-24 devnull if(folded){
735 8a3b2ceb 2004-04-24 devnull ra -= 180*c;
736 8a3b2ceb 2004-04-24 devnull if(ra < 0)
737 8a3b2ceb 2004-04-24 devnull ra += 360*c;
739 8a3b2ceb 2004-04-24 devnull if(textlevel){
740 8a3b2ceb 2004-04-24 devnull name = nameof(r);
741 8a3b2ceb 2004-04-24 devnull if(name==nil && textlevel>1 && r->type==SAO){
742 8a3b2ceb 2004-04-24 devnull snprint(buf, sizeof buf, "SAO%ld", r->index);
743 8a3b2ceb 2004-04-24 devnull name = buf;
745 8a3b2ceb 2004-04-24 devnull if(name)
746 8a3b2ceb 2004-04-24 devnull drawname(scr, nogrey? display->white : alphagrey, name, ra, dec);
748 8a3b2ceb 2004-04-24 devnull if(r->type == Planet){
749 e9ad8994 2005-05-01 devnull drawplanet(scr, &r->u.planet, map(ra, dec));
750 8a3b2ceb 2004-04-24 devnull continue;
752 8a3b2ceb 2004-04-24 devnull if(r->type == SAO){
753 e9ad8994 2005-05-01 devnull m = r->u.sao.mag;
754 8a3b2ceb 2004-04-24 devnull if(m == UNKNOWNMAG)
755 e9ad8994 2005-05-01 devnull m = r->u.sao.mpg;
756 8a3b2ceb 2004-04-24 devnull if(m == UNKNOWNMAG)
757 8a3b2ceb 2004-04-24 devnull continue;
758 8a3b2ceb 2004-04-24 devnull m = dsize(m);
759 8a3b2ceb 2004-04-24 devnull if(m < 3)
760 8a3b2ceb 2004-04-24 devnull fillellipse(scr, map(ra, dec), m, m, nogrey? display->white : lightgrey, ZP);
762 8a3b2ceb 2004-04-24 devnull ellipse(scr, map(ra, dec), m+1, m+1, 0, display->black, ZP);
763 8a3b2ceb 2004-04-24 devnull fillellipse(scr, map(ra, dec), m, m, display->white, ZP);
765 8a3b2ceb 2004-04-24 devnull continue;
767 8a3b2ceb 2004-04-24 devnull if(r->type == Abell){
768 8a3b2ceb 2004-04-24 devnull ellipse(scr, addpt(map(ra, dec), Pt(-3, 2)), 2, 1, 0, lightblue, ZP);
769 8a3b2ceb 2004-04-24 devnull ellipse(scr, addpt(map(ra, dec), Pt(3, 2)), 2, 1, 0, lightblue, ZP);
770 8a3b2ceb 2004-04-24 devnull ellipse(scr, addpt(map(ra, dec), Pt(0, -2)), 1, 2, 0, lightblue, ZP);
771 8a3b2ceb 2004-04-24 devnull continue;
773 e9ad8994 2005-05-01 devnull switch(r->u.ngc.type){
774 8a3b2ceb 2004-04-24 devnull case Galaxy:
775 e9ad8994 2005-05-01 devnull j = npixels(r->u.ngc.diam);
776 8a3b2ceb 2004-04-24 devnull if(j < 4)
778 8a3b2ceb 2004-04-24 devnull if(j > 10)
779 8a3b2ceb 2004-04-24 devnull k = j/3;
781 8a3b2ceb 2004-04-24 devnull k = j/2;
782 8a3b2ceb 2004-04-24 devnull ellipse(scr, map(ra, dec), j, k, 0, lightblue, ZP);
785 8a3b2ceb 2004-04-24 devnull case PlanetaryN:
786 8a3b2ceb 2004-04-24 devnull p = map(ra, dec);
787 e9ad8994 2005-05-01 devnull j = npixels(r->u.ngc.diam);
788 8a3b2ceb 2004-04-24 devnull if(j < 3)
790 8a3b2ceb 2004-04-24 devnull ellipse(scr, p, j, j, 0, green, ZP);
791 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x, p.y+j+1), Pt(p.x, p.y+j+4),
792 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, green, ZP);
793 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x, p.y-(j+1)), Pt(p.x, p.y-(j+4)),
794 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, green, ZP);
795 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x+j+1, p.y), Pt(p.x+j+4, p.y),
796 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, green, ZP);
797 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x-(j+1), p.y), Pt(p.x-(j+4), p.y),
798 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, green, ZP);
801 8a3b2ceb 2004-04-24 devnull case DiffuseN:
802 8a3b2ceb 2004-04-24 devnull case NebularCl:
803 8a3b2ceb 2004-04-24 devnull p = map(ra, dec);
804 e9ad8994 2005-05-01 devnull j = npixels(r->u.ngc.diam);
805 8a3b2ceb 2004-04-24 devnull if(j < 4)
807 8a3b2ceb 2004-04-24 devnull r1.min = Pt(p.x-j, p.y-j);
808 8a3b2ceb 2004-04-24 devnull r1.max = Pt(p.x+j, p.y+j);
809 e9ad8994 2005-05-01 devnull if(r->u.ngc.type != DiffuseN)
810 8a3b2ceb 2004-04-24 devnull draw(scr, r1, ocstipple, ocstipple, ZP);
811 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x-j, p.y-j), Pt(p.x+j, p.y-j),
812 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, green, ZP);
813 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x-j, p.y+j), Pt(p.x+j, p.y+j),
814 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, green, ZP);
815 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x-j, p.y-j), Pt(p.x-j, p.y+j),
816 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, green, ZP);
817 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x+j, p.y-j), Pt(p.x+j, p.y+j),
818 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, green, ZP);
821 8a3b2ceb 2004-04-24 devnull case OpenCl:
822 8a3b2ceb 2004-04-24 devnull p = map(ra, dec);
823 e9ad8994 2005-05-01 devnull j = npixels(r->u.ngc.diam);
824 8a3b2ceb 2004-04-24 devnull if(j < 4)
826 8a3b2ceb 2004-04-24 devnull fillellipse(scr, p, j, j, ocstipple, ZP);
829 8a3b2ceb 2004-04-24 devnull case GlobularCl:
830 e9ad8994 2005-05-01 devnull j = npixels(r->u.ngc.diam);
831 8a3b2ceb 2004-04-24 devnull if(j < 4)
833 8a3b2ceb 2004-04-24 devnull p = map(ra, dec);
834 8a3b2ceb 2004-04-24 devnull ellipse(scr, p, j, j, 0, lightgrey, ZP);
835 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x-(j-1), p.y), Pt(p.x+j, p.y),
836 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, lightgrey, ZP);
837 8a3b2ceb 2004-04-24 devnull line(scr, Pt(p.x, p.y-(j-1)), Pt(p.x, p.y+j),
838 8a3b2ceb 2004-04-24 devnull Endsquare, Endsquare, 0, lightgrey, ZP);
843 8a3b2ceb 2004-04-24 devnull flushimage(display, 1);
844 8a3b2ceb 2004-04-24 devnull displayimage(scr);
848 8a3b2ceb 2004-04-24 devnull runcommand(char *command, int p[2])
850 8a3b2ceb 2004-04-24 devnull switch(rfork(RFPROC|RFFDG|RFNOWAIT)){
851 8a3b2ceb 2004-04-24 devnull case -1:
852 8a3b2ceb 2004-04-24 devnull return -1;
853 8a3b2ceb 2004-04-24 devnull default:
856 8a3b2ceb 2004-04-24 devnull dup(p[1], 1);
857 8a3b2ceb 2004-04-24 devnull close(p[0]);
858 8a3b2ceb 2004-04-24 devnull close(p[1]);
859 8a3b2ceb 2004-04-24 devnull execlp("rc", "rc", "-c", command, nil);
860 8a3b2ceb 2004-04-24 devnull fprint(2, "can't exec %s: %r", command);
861 8a3b2ceb 2004-04-24 devnull exits(nil);
863 8a3b2ceb 2004-04-24 devnull return 1;
867 8a3b2ceb 2004-04-24 devnull parseplanet(char *line, Planetrec *p)
869 8a3b2ceb 2004-04-24 devnull char *fld[6];
870 8a3b2ceb 2004-04-24 devnull int i, nfld;
871 8a3b2ceb 2004-04-24 devnull char *s;
873 8a3b2ceb 2004-04-24 devnull if(line[0] == '\0')
875 8a3b2ceb 2004-04-24 devnull line[10] = '\0'; /* terminate name */
876 8a3b2ceb 2004-04-24 devnull s = strrchr(line, ' ');
877 8a3b2ceb 2004-04-24 devnull if(s == nil)
878 8a3b2ceb 2004-04-24 devnull s = line;
881 8a3b2ceb 2004-04-24 devnull strcpy(p->name, s);
882 8a3b2ceb 2004-04-24 devnull for(i=0; s[i]!='\0'; i++)
883 8a3b2ceb 2004-04-24 devnull p->name[i] = tolower(s[i]);
884 8a3b2ceb 2004-04-24 devnull p->name[i] = '\0';
885 8a3b2ceb 2004-04-24 devnull nfld = getfields(line+11, fld, nelem(fld), 1, " ");
886 8a3b2ceb 2004-04-24 devnull p->ra = dangle(getra(fld[0]));
887 8a3b2ceb 2004-04-24 devnull p->dec = dangle(getra(fld[1]));
888 8a3b2ceb 2004-04-24 devnull p->az = atof(fld[2])*MILLIARCSEC;
889 8a3b2ceb 2004-04-24 devnull p->alt = atof(fld[3])*MILLIARCSEC;
890 8a3b2ceb 2004-04-24 devnull p->semidiam = atof(fld[4])*1000;
891 8a3b2ceb 2004-04-24 devnull if(nfld > 5)
892 8a3b2ceb 2004-04-24 devnull p->phase = atof(fld[5]);
894 8a3b2ceb 2004-04-24 devnull p->phase = 0;
898 8a3b2ceb 2004-04-24 devnull astro(char *flags, int initial)
900 8a3b2ceb 2004-04-24 devnull int p[2];
901 8a3b2ceb 2004-04-24 devnull int i, n, np;
902 8a3b2ceb 2004-04-24 devnull char cmd[256], buf[4096], *lines[20], *fld[10];
904 8a3b2ceb 2004-04-24 devnull snprint(cmd, sizeof cmd, "astro -p %s", flags);
905 8a3b2ceb 2004-04-24 devnull if(pipe(p) < 0){
906 8a3b2ceb 2004-04-24 devnull fprint(2, "can't pipe: %r\n");
909 8a3b2ceb 2004-04-24 devnull if(runcommand(cmd, p) < 0){
910 8a3b2ceb 2004-04-24 devnull close(p[0]);
911 8a3b2ceb 2004-04-24 devnull close(p[1]);
912 8a3b2ceb 2004-04-24 devnull fprint(2, "can't run astro: %r");
915 8a3b2ceb 2004-04-24 devnull close(p[1]);
916 8a3b2ceb 2004-04-24 devnull n = readn(p[0], buf, sizeof buf-1);
917 8a3b2ceb 2004-04-24 devnull if(n <= 0){
918 8a3b2ceb 2004-04-24 devnull fprint(2, "no data from astro\n");
921 8a3b2ceb 2004-04-24 devnull if(!initial)
922 8a3b2ceb 2004-04-24 devnull Bwrite(&bout, buf, n);
923 8a3b2ceb 2004-04-24 devnull buf[n] = '\0';
924 8a3b2ceb 2004-04-24 devnull np = getfields(buf, lines, nelem(lines), 0, "\n");
925 8a3b2ceb 2004-04-24 devnull if(np <= 1){
926 8a3b2ceb 2004-04-24 devnull fprint(2, "astro: not enough output\n");
929 8a3b2ceb 2004-04-24 devnull Bprint(&bout, "%s\n", lines[0]);
930 8a3b2ceb 2004-04-24 devnull Bflush(&bout);
931 8a3b2ceb 2004-04-24 devnull /* get latitude and longitude */
932 8a3b2ceb 2004-04-24 devnull if(getfields(lines[0], fld, nelem(fld), 1, " ") < 8)
933 8a3b2ceb 2004-04-24 devnull fprint(2, "astro: can't read longitude: too few fields\n");
935 8a3b2ceb 2004-04-24 devnull mysid = getra(fld[5])*180./PI;
936 8a3b2ceb 2004-04-24 devnull mylat = getra(fld[6])*180./PI;
937 8a3b2ceb 2004-04-24 devnull mylon = getra(fld[7])*180./PI;
940 8a3b2ceb 2004-04-24 devnull * Each time we run astro, we generate a new planet list
941 8a3b2ceb 2004-04-24 devnull * so multiple appearances of a planet may exist as we plot
942 8a3b2ceb 2004-04-24 devnull * its motion over time.
944 8a3b2ceb 2004-04-24 devnull planet = malloc(NPlanet*sizeof planet[0]);
945 8a3b2ceb 2004-04-24 devnull if(planet == nil){
946 8a3b2ceb 2004-04-24 devnull fprint(2, "astro: malloc failed: %r\n");
947 8a3b2ceb 2004-04-24 devnull exits("malloc");
949 8a3b2ceb 2004-04-24 devnull memset(planet, 0, NPlanet*sizeof planet[0]);
950 8a3b2ceb 2004-04-24 devnull for(i=1; i<np; i++)
951 8a3b2ceb 2004-04-24 devnull parseplanet(lines[i], &planet[i-1]);
