commit - 3f8c70e97c2eb85992424439af56a4dd6412b8c6
commit + 8a3b2ceb0ff632c47e1516d3ffef8572dc8eb974
blob - /dev/null
blob + d9bf8db805482943fa01379ad0482fe0e13cee5a (mode 644)
--- /dev/null
+++ man/man1/scat.1
+.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
+wget -O- http://pdos.lcs.mit.edu/~rsc/scat.tgz | gunzip | tar xf -
blob - 258ca9a64a4314491d8829b8c7489cacdabafab8
blob + 7ebbe368406352a51d9bab661bdd1c4ebff5fabf
--- sky/here
+++ sky/here
-41.0343 73.3936 200
+40.6843333333 74.3996666667 150
blob - 85b20c37f05999acbababc61904b1971e240e2e9
blob + 0f27676d33cf13060cea967946347b4d548e7b41
--- src/cmd/rc/simple.c
+++ src/cmd/rc/simple.c
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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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
+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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include <draw.h>
+#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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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<aend; a++) {
+ if(px = *a) {
+ if(mask == 0) {
+ mask = 0x80;
+ bits = Bgetc(infile);
+ }
+ if(mask & bits)
+ *a = -px;
+ mask >>= 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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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); i++) {
+ j++;
+ if(j >= 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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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<<log2n))
+ log2n++;
+
+ /*
+ * get temporary storage for shuffling elements
+ */
+ tmp = (Pix*)malloc(((nmax+1)/2) * sizeof(*tmp));
+ if(tmp == nil) {
+ fprint(2, "hinv: insufficient memory\n");
+ exits("memory");
+ }
+
+ /*
+ * do log2n expansions
+ *
+ * We're indexing a as a 2-D array with dimensions (nx,ny).
+ */
+ shift = 1;
+ nxtop = 1;
+ nytop = 1;
+ nxf = nx;
+ nyf = ny;
+ c = 1<<log2n;
+ for(k = log2n-1; k>=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<nxtop; i++)
+ unshuffle1(&a[ny*i], nytop, tmp);
+ for(j = 0; j<nytop; j++)
+ unshuffle(&a[j], nxtop, ny, tmp);
+ oddx = nxtop & 1;
+ oddy = nytop & 1;
+ for(i = 0; i<nxtop-oddx; i += 2) {
+ s00 = ny*i; /* s00 is index of a[i,j] */
+ s10 = s00+ny; /* s10 is index of a[i+1,j] */
+ for(j = 0; j<nytop-oddy; j += 2) {
+ /*
+ * Multiply h0,hx,hy,hc by 2 (1 the last time through).
+ */
+ h0 = a[s00 ] << shift;
+ hx = a[s10 ] << shift;
+ hy = a[s00+1] << shift;
+ hc = a[s10+1] << shift;
+
+ /*
+ * Divide sums by 4 (shift right 2 bits).
+ * Add 1 to round -- note that these values are always
+ * positive so we don't need to do anything special
+ * for rounding negative numbers.
+ */
+ a[s10+1] = (h0 + hx + hy + hc + 2) >> 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<nytop-oddy; j += 2) {
+ h0 = a[s00 ] << shift;
+ hy = a[s00+1] << shift;
+ a[s00+1] = (h0 + hy + 2) >> 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<n; i++) {
+ *pt = *p1;
+ pt++;
+ p1 += n2;
+ }
+
+ /*
+ * distribute 1st half of array to even elements
+ */
+ p2 = &a[n2xnhalf]; /* pointer to a[i] */
+ p1 = &a[n2xnhalf<<1]; /* pointer to a[2*i] */
+ for(i=nhalf-1; 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<n; i+=2) {
+ *p1 = *pt;
+ p1 += twon2;
+ pt++;
+ }
+}
+
+static
+void
+unshuffle1(Pix *a, int n, Pix *tmp)
+{
+ int i;
+ int nhalf;
+ Pix *p1, *p2, *pt;
+
+ nhalf = (n+1) >> 1;
+
+ /*
+ * copy 2nd half of array to tmp
+ */
+ pt = tmp;
+ p1 = &a[nhalf]; /* pointer to a[i] */
+ for(i=nhalf; i<n; i++) {
+ *pt = *p1;
+ pt++;
+ p1++;
+ }
+
+ /*
+ * distribute 1st half of array to even elements
+ */
+ p2 = &a[nhalf]; /* pointer to a[i] */
+ p1 = &a[nhalf<<1]; /* pointer to a[2*i] */
+ for(i=nhalf-1; 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<n; i+=2) {
+ *p1 = *pt;
+ p1 += 2;
+ pt++;
+ }
+}
blob - /dev/null
blob + e61f01068e41f79eefdac9324dd2e9e5818318ec (mode 644)
--- /dev/null
+++ src/cmd/scat/image.c
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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; i<nplate; i++) {
+ pp = &plate[i];
+ d = dist(ra, dec, pp->ra, 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; j<ip->ny; 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; i<ip->nx; 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
+<$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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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
+#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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include <draw.h>
+#include <event.h>
+#include <ctype.h>
+#include "map.h"
+#undef RAD
+#undef TWOPI
+#include "sky.h"
+
+ /* convert to milliarcsec */
+static long c = MILLIARCSEC; /* 1 degree */
+static long m5 = 1250*60*60; /* 5 minutes of ra */
+
+DAngle ramin;
+DAngle ramax;
+DAngle decmin;
+DAngle decmax;
+int folded;
+
+Image *grey;
+Image *alphagrey;
+Image *green;
+Image *lightblue;
+Image *lightgrey;
+Image *ocstipple;
+Image *suncolor;
+Image *mooncolor;
+Image *shadowcolor;
+Image *mercurycolor;
+Image *venuscolor;
+Image *marscolor;
+Image *jupitercolor;
+Image *saturncolor;
+Image *uranuscolor;
+Image *neptunecolor;
+Image *plutocolor;
+Image *cometcolor;
+
+Planetrec *planet;
+
+long mapx0, mapy0;
+long mapra, mapdec;
+double mylat, mylon, mysid;
+double mapscale;
+double maps;
+int (*projection)(struct place*, double*, double*);
+
+char *fontname = "/lib/font/bit/lucida/unicode.6.font";
+
+/* types Coord and Loc correspond to types in map(3) thus:
+ Coord == struct coord;
+ Loc == struct place, except longitudes are measured
+ positive east for Loc and west for struct place
+*/
+
+typedef struct Xyz Xyz;
+typedef struct coord Coord;
+typedef struct Loc Loc;
+
+struct Xyz{
+ double x,y,z;
+};
+
+struct Loc{
+ Coord nlat, elon;
+};
+
+Xyz north = { 0, 0, 1 };
+
+int
+plotopen(void)
+{
+ if(display != nil)
+ return 1;
+ 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; i<nrec-nlast; i++,r++)
+ if(r->type == Planet)
+ if(name==nil || strcmp(r->planet.name, name)==0){
+ rr = *r;
+ memmove(rec+i, rec+i+1, (nrec-i-1)*sizeof(Record));
+ rec[nrec-1] = rr;
+ --i;
+ --r;
+ nlast++;
+ }
+}
+
+int
+bbox(long extrara, long extradec, int quantize)
+{
+ int i;
+ Record *r;
+ int ra, dec;
+ int rah, ram, d1, d2;
+ double r0;
+
+ ramin = 0x7FFFFFFF;
+ ramax = -0x7FFFFFFF;
+ decmin = 0x7FFFFFFF;
+ decmax = -0x7FFFFFFF;
+
+ for(i=0,r=rec; i<nrec; i++,r++){
+ if(r->type == Patch){
+ radec(r->index, &rah, &ram, &dec);
+ ra = 15*rah+ram/4;
+ r0 = c/cos(dec*PI/180);
+ 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; j<nelem(pts); j++){
+ /* use double to avoid overflow */
+ v = (double)j / (double)(nelem(pts)-1);
+ v = decmin + v*(decmax-decmin);
+ pts[j] = map(x, v);
+ }
+ bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
+ ra = x;
+ if(folded){
+ ra -= 180*c;
+ if(ra < 0)
+ ra += 360*c;
+ }
+ p = pts[0];
+ p.x -= stringwidth(font, hm5(angle(ra)))/2;
+ string(scr, p, GREY, ZP, font, hm5(angle(ra)));
+ p = pts[nelem(pts)-1];
+ p.x -= stringwidth(font, hm5(angle(ra)))/2;
+ p.y -= font->height;
+ 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; j<nelem(pts); j++){
+ /* use double to avoid overflow */
+ v = (double)j / (double)(nelem(pts)-1);
+ v = ramin + v*(ramax-ramin);
+ pts[j] = map(v, y);
+ }
+ bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
+ p = pts[0];
+ p.x += 3;
+ p.y -= font->height/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; i<nrec; i++,r++){
+ dec = r->ngc.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<np; i++)
+ parseplanet(lines[i], &planet[i-1]);
+}
blob - /dev/null
blob + 5cce14d7cff8984cf1bb379014ed6b242e9dbe26 (mode 644)
--- /dev/null
+++ src/cmd/scat/posn.c
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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<<log2n))
+ log2n++;
+
+ /*
+ * allocate scratch array for working space
+ */
+ nqx2 = (nqx+1)/2;
+ nqy2 = (nqy+1)/2;
+ scratch = (uchar*)malloc(nqx2*nqy2);
+ if(scratch == nil) {
+ fprint(2, "qtree_decode: insufficient memory\n");
+ exits("memory");
+ }
+
+ /*
+ * now decode each bit plane, starting at the top
+ * A is assumed to be initialized to zero
+ */
+ for(bit = nbitplanes-1; bit >= 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<<log2n;
+ for(k = 1; k<log2n; k++) {
+ /*
+ * this somewhat cryptic code generates the sequence
+ * n[k-1] = (n[k]+1)/2 where n[log2n]=nqx or nqy
+ */
+ c = c>>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<nx-1; i += 2) {
+ s00 = n*i; /* s00 is index of b[i,j] */
+ s10 = s00+n; /* s10 is index of b[i+1,j] */
+ for(j = 0; j<ny-1; j += 2) {
+ b[s10+1] = b[s00] & 1;
+ b[s10 ] = (b[s00]>>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<ny-1; j += 2) {
+ b[s00+1] = (b[s00]>>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<nx-1; i+=2) {
+ s00 = &b[n*i]; /* s00 is index of b[i,j] */
+ s10 = s00+n; /* s10 is index of b[i+1,j] */
+ for(j=0; j<ny-1; j+=2) {
+ px = *a++;
+ s10[1] |= ( px & 1) << bit;
+ s10[0] |= ((px>>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<ny-1; j+=2) {
+ px = *a++;
+ s00[1] |= ((px>>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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include <draw.h>
+#include <event.h>
+#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; i<NName; i++){
+ l = Brdline(&b, '\n');
+ if(l == 0)
+ break;
+ p = strchr(l, '\t');
+ if(p == 0){
+ Badformat:
+ Bprint(&bout, "warning: name.scat bad format; line %d\n", i+1);
+ break;
+ }
+ *p++ = 0;
+ strcpy(name[i].name, l);
+ if(strncmp(p, "ngc", 3) == 0)
+ name[i].ngc = atoi(p+3);
+ else if(strncmp(p, "ic", 2) == 0)
+ name[i].ngc = atoi(p+2)+NNGC;
+ else if(strncmp(p, "sao", 3) == 0)
+ name[i].sao = atoi(p+3);
+ else if(strncmp(p, "abell", 5) == 0)
+ name[i].abell = atoi(p+5);
+ else
+ goto Badformat;
+ }
+ if(i == NName)
+ Bprint(&bout, "warning: too many names in name.scat (max %d); extra ignored\n", NName);
+ close(namedb);
+
+ bayerdb = eopen("bayer");
+ Eread(bayerdb, "bayer", bayer, sizeof bayer);
+ close(bayerdb);
+ for(i=0; i<NBayer; i++)
+ bayer[i].sao = Long(&bayer[i].sao);
+ }
+}
+
+void
+saoopen(void)
+{
+ if(saodb == 0){
+ nameopen();
+ saodb = eopen("sao");
+ }
+}
+
+void
+ngcopen(void)
+{
+ if(ngcdb == 0){
+ nameopen();
+ ngcdb = eopen("ngc2000");
+ ngctypedb = eopen("ngc2000type");
+ Eread(ngctypedb, "ngctype", ngctype, sizeof ngctype);
+ close(ngctypedb);
+ }
+}
+
+void
+abellopen(void)
+{
+ /* nothing extra to do with abell: it's directly indexed by number */
+ if(abelldb == 0)
+ abelldb = eopen("abell");
+}
+
+void
+patchopen(void)
+{
+ Biobuf *b;
+ long l, m;
+ char buf[100];
+
+ if(patchdb == 0){
+ patchdb = eopen("patch");
+ sprint(buf, "%s/patchindex.scat", dir);
+ b = Bopen(buf, OREAD);
+ if(b == 0){
+ fprint(2, "can't open %s\n", buf);
+ exits("open");
+ }
+ for(m=0,l=0; l<=Npatch; l++)
+ patchaddr[l] = m += Bgetc(b)*4;
+ Bterm(b);
+ }
+}
+
+void
+mopen(void)
+{
+ int i;
+
+ if(mindexdb == 0){
+ mindexdb = eopen("mindex");
+ Eread(mindexdb, "mindex", mindex, sizeof mindex);
+ close(mindexdb);
+ for(i=0; i<NMrec; i++)
+ mindex[i].ngc = Short(&mindex[i].ngc);
+ }
+}
+
+void
+constelopen(void)
+{
+ int i;
+
+ if(condb == 0){
+ condb = eopen("con");
+ conindexdb = eopen("conindex");
+ Eread(conindexdb, "conindex", conindex, sizeof conindex);
+ close(conindexdb);
+ for(i=0; i<Ncon+1; i++)
+ conindex[i] = Short((short*)&conindex[i]);
+ }
+}
+
+void
+lowercase(char *s)
+{
+ for(; *s; s++)
+ if('A'<=*s && *s<='Z')
+ *s += 'a'-'A';
+}
+
+int
+loadngc(long index)
+{
+ static int failed;
+ long j;
+
+ ngcopen();
+ j = (index-1)*sizeof(NGCrec);
+ grow();
+ cur->type = 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; i<cur->patch.nkey; i++)
+ cur->patch.key[i] = Long(&cur->patch.key[i]);
+ return 1;
+}
+
+int
+loadtype(int t)
+{
+ int i;
+
+ ngcopen();
+ for(i=0; i<NNGCrec; i++)
+ if(t == (ngctype[i])){
+ grow();
+ cur->type = 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; i<norec; i++,or++){
+ switch(or->type){
+ 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; j<or->patch.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<NMrec; i++)
+ if(mindex[i].ngc == index)
+ return 1;
+ return 0;
+}
+
+char*
+alpha(char *s, char *t)
+{
+ int n;
+
+ n = strlen(t);
+ if(strncmp(s, t, n)==0 && (s[n]<'a' || 'z'<s[n]))
+ return skipbl(s+n);
+ return 0;
+
+}
+
+char*
+text(char *s, char *t)
+{
+ int n;
+
+ n = strlen(t);
+ if(strncmp(s, t, n)==0 && (s[n]==0 || s[n]==' ' || s[n]=='\t'))
+ return skipbl(s+n);
+ return 0;
+
+}
+
+int
+cull(char *s, int keep, int dobbox)
+{
+ int i, j, nobj, keepthis;
+ Record *or;
+ char *t;
+ int dogrtr, doless, dom, dosao, dongc, doabell;
+ int mgrtr, mless;
+ char obj[100];
+
+ memset(obj, 0, sizeof(obj));
+ nobj = 0;
+ dogrtr = 0;
+ doless = 0;
+ dom = 0;
+ dongc = 0;
+ dosao = 0;
+ doabell = 0;
+ mgrtr = mless= 0;
+ if(dobbox)
+ goto Cull;
+ for(;;){
+ if(s[0] == '>'){
+ 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; i<norec; i++,or++){
+ keepthis = !keep;
+ if(dobbox && inbbox(or->ngc.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; j<nobj; j++)
+ if(or->type==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; i<nrec; i++,r++){
+ /* may have multiple instances of a planet in the scene */
+ if(r->type==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; i<n; i++)
+ *t++ = *s++;
+ Cont:;
+ }
+ *t = 0;
+ return greekbuf;
+}
+
+char*
+fromgreek(char *s)
+{
+ char *t;
+ int i, n;
+ Rune r;
+
+ t = greekbuf;
+ while(*s){
+ n = chartorune(&r, s);
+ for(i=1; i<=24; i++){
+ if(r == greeklet[i]){
+ strcpy(t, greek[i]);
+ t += strlen(greek[i]);
+ s += n;
+ goto Cont;
+ }
+ }
+ for(i=0; i<n; i++)
+ *t++ = *s++;
+ Cont:;
+ }
+ *t = 0;
+ return greekbuf;
+}
+
+#ifdef OLD
+/*
+ * Old version
+ */
+int
+coords(int deg)
+{
+ int i;
+ int x, y;
+ Record *or;
+ long dec, ra, ndec, nra;
+ int rdeg;
+
+ flatten();
+ copy();
+ reset();
+ deg *= 2;
+ for(i=0,or=orec; i<norec; i++,or++){
+ if(or->type == 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; i<norec; i++,or++)
+ if(or->type == 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; i<nrec; i++,r++){
+ if(r->type == 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; i<nrec; i++,r++)
+ prrec(r);
+ return;
+ }
+
+ if(t = alpha(s, "add")){
+ lookup(t, 0);
+ return;
+ }
+
+ if(t = alpha(s, "sao")){
+ n = strtoul(t, &u, 10);
+ if(n<=0 || n>NSAO)
+ 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].m<n; j++)
+ ;
+ if(doreset)
+ reset();
+ while(mindex[j].m == n){
+ if(mindex[j].ngc){
+ grow();
+ cur->type = 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; k<j; k++){
+ grow();
+ cur->type = 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<NName; i++)
+ if(equal(starts, name[i].name)){
+ grow();
+ if(name[i].sao){
+ rec[j].type = NamedSAO;
+ rec[j].index = name[i].sao;
+ }
+ if(name[i].ngc){
+ rec[j].type = NamedNGC;
+ rec[j].index = name[i].ngc;
+ }
+ if(name[i].abell){
+ rec[j].type = NamedAbell;
+ rec[j].index = name[i].abell;
+ }
+ strcpy(rec[j].named.name, name[i].name);
+ j++;
+ }
+ if(parsename(starts))
+ for(i=0; i<NBayer; i++)
+ if(bayer[i].name[0]==parsed[0] &&
+ (bayer[i].name[1]==parsed[1] || parsed[1]==0) &&
+ bayer[i].name[2]==parsed[2]){
+ grow();
+ rec[j].type = NamedSAO;
+ rec[j].index = bayer[i].sao;
+ strncpy(rec[j].named.name, starts, sizeof(rec[j].named.name));
+ j++;
+ }
+ if(j == 0){
+ *s = '"';
+ goto NotFound;
+ }
+ break;
+
+ case '0': case '1': case '2': case '3': case '4':
+ case '5': case '6': case '7': case '8': case '9':
+ strtoul(s, &t, 10);
+ if(*t != 'h'){
+ BadCoords:
+ fprint(2, "bad coordinates %s\n", inputline);
+ break;
+ }
+ ra = DEG(getra(s));
+ while(*s && *s!=' ' && *s!='\t')
+ s++;
+ rah = ra/15;
+ ra = ra-rah*15;
+ ram = ra*4;
+ deg = strtol(s, &t, 10);
+ if(t == s)
+ goto BadCoords;
+ /* degree sign etc. is optional */
+ chartorune(&c, t);
+ if(c == 0xb0)
+ deg = DEG(getra(s));
+ if(doreset)
+ reset();
+ if(abs(deg)>=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; i<nrec; i++)
+ prrec(rec+i);
+ else
+ Bprint(&bout, "%ld items\n", nrec);
+ return;
+
+ NotFound:
+ fprint(2, "%s not found\n", inputline);
+ return;
+}
+
+char *ngctypes[] =
+{
+[Galaxy] "Gx",
+[PlanetaryN] "Pl",
+[OpenCl] "OC",
+[GlobularCl] "Gb",
+[DiffuseN] "Nb",
+[NebularCl] "C+N",
+[Asterism] "Ast",
+[Knot] "Kt",
+[Triple] "***",
+[Double] "D*",
+[Single] "*",
+[Uncertain] "?",
+[Nonexistent] "-",
+[Unknown] " ",
+[PlateDefect] "PD",
+};
+
+char*
+ngcstring(int d)
+{
+ if(d<Galaxy || d>PlateDefect)
+ return "can't happen";
+ return ngctypes[d];
+}
+
+short descindex[NINDEX];
+
+void
+printnames(Record *r)
+{
+ int i, ok, done;
+
+ done = 0;
+ for(i=0; i<NName; i++){ /* stupid linear search! */
+ ok = 0;
+ if(r->type==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; i<r->patch.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
+#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
+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
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#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;
+}
