#include <u.h>#include <libc.h>#include <bio.h>#include <draw.h>#include <event.h>#include <ctype.h>#include "map.h"#undef	RAD#undef	TWOPI#include "sky.h" 	/* convert to milliarcsec */static long	c = MILLIARCSEC;	/* 1 degree */static long	m5 = 1250*60*60;	/* 5 minutes of ra */ DAngle	ramin;DAngle	ramax;DAngle	decmin;DAngle	decmax;int		folded; Image	*grey;Image	*alphagrey;Image	*green;Image	*lightblue;Image	*lightgrey;Image	*ocstipple;Image	*suncolor;Image	*mooncolor;Image	*shadowcolor;Image	*mercurycolor;Image	*venuscolor;Image	*marscolor;Image	*jupitercolor;Image	*saturncolor;Image	*uranuscolor;Image	*neptunecolor;Image	*plutocolor;Image	*cometcolor; Planetrec	*planet; long	mapx0, mapy0;long	mapra, mapdec;double	mylat, mylon, mysid;double	mapscale;double	maps;int (*projection)(struct place*, double*, double*); char *fontname = "/lib/font/bit/lucida/unicode.6.font"; /* types Coord and Loc correspond to types in map(3) thus:   Coord == struct coord;   Loc == struct place, except longitudes are measured     positive east for Loc and west for struct place*/ typedef struct Xyz Xyz;typedef struct coord Coord;typedef struct Loc Loc; struct Xyz{	double x,y,z;}; struct Loc{	Coord nlat, elon;}; Xyz north = { 0, 0, 1 }; intplotopen(void){	if(display != nil)		return 1;	display = initdisplay(nil, nil, drawerror);	if(display == nil){		fprint(2, "initdisplay failed: %r\n");		return -1;	}	grey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x777777FF);	lightgrey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0xAAAAAAFF);	alphagrey = allocimage(display, Rect(0, 0, 2, 2), RGBA32, 1, 0x77777777);	green = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x00AA00FF);	lightblue = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x009EEEFF);	ocstipple = allocimage(display, Rect(0, 0, 2, 2), CMAP8, 1, 0xAAAAAAFF);	draw(ocstipple, Rect(0, 0, 1, 1), display->black, nil, ZP);	draw(ocstipple, Rect(1, 1, 2, 2), display->black, nil, ZP); 	suncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFF77FF);	mooncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAAAFF);	shadowcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x00000055);	mercurycolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFAAAAFF);	venuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);	marscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFF5555FF);	jupitercolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFFAAFF);	saturncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAFFAAFF);	uranuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77DDDDFF);	neptunecolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77FF77FF);	plutocolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x7777FFFF);	cometcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);	font = openfont(display, fontname);	if(font == nil)		fprint(2, "warning: no font %s: %r\n", fontname);	return 1;} /* * Function heavens() for setting up observer-eye-view * sky charts, plus subsidiary functions. * Written by Doug McIlroy. */ /* heavens(nlato, wlono, nlatc, wlonc) sets up a map(3)-style   coordinate change (by calling orient()) and returns a   projection function heavens for calculating a star map   centered on nlatc,wlonc and rotated so it will appear   upright as seen by a ground observer at nlato,wlono.   all coordinates (north latitude and west longitude)   are in degrees.*/ Coordmkcoord(double degrees){	Coord c; 	c.l = degrees*PI/180;	c.s = sin(c.l);	c.c = cos(c.l);	return c;} Xyzptov(struct place p){	Xyz v; 	v.z = p.nlat.s;	v.x = p.nlat.c * p.wlon.c;	v.y = -p.nlat.c * p.wlon.s;	return v;} doubledot(Xyz a, Xyz b){	return a.x*b.x + a.y*b.y + a.z*b.z;} Xyzcross(Xyz a, Xyz b){	Xyz v; 	v.x = a.y*b.z - a.z*b.y;	v.y = a.z*b.x - a.x*b.z;	v.z = a.x*b.y - a.y*b.x;	return v;} doublelen(Xyz a){	return sqrt(dot(a, a));} /* An azimuth vector from a to b is a tangent to   the sphere at a pointing along the (shorter)   great-circle direction to b.  It lies in the   plane of the vectors a and b, is perpendicular   to a, and has a positive compoent along b,   provided a and b span a 2-space.  Otherwise   it is null.  (aXb)Xa, where X denotes cross   product, is such a vector. */ Xyzazvec(Xyz a, Xyz b){	return cross(cross(a,b), a);} /* Find the azimuth of point b as seen   from point a, given that a is distinct   from either pole and from b */ doubleazimuth(Xyz a, Xyz b){	Xyz aton = azvec(a,north);	Xyz atob = azvec(a,b);	double lenaton = len(aton);	double lenatob = len(atob);	double az = acos(dot(aton,atob)/(lenaton*lenatob)); 	if(dot(b,cross(north,a)) < 0)		az = -az;	return az;} /* Find the rotation (3rd argument of orient() in the   map projection package) for the map described   below.  The return is radians; it must be converted   to degrees for orient().*/ doublerot(struct place center, struct place zenith){	Xyz cen = ptov(center);	Xyz zen = ptov(zenith); 	if(cen.z > 1-FUZZ) 	/* center at N pole */		return PI + zenith.wlon.l;	if(cen.z < FUZZ-1)		/* at S pole */		return -zenith.wlon.l;	if(fabs(dot(cen,zen)) > 1-FUZZ)	/* at zenith */		return 0;	return azimuth(cen, zen);} int (*heavens(double zlatdeg, double zlondeg, double clatdeg, double clondeg))(Loc*, double*, double*){	struct place center;	struct place zenith; 	center.nlat = mkcoord(clatdeg);	center.wlon = mkcoord(clondeg);	zenith.nlat = mkcoord(zlatdeg);	zenith.wlon = mkcoord(zlondeg);	projection = stereographic();	orient(clatdeg, clondeg, rot(center, zenith)*180/PI);	return projection;} intmaptoxy(long ra, long dec, double *x, double *y){	double lat, lon;	struct place pl; 	lat = angle(dec);	lon = angle(ra);	pl.nlat.l = lat;	pl.nlat.s = sin(lat);	pl.nlat.c = cos(lat);	pl.wlon.l = lon;	pl.wlon.s = sin(lon);	pl.wlon.c = cos(lon);	normalize(&pl);	return projection(&pl, x, y);} /* end of 'heavens' section */ intsetmap(long ramin, long ramax, long decmin, long decmax, Rectangle r, int zenithup){	int c;	double minx, maxx, miny, maxy; 	c = 1000*60*60;	mapra = ramax/2+ramin/2;	mapdec = decmax/2+decmin/2; 	if(zenithup)		projection = heavens(mylat, mysid, (double)mapdec/c, (double)mapra/c);	else{		orient((double)mapdec/c, (double)mapra/c, 0.);		projection = stereographic();	}	mapx0 = (r.max.x+r.min.x)/2;	mapy0 = (r.max.y+r.min.y)/2;	maps = ((double)Dy(r))/(double)(decmax-decmin);	if(maptoxy(ramin, decmin, &minx, &miny) < 0)		return -1;	if(maptoxy(ramax, decmax, &maxx, &maxy) < 0)		return -1;	/*	 * It's tricky to get the scale right.  This noble attempt scales	 * to fit the lengths of the sides of the rectangle.	 */	mapscale = Dx(r)/(minx-maxx);	if(mapscale > Dy(r)/(maxy-miny))		mapscale = Dy(r)/(maxy-miny);	/*	 * near the pole It's not a rectangle, though, so this scales	 * to fit the corners of the trapezoid (a triangle at the pole).	 */	mapscale *= (cos(angle(mapdec))+1.)/2;	if(maxy  < miny){		/* upside down, between zenith and pole */		fprint(2, "reverse plot\n");		mapscale = -mapscale;	}	return 1;} Pointmap(long ra, long dec){	double x, y;	Point p; 	if(maptoxy(ra, dec, &x, &y) > 0){		p.x = mapscale*x + mapx0;		p.y = mapscale*-y + mapy0;	}else{		p.x = -100;		p.y = -100;	}	return p;} intdsize(int mag)	/* mag is 10*magnitude; return disc size */{	double d; 	mag += 25;	/* make mags all positive; sirius is -1.6m */	d = (130-mag)/10;	/* if plate scale is huge, adjust */	if(maps < 100.0/MILLIARCSEC)		d *= .71;	if(maps < 50.0/MILLIARCSEC)		d *= .71;	return d;} voiddrawname(Image *scr, Image *col, char *s, int ra, int dec){	Point p; 	if(font == nil)		return;	p = addpt(map(ra, dec), Pt(4, -1));	/* font has huge ascent */	string(scr, p, col, ZP, font, s);} intnpixels(DAngle diam){	Point p0, p1; 	diam = DEG(angle(diam)*MILLIARCSEC);	/* diam in milliarcsec */	diam /= 2;	/* radius */	/* convert to plate scale */	/* BUG: need +100 because we sometimes crash in library if we plot the exact center */	p0 = map(mapra+100, mapdec);	p1 = map(mapra+100+diam, mapdec);	return hypot(p0.x-p1.x, p0.y-p1.y);} voiddrawdisc(Image *scr, DAngle semidiam, int ring, Image *color, Point pt, char *s){	int d; 	d = npixels(semidiam*2);	if(d < 5)		d = 5;	fillellipse(scr, pt, d, d, color, ZP);	if(ring == 1)		ellipse(scr, pt, 5*d/2, d/2, 0, color, ZP);	if(ring == 2)		ellipse(scr, pt, d, d, 0, grey, ZP);	if(s){		d = stringwidth(font, s);		pt.x -= d/2;		pt.y -= font->height/2 + 1;		string(scr, pt, display->black, ZP, font, s);	}} voiddrawplanet(Image *scr, Planetrec *p, Point pt){	if(strcmp(p->name, "sun") == 0){		drawdisc(scr, p->semidiam, 0, suncolor, pt, nil);		return;	}	if(strcmp(p->name, "moon") == 0){		drawdisc(scr, p->semidiam, 0, mooncolor, pt, nil);		return;	}	if(strcmp(p->name, "shadow") == 0){		drawdisc(scr, p->semidiam, 2, shadowcolor, pt, nil);		return;	}	if(strcmp(p->name, "mercury") == 0){		drawdisc(scr, p->semidiam, 0, mercurycolor, pt, "m");		return;	}	if(strcmp(p->name, "venus") == 0){		drawdisc(scr, p->semidiam, 0, venuscolor, pt, "v");		return;	}	if(strcmp(p->name, "mars") == 0){		drawdisc(scr, p->semidiam, 0, marscolor, pt, "M");		return;	}	if(strcmp(p->name, "jupiter") == 0){		drawdisc(scr, p->semidiam, 0, jupitercolor, pt, "J");		return;	}	if(strcmp(p->name, "saturn") == 0){		drawdisc(scr, p->semidiam, 1, saturncolor, pt, "S");				return;	}	if(strcmp(p->name, "uranus") == 0){		drawdisc(scr, p->semidiam, 0, uranuscolor, pt, "U");				return;	}	if(strcmp(p->name, "neptune") == 0){		drawdisc(scr, p->semidiam, 0, neptunecolor, pt, "N");				return;	}	if(strcmp(p->name, "pluto") == 0){		drawdisc(scr, p->semidiam, 0, plutocolor, pt, "P");				return;	}	if(strcmp(p->name, "comet") == 0){		drawdisc(scr, p->semidiam, 0, cometcolor, pt, "C");		return;	} 	pt.x -= stringwidth(font, p->name)/2;	pt.y -= font->height/2;	string(scr, pt, grey, ZP, font, p->name);} voidtolast(char *name){	int i, nlast;	Record *r, rr; 	/* stop early to simplify inner loop adjustment */	nlast = 0;	for(i=0,r=rec; i<nrec-nlast; i++,r++)		if(r->type == Planet)			if(name==nil || strcmp(r->planet.name, name)==0){				rr = *r;				memmove(rec+i, rec+i+1, (nrec-i-1)*sizeof(Record));				rec[nrec-1] = rr;				--i;				--r;				nlast++;			}} intbbox(long extrara, long extradec, int quantize){	int i;	Record *r;	int ra, dec;	int rah, ram, d1, d2;	double r0; 	ramin = 0x7FFFFFFF;	ramax = -0x7FFFFFFF;	decmin = 0x7FFFFFFF;	decmax = -0x7FFFFFFF; 	for(i=0,r=rec; i<nrec; i++,r++){		if(r->type == Patch){			radec(r->index, &rah, &ram, &dec);			ra = 15*rah+ram/4;			r0 = c/cos(dec*PI/180);			ra *= c;			dec *= c;			if(dec == 0)				d1 = c, d2 = c;			else if(dec < 0)				d1 = c, d2 = 0;			else				d1 = 0, d2 = c;		}else if(r->type==SAO || r->type==NGC || r->type==Planet || r->type==Abell){			ra = r->ngc.ra;			dec = r->ngc.dec;			d1 = 0, d2 = 0, r0 = 0;		}else			continue;		if(dec+d2+extradec > decmax)			decmax = dec+d2+extradec;		if(dec-d1-extradec < decmin)			decmin = dec-d1-extradec;		if(folded){			ra -= 180*c;			if(ra < 0)				ra += 360*c;		}

Lines 1–500 of 953. Load more