plot.c

这是一个同样来自贝尔实验室的和UNIX有着渊源的操作系统, 其简洁的设计和实现易于我们学习和理解

#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 };

int
plotopen(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.
*/

Coord
mkcoord(double degrees)
{
	Coord c;

	c.l = degrees*PI/180;
	c.s = sin(c.l);
	c.c = cos(c.l);
	return c;
}

Xyz
ptov(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;
}

double
dot(Xyz a, Xyz b)
{
	return a.x*b.x + a.y*b.y + a.z*b.z;
}

Xyz
cross(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;
}

double
len(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. */

Xyz
azvec(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 */

double
azimuth(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().
*/

double
rot(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;
}

int
maptoxy(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 */

int
setmap(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;
}

Point
map(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;
}

int
dsize(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;
}

void
drawname(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);
}

int
npixels(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);
}

void
drawdisc(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);
	}
}

void
drawplanet(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);
}

void
tolast(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++;
			}
}

int
bbox(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);