commit 8a3b2ceb0ff632c47e1516d3ffef8572dc8eb974 from: rsc date: Sat Apr 24 17:05:43 2004 UTC Add scat. Temporary fix to rc r.e. note groups. commit - 3f8c70e97c2eb85992424439af56a4dd6412b8c6 commit + 8a3b2ceb0ff632c47e1516d3ffef8572dc8eb974 blob - /dev/null blob + d9bf8db805482943fa01379ad0482fe0e13cee5a (mode 644) --- /dev/null +++ man/man1/scat.1 @@ -0,0 +1,335 @@ +.TH SCAT 7 +.SH NAME +scat \- sky catalogue and Digitized Sky Survey +.SH SYNOPSIS +.B scat +.SH DESCRIPTION +.I Scat +looks up items in catalogues of objects +outside the solar system +and implements database-like manipulations +on sets of such objects. +It also provides an interface to +.IR astro (7) +to plot the locations of solar system objects. +Finally, it displays images from the +Space Telescope Science Institute's +Digitized Sky Survey, keyed to the catalogues. +.PP +Items are read, one per line, from the standard input +and looked up in the catalogs. +Input is case-insensitive. +The result of the lookup becomes the set of objects available +to the database commands. +After each lookup or command, if more than two objects are +in the set, +.I scat +prints how many objects are in the set; otherwise it +prints the objects' +descriptions or cross-index listings (suitable for input to +.IR scat ). +An item is in one of the following formats: +.TP +.B ngc1234 +Number 1234 in the New General Catalogue of +Nonstellar Objects, NGC2000.0. +The output identifies the type +.RB( Gx =galaxy, +.BR Pl =planetary +nebula, +.BR OC =open +cluster, +.BR Gb =globular +cluster, +.BR Nb =bright +nebula, +.BR C+N =cluster +associated with nebulosity, +.BR Ast =asterism, +.BR Kt =knot +or nebulous region in a galaxy, +.BR *** =triple +star, +.BR D* =double +star, +.BR ? =uncertain, +.BR - =nonexistent, +.BR PD =plate +defect, and +(blank)=unverified or unknown), +its position in 2000.0 coordinates, +its size in minutes of arc, a brief description, and popular names. +.TP +.B ic1234 +Like NGC references, but from the Index Catalog. +.TP +.B sao12345 +Number 12345 in the Smithsonian Astrophysical Star Catalogue. +Output identifies the visual and photographic magnitudes, +2000.0 coordinates, proper motion, spectral type, multiplicity and variability +class, and HD number. +.TP +.B m4 +Catalog number 4 in Messier's catalog. +The output is the NGC number. +.TP +.B abell1701 +Catalog number 1701 in the Abell and Zwicky +catalog of clusters of galaxies. +Output identifies the magnitude of the tenth brightest member of the cluster, +radius of the cluster in degrees, its distance in megaparsecs, +2000.0 coordinates, galactic latitude and longitude, +magnitude range of the cluster (the `distance group'), +number of members (the `richness group'), population +per square degree, and popular names. +.TP +.B planetarynebula +The set of NGC objects of the specified type. +The type may be a compact NGC code or a full name, as above, with no blank. +.TP +\fL"α umi"\fP +Names are provided in double quotes. +Known names are the Greek +letter designations, proper names such as Betelgeuse, bright variable stars, +and some proper names of stars, NGC objects, and Abell clusters. +Greek letters may be spelled out, e.g. +.BR alpha . +Constellation names must be the three-letter abbreviations. +The output +is the SAO number. +For non-Greek names, catalog numbers and names are listed for all objects with +names for which the given name is a prefix. +.TP +.B 12h34m -16 +Coordinates in the sky are translated to the nearest `patch', +approximately one square degree of sky. +The output is the coordinates identifying the patch, +the constellations touching the patch, and the Abell, NGC, and SAO +objects in the patch. +The program prints sky positions in several formats corresponding +to different precisions; any output format is understood as input. +.TP +.B umi +All the patches in the named constellation. +.TP +.B mars +The planets are identified by their names. +The names +.B shadow +and +.B comet +refer to the earth's penumbra at lunar distance and the comet installed in the current +.IR astro (7). +The output is the planet's name, right ascension and declination, azimuth and altitude, and phase +for the moon and sun, as shown by +.BR astro . +The positions are current at the start of +.I scat 's +execution; see the +.B astro +command in the next section for more information. +.PP +The commands are: +.TF print +.TP +.BI add " item" +Add the named item to the set. +.TP +.BI keep " class ..." +Flatten the set and cull it, keeping only the specified classes. +The classes may be specific NGC types, +all stars +.RB ( sao ), +all NGC objects +.RB ( ngc ), +all M objects +.RB ( m ), +all Abell clusters +.RB ( abell ), +or a specified brightness range. +Brightness ranges are specified by a leading +.B > +or +.B < +followed by a magnitude. +Remember that brighter objects have lesser magnitudes. +.TP +.BI drop " class ..." +Complement to +.BR keep . +.TP +.BI flat +Some items such as patches represents sets of items. +.I Flat +flattens the set so +.I scat +holds all the information available for the objects in the set. +.TP +.BI print +Print the contents of the set. If the information seems meager, try +flattening the set. +.TP +.BI expand " n" +Flatten the set, +expand the area of the sky covered by the set to be +.I n +degrees wider, and collect all the objects in that area. +If +.I n +is zero, +.I expand +collects all objects in the patches that cover the current set. +.TP +.BI astro " option" +Run +.IR astro (7) +with the specified +.I options +(to which will be appended +.BR -p ), +to discover the positions of the planets. +.BR Astro 's +.B -d +and +.B -l +options can be used to set the time and place; by default, it's right now at the coordinates in +.BR /lib/sky/here . +Running +.B astro +does not change the positions of planets already in the display set, +so +.B astro +may be run multiple times, executing e.g. +.B "add mars" +each time, to plot a series of planetary positions. +.TP +.BI plot " option" +Expand and plot the set in a new window on the screen. +Symbols for NGC objects are as in Sky Atlas 2000.0, except that open clusters +are shown as stippled disks rather than circles. +Abell clusters are plotted as a triangle of ellipses. +The planets are drawn as disks of representative color with the first letter of the name +in the disk (lower case for inferior planets; upper case for superior); +the sun, moon, and earth's shadow are unlabeled disks. +Objects larger than a few pixels are plotted to scale; however, +.I scat +does not have the information necessary to show the correct orientation for galaxies. +.IP +The option +.B nogrid +suppresses the lines of declination and right ascension. +By default, +.I scat +labels NGC objects, Abell clusters, and bright stars; option +.B nolabel +suppresses these while +.B alllabel +labels stars with their SAO number as well. +The default size is 512×512; options +.B dx +.I n +and +.BR dy +.I n +set the +.I x +and +.I y +extent. +The option +.B zenithup +orients the map so it appears as it would in the sky at the time and +location used by the +.B astro +command +.RI ( q.v. ). +.IP +The output is designed to look best on an LCD display. +CRTs have trouble with the thin, grey lines and dim stars. +The option +.B nogrey +uses white instead of grey for these details, improving visibility +at the cost of legibility when plotting on CRTs. +.TP +.B "plate \f1[[\f2ra dec\f1] \f2rasize\f1 [\f2decsize\f1]]" +Display the section of the Digitized Sky Survey (plate scale +approximately 1.7 arcseconds per pixel) centered on the +given right ascension and declination or, if no position is specified, the +current set of objects. The maximum area that will be displayed +is one degree on a side. The horizontal and vertical sizes may +be specified in the usual notation for angles. +If the second size is omitted, a square region is displayed. +If no size is specified, the size is sufficient to display the centers +of all the +objects in the current set. If a single object is in the set, the +500×500 pixel block from the survey containing the center +of the object is displayed. +The survey is stored in the CD-ROM juke box; run +.B 9fs +.B juke +before running +.IR scat . +.TP +.BI gamma " value" +Set the gamma for converting plates to images. Default is \-1.0. +Negative values display white stars, positive black. +The images look best on displays with depth 8 or greater. +.I Scat +does not change the hardware color map, which +should be set externally to a grey scale; try the command +.B getmap gamma +(see +.IR getmap (9.1)) +on an 8-bit color-mapped display. +.PD +.SH EXAMPLES +Plot the Messier objects and naked-eye stars in Orion. +.EX + ori + keep m <6 + plot nogrid +.EE +.PP +Draw a finder chart for Uranus: +.EX + uranus + expand 5 + plot +.EE +.PP +Show a partial lunar eclipse: +.EX + astro -d + 2000 07 16 12 45 + moon + add shadow + expand 2 + plot +.EE +.PP +Draw a map of the Pleiades. +.EX + "alcyone" + expand 1 + plot +.EE +.PP +Show a pretty galaxy. +.EX + ngc1300 + plate 10' +.EE +.SH FILES +.B /lib/sky/*.scat +.SH SOURCE +.B /sys/src/cmd/scat +.SH SEE ALSO +.IR astro (7) +.br +.B /lib/sky/constelnames\ \ +the three-letter abbreviations of the constellation names. +.PP +The data was provided by the Astronomical Data Center at the NASA Goddard +Space Flight Center, except for NGC2000.0, which is Copyright © 1988, Sky +Publishing Corporation, used (but not distributed) by permission. The Digitized Sky Survey, 102 +CD-ROMs, is not distributed with the system. blob - /dev/null blob + 7bcdd971dde016528cf284d3eba5c3578c2d250b (mode 644) --- /dev/null +++ sky/README @@ -0,0 +1 @@ +wget -O- http://pdos.lcs.mit.edu/~rsc/scat.tgz | gunzip | tar xf - blob - 258ca9a64a4314491d8829b8c7489cacdabafab8 blob + 7ebbe368406352a51d9bab661bdd1c4ebff5fabf --- sky/here +++ sky/here @@ -1 +1 @@ -41.0343 73.3936 200 +40.6843333333 74.3996666667 150 blob - 85b20c37f05999acbababc61904b1971e240e2e9 blob + 0f27676d33cf13060cea967946347b4d548e7b41 --- src/cmd/rc/simple.c +++ src/cmd/rc/simple.c @@ -63,7 +63,7 @@ void Xsimple(void){ Xerror("try again"); return; case 0: - rfork(RFNOTEG); + // rfork(RFNOTEG); pushword("exec"); execexec(); strcpy(buf, "can't exec: "); blob - /dev/null blob + b4bd286baa51490af3e0ad91ce0582824a942245 (mode 644) --- /dev/null +++ src/cmd/scat/bitinput.c @@ -0,0 +1,76 @@ +#include +#include +#include +#include "sky.h" + +static int hufvals[] = { + 1, 1, 1, 1, 1, 1, 1, 1, + 2, 2, 2, 2, 2, 2, 2, 2, + 4, 4, 4, 4, 4, 4, 4, 4, + 8, 8, 8, 8, 8, 8, 8, 8, + 3, 3, 3, 3, 5, 5, 5, 5, + 10, 10, 10, 10, 12, 12, 12, 12, + 15, 15, 15, 15, 6, 6, 7, 7, + 9, 9, 11, 11, 13, 13, 0, 14, +}; + +static int huflens[] = { + 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, + 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 6, 6, +}; + +static int buffer; +static int bits_to_go; /* Number of bits still in buffer */ + +void +start_inputing_bits(void) +{ + bits_to_go = 0; +} + +int +input_huffman(Biobuf *infile) +{ + int c; + + if(bits_to_go < 6) { + c = Bgetc(infile); + if(c < 0) { + fprint(2, "input_huffman: unexpected EOF\n"); + exits("format"); + } + buffer = (buffer<<8) | c; + bits_to_go += 8; + } + c = (buffer >> (bits_to_go-6)) & 0x3f; + bits_to_go -= huflens[c]; + return hufvals[c]; +} + +int +input_nybble(Biobuf *infile) +{ + int c; + + if(bits_to_go < 4) { + c = Bgetc(infile); + if(c < 0){ + fprint(2, "input_nybble: unexpected EOF\n"); + exits("format"); + } + buffer = (buffer<<8) | c; + bits_to_go += 8; + } + + /* + * pick off the first 4 bits + */ + bits_to_go -= 4; + return (buffer>>bits_to_go) & 0x0f; +} blob - /dev/null blob + 968229cccfd7c003a7ff6d6e83e394469b931bb3 (mode 644) --- /dev/null +++ src/cmd/scat/desc.c @@ -0,0 +1,327 @@ +char *desctab[][2]={ + "!!!", "(magnificent or otherwise interesting object)", + "!!", "(superb)", + "!", "(remarkable)", + "1st", "first", + "2nd", "second", + "3rd", "third", + "4th", "fourth", + "5th", "fifth", + "6th", "sixth", + "7th", "seventh", + "8th", "eighth", + "9th", "ninth", + "B", "bright", + "Borealis", "Borealis", + "C", "compressed", + "Car", "Car", + "Cas", "Cas", + "Cen", "Cen", + "Cl", "cluster", + "Cyg", "Cyg", + "D", "double", + "Dra", "Dra", + "Dumbbell", "Dumbbell", + "E", "extended", + "F", "faint", + "L", "large", + "Lyra", "Lyra", + "M", "(in the) middle", + "Merope", "Merope", + "Milky", "Milky", + "Mon", "Mon", + "N", "(to a) nucleus", + "Neb", "Nebula", + "Nor", "Nor", + "Nucl", "nucleus", + "Nuclei", "nuclei", + "P", "poor in stars", + "PN", "planetary nebula", + "Polarissima", "Polarissima", + "Praesepe", "Praesepe", + "Psc", "Psc", + "R", "round", + "RA", "RA", + "RR", "exactly round", + "Ri", "rich in stars", + "S", "small", + "Sco", "Sco", + "S*", "small star", + "Way", "Way", + "ab", "about", + "about", "about", + "alm", "almost", + "alpha", "α", + "am", "among", + "annul", "annular", + "att", "attached", + "b", "brighter", + "beautiful", "beautiful", + "bet", "between", + "beta", "β", + "bf", "brightest to f side", + "biN", "binuclear", + "bifid", "bifid", + "bifurcated", "bifurcated", + "black", "black", + "blue", "blue", + "bn", "brightest to n side", + "border", "border", + "bp", "brightest to p side", + "branch", "branch", + "branched", "branched", + "branches", "branches", + "bright", "bright", + "brighter", "brighter", + "brightest", "brightest", + "brightness", "brightness", + "brush", "brush", + "bs", "brightest to s side", + "but", "but", + "by", "by", + "c", "considerably", + "centre", "centre", + "certain", "certain", + "chev", "chevelure", + "chief", "chief", + "chiefly", "chiefly", + "circle", "circle", + "close", "close", + "cloud", "cloud", + "cluster", "cluster", + "clusters", "clusters", + "co", "coarse, coarsely", + "com", "cometic", + "comet", "comet", + "cometary", "cometary", + "comp", "companion", + "condens", "condensations", + "condensed", "condensed", + "conn", "connected", + "connected", "connected", + "connecting", "connecting", + "cont", "in contact", + "corner", "corner", + "curved", "curved", + "d", "diameter", + "decl", "declination", + "def", "defined", + "defect", "defect", + "deg", "°", + "delta", "δ", + "dense", "dense", + "densest", "densest", + "descr", "description", + "description", "description", + "dif", "diffuse", + "different", "different", + "diffic", "difficult", + "difficult", "difficult", + "diffuse", "diffuse", + "diffused", "diffused", + "disc", "disc", + "dist", "distant", + "distant", "distant", + "distinct", "distinct", + "double", "double", + "doubtful", "doubtful", + "dozen", "dozen", + "e", "extremely", + "each", "each", + "edge", "edge", + "ee", "most extremely", + "ellipt", "elliptical", + "elliptic", "elliptical", + "end", "end", + "ends", "ends", + "epsilon", "ε", + "equilateral", "equilateral", + "er", "easily resolvable", + "eta", "η", + "evidently", "evidently", + "exc", "excentric", + "excen", "excentric", + "excent", "excentric", + "excentric", "excentric", + "extends", "extends", + "f", "following", + "faint", "faint", + "fainter", "fainter", + "faintest", "faintest", + "falcate", "falcate", + "fan", "fan", + "farther", "farther", + "field", "field", + "fine", "fine", + "forming", "forming", + "forms", "forms", + "found", "found", + "from", "from", + "g", "gradually", + "gamma", "γ", + "gaseous", "gaseous", + "gl", "gradually a little", + "glob. cl.", "globular cluster", + "gr", "group", + "great", "great", + "greater", "greater", + "group", "group", + "groups", "groups", + "i", "irregular", + "iF", "irregular figure", + "if", "if", + "in", "in", + "indistinct", "indistinct", + "incl", "including", + "inv", "involved", + "iota", "ι", + "irr", "irregular", + "is", "is", + "it", "it", + "kappa", "κ", + "l", "little, long", + "lC", "little compressed", + "lE", "little extended", + "lambda", "λ", + "larger", "larger", + "last", "last", + "lb", "little brighter", + "least", "least", + "like", "like", + "line", "in a line", + "little", "little", + "long", "long", + "looks", "looks", + "looped", "looped", + "m", "magnitude", + "mE", "much extended", + "mag", "mag", + "makes", "makes", + "many", "many", + "mb", "much brighter", + "more", "more", + "mottled", "mottled", + "mu", "μ", + "mult", "multiple", + "n", "north", + "narrow", "narrow", + "near", "near", + "nearly", "nearly", + "neb", "nebula", + "nebs", "nebulous", + "nebula", "nebula", + "nebulosity", "nebulosity", + "nebulous", "nebulous", + "neby", "nebulosity", + "nf", "north following", + "no", "no", + "nonexistent", "nonexistent", + "not", "not", + "np", "north preceding", + "nr", "near", + "ns", "north-south", + "nu", "ν", + "omega", "ω", + "p", "preceding", + "pB", "pretty bright", + "pC", "pretty compressed", + "pF", "pretty faint", + "pL", "pretty large", + "pR", "pretty round", + "pS", "pretty small", + "parallel", "parallel", + "part", "part", + "partly", "partly", + "patch", "patch", + "patches", "patches", + "perhaps", "perhaps", + "perpendicular", "perpendicular", + "pf", "preceding-following", + "pg", "pretty gradually", + "photosphere", "photosphere", + "pi", "π", + "place", "place", + "plate", "plate", + "plan", "planetary nebula", + "pointed", "pointed", + "portion", "portion", + "pos", "position angle", + "possibly", "possibly", + "prob", "probably", + "probably", "probably", + "ps", "pretty suddenly", + "r", "mottled", + "requires", "requires", + "resolved", "resolved", + "rho", "ρ", + "ring", "ring", + "rough", "rough", + "rr", "some stars seen", + "rrr", "clearly consisting of stars", + "ruby", "ruby", + "s", "south", + "same", "same", + "sb", "suddenly brighter", + "sc", "scattered", + "second", "second", + "seems", "seems", + "seen", "seen", + "segment", "segment", + "semi", "semi", + "sev", "several", + "several", "several", + "sf", "south following", + "shape", "shape", + "shaped", "shaped", + "sharp", "sharp", + "sigma", "σ", + "sl", "suddenly a little", + "slightly", "slightly", + "small", "small", + "south", "south", + "sp", "south preceding", + "spectrum", "spectrum", + "spindle", "spindle", + "spir", "spiral", + "spiral", "spiral", + "st 9...", "stars of mag. 9 and fainter", + "st 9...13", "stars of mag. 9 to 13", + "st", "stars", + "stell", "stellar, pointlike", + "stellar", "stellar", + "straight", "straight", + "streak", "streak", + "strongly", "strongly", + "surrounded", "surrounded", + "surrounds", "surrounds", + "susp", "suspected", + "suspected", "suspected", + "tau", "τ", + "theta", "θ", + "trap", "trapezium", + "trapez", "trapezium", + "trapezium", "trapezium", + "triN", "trinuclear", + "v", "very", + "var", "variable", + "variable", "variable", + "verification", "verification", + "verified", "verified", + "very", "very", + "vl", "very little", + "vm", "very much", + "vs", "very suddenly", + "vv", "very very", + "zeta", "ζ", + 0, 0, +}; + +/*& + "ST9", "stars from the 9th magnitude downward", + "ST9...13", "stars from 9th to 13th magnitude", + "()", "items questioned by Dreyer enclosed in parentheses", + *10 star of 10th mag + "*", "a star (or stars)", + "**", "double star", + "***", "triple star", +*/ blob - /dev/null blob + 11642e59902daa61024f5cd1db5ffa58a61d181a (mode 644) --- /dev/null +++ src/cmd/scat/display.c @@ -0,0 +1,88 @@ +#include +#include +#include +#include +#include "sky.h" + +void +displaypic(Picture *pic) +{ + int p[2]; + int i, n; + uchar *a; + + + if(pipe(p) < 0){ + fprint(2, "pipe failed: %r\n"); + return; + } + switch(rfork(RFPROC|RFFDG)){ + case -1: + fprint(2, "fork failed: %r\n"); + return; + + case 0: + close(p[1]); + dup(p[0], 0); + close(p[0]); + // execl("/bin/page", "page", "-w", 0); + execlp("img", "img", 0); + fprint(2, "exec failed: %r\n"); + exits("exec"); + + default: + close(p[0]); + fprint(p[1], "%11s %11d %11d %11d %11d ", + "k8", pic->minx, pic->miny, pic->maxx, pic->maxy); + n = (pic->maxx-pic->minx)*(pic->maxy-pic->miny); + /* release the memory as we hand it off; this could be a big piece of data */ + a = pic->data; + while(n > 0){ + i = 8192 - (((int)a)&8191); + if(i > n) + i = n; + if(write(p[1], a, i)!=i) + fprint(2, "write error: %r\n"); + // if(i == 8192) /* page aligned */ + // segfree(a, i); + n -= i; + a += i; + } + free(pic->data); + free(pic); + close(p[1]); + break; + } +} + +void +displayimage(Image *im) +{ + int p[2]; + + if(pipe(p) < 0){ + fprint(2, "pipe failed: %r\n"); + return; + } + switch(rfork(RFPROC|RFFDG)){ + case -1: + fprint(2, "fork failed: %r\n"); + return; + + case 0: + close(p[1]); + dup(p[0], 0); + close(p[0]); + execlp("img", "img", 0); + // execl("/bin/page", "page", "-w", 0); + fprint(2, "exec failed: %r\n"); + exits("exec"); + + default: + close(p[0]); + writeimage(p[1], im, 0); + freeimage(im); + close(p[1]); + break; + } +} blob - /dev/null blob + 06f0135e0873aa46dbeb8b8a6ae6d1784cab2498 (mode 644) --- /dev/null +++ src/cmd/scat/dssread.c @@ -0,0 +1,113 @@ +#include +#include +#include +#include "sky.h" + +static void dodecode(Biobuf*, Pix*, int, int, uchar*); +static long getlong(uchar*); +int debug; + +Img* +dssread(char *file) +{ + int nx, ny, scale, sumall; + Pix *p, *pend; + uchar buf[21]; + Biobuf *bp; + Img *ip; + + if(debug) + Bprint(&bout, "reading %s\n", file); + bp = Bopen(file, OREAD); + if(bp == 0) + return 0; + if(Bread(bp, buf, sizeof(buf)) != sizeof(buf) || + buf[0] != 0xdd || buf[1] != 0x99){ + werrstr("bad format"); + return 0; + } + nx = getlong(buf+2); + ny = getlong(buf+6); + scale = getlong(buf+10); + sumall = getlong(buf+14); + if(debug) + fprint(2, "%s: nx=%d, ny=%d, scale=%d, sumall=%d, nbitplanes=%d,%d,%d\n", + file, nx, ny, scale, sumall, buf[18], buf[19], buf[20]); + ip = malloc(sizeof(Img) + (nx*ny-1)*sizeof(int)); + if(ip == 0){ + Bterm(bp); + werrstr("no memory"); + return 0; + } + ip->nx = nx; + ip->ny = ny; + dodecode(bp, ip->a, nx, ny, buf+18); + ip->a[0] = sumall; /* sum of all pixels */ + Bterm(bp); + if(scale > 1){ + p = ip->a; + pend = &ip->a[nx*ny]; + while(p < pend) + *p++ *= scale; + } + hinv(ip->a, nx, ny); + return ip; +} + +static +void +dodecode(Biobuf *infile, Pix *a, int nx, int ny, uchar *nbitplanes) +{ + int nel, nx2, ny2, bits, mask; + Pix *aend, px; + + nel = nx*ny; + nx2 = (nx+1)/2; + ny2 = (ny+1)/2; + memset(a, 0, nel*sizeof(*a)); + + /* + * Initialize bit input + */ + start_inputing_bits(); + + /* + * read bit planes for each quadrant + */ + qtree_decode(infile, &a[0], ny, nx2, ny2, nbitplanes[0]); + qtree_decode(infile, &a[ny2], ny, nx2, ny/2, nbitplanes[1]); + qtree_decode(infile, &a[ny*nx2], ny, nx/2, ny2, nbitplanes[1]); + qtree_decode(infile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]); + + /* + * make sure there is an EOF symbol (nybble=0) at end + */ + if(input_nybble(infile) != 0){ + fprint(2, "dodecode: bad bit plane values\n"); + exits("format"); + } + + /* + * get the sign bits + */ + aend = &a[nel]; + mask = 0; + bits = 0;; + for(; a>= 1; + } + } +} + +static +long getlong(uchar *p) +{ + return (p[0]<<24) | (p[1]<<16) | (p[2]<<8) | p[3]; +} blob - /dev/null blob + 519d98f9aa8fa081d4c40739bc2427ddea137f89 (mode 644) --- /dev/null +++ src/cmd/scat/header.c @@ -0,0 +1,247 @@ +#include +#include +#include +#include "sky.h" + +struct +{ + char name[9]; + char offset; +} Hproto[] = +{ + "ppo1", Pppo1, + "ppo2", Pppo2, + "ppo3", Pppo3, + "ppo4", Pppo4, + "ppo5", Pppo5, + "ppo6", Pppo6, + + "amdx1", Pamdx1, + "amdx2", Pamdx2, + "amdx3", Pamdx3, + "amdx4", Pamdx4, + "amdx5", Pamdx5, + "amdx6", Pamdx6, + "amdx7", Pamdx7, + "amdx8", Pamdx8, + "amdx9", Pamdx9, + "amdx10", Pamdx10, + "amdx11", Pamdx11, + "amdx12", Pamdx12, + "amdx13", Pamdx13, + "amdx14", Pamdx14, + "amdx15", Pamdx15, + "amdx16", Pamdx16, + "amdx17", Pamdx17, + "amdx18", Pamdx18, + "amdx19", Pamdx19, + "amdx20", Pamdx20, + + "amdy1", Pamdy1, + "amdy2", Pamdy2, + "amdy3", Pamdy3, + "amdy4", Pamdy4, + "amdy5", Pamdy5, + "amdy6", Pamdy6, + "amdy7", Pamdy7, + "amdy8", Pamdy8, + "amdy9", Pamdy9, + "amdy10", Pamdy10, + "amdy11", Pamdy11, + "amdy12", Pamdy12, + "amdy13", Pamdy13, + "amdy14", Pamdy14, + "amdy15", Pamdy15, + "amdy16", Pamdy16, + "amdy17", Pamdy17, + "amdy18", Pamdy18, + "amdy19", Pamdy19, + "amdy20", Pamdy20, + + "pltscale", Ppltscale, + "xpixelsz", Pxpixelsz, + "ypixelsz", Pypixelsz, + + "pltrah", Ppltrah, + "pltram", Ppltram, + "pltras", Ppltras, + "pltdecd", Ppltdecd, + "pltdecm", Ppltdecm, + "pltdecs", Ppltdecs, + +}; + +Header* +getheader(char *rgn) +{ + char rec[81], name[81], value[81]; + char *p; + Biobuf *bin; + Header hd, *h; + int i, j, decsn, dss; + + dss = 0; + sprint(rec, "/lib/sky/dssheaders/%s.hhh", rgn); + bin = Bopen(unsharp(rec), OREAD); +/* + if(bin == 0) { + dss = 102; + sprint(rec, "/n/juke/dss/dss.102/headers/%s.hhh", rgn); + bin = Bopen(rec, OREAD); + } + if(bin == 0) { + dss = 61; + sprint(rec, "/n/juke/dss/dss.061/headers/%s.hhh", rgn); + bin = Bopen(rec, OREAD); + } +*/ + if(bin == 0) { + fprint(2, "cannot open %s\n", rgn); + exits("file"); + } + if(debug) + Bprint(&bout, "reading %s\n", rec); + if(dss) + Bprint(&bout, "warning: reading %s from jukebox\n", rec); + + memset(&hd, 0, sizeof(hd)); + j = 0; + decsn = 0; + for(;;) { + if(dss) { + if(Bread(bin, rec, 80) != 80) + break; + rec[80] = 0; + } else { + p = Brdline(bin, '\n'); + if(p == 0) + break; + p[Blinelen(bin)-1] = 0; + strcpy(rec, p); + } + + p = strchr(rec, '/'); + if(p) + *p = 0; + p = strchr(rec, '='); + if(p == 0) + continue; + *p++ = 0; + if(getword(name, rec) == 0) + continue; + if(getword(value, p) == 0) + continue; + if(strcmp(name, "pltdecsn") == 0) { + if(strchr(value, '-')) + decsn = 1; + continue; + } + for(i=0; i= nelem(Hproto)) + j = 0; + if(strcmp(name, Hproto[j].name) == 0) { + hd.param[(uchar)Hproto[j].offset] = atof(value); + break; + } + } + } + Bterm(bin); + + hd.param[Ppltra] = RAD(hd.param[Ppltrah]*15 + + hd.param[Ppltram]/4 + hd.param[Ppltras]/240); + hd.param[Ppltdec] = RAD(hd.param[Ppltdecd] + + hd.param[Ppltdecm]/60 + hd.param[Ppltdecs]/3600); + if(decsn) + hd.param[Ppltdec] = -hd.param[Ppltdec]; + hd.amdflag = 0; + for(i=Pamdx1; i<=Pamdx20; i++) + if(hd.param[i] != 0) { + hd.amdflag = 1; + break; + } + h = malloc(sizeof(*h)); + *h = hd; + return h; +} + +void +getplates(void) +{ + char rec[81], *q; + Plate *p; + Biobuf *bin; + int c, i, dss; + + dss = 0; + sprint(rec, "/lib/sky/dssheaders/lo_comp.lis"); + bin = Bopen(rec, OREAD); + if(bin == 0) { + dss = 102; + sprint(rec, "%s/headers/lo_comp.lis", dssmount(dss)); + bin = Bopen(rec, OREAD); + } + if(bin == 0) { + dss = 61; + sprint(rec, "%s/headers/lo_comp.lis", dssmount(dss)); + bin = Bopen(rec, OREAD); + } + if(bin == 0) { + fprint(2, "can't open lo_comp.lis; try 9fs juke\n"); + exits("open"); + } + if(debug) + Bprint(&bout, "reading %s\n", rec); + if(dss) + Bprint(&bout, "warning: reading %s from jukebox\n", rec); + for(nplate=0;;) { + if(dss) { + if(Bread(bin, rec, 80) != 80) + break; + rec[80] = 0; + } else { + q = Brdline(bin, '\n'); + if(q == 0) + break; + q[Blinelen(bin)-1] = 0; + strcpy(rec, q); + } + if(rec[0] == '#') + continue; + if(nplate < nelem(plate)) { + p = &plate[nplate]; + memmove(p->rgn, rec+0, 5); + if(p->rgn[4] == ' ') + p->rgn[4] = 0; + for(i=0; c=p->rgn[i]; i++) + if(c >= 'A' && c <= 'Z') + p->rgn[i] += 'a'-'A'; + p->ra = RAD(atof(rec+12)*15 + + atof(rec+15)/4 + + atof(rec+18)/240); + p->dec = RAD(atof(rec+26) + + atof(rec+29)/60 + + atof(rec+32)/3600); + if(rec[25] == '-') + p->dec = -p->dec; + p->disk = atoi(rec+53); + if(p->disk == 0) + continue; + } + nplate++; + } + Bterm(bin); + + if(nplate >= nelem(plate)) + fprint(2, "nplate too small %d %d\n", nelem(plate), nplate); + if(debug) + Bprint(&bout, "%d plates\n", nplate); +} + +char* +dssmount(int dskno) +{ + Bprint(&bout, "not mounting dss\n"); + return "/n/dss"; +} + blob - /dev/null blob + c76b00ea127641ca5e77d657d062c9b1e62dd0d0 (mode 644) --- /dev/null +++ src/cmd/scat/hinv.c @@ -0,0 +1,231 @@ +#include +#include +#include +#include "sky.h" + +static void unshuffle(Pix*, int, int, Pix*); +static void unshuffle1(Pix*, int, Pix*); + +void +hinv(Pix *a, int nx, int ny) +{ + int nmax, log2n, h0, hx, hy, hc, i, j, k; + int nxtop, nytop, nxf, nyf, c; + int oddx, oddy; + int shift; + int s10, s00; + Pix *tmp; + + /* + * log2n is log2 of max(nx, ny) rounded up to next power of 2 + */ + nmax = ny; + if(nx > nmax) + nmax = nx; + log2n = log(nmax)/LN2 + 0.5; + if(nmax > (1<=0; k--) { + /* + * this somewhat cryptic code generates the sequence + * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n + */ + c = c>>1; + nxtop = nxtop<<1; + nytop = nytop<<1; + if(nxf <= c) + nxtop--; + else + nxf -= c; + if(nyf <= c) + nytop--; + else + nyf -= c; + + /* + * halve divisors on last pass + */ + if(k == 0) + shift = 0; + + /* + * unshuffle in each dimension to interleave coefficients + */ + for(i = 0; i> 2; + a[s10 ] = (h0 + hx - hy - hc + 2) >> 2; + a[s00+1] = (h0 - hx + hy - hc + 2) >> 2; + a[s00 ] = (h0 - hx - hy + hc + 2) >> 2; + s00 += 2; + s10 += 2; + } + if(oddy) { + /* + * do last element in row if row length is odd + * s00+1, s10+1 are off edge + */ + h0 = a[s00 ] << shift; + hx = a[s10 ] << shift; + a[s10 ] = (h0 + hx + 2) >> 2; + a[s00 ] = (h0 - hx + 2) >> 2; + } + } + if(oddx) { + /* + * do last row if column length is odd + * s10, s10+1 are off edge + */ + s00 = ny*i; + for(j = 0; j> 2; + a[s00 ] = (h0 - hy + 2) >> 2; + s00 += 2; + } + if(oddy) { + /* + * do corner element if both row and column lengths are odd + * s00+1, s10, s10+1 are off edge + */ + h0 = a[s00 ] << shift; + a[s00 ] = (h0 + 2) >> 2; + } + } + } + free(tmp); +} + +static +void +unshuffle(Pix *a, int n, int n2, Pix *tmp) +{ + int i; + int nhalf, twon2, n2xnhalf; + Pix *p1, *p2, *pt; + + twon2 = n2<<1; + nhalf = (n+1)>>1; + n2xnhalf = n2*nhalf; /* pointer to a[i] */ + + /* + * copy 2nd half of array to tmp + */ + pt = tmp; + p1 = &a[n2xnhalf]; /* pointer to a[i] */ + for(i=nhalf; i=0; i--) { + p1 -= twon2; + p2 -= n2; + *p1 = *p2; + } + + /* + * now distribute 2nd half of array (in tmp) to odd elements + */ + pt = tmp; + p1 = &a[n2]; /* pointer to a[i] */ + for(i=1; i> 1; + + /* + * copy 2nd half of array to tmp + */ + pt = tmp; + p1 = &a[nhalf]; /* pointer to a[i] */ + for(i=nhalf; i=0; i--) { + p1 -= 2; + p2--; + *p1 = *p2; + } + + /* + * now distribute 2nd half of array (in tmp) to odd elements + */ + pt = tmp; + p1 = &a[1]; /* pointer to a[i] */ + for(i=1; i +#include +#include +#include "sky.h" + +char rad28[] = "0123456789abcdefghijklmnopqr"; + +Picture* +image(Angle ra, Angle dec, Angle wid, Angle hig) +{ + Pix *p; + uchar *b, *up; + int i, j, sx, sy, x, y; + char file[50]; + Picture *pic; + Img* ip; + int lowx, lowy, higx, higy; + int slowx, slowy, shigx, shigy; + Header *h; + Angle d, bd; + Plate *pp, *bp; + + if(gam.gamma == 0) + gam.gamma = -1; + if(gam.max == gam.min) { + gam.max = 17600; + gam.min = 2500; + } + gam.absgamma = gam.gamma; + gam.neg = 0; + if(gam.absgamma < 0) { + gam.absgamma = -gam.absgamma; + gam.neg = 1; + } + gam.mult1 = 1. / (gam.max - gam.min); + gam.mult2 = 255. * gam.mult1; + + if(nplate == 0) + getplates(); + + bp = 0; + bd = 0; + for(i=0; ira, pp->dec); + if(bp == 0 || d < bd) { + bp = pp; + bd = d; + } + } + + if(debug) + Bprint(&bout, "best plate: %s %s disk %d %s\n", + hms(bp->ra), dms(bp->dec), + bp->disk, bp->rgn); + + h = getheader(bp->rgn); + xypos(h, ra, dec, 0, 0); + if(wid <= 0 || hig <= 0) { + lowx = h->x; + lowy = h->y; + lowx = (lowx/500) * 500; + lowy = (lowy/500) * 500; + higx = lowx + 500; + higy = lowy + 500; + } else { + lowx = h->x - wid*ARCSECONDS_PER_RADIAN*1000 / + (h->param[Pxpixelsz]*h->param[Ppltscale]*2); + lowy = h->y - hig*ARCSECONDS_PER_RADIAN*1000 / + (h->param[Pypixelsz]*h->param[Ppltscale]*2); + higx = h->x + wid*ARCSECONDS_PER_RADIAN*1000 / + (h->param[Pxpixelsz]*h->param[Ppltscale]*2); + higy = h->y + hig*ARCSECONDS_PER_RADIAN*1000 / + (h->param[Pypixelsz]*h->param[Ppltscale]*2); + } + free(h); + + if(lowx < 0) lowx = 0; + if(higx < 0) higx = 0; + if(lowy < 0) lowy = 0; + if(higy < 0) higy = 0; + if(lowx > 14000) lowx = 14000; + if(higx > 14000) higx = 14000; + if(lowy > 14000) lowy = 14000; + if(higy > 14000) higy = 14000; + + if(debug) + Bprint(&bout, "xy on plate: %d,%d %d,%d\n", + lowx,lowy, higx, higy); + + if(lowx >= higx || lowy >=higy) { + Bprint(&bout, "no image found\n"); + return 0; + } + + b = malloc((higx-lowx)*(higy-lowy)*sizeof(*b)); + if(b == 0) { + emalloc: + fprint(2, "malloc error\n"); + return 0; + } + memset(b, 0, (higx-lowx)*(higy-lowy)*sizeof(*b)); + + slowx = lowx/500; + shigx = (higx-1)/500; + slowy = lowy/500; + shigy = (higy-1)/500; + + for(sx=slowx; sx<=shigx; sx++) + for(sy=slowy; sy<=shigy; sy++) { + if(sx < 0 || sx >= nelem(rad28) || sy < 0 || sy >= nelem(rad28)) { + fprint(2, "bad subplate %d %d\n", sy, sx); + free(b); + return 0; + } + sprint(file, "%s/%s/%s.%c%c", + dssmount(bp->disk), + bp->rgn, bp->rgn, + rad28[sy], + rad28[sx]); + + ip = dssread(file); + if(ip == 0) { + fprint(2, "can't read %s: %r\n", file); + free(b); + return 0; + } + + x = sx*500; + y = sy*500; + for(j=0; jny; j++) { + if(y+j < lowy || y+j >= higy) + continue; + p = &ip->a[j*ip->ny]; + up = b + (higy - (y+j+1))*(higx-lowx) + (x - lowx); + for(i=0; inx; i++) { + if(x+i >= lowx && x+i < higx) + *up = dogamma(*p); + up++; + p += 1; + } + } + free(ip); + } + + pic = malloc(sizeof(Picture)); + if(pic == 0){ + free(b); + goto emalloc; + } + pic->minx = lowx; + pic->miny = lowy; + pic->maxx = higx; + pic->maxy = higy; + pic->data = b; + strcpy(pic->name, bp->rgn); + return pic; +} blob - /dev/null blob + c0258a838704510f3ce1860f7bdb33fd996b45d2 (mode 644) --- /dev/null +++ src/cmd/scat/mkfile @@ -0,0 +1,31 @@ +<$PLAN9/src/mkhdr + +TARG=scat +OFILES=scat.$O\ + bitinput.$O\ + desc.$O\ + display.$O\ + dssread.$O\ + header.$O\ + hinv.$O\ + image.$O\ + patch.$O\ + plot.$O\ + posn.$O\ + prose.$O\ + qtree.$O\ + util.$O\ + +HFILES=sky.h +CFLAGS=$CFLAGS -I../map +LDFLAGS=$LDFLAGS -L$X11/lib -lX11 + +SHORTLIB=draw bio 9 +LIB=../map/libmap/libmap.a +<$PLAN9/src/mkone + +scat.$O: strings.c + +$LIB:V: + cd ../map/libmap; mk + blob - /dev/null blob + 27d04edbbf02a18ed0fe7b91d7c61d60a1828514 (mode 644) --- /dev/null +++ src/cmd/scat/patch.c @@ -0,0 +1,101 @@ +#include +#include +#include +#include "sky.h" + +/* + * dec varies from -89 to 89, inclusive. + * ra varies depending on dec; each patch is about 1 square degree. + * + * Northern hemisphere (0<=dec<=89): + * from 0<=dec<=59, ra is every 4m, 360 values + * from 60<=dec<=69, ra is every 8m, 180 values + * from 70<=dec<=79, ra is every 12m, 120 values + * from 80<=dec<=84, ra is every 24m, 60 values + * at dec=85 and 86, ra is every 48m, 30 values + * at dec=87, ra is every 60m, 24 values + * at dec=88, ra is every 120m, 12 values + * at dec=89, ra is 12h, 1 value + * + * Total number of patches in northern hemisphere is therefore: + * 360*60+180*10+120*10+60*5+30*2+24*1+12*1+1 = 24997 + * Total number of patches is therefore + * 2*24997-360 = 49634 (dec=0 has been counted twice) + * (cf. 41253 square degrees in the sky) + */ + +void +radec(int p, int *rah, int *ram, int *deg) +{ + *deg = (p&255)-90; + p >>= 8; + *rah = p/15; + *ram = (p%15)*4; + if(*deg<0) + (*deg)++; +} + +long +patcha(Angle ra, Angle dec) +{ + ra = DEG(ra); + dec = DEG(dec); + if(dec >= 0) + return patch(floor(ra/15), ((int)floor(ra*4))%60, floor(dec)); + dec = -dec; + return patch(floor(ra/15), ((int)floor(ra*4))%60, -floor(dec)); +} + +#define round scatround + +char round[91]={ /* extra 0 is to offset the array */ + /* 0 */ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + /* 10 */ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + /* 20 */ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + /* 30 */ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + /* 40 */ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + /* 50 */ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + /* 60 */ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + /* 70 */ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + /* 80 */ 6, 6, 6, 6, 6, 12, 12, 15, 30, -1, + /* 90 */ +}; + +long +patch(int rah, int ram, int deg) +{ + int ra, dec; + + /* + * patches go from lower limit <= position < upper limit. + * therefore dec ' " can be ignored; always inc. dec degrees. + * the computed angle is then the upper limit (ignoring sign). + * when done, +ve values are shifted down so 90 (0 degrees) is a value; + */ + if(rah<0 || rah>=24 || ram<0 || abs(deg)>=90){ + fprint(2, "scat: patch: bad ra or dec %dh%dm %d\n", rah, ram, deg); + abort(); + } + if(deg < 0) + deg--; + else if(deg < 90) + deg++; + dec = deg+90; + deg = abs(deg); + if(deg<1 || deg>90){ + fprint(2, "scat: patch: panic %dh%dm %d\n", rah, ram, deg); + abort(); + } + if(deg == 90) + ra = 180; + else{ + ra = 15*rah+ram/4; + ra -= ra%round[deg]; + } + /* close the hole at 0 */ + if(dec > 90) + --dec; + if(ra >= 360) + ra -= 360; + return (ra<<8)|dec; +} blob - /dev/null blob + 4d10bb5b24134ee8ebf8fcfba3d91647acae1fc0 (mode 644) --- /dev/null +++ src/cmd/scat/plate.h @@ -0,0 +1,138 @@ +#define RAD(x) ((x)*PI_180) +#define DEG(x) ((x)/PI_180) +#define ARCSECONDS_PER_RADIAN (DEG(1)*3600) +#define input_nybble(infile) input_nbits(infile,4) + +typedef float Angle; /* in radians */ + +enum +{ + /* + * parameters for plate + */ + Pppo1 = 0, + Pppo2, + Pppo3, + Pppo4, + Pppo5, + Pppo6, + Pamdx1, + Pamdx2, + Pamdx3, + Pamdx4, + Pamdx5, + Pamdx6, + Pamdx7, + Pamdx8, + Pamdx9, + Pamdx10, + Pamdx11, + Pamdx12, + Pamdx13, + Pamdx14, + Pamdx15, + Pamdx16, + Pamdx17, + Pamdx18, + Pamdx19, + Pamdx20, + Pamdy1, + Pamdy2, + Pamdy3, + Pamdy4, + Pamdy5, + Pamdy6, + Pamdy7, + Pamdy8, + Pamdy9, + Pamdy10, + Pamdy11, + Pamdy12, + Pamdy13, + Pamdy14, + Pamdy15, + Pamdy16, + Pamdy17, + Pamdy18, + Pamdy19, + Pamdy20, + Ppltscale, + Pxpixelsz, + Pypixelsz, + Ppltra, + Ppltrah, + Ppltram, + Ppltras, + Ppltdec, + Ppltdecd, + Ppltdecm, + Ppltdecs, + Pnparam, +}; + +typedef struct Plate Plate; +struct Plate +{ + char rgn[7]; + char disk; + Angle ra; + Angle dec; +}; + +typedef struct Header Header; +struct Header +{ + float param[Pnparam]; + int amdflag; + + float x; + float y; + float xi; + float eta; +}; +typedef long Type; + +typedef struct Image Image; +struct Image +{ + int nx; + int ny; /* ny is the fast-varying dimension */ + Type a[1]; +}; + +int nplate; +Plate plate[2000]; /* needs to go to 2000 when the north comes */ +double PI_180; +double TWOPI; +int debug; +struct +{ + float min; + float max; + float del; + double gamma; + int neg; +} gam; + +char* hms(Angle); +char* dms(Angle); +double xsqrt(double); +Angle dist(Angle, Angle, Angle, Angle); +Header* getheader(char*); +char* getword(char*, char*); +void amdinv(Header*, Angle, Angle, float, float); +void ppoinv(Header*, Angle, Angle); +void xypos(Header*, Angle, Angle, float, float); +void traneqstd(Header*, Angle, Angle); +Angle getra(char*); +Angle getdec(char*); +void getplates(void); + +Image* dssread(char*); +void hinv(Type*, int, int); +int input_bit(Biobuf*); +int input_nbits(Biobuf*, int); +void qtree_decode(Biobuf*, Type*, int, int, int, int); +void start_inputing_bits(void); +Bitmap* image(Angle, Angle, Angle, Angle); +int dogamma(int); blob - /dev/null blob + e87985ff8b219b434ec55b869aba9f3cef5bfc55 (mode 644) --- /dev/null +++ src/cmd/scat/plot.c @@ -0,0 +1,952 @@ +#include +#include +#include +#include +#include +#include +#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; + if(initdraw(drawerror, nil, nil) < 0){ + 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))(struct place*, 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; itype == 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; itype == 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; + } + if(ra+r0+extrara > ramax) + ramax = ra+r0+extrara; + if(ra-extrara < ramin) + ramin = ra-extrara; + } + + if(decmax > 90*c) + decmax = 90*c; + if(decmin < -90*c) + decmin = -90*c; + if(ramin < 0) + ramin += 360*c; + if(ramax >= 360*c) + ramax -= 360*c; + + if(quantize){ + /* quantize to degree boundaries */ + ramin -= ramin%m5; + if(ramax%m5 != 0) + ramax += m5-(ramax%m5); + if(decmin > 0) + decmin -= decmin%c; + else + decmin -= c - (-decmin)%c; + if(decmax > 0){ + if(decmax%c != 0) + decmax += c-(decmax%c); + }else if(decmax < 0){ + if(decmax%c != 0) + decmax += ((-decmax)%c); + } + } + + if(folded){ + if(ramax-ramin > 270*c){ + fprint(2, "ra range too wide %ld°\n", (ramax-ramin)/c); + return -1; + } + }else if(ramax-ramin > 270*c){ + folded = 1; + return bbox(0, 0, quantize); + } + + return 1; +} + +int +inbbox(DAngle ra, DAngle dec) +{ + int min; + + if(dec < decmin) + return 0; + if(dec > decmax) + return 0; + min = ramin; + if(ramin > ramax){ /* straddling 0h0m */ + min -= 360*c; + if(ra > 180*c) + ra -= 360*c; + } + if(ra < min) + return 0; + if(ra > ramax) + return 0; + return 1; +} + +int +gridra(long mapdec) +{ + mapdec = abs(mapdec)+c/2; + mapdec /= c; + if(mapdec < 60) + return m5; + if(mapdec < 80) + return 2*m5; + if(mapdec < 85) + return 4*m5; + return 8*m5; +} + +#define GREY (nogrey? display->white : grey) + +void +plot(char *flags) +{ + int i, j, k; + char *t; + long x, y; + int ra, dec; + int m; + Point p, pts[10]; + Record *r; + Rectangle rect, r1; + int dx, dy, nogrid, textlevel, nogrey, zenithup; + Image *scr; + char *name, buf[32]; + double v; + + if(plotopen() < 0) + return; + nogrid = 0; + nogrey = 0; + textlevel = 1; + dx = 512; + dy = 512; + zenithup = 0; + for(;;){ + if(t = alpha(flags, "nogrid")){ + nogrid = 1; + flags = t; + continue; + } + if((t = alpha(flags, "zenith")) || (t = alpha(flags, "zenithup")) ){ + zenithup = 1; + flags = t; + continue; + } + if((t = alpha(flags, "notext")) || (t = alpha(flags, "nolabel")) ){ + textlevel = 0; + flags = t; + continue; + } + if((t = alpha(flags, "alltext")) || (t = alpha(flags, "alllabel")) ){ + textlevel = 2; + flags = t; + continue; + } + if(t = alpha(flags, "dx")){ + dx = strtol(t, &t, 0); + if(dx < 100){ + fprint(2, "dx %d too small (min 100) in plot\n", dx); + return; + } + flags = skipbl(t); + continue; + } + if(t = alpha(flags, "dy")){ + dy = strtol(t, &t, 0); + if(dy < 100){ + fprint(2, "dy %d too small (min 100) in plot\n", dy); + return; + } + flags = skipbl(t); + continue; + } + if((t = alpha(flags, "nogrey")) || (t = alpha(flags, "nogray"))){ + nogrey = 1; + flags = skipbl(t); + continue; + } + if(*flags){ + fprint(2, "syntax error in plot\n"); + return; + } + break; + } + flatten(); + folded = 0; + + if(bbox(0, 0, 1) < 0) + return; + if(ramax-ramin<100 || decmax-decmin<100){ + fprint(2, "plot too small\n"); + return; + } + + scr = allocimage(display, Rect(0, 0, dx, dy), RGB24, 0, DBlack); + if(scr == nil){ + fprint(2, "can't allocate image: %r\n"); + return; + } + rect = scr->r; + rect.min.x += 16; + rect = insetrect(rect, 40); + if(setmap(ramin, ramax, decmin, decmax, rect, zenithup) < 0){ + fprint(2, "can't set up map coordinates\n"); + return; + } + if(!nogrid){ + for(x=ramin; ; ){ + for(j=0; jheight; + string(scr, p, GREY, ZP, font, hm5(angle(ra))); + if(x == ramax) + break; + x += gridra(mapdec); + if(x > ramax) + x = ramax; + } + for(y=decmin; y<=decmax; y+=c){ + for(j=0; jheight/2; + string(scr, p, GREY, ZP, font, deg(angle(y))); + p = pts[nelem(pts)-1]; + p.x -= 3+stringwidth(font, deg(angle(y))); + p.y -= font->height/2; + string(scr, p, GREY, ZP, font, deg(angle(y))); + } + } + /* reorder to get planets in front of stars */ + tolast(nil); + tolast("moon"); /* moon is in front of everything... */ + tolast("shadow"); /* ... except the shadow */ + + for(i=0,r=rec; ingc.dec; + ra = r->ngc.ra; + if(folded){ + ra -= 180*c; + if(ra < 0) + ra += 360*c; + } + if(textlevel){ + name = nameof(r); + if(name==nil && textlevel>1 && r->type==SAO){ + snprint(buf, sizeof buf, "SAO%ld", r->index); + name = buf; + } + if(name) + drawname(scr, nogrey? display->white : alphagrey, name, ra, dec); + } + if(r->type == Planet){ + drawplanet(scr, &r->planet, map(ra, dec)); + continue; + } + if(r->type == SAO){ + m = r->sao.mag; + if(m == UNKNOWNMAG) + m = r->sao.mpg; + if(m == UNKNOWNMAG) + continue; + m = dsize(m); + if(m < 3) + fillellipse(scr, map(ra, dec), m, m, nogrey? display->white : lightgrey, ZP); + else{ + ellipse(scr, map(ra, dec), m+1, m+1, 0, display->black, ZP); + fillellipse(scr, map(ra, dec), m, m, display->white, ZP); + } + continue; + } + if(r->type == Abell){ + ellipse(scr, addpt(map(ra, dec), Pt(-3, 2)), 2, 1, 0, lightblue, ZP); + ellipse(scr, addpt(map(ra, dec), Pt(3, 2)), 2, 1, 0, lightblue, ZP); + ellipse(scr, addpt(map(ra, dec), Pt(0, -2)), 1, 2, 0, lightblue, ZP); + continue; + } + switch(r->ngc.type){ + case Galaxy: + j = npixels(r->ngc.diam); + if(j < 4) + j = 4; + if(j > 10) + k = j/3; + else + k = j/2; + ellipse(scr, map(ra, dec), j, k, 0, lightblue, ZP); + break; + + case PlanetaryN: + p = map(ra, dec); + j = npixels(r->ngc.diam); + if(j < 3) + j = 3; + ellipse(scr, p, j, j, 0, green, ZP); + line(scr, Pt(p.x, p.y+j+1), Pt(p.x, p.y+j+4), + Endsquare, Endsquare, 0, green, ZP); + line(scr, Pt(p.x, p.y-(j+1)), Pt(p.x, p.y-(j+4)), + Endsquare, Endsquare, 0, green, ZP); + line(scr, Pt(p.x+j+1, p.y), Pt(p.x+j+4, p.y), + Endsquare, Endsquare, 0, green, ZP); + line(scr, Pt(p.x-(j+1), p.y), Pt(p.x-(j+4), p.y), + Endsquare, Endsquare, 0, green, ZP); + break; + + case DiffuseN: + case NebularCl: + p = map(ra, dec); + j = npixels(r->ngc.diam); + if(j < 4) + j = 4; + r1.min = Pt(p.x-j, p.y-j); + r1.max = Pt(p.x+j, p.y+j); + if(r->ngc.type != DiffuseN) + draw(scr, r1, ocstipple, ocstipple, ZP); + line(scr, Pt(p.x-j, p.y-j), Pt(p.x+j, p.y-j), + Endsquare, Endsquare, 0, green, ZP); + line(scr, Pt(p.x-j, p.y+j), Pt(p.x+j, p.y+j), + Endsquare, Endsquare, 0, green, ZP); + line(scr, Pt(p.x-j, p.y-j), Pt(p.x-j, p.y+j), + Endsquare, Endsquare, 0, green, ZP); + line(scr, Pt(p.x+j, p.y-j), Pt(p.x+j, p.y+j), + Endsquare, Endsquare, 0, green, ZP); + break; + + case OpenCl: + p = map(ra, dec); + j = npixels(r->ngc.diam); + if(j < 4) + j = 4; + fillellipse(scr, p, j, j, ocstipple, ZP); + break; + + case GlobularCl: + j = npixels(r->ngc.diam); + if(j < 4) + j = 4; + p = map(ra, dec); + ellipse(scr, p, j, j, 0, lightgrey, ZP); + line(scr, Pt(p.x-(j-1), p.y), Pt(p.x+j, p.y), + Endsquare, Endsquare, 0, lightgrey, ZP); + line(scr, Pt(p.x, p.y-(j-1)), Pt(p.x, p.y+j), + Endsquare, Endsquare, 0, lightgrey, ZP); + break; + + } + } + flushimage(display, 1); + displayimage(scr); +} + +int +runcommand(char *command, int p[2]) +{ + switch(rfork(RFPROC|RFFDG|RFNOWAIT)){ + case -1: + return -1; + default: + break; + case 0: + dup(p[1], 1); + close(p[0]); + close(p[1]); + execlp("rc", "rc", "-c", command, nil); + fprint(2, "can't exec %s: %r", command); + exits(nil); + } + return 1; +} + +void +parseplanet(char *line, Planetrec *p) +{ + char *fld[6]; + int i, nfld; + char *s; + + if(line[0] == '\0') + return; + line[10] = '\0'; /* terminate name */ + s = strrchr(line, ' '); + if(s == nil) + s = line; + else + s++; + strcpy(p->name, s); + for(i=0; s[i]!='\0'; i++) + p->name[i] = tolower(s[i]); + p->name[i] = '\0'; + nfld = getfields(line+11, fld, nelem(fld), 1, " "); + p->ra = dangle(getra(fld[0])); + p->dec = dangle(getra(fld[1])); + p->az = atof(fld[2])*MILLIARCSEC; + p->alt = atof(fld[3])*MILLIARCSEC; + p->semidiam = atof(fld[4])*1000; + if(nfld > 5) + p->phase = atof(fld[5]); + else + p->phase = 0; +} + +void +astro(char *flags, int initial) +{ + int p[2]; + int i, n, np; + char cmd[256], buf[4096], *lines[20], *fld[10]; + + snprint(cmd, sizeof cmd, "astro -p %s", flags); + if(pipe(p) < 0){ + fprint(2, "can't pipe: %r\n"); + return; + } + if(runcommand(cmd, p) < 0){ + close(p[0]); + close(p[1]); + fprint(2, "can't run astro: %r"); + return; + } + close(p[1]); + n = readn(p[0], buf, sizeof buf-1); + if(n <= 0){ + fprint(2, "no data from astro\n"); + return; + } + if(!initial) + Bwrite(&bout, buf, n); + buf[n] = '\0'; + np = getfields(buf, lines, nelem(lines), 0, "\n"); + if(np <= 1){ + fprint(2, "astro: not enough output\n"); + return; + } + Bprint(&bout, "%s\n", lines[0]); + Bflush(&bout); + /* get latitude and longitude */ + if(getfields(lines[0], fld, nelem(fld), 1, " ") < 8) + fprint(2, "astro: can't read longitude: too few fields\n"); + else{ + mysid = getra(fld[5])*180./PI; + mylat = getra(fld[6])*180./PI; + mylon = getra(fld[7])*180./PI; + } + /* + * Each time we run astro, we generate a new planet list + * so multiple appearances of a planet may exist as we plot + * its motion over time. + */ + planet = malloc(NPlanet*sizeof planet[0]); + if(planet == nil){ + fprint(2, "astro: malloc failed: %r\n"); + exits("malloc"); + } + memset(planet, 0, NPlanet*sizeof planet[0]); + for(i=1; i +#include +#include +#include "sky.h" + +void +amdinv(Header *h, Angle ra, Angle dec, float mag, float col) +{ + int i, max_iterations; + float tolerance; + float object_x, object_y, delta_x, delta_y, f, fx, fy, g, gx, gy; + float x1, x2, x3, x4; + float y1, y2, y3, y4; + + /* + * Initialize + */ + max_iterations = 50; + tolerance = 0.0000005; + + /* + * Convert RA and Dec to St.coords + */ + traneqstd(h, ra, dec); + + /* + * Set initial value for x,y + */ + object_x = h->xi/h->param[Ppltscale]; + object_y = h->eta/h->param[Ppltscale]; + + /* + * Iterate by Newtons method + */ + for(i = 0; i < max_iterations; i++) { + /* + * X plate model + */ + x1 = object_x; + x2 = x1 * object_x; + x3 = x1 * object_x; + x4 = x1 * object_x; + + y1 = object_y; + y2 = y1 * object_y; + y3 = y1 * object_y; + y4 = y1 * object_y; + + f = + h->param[Pamdx1] * x1 + + h->param[Pamdx2] * y1 + + h->param[Pamdx3] + + h->param[Pamdx4] * x2 + + h->param[Pamdx5] * x1*y1 + + h->param[Pamdx6] * y2 + + h->param[Pamdx7] * (x2+y2) + + h->param[Pamdx8] * x2*x1 + + h->param[Pamdx9] * x2*y1 + + h->param[Pamdx10] * x1*y2 + + h->param[Pamdx11] * y3 + + h->param[Pamdx12] * x1* (x2+y2) + + h->param[Pamdx13] * x1 * (x2+y2) * (x2+y2) + + h->param[Pamdx14] * mag + + h->param[Pamdx15] * mag*mag + + h->param[Pamdx16] * mag*mag*mag + + h->param[Pamdx17] * mag*x1 + + h->param[Pamdx18] * mag * (x2+y2) + + h->param[Pamdx19] * mag*x1 * (x2+y2) + + h->param[Pamdx20] * col; + + /* + * Derivative of X model wrt x + */ + fx = + h->param[Pamdx1] + + h->param[Pamdx4] * 2*x1 + + h->param[Pamdx5] * y1 + + h->param[Pamdx7] * 2*x1 + + h->param[Pamdx8] * 3*x2 + + h->param[Pamdx9] * 2*x1*y1 + + h->param[Pamdx10] * y2 + + h->param[Pamdx12] * (3*x2+y2) + + h->param[Pamdx13] * (5*x4 + 6*x2*y2 + y4) + + h->param[Pamdx17] * mag + + h->param[Pamdx18] * mag*2*x1 + + h->param[Pamdx19] * mag*(3*x2+y2); + + /* + * Derivative of X model wrt y + */ + fy = + h->param[Pamdx2] + + h->param[Pamdx5] * x1 + + h->param[Pamdx6] * 2*y1 + + h->param[Pamdx7] * 2*y1 + + h->param[Pamdx9] * x2 + + h->param[Pamdx10] * x1*2*y1 + + h->param[Pamdx11] * 3*y2 + + h->param[Pamdx12] * 2*x1*y1 + + h->param[Pamdx13] * 4*x1*y1* (x2+y2) + + h->param[Pamdx18] * mag*2*y1 + + h->param[Pamdx19] * mag*2*x1*y1; + /* + * Y plate model + */ + g = + h->param[Pamdy1] * y1 + + h->param[Pamdy2] * x1 + + h->param[Pamdy3] + + h->param[Pamdy4] * y2 + + h->param[Pamdy5] * y1*x1 + + h->param[Pamdy6] * x2 + + h->param[Pamdy7] * (x2+y2) + + h->param[Pamdy8] * y3 + + h->param[Pamdy9] * y2*x1 + + h->param[Pamdy10] * y1*x3 + + h->param[Pamdy11] * x3 + + h->param[Pamdy12] * y1 * (x2+y2) + + h->param[Pamdy13] * y1 * (x2+y2) * (x2+y2) + + h->param[Pamdy14] * mag + + h->param[Pamdy15] * mag*mag + + h->param[Pamdy16] * mag*mag*mag + + h->param[Pamdy17] * mag*y1 + + h->param[Pamdy18] * mag * (x2+y2) + + h->param[Pamdy19] * mag*y1 * (x2+y2) + + h->param[Pamdy20] * col; + + /* + * Derivative of Y model wrt x + */ + gx = + h->param[Pamdy2] + + h->param[Pamdy5] * y1 + + h->param[Pamdy6] * 2*x1 + + h->param[Pamdy7] * 2*x1 + + h->param[Pamdy9] * y2 + + h->param[Pamdy10] * y1*2*x1 + + h->param[Pamdy11] * 3*x2 + + h->param[Pamdy12] * 2*x1*y1 + + h->param[Pamdy13] * 4*x1*y1 * (x2+y2) + + h->param[Pamdy18] * mag*2*x1 + + h->param[Pamdy19] * mag*y1*2*x1; + + /* + * Derivative of Y model wrt y + */ + gy = + h->param[Pamdy1] + + h->param[Pamdy4] * 2*y1 + + h->param[Pamdy5] * x1 + + h->param[Pamdy7] * 2*y1 + + h->param[Pamdy8] * 3*y2 + + h->param[Pamdy9] * 2*y1*x1 + + h->param[Pamdy10] * x2 + + h->param[Pamdy12] * 3*y2 + + h->param[Pamdy13] * (5*y4 + 6*x2*y2 + x4) + + h->param[Pamdy17] * mag + + h->param[Pamdy18] * mag*2*y1 + + h->param[Pamdy19] * mag*(x2 + 3*y2); + + f = f - h->xi; + g = g - h->eta; + delta_x = (-f*gy+g*fy) / (fx*gy-fy*gx); + delta_y = (-g*fx+f*gx) / (fx*gy-fy*gx); + object_x = object_x + delta_x; + object_y = object_y + delta_y; + if((fabs(delta_x) < tolerance) && (fabs(delta_y) < tolerance)) + break; + } + + /* + * Convert mm from plate center to pixels + */ + h->x = (h->param[Pppo3]-object_x*1000.0)/h->param[Pxpixelsz]; + h->y = (h->param[Pppo6]+object_y*1000.0)/h->param[Pypixelsz]; +} + +void +ppoinv(Header *h, Angle ra, Angle dec) +{ + + /* + * Convert RA and Dec to standard coords. + */ + traneqstd(h, ra, dec); + + /* + * Convert st.coords from arcsec to radians + */ + h->xi /= ARCSECONDS_PER_RADIAN; + h->eta /= ARCSECONDS_PER_RADIAN; + + /* + * Compute PDS coordinates from solution + */ + h->x = + h->param[Pppo1]*h->xi + + h->param[Pppo2]*h->eta + + h->param[Pppo3]; + h->y = + h->param[Pppo4]*h->xi + + h->param[Pppo5]*h->eta + + h->param[Pppo6]; + /* + * Convert x,y from microns to pixels + */ + h->x /= h->param[Pxpixelsz]; + h->y /= h->param[Pypixelsz]; + +} + +void +traneqstd(Header *h, Angle object_ra, Angle object_dec) +{ + float div; + + /* + * Find divisor + */ + div = + (sin(object_dec)*sin(h->param[Ppltdec]) + + cos(object_dec)*cos(h->param[Ppltdec]) * + cos(object_ra - h->param[Ppltra])); + + /* + * Compute standard coords and convert to arcsec + */ + h->xi = cos(object_dec) * + sin(object_ra - h->param[Ppltra]) * + ARCSECONDS_PER_RADIAN/div; + + h->eta = (sin(object_dec)*cos(h->param[Ppltdec])- + cos(object_dec)*sin(h->param[Ppltdec])* + cos(object_ra - h->param[Ppltra]))* + ARCSECONDS_PER_RADIAN/div; + +} + +void +xypos(Header *h, Angle ra, Angle dec, float mag, float col) +{ + if (h->amdflag) { + amdinv(h, ra, dec, mag, col); + } else { + ppoinv(h, ra, dec); + } +} blob - /dev/null blob + 17a09750f7f7eb37a4e956dffa47d5c11957c804 (mode 644) --- /dev/null +++ src/cmd/scat/prose.c @@ -0,0 +1,168 @@ +#include +#include +#include +#include "sky.h" + +extern Biobuf bout; + +char* +append(char *p, char *s) +{ + while(*s) + *p++ = *s++; + return p; +} + +int +matchlen(char *a, char *b) +{ + int n; + + for(n=0; *a==*b; a++, b++, n++) + if(*a == 0) + return n; + if(*a == 0) + return n; + return 0; +} + +char* +prose(char *s, char *desc[][2], short index[]) +{ + static char buf[512]; + char *p=buf; + int i, j, k, max; + + j = 0; + while(*s){ + if(p >= buf+sizeof buf) + abort(); + if(*s == ' '){ + if(p>buf && p[-1]!=' ') + *p++ = ' '; + s++; + continue; + } + if(*s == ','){ + *p++ = ';', s++; + continue; + } + if(s[0]=='M' && '0'<=s[1] && s[1]<='9'){ /* Messier tag */ + *p++ = *s++; + continue; /* below will copy the number */ + } + if((i=index[(uchar)*s]) == -1){ + Dup: + switch(*s){ + default: + while(*s && *s!=',' && *s!=' ') + *p++=*s++; + break; + + case '0': case '1': case '2': case '3': case '4': + case '5': case '6': case '7': case '8': case '9': + while('0'<=*s && *s<='9') + *p++ = *s++; + if(*s=='\'' || *s=='s') + *p++ = *s++; + break; + + case '(': case ')': + case '\'': case '"': + case '&': case '-': case '+': + *p++ = *s++; + break; + + case '*': + if('0'<=s[1] && s[1]<='9'){ + int flag=0; + s++; + Pnumber: + while('0'<=*s && *s<='9') + *p++=*s++; + if(s[0] == '-'){ + *p++ = *s++; + flag++; + goto Pnumber; + } + if(s[0]==',' && s[1]==' ' && '0'<=s[2] && s[2]<='9'){ + *p++ = *s++; + s++; /* skip blank */ + flag++; + goto Pnumber; + } + if(s[0] == '.'){ + if(s[1]=='.' && s[2]=='.'){ + *p++ = '-'; + s += 3; + flag++; + goto Pnumber; + } + *p++ = *s++; + goto Pnumber; + } + p = append(p, "m star"); + if(flag) + *p++ = 's'; + *p++ = ' '; + break; + } + if(s[1] == '*'){ + if(s[2] == '*'){ + p = append(p, "triple star "); + s += 3; + }else{ + p = append(p, "double star "); + s += 2; + } + break; + } + p = append(p, "star "); + s++; + break; + } + continue; + } + for(max=-1; desc[i][0] && desc[i][0][0]==*s; i++){ + k = matchlen(desc[i][0], s); + if(k > max) + max = k, j = i; + } + if(max == 0) + goto Dup; + s += max; + for(k=0; desc[j][1][k]; k++) + *p++=desc[j][1][k]; + if(*s == ' ') + *p++ = *s++; + else if(*s == ',') + *p++ = ';', s++; + else + *p++ = ' '; + } + *p = 0; + return buf; +} + +void +prdesc(char *s, char *desc[][2], short index[]) +{ + int c, j; + + if(index[0] == 0){ + index[0] = 1; + for(c=1, j=0; c<128; c++) + if(desc[j][0]==0 || desc[j][0][0]>c) + index[c] = -1; + else if(desc[j][0][0] == c){ + index[c] = j; + while(desc[j][0] && desc[j][0][0] == c) + j++; + if(j >= NINDEX){ + fprint(2, "scat: internal error: too many prose entries\n"); + exits("NINDEX"); + } + } + } + Bprint(&bout, "\t%s [%s]\n", prose(s, desc, index), s); +} blob - /dev/null blob + a3549c65ef45be1266a682310c5dfbc47ea8866b (mode 644) --- /dev/null +++ src/cmd/scat/qtree.c @@ -0,0 +1,278 @@ +#include +#include +#include +#include "sky.h" + +static void qtree_expand(Biobuf*, uchar*, int, int, uchar*); +static void qtree_copy(uchar*, int, int, uchar*, int); +static void qtree_bitins(uchar*, int, int, Pix*, int, int); +static void read_bdirect(Biobuf*, Pix*, int, int, int, uchar*, int); + +void +qtree_decode(Biobuf *infile, Pix *a, int n, int nqx, int nqy, int nbitplanes) +{ + int log2n, k, bit, b, nqmax; + int nx,ny,nfx,nfy,c; + int nqx2, nqy2; + unsigned char *scratch; + + /* + * log2n is log2 of max(nqx,nqy) rounded up to next power of 2 + */ + nqmax = nqy; + if(nqx > nqmax) + nqmax = nqx; + log2n = log(nqmax)/LN2+0.5; + if (nqmax > (1<= 0; bit--) { + /* + * Was bitplane was quadtree-coded or written directly? + */ + b = input_nybble(infile); + if(b == 0) { + /* + * bit map was written directly + */ + read_bdirect(infile, a, n, nqx, nqy, scratch, bit); + } else + if(b != 0xf) { + fprint(2, "qtree_decode: bad format code %x\n",b); + exits("format"); + } else { + /* + * bitmap was quadtree-coded, do log2n expansions + * read first code + */ + + scratch[0] = input_huffman(infile); + + /* + * do log2n expansions, reading codes from file as necessary + */ + nx = 1; + ny = 1; + nfx = nqx; + nfy = nqy; + c = 1<>1; + nx = nx<<1; + ny = ny<<1; + if(nfx <= c) + nx--; + else + nfx -= c; + if(nfy <= c) + ny--; + else + nfy -= c; + qtree_expand(infile, scratch, nx, ny, scratch); + } + + /* + * copy last set of 4-bit codes to bitplane bit of array a + */ + qtree_bitins(scratch, nqx, nqy, a, n, bit); + } + } + free(scratch); +} + +/* + * do one quadtree expansion step on array a[(nqx+1)/2,(nqy+1)/2] + * results put into b[nqx,nqy] (which may be the same as a) + */ +static +void +qtree_expand(Biobuf *infile, uchar *a, int nx, int ny, uchar *b) +{ + uchar *b1; + + /* + * first copy a to b, expanding each 4-bit value + */ + qtree_copy(a, nx, ny, b, ny); + + /* + * now read new 4-bit values into b for each non-zero element + */ + b1 = &b[nx*ny]; + while(b1 > b) { + b1--; + if(*b1 != 0) + *b1 = input_huffman(infile); + } +} + +/* + * copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding + * each value to 2x2 pixels + * a,b may be same array + */ +static +void +qtree_copy(uchar *a, int nx, int ny, uchar *b, int n) +{ + int i, j, k, nx2, ny2; + int s00, s10; + + /* + * first copy 4-bit values to b + * start at end in case a,b are same array + */ + nx2 = (nx+1)/2; + ny2 = (ny+1)/2; + k = ny2*(nx2-1) + ny2-1; /* k is index of a[i,j] */ + for (i = nx2-1; i >= 0; i--) { + s00 = 2*(n*i+ny2-1); /* s00 is index of b[2*i,2*j] */ + for (j = ny2-1; j >= 0; j--) { + b[s00] = a[k]; + k -= 1; + s00 -= 2; + } + } + + /* + * now expand each 2x2 block + */ + for(i = 0; i>1) & 1; + b[s00+1] = (b[s00]>>2) & 1; + b[s00 ] = (b[s00]>>3) & 1; + s00 += 2; + s10 += 2; + } + if(j < ny) { + /* + * row size is odd, do last element in row + * s00+1, s10+1 are off edge + */ + b[s10 ] = (b[s00]>>1) & 1; + b[s00 ] = (b[s00]>>3) & 1; + } + } + if(i < nx) { + /* + * column size is odd, do last row + * s10, s10+1 are off edge + */ + s00 = n*i; + for (j = 0; j>2) & 1; + b[s00 ] = (b[s00]>>3) & 1; + s00 += 2; + } + if(j < ny) { + /* + * both row and column size are odd, do corner element + * s00+1, s10, s10+1 are off edge + */ + b[s00 ] = (b[s00]>>3) & 1; + } + } +} + +/* + * Copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding + * each value to 2x2 pixels and inserting into bitplane BIT of B. + * A,B may NOT be same array (it wouldn't make sense to be inserting + * bits into the same array anyway.) + */ +static +void +qtree_bitins(uchar *a, int nx, int ny, Pix *b, int n, int bit) +{ + int i, j; + Pix *s00, *s10; + Pix px; + + /* + * expand each 2x2 block + */ + for(i=0; i>1) & 1) << bit; + s00[1] |= ((px>>2) & 1) << bit; + s00[0] |= ((px>>3) & 1) << bit; + s00 += 2; + s10 += 2; + } + if(j < ny) { + /* + * row size is odd, do last element in row + * s00+1, s10+1 are off edge + */ + px = *a++; + s10[0] |= ((px>>1) & 1) << bit; + s00[0] |= ((px>>3) & 1) << bit; + } + } + if(i < nx) { + /* + * column size is odd, do last row + * s10, s10+1 are off edge + */ + s00 = &b[n*i]; + for(j=0; j>2) & 1) << bit; + s00[0] |= ((px>>3) & 1) << bit; + s00 += 2; + } + if(j < ny) { + /* + * both row and column size are odd, do corner element + * s00+1, s10, s10+1 are off edge + */ + s00[0] |= ((*a>>3) & 1) << bit; + } + } +} + +static +void +read_bdirect(Biobuf *infile, Pix *a, int n, int nqx, int nqy, uchar *scratch, int bit) +{ + int i; + + /* + * read bit image packed 4 pixels/nybble + */ + for(i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) { + scratch[i] = input_nybble(infile); + } + + /* + * insert in bitplane BIT of image A + */ + qtree_bitins(scratch, nqx, nqy, a, n, bit); +} blob - /dev/null blob + 9579706aab71f1388d12277b0e958a83c6dfadba (mode 644) --- /dev/null +++ src/cmd/scat/scat.c @@ -0,0 +1,1671 @@ +#include +#include +#include +#include +#include +#include "sky.h" +#include "strings.c" + +enum +{ + NNGC=7840, /* number of NGC numbers [1..NNGC] */ + NIC = 5386, /* number of IC numbers */ + NNGCrec=NNGC+NIC, /* number of records in the NGC catalog (including IC's, starting at NNGC */ + NMrec=122, /* number of M records */ + NM=110, /* number of M numbers */ + NAbell=2712, /* number of records in the Abell catalog */ + NName=1000, /* number of prose names; estimated maximum (read from editable text file) */ + NBayer=1517, /* number of bayer entries */ + NSAO=258998, /* number of SAO stars */ + MAXcon=1932, /* maximum number of patches in a constellation */ + Ncon=88, /* number of constellations */ + Npatch=92053, /* highest patch number */ +}; + +char ngctype[NNGCrec]; +Mindexrec mindex[NMrec]; +Namerec name[NName]; +Bayerec bayer[NBayer]; +long con[MAXcon]; +ushort conindex[Ncon+1]; +long patchaddr[Npatch+1]; + +Record *rec; +Record *orec; +Record *cur; + +char *dir; +int saodb; +int ngcdb; +int abelldb; +int ngctypedb; +int mindexdb; +int namedb; +int bayerdb; +int condb; +int conindexdb; +int patchdb; +char parsed[3]; +long nrec; +long nreca; +long norec; +long noreca; + +Biobuf bin; +Biobuf bout; + +void +main(int argc, char *argv[]) +{ + char *line; + + dir = unsharp(DIR); + Binit(&bin, 0, OREAD); + Binit(&bout, 1, OWRITE); + if(argc != 1) + dir = argv[1]; + astro("", 1); + while(line = Brdline(&bin, '\n')){ + line[Blinelen(&bin)-1] = 0; + lookup(line, 1); + Bflush(&bout); + } + if(display != nil){ + closedisplay(display); + /* automatic refresh of rio window is triggered by mouse */ + // close(open("/dev/mouse", OREAD)); + } + return; +} + +void +reset(void) +{ + nrec = 0; + cur = rec; +} + +void +grow(void) +{ + nrec++; + if(nreca < nrec){ + nreca = nrec+50; + rec = realloc(rec, nreca*sizeof(Record)); + if(rec == 0){ + fprint(2, "scat: realloc fails\n"); + exits("realloc"); + } + } + cur = rec+nrec-1; +} + +void +copy(void) +{ + if(noreca < nreca){ + noreca = nreca; + orec = realloc(orec, nreca*sizeof(Record)); + if(orec == 0){ + fprint(2, "scat: realloc fails\n"); + exits("realloc"); + } + } + memmove(orec, rec, nrec*sizeof(Record)); + norec = nrec; +} + +int +eopen(char *s) +{ + char buf[128]; + int f; + + sprint(buf, "%s/%s.scat", dir, s); + f = open(buf, 0); + if(f<0){ + fprint(2, "scat: can't open %s\n", buf); + exits("open"); + } + return f; +} + + +void +Eread(int f, char *name, void *addr, long n) +{ + if(read(f, addr, n) != n){ /* BUG! */ + fprint(2, "scat: read error on %s\n", name); + exits("read"); + } +} + +char* +skipbl(char *s) +{ + while(*s!=0 && (*s==' ' || *s=='\t')) + s++; + return s; +} + +char* +skipstr(char *s, char *t) +{ + while(*s && *s==*t) + s++, t++; + return skipbl(s); +} + +/* produce little-endian long at address l */ +long +Long(long *l) +{ + uchar *p; + + p = (uchar*)l; + return (long)p[0]|((long)p[1]<<8)|((long)p[2]<<16)|((long)p[3]<<24); +} + +/* produce little-endian long at address l */ +int +Short(short *s) +{ + uchar *p; + + p = (uchar*)s; + return p[0]|(p[1]<<8); +} + +void +nameopen(void) +{ + Biobuf b; + int i; + char *l, *p; + + if(namedb == 0){ + namedb = eopen("name"); + Binit(&b, namedb, OREAD); + for(i=0; itype = NGC; + cur->index = index; + seek(ngcdb, j, 0); + /* special case: NGC data may not be available */ + if(read(ngcdb, &cur->ngc, sizeof(NGCrec)) != sizeof(NGCrec)){ + if(!failed){ + fprint(2, "scat: NGC database not available\n"); + failed++; + } + cur->type = NONGC; + cur->ngc.ngc = 0; + cur->ngc.ra = 0; + cur->ngc.dec = 0; + cur->ngc.diam = 0; + cur->ngc.mag = 0; + return 0; + } + cur->ngc.ngc = Short(&cur->ngc.ngc); + cur->ngc.ra = Long(&cur->ngc.ra); + cur->ngc.dec = Long(&cur->ngc.dec); + cur->ngc.diam = Long(&cur->ngc.diam); + cur->ngc.mag = Short(&cur->ngc.mag); + return 1; +} + +int +loadabell(long index) +{ + long j; + + abellopen(); + j = index-1; + grow(); + cur->type = Abell; + cur->index = index; + seek(abelldb, j*sizeof(Abellrec), 0); + Eread(abelldb, "abell", &cur->abell, sizeof(Abellrec)); + cur->abell.abell = Short(&cur->abell.abell); + if(cur->abell.abell != index){ + fprint(2, "bad format in abell catalog\n"); + exits("abell"); + } + cur->abell.ra = Long(&cur->abell.ra); + cur->abell.dec = Long(&cur->abell.dec); + cur->abell.glat = Long(&cur->abell.glat); + cur->abell.glong = Long(&cur->abell.glong); + cur->abell.rad = Long(&cur->abell.rad); + cur->abell.mag10 = Short(&cur->abell.mag10); + cur->abell.pop = Short(&cur->abell.pop); + cur->abell.dist = Short(&cur->abell.dist); + return 1; +} + +int +loadsao(int index) +{ + if(index<=0 || index>NSAO) + return 0; + saoopen(); + grow(); + cur->type = SAO; + cur->index = index; + seek(saodb, (index-1)*sizeof(SAOrec), 0); + Eread(saodb, "sao", &cur->sao, sizeof(SAOrec)); + cur->sao.ra = Long(&cur->sao.ra); + cur->sao.dec = Long(&cur->sao.dec); + cur->sao.dra = Long(&cur->sao.dra); + cur->sao.ddec = Long(&cur->sao.ddec); + cur->sao.mag = Short(&cur->sao.mag); + cur->sao.mpg = Short(&cur->sao.mpg); + cur->sao.hd = Long(&cur->sao.hd); + return 1; +} + +int +loadplanet(int index, Record *r) +{ + if(index<0 || index>NPlanet || planet[index].name[0]=='\0') + return 0; + grow(); + cur->type = Planet; + cur->index = index; + /* check whether to take new or existing record */ + if(r == nil) + memmove(&cur->planet, &planet[index], sizeof(Planetrec)); + else + memmove(&cur->planet, &r->planet, sizeof(Planetrec)); + return 1; +} + +int +loadpatch(long index) +{ + int i; + + patchopen(); + if(index<=0 || index>Npatch) + return 0; + grow(); + cur->type = Patch; + cur->index = index; + seek(patchdb, patchaddr[index-1], 0); + cur->patch.nkey = (patchaddr[index]-patchaddr[index-1])/4; + Eread(patchdb, "patch", cur->patch.key, cur->patch.nkey*4); + for(i=0; ipatch.nkey; i++) + cur->patch.key[i] = Long(&cur->patch.key[i]); + return 1; +} + +int +loadtype(int t) +{ + int i; + + ngcopen(); + for(i=0; itype = NGCN; + cur->index = i+1; + } + return 1; +} + +void +flatten(void) +{ + int i, j, notflat; + Record *or; + long key; + + loop: + copy(); + reset(); + notflat = 0; + for(i=0,or=orec; itype){ + default: + fprint(2, "bad type %d in flatten\n", or->type); + break; + + case NONGC: + break; + + case Planet: + case Abell: + case NGC: + case SAO: + grow(); + memmove(cur, or, sizeof(Record)); + break; + + case NGCN: + if(loadngc(or->index)) + notflat = 1; + break; + + case NamedSAO: + loadsao(or->index); + notflat = 1; + break; + + case NamedNGC: + if(loadngc(or->index)) + notflat = 1; + break; + + case NamedAbell: + loadabell(or->index); + notflat = 1; + break; + + case PatchC: + loadpatch(or->index); + notflat = 1; + break; + + case Patch: + for(j=1; jpatch.nkey; j++){ + key = or->patch.key[j]; + if((key&0x3F) == SAO) + loadsao((key>>8)&0xFFFFFF); + else if((key&0x3F) == Abell) + loadabell((key>>8)&0xFFFFFF); + else + loadngc((key>>16)&0xFFFF); + } + break; + } + } + if(notflat) + goto loop; +} + +int +ism(int index) +{ + int i; + + for(i=0; i'){ + dogrtr = 1; + mgrtr = 10 * strtod(s+1, &t); + if(mgrtr==0 && t==s+1){ + fprint(2, "bad magnitude\n"); + return 0; + } + s = skipbl(t); + continue; + } + if(s[0] == '<'){ + doless = 1; + mless = 10 * strtod(s+1, &t); + if(mless==0 && t==s+1){ + fprint(2, "bad magnitude\n"); + return 0; + } + s = skipbl(t); + continue; + } + if(t = text(s, "m")){ + dom = 1; + s = t; + continue; + } + if(t = text(s, "sao")){ + dosao = 1; + s = t; + continue; + } + if(t = text(s, "ngc")){ + dongc = 1; + s = t; + continue; + } + if(t = text(s, "abell")){ + doabell = 1; + s = t; + continue; + } + for(i=0; names[i].name; i++) + if(t = alpha(s, names[i].name)){ + if(nobj > 100){ + fprint(2, "too many object types\n"); + return 0; + } + obj[nobj++] = names[i].type; + s = t; + goto Continue; + } + break; + Continue:; + } + if(*s){ + fprint(2, "syntax error in object list\n"); + return 0; + } + + Cull: + flatten(); + copy(); + reset(); + if(dom) + mopen(); + if(dosao) + saoopen(); + if(dongc || nobj) + ngcopen(); + if(doabell) + abellopen(); + for(i=0,or=orec; ingc.ra, or->ngc.dec)) + keepthis = keep; + if(doless && or->ngc.mag <= mless) + keepthis = keep; + if(dogrtr && or->ngc.mag >= mgrtr) + keepthis = keep; + if(dom && (or->type==NGC && ism(or->ngc.ngc))) + keepthis = keep; + if(dongc && or->type==NGC) + keepthis = keep; + if(doabell && or->type==Abell) + keepthis = keep; + if(dosao && or->type==SAO) + keepthis = keep; + for(j=0; jtype==NGC && or->ngc.type==obj[j]) + keepthis = keep; + if(keepthis){ + grow(); + memmove(cur, or, sizeof(Record)); + } + } + return 1; +} + +int +compar(const void *va, const void *vb) +{ + Record *a=(Record*)va, *b=(Record*)vb; + + if(a->type == b->type) + return a->index - b->index; + return a->type - b->type; +} + +void +sort(void) +{ + int i; + Record *r, *s; + + if(nrec == 0) + return; + qsort(rec, nrec, sizeof(Record), compar); + r = rec+1; + s = rec; + for(i=1; itype==s->type && r->index==s->index && r->type!=Planet) + continue; + memmove(++s, r, sizeof(Record)); + } + nrec = (s+1)-rec; +} + +char greekbuf[128]; + +char* +togreek(char *s) +{ + char *t; + int i, n; + Rune r; + + t = greekbuf; + while(*s){ + for(i=1; i<=24; i++){ + n = strlen(greek[i]); + if(strncmp(s, greek[i], n)==0 && (s[n]==' ' || s[n]=='\t')){ + s += n; + t += runetochar(t, &greeklet[i]); + goto Cont; + } + } + n = chartorune(&r, s); + for(i=0; itype == Planet) /* must keep it here */ + loadplanet(or->index, or); + dec = or->ngc.dec/MILLIARCSEC; + ra = or->ngc.ra/MILLIARCSEC; + rdeg = deg/cos((dec*PI)/180); + for(y=-deg; y<=+deg; y++){ + ndec = dec*2+y; + if(ndec/2>=90 || ndec/2<=-90) + continue; + /* fp errors hurt here, so we round 1' to the pole */ + if(ndec >= 0) + ndec = ndec*500*60*60 + 60000; + else + ndec = ndec*500*60*60 - 60000; + for(x=-rdeg; x<=+rdeg; x++){ + nra = ra*2+x; + if(nra/2 < 0) + nra += 360*2; + if(nra/2 >= 360) + nra -= 360*2; + /* fp errors hurt here, so we round up 1' */ + nra = nra/2*MILLIARCSEC + 60000; + loadpatch(patcha(angle(nra), angle(ndec))); + } + } + } + sort(); + return 1; +} +#endif + +/* + * New version attempts to match the boundaries of the plot better. + */ +int +coords(int deg) +{ + int i; + int x, y, xx; + Record *or; + long min, circle; + double factor; + + flatten(); + circle = 360*MILLIARCSEC; + deg *= MILLIARCSEC; + + /* find center */ + folded = 0; + bbox(0, 0, 0); + /* now expand */ + factor = cos(angle((decmax+decmin)/2)); + if(factor < .2) + factor = .2; + factor = floor(1/factor); + folded = 0; + bbox(factor*deg, deg, 1); + Bprint(&bout, "%s to ", hms(angle(ramin))); + Bprint(&bout, "%s\n", hms(angle(ramax))); + Bprint(&bout, "%s to ", dms(angle(decmin))); + Bprint(&bout, "%s\n", dms(angle(decmax))); + copy(); + reset(); + for(i=0,or=orec; itype == Planet) /* must keep it here */ + loadplanet(or->index, or); + min = ramin; + if(ramin > ramax) + min -= circle; + for(x=min; x<=ramax; x+=250*60*60){ + xx = x; + if(xx < 0) + xx += circle; + for(y=decmin; y<=decmax; y+=250*60*60) + if(-circle/4 < y && y < circle/4) + loadpatch(patcha(angle(xx), angle(y))); + } + sort(); + cull(nil, 1, 1); + return 1; +} + +void +pplate(char *flags) +{ + int i; + long c; + int na, rah, ram, d1, d2; + double r0; + int ra, dec; + long ramin, ramax, decmin, decmax; /* all in degrees */ + Record *r; + int folded; + Angle racenter, deccenter, rasize, decsize, a[4]; + Picture *pic; + + rasize = -1.0; + decsize = -1.0; + na = 0; + for(;;){ + while(*flags==' ') + flags++; + if(('0'<=*flags && *flags<='9') || *flags=='+' || *flags=='-'){ + if(na >= 3) + goto err; + a[na++] = getra(flags); + while(*flags && *flags!=' ') + flags++; + continue; + } + if(*flags){ + err: + Bprint(&bout, "syntax error in plate\n"); + return; + } + break; + } + switch(na){ + case 0: + break; + case 1: + rasize = a[0]; + decsize = rasize; + break; + case 2: + rasize = a[0]; + decsize = a[1]; + break; + case 3: + case 4: + racenter = a[0]; + deccenter = a[1]; + rasize = a[2]; + if(na == 4) + decsize = a[3]; + else + decsize = rasize; + if(rasize<0.0 || decsize<0.0){ + Bprint(&bout, "negative sizes\n"); + return; + } + goto done; + } + folded = 0; + /* convert to milliarcsec */ + c = 1000*60*60; + Again: + if(nrec == 0){ + Bprint(&bout, "empty\n"); + return; + } + ramin = 0x7FFFFFFF; + ramax = -0x7FFFFFFF; + decmin = 0x7FFFFFFF; + decmax = -0x7FFFFFFF; + for(r=rec,i=0; itype == Patch){ + radec(r->index, &rah, &ram, &dec); + ra = 15*rah+ram/4; + r0 = c/cos(RAD(dec)); + 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==Abell){ + ra = r->ngc.ra; + dec = r->ngc.dec; + d1 = 0, d2 = 0, r0 = 0; + }else if(r->type==NGCN){ + loadngc(r->index); + continue; + }else if(r->type==NamedSAO){ + loadsao(r->index); + continue; + }else if(r->type==NamedNGC){ + loadngc(r->index); + continue; + }else if(r->type==NamedAbell){ + loadabell(r->index); + continue; + }else + continue; + if(dec+d2 > decmax) + decmax = dec+d2; + if(dec-d1 < decmin) + decmin = dec-d1; + if(folded){ + ra -= 180*c; + if(ra < 0) + ra += 360*c; + } + if(ra+r0 > ramax) + ramax = ra+r0; + if(ra < ramin) + ramin = ra; + } + if(!folded && ramax-ramin>270*c){ + folded = 1; + goto Again; + } + racenter = angle(ramin+(ramax-ramin)/2); + deccenter = angle(decmin+(decmax-decmin)/2); + if(rasize<0 || decsize<0){ + rasize = angle(ramax-ramin)*cos(deccenter); + decsize = angle(decmax-decmin); + } + done: + if(DEG(rasize)>1.1 || DEG(decsize)>1.1){ + Bprint(&bout, "plate too big: %s", ms(rasize)); + Bprint(&bout, " x %s\n", ms(decsize)); + Bprint(&bout, "trimming to 30'x30'\n"); + rasize = RAD(0.5); + decsize = RAD(0.5); + } + Bprint(&bout, "%s %s ", hms(racenter), dms(deccenter)); + Bprint(&bout, "%s", ms(rasize)); + Bprint(&bout, " x %s\n", ms(decsize)); + Bflush(&bout); + flatten(); + pic = image(racenter, deccenter, rasize, decsize); + if(pic == 0) + return; + Bprint(&bout, "plate %s locn %d %d %d %d\n", pic->name, pic->minx, pic->miny, pic->maxx, pic->maxy); + Bflush(&bout); + displaypic(pic); +} + +void +lookup(char *s, int doreset) +{ + int i, j, k; + int rah, ram, deg; + char *starts, *inputline=s, *t, *u; + Record *r; + Rune c; + long n; + double x; + Angle ra; + + lowercase(s); + s = skipbl(s); + + if(*s == 0) + goto Print; + + if(t = alpha(s, "flat")){ + if(*t){ + fprint(2, "flat takes no arguments\n"); + return; + } + if(nrec == 0){ + fprint(2, "no records\n"); + return; + } + flatten(); + goto Print; + } + + if(t = alpha(s, "print")){ + if(*t){ + fprint(2, "print takes no arguments\n"); + return; + } + for(i=0,r=rec; iNSAO) + goto NotFound; + t = skipbl(u); + if(*t){ + fprint(2, "syntax error in sao\n"); + return; + } + if(doreset) + reset(); + if(!loadsao(n)) + goto NotFound; + goto Print; + } + + if(t = alpha(s, "ngc")){ + n = strtoul(t, &u, 10); + if(n<=0 || n>NNGC) + goto NotFound; + t = skipbl(u); + if(*t){ + fprint(2, "syntax error in ngc\n"); + return; + } + if(doreset) + reset(); + if(!loadngc(n)) + goto NotFound; + goto Print; + } + + if(t = alpha(s, "ic")){ + n = strtoul(t, &u, 10); + if(n<=0 || n>NIC) + goto NotFound; + t = skipbl(u); + if(*t){ + fprint(2, "syntax error in ic\n"); + return; + } + if(doreset) + reset(); + if(!loadngc(n+NNGC)) + goto NotFound; + goto Print; + } + + if(t = alpha(s, "abell")){ + n = strtoul(t, &u, 10); + if(n<=0 || n>NAbell) + goto NotFound; + if(doreset) + reset(); + if(!loadabell(n)) + goto NotFound; + goto Print; + } + + if(t = alpha(s, "m")){ + n = strtoul(t, &u, 10); + if(n<=0 || n>NM) + goto NotFound; + mopen(); + for(j=n-1; mindex[j].mtype = NGCN; + cur->index = mindex[j].ngc; + } + j++; + } + goto Print; + } + + for(i=1; i<=Ncon; i++) + if(t = alpha(s, constel[i])){ + if(*t){ + fprint(2, "syntax error in constellation\n"); + return; + } + constelopen(); + seek(condb, 4L*conindex[i-1], 0); + j = conindex[i]-conindex[i-1]; + Eread(condb, "con", con, 4*j); + if(doreset) + reset(); + for(k=0; ktype = PatchC; + cur->index = Long(&con[k]); + } + goto Print; + } + + if(t = alpha(s, "expand")){ + n = 0; + if(*t){ + if(*t<'0' && '9'<*t){ + Expanderr: + fprint(2, "syntax error in expand\n"); + return; + } + n = strtoul(t, &u, 10); + t = skipbl(u); + if(*t) + goto Expanderr; + } + coords(n); + goto Print; + } + + if(t = alpha(s, "plot")){ + if(nrec == 0){ + Bprint(&bout, "empty\n"); + return; + } + plot(t); + return; + } + + if(t = alpha(s, "astro")){ + astro(t, 0); + return; + } + + if(t = alpha(s, "plate")){ + pplate(t); + return; + } + + if(t = alpha(s, "gamma")){ + while(*t==' ') + t++; + u = t; + x = strtod(t, &u); + if(u > t) + gam.gamma = x; + Bprint(&bout, "%.2f\n", gam.gamma); + return; + } + + if(t = alpha(s, "keep")){ + if(!cull(t, 1, 0)) + return; + goto Print; + } + + if(t = alpha(s, "drop")){ + if(!cull(t, 0, 0)) + return; + goto Print; + } + + for(i=0; planet[i].name[0]; i++){ + if(t = alpha(s, planet[i].name)){ + if(doreset) + reset(); + loadplanet(i, nil); + goto Print; + } + } + + for(i=0; names[i].name; i++){ + if(t = alpha(s, names[i].name)){ + if(*t){ + fprint(2, "syntax error in type\n"); + return; + } + if(doreset) + reset(); + loadtype(names[i].type); + goto Print; + } + } + + switch(s[0]){ + case '"': + starts = ++s; + while(*s != '"') + if(*s++ == 0){ + fprint(2, "bad star name\n"); + return; + } + *s = 0; + if(doreset) + reset(); + j = nrec; + saoopen(); + starts = fromgreek(starts); + for(i=0; i=90 || rah>=24) + goto BadCoords; + if(!loadpatch(patch(rah, ram, deg))) + goto NotFound; + break; + + default: + fprint(2, "unknown command %s\n", inputline); + return; + } + + Print: + if(nrec == 0) + Bprint(&bout, "empty\n"); + else if(nrec <= 2) + for(i=0; iPlateDefect) + return "can't happen"; + return ngctypes[d]; +} + +short descindex[NINDEX]; + +void +printnames(Record *r) +{ + int i, ok, done; + + done = 0; + for(i=0; itype==SAO && r->index==name[i].sao) + ok = 1; + if(r->type==NGC && r->ngc.ngc==name[i].ngc) + ok = 1; + if(r->type==Abell && r->abell.abell==name[i].abell) + ok = 1; + if(ok){ + if(done++ == 0) + Bprint(&bout, "\t"); + Bprint(&bout, " \"%s\"", togreek(name[i].name)); + } + } + if(done) + Bprint(&bout, "\n"); +} + +int +equal(char *s1, char *s2) +{ + int c; + + while(*s1){ + if(*s1==' '){ + while(*s1==' ') + s1++; + continue; + } + while(*s2==' ') + s2++; + c=*s2; + if('A'<=*s2 && *s2<='Z') + c^=' '; + if(*s1!=c) + return 0; + s1++, s2++; + } + return 1; +} + +int +parsename(char *s) +{ + char *blank; + int i; + + blank = strchr(s, ' '); + if(blank==0 || strchr(blank+1, ' ') || strlen(blank+1)!=3) + return 0; + blank++; + parsed[0] = parsed[1] = parsed[2] = 0; + if('0'<=s[0] && s[0]<='9'){ + i = atoi(s); + parsed[0] = i; + if(i > 100) + return 0; + }else{ + for(i=1; i<=24; i++) + if(strncmp(greek[i], s, strlen(greek[i]))==0){ + parsed[0]=100+i; + goto out; + } + return 0; + out: + if('0'<=s[strlen(greek[i])] && s[strlen(greek[i])]<='9') + parsed[1]=s[strlen(greek[i])]-'0'; + } + for(i=1; i<=88; i++) + if(strcmp(constel[i], blank)==0){ + parsed[2] = i; + return 1; + } + return 0; +} + +char* +dist_grp(int dg) +{ + switch(dg){ + default: + return "unknown"; + case 1: + return "13.3-14.0"; + case 2: + return "14.1-14.8"; + case 3: + return "14.9-15.6"; + case 4: + return "15.7-16.4"; + case 5: + return "16.5-17.2"; + case 6: + return "17.3-18.0"; + case 7: + return ">18.0"; + } +} + +char* +rich_grp(int dg) +{ + switch(dg){ + default: + return "unknown"; + case 0: + return "30-40"; + case 1: + return "50-79"; + case 2: + return "80-129"; + case 3: + return "130-199"; + case 4: + return "200-299"; + case 5: + return ">=300"; + } +} + +char* +nameof(Record *r) +{ + NGCrec *n; + SAOrec *s; + Abellrec *a; + static char buf[128]; + int i; + + switch(r->type){ + default: + return nil; + case SAO: + s = &r->sao; + if(s->name[0] == 0) + return nil; + if(s->name[0] >= 100){ + i = snprint(buf, sizeof buf, "%C", greeklet[s->name[0]-100]); + if(s->name[1]) + i += snprint(buf+i, sizeof buf-i, "%d", s->name[1]); + }else + i = snprint(buf, sizeof buf, " %d", s->name[0]); + snprint(buf+i, sizeof buf-i, " %s", constel[(uchar)s->name[2]]); + break; + case NGC: + n = &r->ngc; + if(n->type >= Uncertain) + return nil; + if(n->ngc <= NNGC) + snprint(buf, sizeof buf, "NGC%4d ", n->ngc); + else + snprint(buf, sizeof buf, "IC%4d ", n->ngc-NNGC); + break; + case Abell: + a = &r->abell; + snprint(buf, sizeof buf, "Abell%4d", a->abell); + break; + } + return buf; +} + +void +prrec(Record *r) +{ + NGCrec *n; + SAOrec *s; + Abellrec *a; + Planetrec *p; + int i, rah, ram, dec, nn; + long key; + + if(r) switch(r->type){ + default: + fprint(2, "can't prrec type %d\n", r->type); + exits("type"); + + case Planet: + p = &r->planet; + Bprint(&bout, "%s", p->name); + Bprint(&bout, "\t%s %s", + hms(angle(p->ra)), + dms(angle(p->dec))); + Bprint(&bout, " %3.2f° %3.2f°", + p->az/(double)MILLIARCSEC, p->alt/(double)MILLIARCSEC); + Bprint(&bout, " %s", + ms(angle(p->semidiam))); + if(r->index <= 1) + Bprint(&bout, " %g", p->phase); + Bprint(&bout, "\n"); + break; + + case NGC: + n = &r->ngc; + if(n->ngc <= NNGC) + Bprint(&bout, "NGC%4d ", n->ngc); + else + Bprint(&bout, "IC%4d ", n->ngc-NNGC); + Bprint(&bout, "%s ", ngcstring(n->type)); + if(n->mag == UNKNOWNMAG) + Bprint(&bout, "----"); + else + Bprint(&bout, "%.1f%c", n->mag/10.0, n->magtype); + Bprint(&bout, "\t%s %s\t%c%.1f'\n", + hm(angle(n->ra)), + dm(angle(n->dec)), + n->diamlim, + DEG(angle(n->diam))*60.); + prdesc(n->desc, desctab, descindex); + printnames(r); + break; + + case Abell: + a = &r->abell; + Bprint(&bout, "Abell%4d %.1f %.2f° %dMpc", a->abell, a->mag10/10.0, + DEG(angle(a->rad)), a->dist); + Bprint(&bout, "\t%s %s\t%.2f %.2f\n", + hm(angle(a->ra)), + dm(angle(a->dec)), + DEG(angle(a->glat)), + DEG(angle(a->glong))); + Bprint(&bout, "\tdist grp: %s rich grp: %s %d galaxies/°²\n", + dist_grp(a->distgrp), + rich_grp(a->richgrp), + a->pop); + printnames(r); + break; + + case SAO: + s = &r->sao; + Bprint(&bout, "SAO%6ld ", r->index); + if(s->mag==UNKNOWNMAG) + Bprint(&bout, "---"); + else + Bprint(&bout, "%.1f", s->mag/10.0); + if(s->mpg==UNKNOWNMAG) + Bprint(&bout, ",---"); + else + Bprint(&bout, ",%.1f", s->mpg/10.0); + Bprint(&bout, " %s %s %.4fs %.3f\"", + hms(angle(s->ra)), + dms(angle(s->dec)), + DEG(angle(s->dra))*(4*60), + DEG(angle(s->ddec))*(60*60)); + Bprint(&bout, " %.3s %c %.2s %ld %d", + s->spec, s->code, s->compid, s->hd, s->hdcode); + if(s->name[0]) + Bprint(&bout, " \"%s\"", nameof(r)); + Bprint(&bout, "\n"); + printnames(r); + break; + + case Patch: + radec(r->index, &rah, &ram, &dec); + Bprint(&bout, "%dh%dm %d°", rah, ram, dec); + key = r->patch.key[0]; + Bprint(&bout, " %s", constel[key&0xFF]); + if((key>>=8) & 0xFF) + Bprint(&bout, " %s", constel[key&0xFF]); + if((key>>=8) & 0xFF) + Bprint(&bout, " %s", constel[key&0xFF]); + if((key>>=8) & 0xFF) + Bprint(&bout, " %s", constel[key&0xFF]); + for(i=1; ipatch.nkey; i++){ + key = r->patch.key[i]; + switch(key&0x3F){ + case SAO: + Bprint(&bout, " SAO%ld", (key>>8)&0xFFFFFF); + break; + case Abell: + Bprint(&bout, " Abell%ld", (key>>8)&0xFFFFFF); + break; + default: /* NGC */ + nn = (key>>16)&0xFFFF; + if(nn > NNGC) + Bprint(&bout, " IC%d", nn-NNGC); + else + Bprint(&bout, " NGC%d", nn); + Bprint(&bout, "(%s)", ngcstring(key&0x3F)); + break; + } + } + Bprint(&bout, "\n"); + break; + + case NGCN: + if(r->index <= NNGC) + Bprint(&bout, "NGC%ld\n", r->index); + else + Bprint(&bout, "IC%ld\n", r->index-NNGC); + break; + + case NamedSAO: + Bprint(&bout, "SAO%ld \"%s\"\n", r->index, togreek(r->named.name)); + break; + + case NamedNGC: + if(r->index <= NNGC) + Bprint(&bout, "NGC%ld \"%s\"\n", r->index, togreek(r->named.name)); + else + Bprint(&bout, "IC%ld \"%s\"\n", r->index-NNGC, togreek(r->named.name)); + break; + + case NamedAbell: + Bprint(&bout, "Abell%ld \"%s\"\n", r->index, togreek(r->named.name)); + break; + + case PatchC: + radec(r->index, &rah, &ram, &dec); + Bprint(&bout, "%dh%dm %d\n", rah, ram, dec); + break; + } +} blob - /dev/null blob + 420a2a9a3ae39abdace6cd68c796ed71b4459114 (mode 644) --- /dev/null +++ src/cmd/scat/sky.h @@ -0,0 +1,413 @@ +#define DIR "#9/sky" +/* + * This code reflects many years of changes. There remain residues + * of prior implementations. + * + * Keys: + * 32 bits long. High 26 bits are encoded as described below. + * Low 6 bits are types: + * + * Patch is ~ one square degree of sky. It points to an otherwise + * anonymous list of Catalog keys. The 0th key is special: + * it contains up to 4 constellation identifiers. + * Catalogs (SAO,NGC,M,...) are: + * 31.........8|76|543210 + * catalog # |BB|catalog name + * BB is two bits of brightness: + * 00 -inf < m <= 7 + * 01 7 < m <= 10 + * 10 10 < m <= 13 + * 11 13 < m < inf + * The BB field is a dreg, and correct only for SAO and NGC. + * IC(n) is just NGC(n+7840) + * Others should be self-explanatory. + * + * Records: + * + * Star is an SAOrec + * Galaxy, PlanetaryN, OpenCl, GlobularCl, DiffuseN, etc., are NGCrecs. + * Abell is an Abellrec + * The Namedxxx records hold a name and a catalog entry; they result from + * name lookups. + */ + +typedef enum +{ + Planet, + Patch, + SAO, + NGC, + M, + Constel_deprecated, + Nonstar_deprecated, + NamedSAO, + NamedNGC, + NamedAbell, + Abell, + /* NGC types */ + Galaxy, + PlanetaryN, + OpenCl, + GlobularCl, + DiffuseN, + NebularCl, + Asterism, + Knot, + Triple, + Double, + Single, + Uncertain, + Nonexistent, + Unknown, + PlateDefect, + /* internal */ + NGCN, + PatchC, + NONGC, +}Type; + +enum +{ + /* + * parameters for plate + */ + Pppo1 = 0, + Pppo2, + Pppo3, + Pppo4, + Pppo5, + Pppo6, + Pamdx1, + Pamdx2, + Pamdx3, + Pamdx4, + Pamdx5, + Pamdx6, + Pamdx7, + Pamdx8, + Pamdx9, + Pamdx10, + Pamdx11, + Pamdx12, + Pamdx13, + Pamdx14, + Pamdx15, + Pamdx16, + Pamdx17, + Pamdx18, + Pamdx19, + Pamdx20, + Pamdy1, + Pamdy2, + Pamdy3, + Pamdy4, + Pamdy5, + Pamdy6, + Pamdy7, + Pamdy8, + Pamdy9, + Pamdy10, + Pamdy11, + Pamdy12, + Pamdy13, + Pamdy14, + Pamdy15, + Pamdy16, + Pamdy17, + Pamdy18, + Pamdy19, + Pamdy20, + Ppltscale, + Pxpixelsz, + Pypixelsz, + Ppltra, + Ppltrah, + Ppltram, + Ppltras, + Ppltdec, + Ppltdecd, + Ppltdecm, + Ppltdecs, + Pnparam, +}; + +#define UNKNOWNMAG 32767 +#define NPlanet 20 + +typedef float Angle; /* in radians */ +typedef long DAngle; /* on disk: in units of milliarcsec */ +typedef short Mag; /* multiplied by 10 */ +typedef long Key; /* known to be 4 bytes, unfortunately */ + +/* + * All integers are stored in little-endian order. + */ +typedef struct NGCrec NGCrec; +struct NGCrec{ + DAngle ra; + DAngle dec; + DAngle dummy1; /* compatibility with old RNGC version */ + DAngle diam; + Mag mag; + short ngc; /* if >NNGC, IC number is ngc-NNGC */ + char diamlim; + char type; + char magtype; + char dummy2; + char desc[52]; /* 0-terminated Dreyer description */ +}; + +typedef struct Abellrec Abellrec; +struct Abellrec{ + DAngle ra; + DAngle dec; + DAngle glat; + DAngle glong; + Mag mag10; /* mag of 10th brightest cluster member; in same place as ngc.mag*/ + short abell; + DAngle rad; + short pop; + short dist; + char distgrp; + char richgrp; + char flag; + char pad; +}; + +typedef struct Planetrec Planetrec; +struct Planetrec{ + DAngle ra; + DAngle dec; + DAngle az; + DAngle alt; + DAngle semidiam; + double phase; + char name[16]; +}; + +/* + * Star names: 0,0==unused. Numbers are name[0]=1,..,99. + * Greek letters are alpha=101, etc. + * Constellations are alphabetical order by abbreviation, and=1, etc. + */ +typedef struct SAOrec SAOrec; +struct SAOrec{ + DAngle ra; + DAngle dec; + DAngle dra; + DAngle ddec; + Mag mag; /* visual */ + Mag mpg; + char spec[3]; + char code; + char compid[2]; + char hdcode; + char pad1; + long hd; /* HD catalog number */ + char name[3]; /* name[0]=alpha name[1]=2 name[3]=ori */ + char nname; /* number of prose names */ + /* 36 bytes to here */ +}; + +typedef struct Mindexrec Mindexrec; +struct Mindexrec{ /* code knows the bit patterns in here; this is a long */ + char m; /* M number */ + char dummy; + short ngc; +}; + +typedef struct Bayerec Bayerec; +struct Bayerec{ + long sao; + char name[3]; + char pad; +}; + +/* + * Internal form + */ + +typedef struct Namedrec Namedrec; +struct Namedrec{ + char name[36]; +}; + +typedef struct Namerec Namerec; +struct Namerec{ + long sao; + long ngc; + long abell; + char name[36]; /* null terminated */ +}; + +typedef struct Patchrec Patchrec; +struct Patchrec{ + int nkey; + long key[60]; +}; + +typedef struct Record Record; +struct Record{ + Type type; + long index; + union{ + SAOrec sao; + NGCrec ngc; + Abellrec abell; + Namedrec named; + Patchrec patch; + Planetrec planet; + /* PatchCrec is empty */ + }; +}; + +typedef struct Name Name; +struct Name{ + char *name; + int type; +}; + +typedef struct Plate Plate; +struct Plate +{ + char rgn[7]; + char disk; + Angle ra; + Angle dec; +}; + +typedef struct Header Header; +struct Header +{ + float param[Pnparam]; + int amdflag; + + float x; + float y; + float xi; + float eta; +}; +typedef long Pix; + +typedef struct Img Img; +struct Img +{ + int nx; + int ny; /* ny is the fast-varying dimension */ + Pix a[1]; +}; + +#define RAD(x) ((x)*PI_180) +#define DEG(x) ((x)/PI_180) +#define ARCSECONDS_PER_RADIAN (DEG(1)*3600) +#define MILLIARCSEC (1000*60*60) + +int nplate; +Plate plate[2000]; /* needs to go to 2000 when the north comes */ +double PI_180; +double TWOPI; +double LN2; +int debug; +struct +{ + float min; + float max; + float gamma; + float absgamma; + float mult1; + float mult2; + int neg; +} gam; + +typedef struct Picture Picture; +struct Picture +{ + int minx; + int miny; + int maxx; + int maxy; + char name[16]; + uchar *data; +}; + +#ifndef _DRAW_H_ +typedef struct Image Image; +#endif + +extern double PI_180; +extern double TWOPI; +extern char *progname; +extern char *desctab[][2]; +extern Name names[]; +extern Record *rec; +extern long nrec; +extern Planetrec *planet; +/* for bbox: */ +extern int folded; +extern DAngle ramin; +extern DAngle ramax; +extern DAngle decmin; +extern DAngle decmax; +extern Biobuf bout; + +extern void saoopen(void); +extern void ngcopen(void); +extern void patchopen(void); +extern void mopen(void); +extern void constelopen(void); +extern void lowercase(char*); +extern void lookup(char*, int); +extern int typetab(int); +extern char*ngcstring(int); +extern char*skip(int, char*); +extern void prrec(Record*); +extern int equal(char*, char*); +extern int parsename(char*); +extern void radec(int, int*, int*, int*); +extern int btag(short); +extern long patcha(Angle, Angle); +extern long patch(int, int, int); +extern char*hms(Angle); +extern char*dms(Angle); +extern char*ms(Angle); +extern char*hm(Angle); +extern char*dm(Angle); +extern char*deg(Angle); +extern char*hm5(Angle); +extern long dangle(Angle); +extern Angle angle(DAngle); +extern void prdesc(char*, char*(*)[2], short*); +extern double xsqrt(double); +extern Angle dist(Angle, Angle, Angle, Angle); +extern Header* getheader(char*); +extern char* getword(char*, char*); +extern void amdinv(Header*, Angle, Angle, float, float); +extern void ppoinv(Header*, Angle, Angle); +extern void xypos(Header*, Angle, Angle, float, float); +extern void traneqstd(Header*, Angle, Angle); +extern Angle getra(char*); +extern Angle getdec(char*); +extern void getplates(void); +extern Img* dssread(char*); +extern void hinv(Pix*, int, int); +extern int input_bit(Biobuf*); +extern int input_nbits(Biobuf*, int); +extern int input_huffman(Biobuf*); +extern int input_nybble(Biobuf*); +extern void qtree_decode(Biobuf*, Pix*, int, int, int, int); +extern void start_inputing_bits(void); +extern Picture* image(Angle, Angle, Angle, Angle); +extern char* dssmount(int); +extern int dogamma(Pix); +extern void displaypic(Picture*); +extern void displayimage(Image*); +extern void plot(char*); +extern void astro(char*, int); +extern char* alpha(char*, char*); +extern char* skipbl(char*); +extern void flatten(void); +extern int bbox(long, long, int); +extern int inbbox(DAngle, DAngle); +extern char* nameof(Record*); + +#define NINDEX 400 blob - /dev/null blob + e60246c708c80361c7f3e5113f84a20cdce34bbc (mode 644) --- /dev/null +++ src/cmd/scat/strings.c @@ -0,0 +1,52 @@ +char *greek[]={ 0, /* 1-indexed */ + "alpha", "beta", "gamma", "delta", "epsilon", "zeta", "eta", "theta", + "iota", "kappa", "lambda", "mu", "nu", "xsi", "omicron", "pi", "rho", + "sigma", "tau", "upsilon", "phi", "chi", "psi", "omega", +}; + +Rune greeklet[]={ 0, + 0x3b1, 0x3b2, 0x3b3, 0x3b4, 0x3b5, 0x3b6, 0x3b7, 0x3b8, 0x3b9, 0x3ba, 0x3bb, + 0x3bc, 0x3bd, 0x3be, 0x3bf, 0x3c0, 0x3c1, 0x3c3, 0x3c4, 0x3c5, 0x3c6, 0x3c7, + 0x3c8, 0x3c9, +}; + +char *constel[]={ 0, /* 1-indexed */ + "and", "ant", "aps", "aql", "aqr", "ara", "ari", "aur", "boo", "cae", + "cam", "cap", "car", "cas", "cen", "cep", "cet", "cha", "cir", "cma", + "cmi", "cnc", "col", "com", "cra", "crb", "crt", "cru", "crv", "cvn", + "cyg", "del", "dor", "dra", "equ", "eri", "for", "gem", "gru", "her", + "hor", "hya", "hyi", "ind", "lac", "leo", "lep", "lib", "lmi", "lup", + "lyn", "lyr", "men", "mic", "mon", "mus", "nor", "oct", "oph", "ori", + "pav", "peg", "per", "phe", "pic", "psa", "psc", "pup", "pyx", "ret", + "scl", "sco", "sct", "ser", "sex", "sge", "sgr", "tau", "tel", "tra", + "tri", "tuc", "uma", "umi", "vel", "vir", "vol", "vul", +}; +Name names[]={ + "gx", Galaxy, + "pl", PlanetaryN, + "oc", OpenCl, + "gb", GlobularCl, + "nb", DiffuseN, + "c+n",NebularCl, + "ast", Asterism, + "kt", Knot, + "***", Triple, + "d*", Double, + "*", Single, + "pd", PlateDefect, + "galaxy", Galaxy, + "planetary", PlanetaryN, + "opencluster", OpenCl, + "globularcluster", GlobularCl, + "nebula", DiffuseN, + "nebularcluster",NebularCl, + "asterism", Asterism, + "knot", Knot, + "triple", Triple, + "double", Double, + "single", Single, + "nonexistent", Nonexistent, + "unknown", Unknown, + "platedefect", PlateDefect, + 0, 0, +}; blob - /dev/null blob + 33dca378b416dcef370d2f749468da764722abce (mode 644) --- /dev/null +++ src/cmd/scat/util.c @@ -0,0 +1,368 @@ +#include +#include +#include +#include "sky.h" + +double PI_180 = 0.0174532925199432957692369; +double TWOPI = 6.2831853071795864769252867665590057683943387987502; +double LN2 = 0.69314718055994530941723212145817656807550013436025; +static double angledangle=(180./PI)*MILLIARCSEC; + +#define rint scatrint + +int +rint(char *p, int n) +{ + int i=0; + + while(*p==' ' && n) + p++, --n; + while(n--) + i=i*10+*p++-'0'; + return i; +} + +DAngle +dangle(Angle angle) +{ + return angle*angledangle; +} + +Angle +angle(DAngle dangle) +{ + return dangle/angledangle; +} + +double +rfloat(char *p, int n) +{ + double i, d=0; + + while(*p==' ' && n) + p++, --n; + if(*p == '+') + return rfloat(p+1, n-1); + if(*p == '-') + return -rfloat(p+1, n-1); + while(*p == ' ' && n) + p++, --n; + if(n == 0) + return 0.0; + while(n-- && *p!='.') + d = d*10+*p++-'0'; + if(n <= 0) + return d; + p++; + i = 1; + while(n--) + d+=(*p++-'0')/(i*=10.); + return d; +} + +int +sign(int c) +{ + if(c=='-') + return -1; + return 1; +} + +char* +hms(Angle a) +{ + static char buf[20]; + double x; + int h, m, s, ts; + + x=DEG(a)/15; + x += 0.5/36000.; /* round up half of 0.1 sec */ + h = floor(x); + x -= h; + x *= 60; + m = floor(x); + x -= m; + x *= 60; + s = floor(x); + x -= s; + ts = 10*x; + sprint(buf, "%dh%.2dm%.2d.%ds", h, m, s, ts); + return buf; +} + +char* +dms(Angle a) +{ + static char buf[20]; + double x; + int sign, d, m, s, ts; + + x = DEG(a); + sign='+'; + if(a<0){ + sign='-'; + x=-x; + } + x += 0.5/36000.; /* round up half of 0.1 arcsecond */ + d = floor(x); + x -= d; + x *= 60; + m = floor(x); + x -= m; + x *= 60; + s = floor(x); + x -= s; + ts = floor(10*x); + sprint(buf, "%c%d°%.2d'%.2d.%d\"", sign, d, m, s, ts); + return buf; +} + +char* +ms(Angle a) +{ + static char buf[20]; + double x; + int d, m, s, ts; + + x = DEG(a); + x += 0.5/36000.; /* round up half of 0.1 arcsecond */ + d = floor(x); + x -= d; + x *= 60; + m = floor(x); + x -= m; + x *= 60; + s = floor(x); + x -= s; + ts = floor(10*x); + if(d != 0) + sprint(buf, "%d°%.2d'%.2d.%d\"", d, m, s, ts); + else + sprint(buf, "%.2d'%.2d.%d\"", m, s, ts); + return buf; +} + +char* +hm(Angle a) +{ + static char buf[20]; + double x; + int h, m, n; + + x = DEG(a)/15; + x += 0.5/600.; /* round up half of tenth of minute */ + h = floor(x); + x -= h; + x *= 60; + m = floor(x); + x -= m; + x *= 10; + n = floor(x); + sprint(buf, "%dh%.2d.%1dm", h, m, n); + return buf; +} + +char* +hm5(Angle a) +{ + static char buf[20]; + double x; + int h, m; + + x = DEG(a)/15; + x += 2.5/60.; /* round up 2.5m */ + h = floor(x); + x -= h; + x *= 60; + m = floor(x); + m -= m % 5; + sprint(buf, "%dh%.2dm", h, m); + return buf; +} + +char* +dm(Angle a) +{ + static char buf[20]; + double x; + int sign, d, m, n; + + x = DEG(a); + sign='+'; + if(a<0){ + sign='-'; + x=-x; + } + x += 0.5/600.; /* round up half of tenth of arcminute */ + d = floor(x); + x -= d; + x *= 60; + m = floor(x); + x -= m; + x *= 10; + n = floor(x); + sprint(buf, "%c%d°%.2d.%.1d'", sign, d, m, n); + return buf; +} + +char* +deg(Angle a) +{ + static char buf[20]; + double x; + int sign, d; + + x = DEG(a); + sign='+'; + if(a<0){ + sign='-'; + x=-x; + } + x += 0.5; /* round up half degree */ + d = floor(x); + sprint(buf, "%c%d°", sign, d); + return buf; +} + +char* +getword(char *ou, char *in) +{ + int c; + + for(;;) { + c = *in++; + if(c == ' ' || c == '\t') + continue; + if(c == 0) + return 0; + break; + } + + if(c == '\'') + for(;;) { + if(c >= 'A' && c <= 'Z') + c += 'a' - 'A'; + *ou++ = c; + c = *in++; + if(c == 0) + return 0; + if(c == '\'') { + *ou = 0; + return in-1; + } + } + for(;;) { + if(c >= 'A' && c <= 'Z') + c += 'a' - 'A'; + *ou++ = c; + c = *in++; + if(c == ' ' || c == '\t' || c == 0) { + *ou = 0; + return in-1; + } + } +} + +/* + * Read formatted angle. Must contain no embedded blanks + */ +Angle +getra(char *p) +{ + Rune r; + char *q; + Angle f, d; + int neg; + + neg = 0; + d = 0; + while(*p == ' ') + p++; + for(;;) { + if(*p == ' ' || *p=='\0') + goto Return; + if(*p == '-') { + neg = 1; + p++; + } + if(*p == '+') { + neg = 0; + p++; + } + q = p; + f = strtod(p, &q); + if(q > p) { + p = q; + } + p += chartorune(&r, p); + switch(r) { + default: + Return: + if(neg) + d = -d; + return RAD(d); + case 'h': + d += f*15; + break; + case 'm': + d += f/4; + break; + case 's': + d += f/240; + break; + case 0xB0: /* ° */ + d += f; + break; + case '\'': + d += f/60; + break; + case '\"': + d += f/3600; + break; + } + } + return 0; +} + +double +xsqrt(double a) +{ + + if(a < 0) + return 0; + return sqrt(a); +} + +Angle +dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2) +{ + double a; + + a = sin(dec1) * sin(dec2) + + cos(dec1) * cos(dec2) * + cos(ra1 - ra2); + a = atan2(xsqrt(1 - a*a), a); + if(a < 0) + a = -a; + return a; +} + +int +dogamma(Pix c) +{ + float f; + + f = c - gam.min; + if(f < 1) + f = 1; + + if(gam.absgamma == 1) + c = f * gam.mult2; + else + c = exp(log(f*gam.mult1) * gam.absgamma) * 255; + if(c > 255) + c = 255; + if(gam.neg) + c = 255-c; + return c; +}