Add scat. Temporary fix to rc r.e. note groups.

This commit is contained in:
rsc 2004-04-24 17:05:43 +00:00
parent 3f8c70e97c
commit 8a3b2ceb0f
22 changed files with 5997 additions and 2 deletions

View file

@ -63,7 +63,7 @@ void Xsimple(void){
Xerror("try again");
return;
case 0:
rfork(RFNOTEG);
// rfork(RFNOTEG);
pushword("exec");
execexec();
strcpy(buf, "can't exec: ");

76
src/cmd/scat/bitinput.c Normal file
View file

@ -0,0 +1,76 @@
#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;
}

327
src/cmd/scat/desc.c Normal file
View file

@ -0,0 +1,327 @@
char *desctab[][2]={
"!!!", "(magnificent or otherwise interesting object)",
"!!", "(superb)",
"!", "(remarkable)",
"1st", "first",
"2nd", "second",
"3rd", "third",
"4th", "fourth",
"5th", "fifth",
"6th", "sixth",
"7th", "seventh",
"8th", "eighth",
"9th", "ninth",
"B", "bright",
"Borealis", "Borealis",
"C", "compressed",
"Car", "Car",
"Cas", "Cas",
"Cen", "Cen",
"Cl", "cluster",
"Cyg", "Cyg",
"D", "double",
"Dra", "Dra",
"Dumbbell", "Dumbbell",
"E", "extended",
"F", "faint",
"L", "large",
"Lyra", "Lyra",
"M", "(in the) middle",
"Merope", "Merope",
"Milky", "Milky",
"Mon", "Mon",
"N", "(to a) nucleus",
"Neb", "Nebula",
"Nor", "Nor",
"Nucl", "nucleus",
"Nuclei", "nuclei",
"P", "poor in stars",
"PN", "planetary nebula",
"Polarissima", "Polarissima",
"Praesepe", "Praesepe",
"Psc", "Psc",
"R", "round",
"RA", "RA",
"RR", "exactly round",
"Ri", "rich in stars",
"S", "small",
"Sco", "Sco",
"S*", "small star",
"Way", "Way",
"ab", "about",
"about", "about",
"alm", "almost",
"alpha", "α",
"am", "among",
"annul", "annular",
"att", "attached",
"b", "brighter",
"beautiful", "beautiful",
"bet", "between",
"beta", "β",
"bf", "brightest to f side",
"biN", "binuclear",
"bifid", "bifid",
"bifurcated", "bifurcated",
"black", "black",
"blue", "blue",
"bn", "brightest to n side",
"border", "border",
"bp", "brightest to p side",
"branch", "branch",
"branched", "branched",
"branches", "branches",
"bright", "bright",
"brighter", "brighter",
"brightest", "brightest",
"brightness", "brightness",
"brush", "brush",
"bs", "brightest to s side",
"but", "but",
"by", "by",
"c", "considerably",
"centre", "centre",
"certain", "certain",
"chev", "chevelure",
"chief", "chief",
"chiefly", "chiefly",
"circle", "circle",
"close", "close",
"cloud", "cloud",
"cluster", "cluster",
"clusters", "clusters",
"co", "coarse, coarsely",
"com", "cometic",
"comet", "comet",
"cometary", "cometary",
"comp", "companion",
"condens", "condensations",
"condensed", "condensed",
"conn", "connected",
"connected", "connected",
"connecting", "connecting",
"cont", "in contact",
"corner", "corner",
"curved", "curved",
"d", "diameter",
"decl", "declination",
"def", "defined",
"defect", "defect",
"deg", "°",
"delta", "δ",
"dense", "dense",
"densest", "densest",
"descr", "description",
"description", "description",
"dif", "diffuse",
"different", "different",
"diffic", "difficult",
"difficult", "difficult",
"diffuse", "diffuse",
"diffused", "diffused",
"disc", "disc",
"dist", "distant",
"distant", "distant",
"distinct", "distinct",
"double", "double",
"doubtful", "doubtful",
"dozen", "dozen",
"e", "extremely",
"each", "each",
"edge", "edge",
"ee", "most extremely",
"ellipt", "elliptical",
"elliptic", "elliptical",
"end", "end",
"ends", "ends",
"epsilon", "ε",
"equilateral", "equilateral",
"er", "easily resolvable",
"eta", "η",
"evidently", "evidently",
"exc", "excentric",
"excen", "excentric",
"excent", "excentric",
"excentric", "excentric",
"extends", "extends",
"f", "following",
"faint", "faint",
"fainter", "fainter",
"faintest", "faintest",
"falcate", "falcate",
"fan", "fan",
"farther", "farther",
"field", "field",
"fine", "fine",
"forming", "forming",
"forms", "forms",
"found", "found",
"from", "from",
"g", "gradually",
"gamma", "γ",
"gaseous", "gaseous",
"gl", "gradually a little",
"glob. cl.", "globular cluster",
"gr", "group",
"great", "great",
"greater", "greater",
"group", "group",
"groups", "groups",
"i", "irregular",
"iF", "irregular figure",
"if", "if",
"in", "in",
"indistinct", "indistinct",
"incl", "including",
"inv", "involved",
"iota", "ι",
"irr", "irregular",
"is", "is",
"it", "it",
"kappa", "κ",
"l", "little, long",
"lC", "little compressed",
"lE", "little extended",
"lambda", "λ",
"larger", "larger",
"last", "last",
"lb", "little brighter",
"least", "least",
"like", "like",
"line", "in a line",
"little", "little",
"long", "long",
"looks", "looks",
"looped", "looped",
"m", "magnitude",
"mE", "much extended",
"mag", "mag",
"makes", "makes",
"many", "many",
"mb", "much brighter",
"more", "more",
"mottled", "mottled",
"mu", "μ",
"mult", "multiple",
"n", "north",
"narrow", "narrow",
"near", "near",
"nearly", "nearly",
"neb", "nebula",
"nebs", "nebulous",
"nebula", "nebula",
"nebulosity", "nebulosity",
"nebulous", "nebulous",
"neby", "nebulosity",
"nf", "north following",
"no", "no",
"nonexistent", "nonexistent",
"not", "not",
"np", "north preceding",
"nr", "near",
"ns", "north-south",
"nu", "ν",
"omega", "ω",
"p", "preceding",
"pB", "pretty bright",
"pC", "pretty compressed",
"pF", "pretty faint",
"pL", "pretty large",
"pR", "pretty round",
"pS", "pretty small",
"parallel", "parallel",
"part", "part",
"partly", "partly",
"patch", "patch",
"patches", "patches",
"perhaps", "perhaps",
"perpendicular", "perpendicular",
"pf", "preceding-following",
"pg", "pretty gradually",
"photosphere", "photosphere",
"pi", "π",
"place", "place",
"plate", "plate",
"plan", "planetary nebula",
"pointed", "pointed",
"portion", "portion",
"pos", "position angle",
"possibly", "possibly",
"prob", "probably",
"probably", "probably",
"ps", "pretty suddenly",
"r", "mottled",
"requires", "requires",
"resolved", "resolved",
"rho", "ρ",
"ring", "ring",
"rough", "rough",
"rr", "some stars seen",
"rrr", "clearly consisting of stars",
"ruby", "ruby",
"s", "south",
"same", "same",
"sb", "suddenly brighter",
"sc", "scattered",
"second", "second",
"seems", "seems",
"seen", "seen",
"segment", "segment",
"semi", "semi",
"sev", "several",
"several", "several",
"sf", "south following",
"shape", "shape",
"shaped", "shaped",
"sharp", "sharp",
"sigma", "σ",
"sl", "suddenly a little",
"slightly", "slightly",
"small", "small",
"south", "south",
"sp", "south preceding",
"spectrum", "spectrum",
"spindle", "spindle",
"spir", "spiral",
"spiral", "spiral",
"st 9...", "stars of mag. 9 and fainter",
"st 9...13", "stars of mag. 9 to 13",
"st", "stars",
"stell", "stellar, pointlike",
"stellar", "stellar",
"straight", "straight",
"streak", "streak",
"strongly", "strongly",
"surrounded", "surrounded",
"surrounds", "surrounds",
"susp", "suspected",
"suspected", "suspected",
"tau", "τ",
"theta", "θ",
"trap", "trapezium",
"trapez", "trapezium",
"trapezium", "trapezium",
"triN", "trinuclear",
"v", "very",
"var", "variable",
"variable", "variable",
"verification", "verification",
"verified", "verified",
"very", "very",
"vl", "very little",
"vm", "very much",
"vs", "very suddenly",
"vv", "very very",
"zeta", "ζ",
0, 0,
};
/*&
"ST9", "stars from the 9th magnitude downward",
"ST9...13", "stars from 9th to 13th magnitude",
"()", "items questioned by Dreyer enclosed in parentheses",
*10 star of 10th mag
"*", "a star (or stars)",
"**", "double star",
"***", "triple star",
*/

88
src/cmd/scat/display.c Normal file
View file

@ -0,0 +1,88 @@
#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;
}
}

113
src/cmd/scat/dssread.c Normal file
View file

@ -0,0 +1,113 @@
#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];
}

247
src/cmd/scat/header.c Normal file
View file

@ -0,0 +1,247 @@
#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";
}

231
src/cmd/scat/hinv.c Normal file
View file

@ -0,0 +1,231 @@
#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++;
}
}

158
src/cmd/scat/image.c Normal file
View file

@ -0,0 +1,158 @@
#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;
}

31
src/cmd/scat/mkfile Normal file
View file

@ -0,0 +1,31 @@
<$PLAN9/src/mkhdr
TARG=scat
OFILES=scat.$O\
bitinput.$O\
desc.$O\
display.$O\
dssread.$O\
header.$O\
hinv.$O\
image.$O\
patch.$O\
plot.$O\
posn.$O\
prose.$O\
qtree.$O\
util.$O\
HFILES=sky.h
CFLAGS=$CFLAGS -I../map
LDFLAGS=$LDFLAGS -L$X11/lib -lX11
SHORTLIB=draw bio 9
LIB=../map/libmap/libmap.a
<$PLAN9/src/mkone
scat.$O: strings.c
$LIB:V:
cd ../map/libmap; mk

101
src/cmd/scat/patch.c Normal file
View file

@ -0,0 +1,101 @@
#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;
}

138
src/cmd/scat/plate.h Normal file
View file

@ -0,0 +1,138 @@
#define RAD(x) ((x)*PI_180)
#define DEG(x) ((x)/PI_180)
#define ARCSECONDS_PER_RADIAN (DEG(1)*3600)
#define input_nybble(infile) input_nbits(infile,4)
typedef float Angle; /* in radians */
enum
{
/*
* parameters for plate
*/
Pppo1 = 0,
Pppo2,
Pppo3,
Pppo4,
Pppo5,
Pppo6,
Pamdx1,
Pamdx2,
Pamdx3,
Pamdx4,
Pamdx5,
Pamdx6,
Pamdx7,
Pamdx8,
Pamdx9,
Pamdx10,
Pamdx11,
Pamdx12,
Pamdx13,
Pamdx14,
Pamdx15,
Pamdx16,
Pamdx17,
Pamdx18,
Pamdx19,
Pamdx20,
Pamdy1,
Pamdy2,
Pamdy3,
Pamdy4,
Pamdy5,
Pamdy6,
Pamdy7,
Pamdy8,
Pamdy9,
Pamdy10,
Pamdy11,
Pamdy12,
Pamdy13,
Pamdy14,
Pamdy15,
Pamdy16,
Pamdy17,
Pamdy18,
Pamdy19,
Pamdy20,
Ppltscale,
Pxpixelsz,
Pypixelsz,
Ppltra,
Ppltrah,
Ppltram,
Ppltras,
Ppltdec,
Ppltdecd,
Ppltdecm,
Ppltdecs,
Pnparam,
};
typedef struct Plate Plate;
struct Plate
{
char rgn[7];
char disk;
Angle ra;
Angle dec;
};
typedef struct Header Header;
struct Header
{
float param[Pnparam];
int amdflag;
float x;
float y;
float xi;
float eta;
};
typedef long Type;
typedef struct Image Image;
struct Image
{
int nx;
int ny; /* ny is the fast-varying dimension */
Type a[1];
};
int nplate;
Plate plate[2000]; /* needs to go to 2000 when the north comes */
double PI_180;
double TWOPI;
int debug;
struct
{
float min;
float max;
float del;
double gamma;
int neg;
} gam;
char* hms(Angle);
char* dms(Angle);
double xsqrt(double);
Angle dist(Angle, Angle, Angle, Angle);
Header* getheader(char*);
char* getword(char*, char*);
void amdinv(Header*, Angle, Angle, float, float);
void ppoinv(Header*, Angle, Angle);
void xypos(Header*, Angle, Angle, float, float);
void traneqstd(Header*, Angle, Angle);
Angle getra(char*);
Angle getdec(char*);
void getplates(void);
Image* dssread(char*);
void hinv(Type*, int, int);
int input_bit(Biobuf*);
int input_nbits(Biobuf*, int);
void qtree_decode(Biobuf*, Type*, int, int, int, int);
void start_inputing_bits(void);
Bitmap* image(Angle, Angle, Angle, Angle);
int dogamma(int);

952
src/cmd/scat/plot.c Normal file
View file

@ -0,0 +1,952 @@
#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]);
}

247
src/cmd/scat/posn.c Normal file
View file

@ -0,0 +1,247 @@
#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);
}
}

168
src/cmd/scat/prose.c Normal file
View file

@ -0,0 +1,168 @@
#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);
}

278
src/cmd/scat/qtree.c Normal file
View file

@ -0,0 +1,278 @@
#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);
}

1671
src/cmd/scat/scat.c Normal file

File diff suppressed because it is too large Load diff

413
src/cmd/scat/sky.h Normal file
View file

@ -0,0 +1,413 @@
#define DIR "#9/sky"
/*
* This code reflects many years of changes. There remain residues
* of prior implementations.
*
* Keys:
* 32 bits long. High 26 bits are encoded as described below.
* Low 6 bits are types:
*
* Patch is ~ one square degree of sky. It points to an otherwise
* anonymous list of Catalog keys. The 0th key is special:
* it contains up to 4 constellation identifiers.
* Catalogs (SAO,NGC,M,...) are:
* 31.........8|76|543210
* catalog # |BB|catalog name
* BB is two bits of brightness:
* 00 -inf < m <= 7
* 01 7 < m <= 10
* 10 10 < m <= 13
* 11 13 < m < inf
* The BB field is a dreg, and correct only for SAO and NGC.
* IC(n) is just NGC(n+7840)
* Others should be self-explanatory.
*
* Records:
*
* Star is an SAOrec
* Galaxy, PlanetaryN, OpenCl, GlobularCl, DiffuseN, etc., are NGCrecs.
* Abell is an Abellrec
* The Namedxxx records hold a name and a catalog entry; they result from
* name lookups.
*/
typedef enum
{
Planet,
Patch,
SAO,
NGC,
M,
Constel_deprecated,
Nonstar_deprecated,
NamedSAO,
NamedNGC,
NamedAbell,
Abell,
/* NGC types */
Galaxy,
PlanetaryN,
OpenCl,
GlobularCl,
DiffuseN,
NebularCl,
Asterism,
Knot,
Triple,
Double,
Single,
Uncertain,
Nonexistent,
Unknown,
PlateDefect,
/* internal */
NGCN,
PatchC,
NONGC,
}Type;
enum
{
/*
* parameters for plate
*/
Pppo1 = 0,
Pppo2,
Pppo3,
Pppo4,
Pppo5,
Pppo6,
Pamdx1,
Pamdx2,
Pamdx3,
Pamdx4,
Pamdx5,
Pamdx6,
Pamdx7,
Pamdx8,
Pamdx9,
Pamdx10,
Pamdx11,
Pamdx12,
Pamdx13,
Pamdx14,
Pamdx15,
Pamdx16,
Pamdx17,
Pamdx18,
Pamdx19,
Pamdx20,
Pamdy1,
Pamdy2,
Pamdy3,
Pamdy4,
Pamdy5,
Pamdy6,
Pamdy7,
Pamdy8,
Pamdy9,
Pamdy10,
Pamdy11,
Pamdy12,
Pamdy13,
Pamdy14,
Pamdy15,
Pamdy16,
Pamdy17,
Pamdy18,
Pamdy19,
Pamdy20,
Ppltscale,
Pxpixelsz,
Pypixelsz,
Ppltra,
Ppltrah,
Ppltram,
Ppltras,
Ppltdec,
Ppltdecd,
Ppltdecm,
Ppltdecs,
Pnparam,
};
#define UNKNOWNMAG 32767
#define NPlanet 20
typedef float Angle; /* in radians */
typedef long DAngle; /* on disk: in units of milliarcsec */
typedef short Mag; /* multiplied by 10 */
typedef long Key; /* known to be 4 bytes, unfortunately */
/*
* All integers are stored in little-endian order.
*/
typedef struct NGCrec NGCrec;
struct NGCrec{
DAngle ra;
DAngle dec;
DAngle dummy1; /* compatibility with old RNGC version */
DAngle diam;
Mag mag;
short ngc; /* if >NNGC, IC number is ngc-NNGC */
char diamlim;
char type;
char magtype;
char dummy2;
char desc[52]; /* 0-terminated Dreyer description */
};
typedef struct Abellrec Abellrec;
struct Abellrec{
DAngle ra;
DAngle dec;
DAngle glat;
DAngle glong;
Mag mag10; /* mag of 10th brightest cluster member; in same place as ngc.mag*/
short abell;
DAngle rad;
short pop;
short dist;
char distgrp;
char richgrp;
char flag;
char pad;
};
typedef struct Planetrec Planetrec;
struct Planetrec{
DAngle ra;
DAngle dec;
DAngle az;
DAngle alt;
DAngle semidiam;
double phase;
char name[16];
};
/*
* Star names: 0,0==unused. Numbers are name[0]=1,..,99.
* Greek letters are alpha=101, etc.
* Constellations are alphabetical order by abbreviation, and=1, etc.
*/
typedef struct SAOrec SAOrec;
struct SAOrec{
DAngle ra;
DAngle dec;
DAngle dra;
DAngle ddec;
Mag mag; /* visual */
Mag mpg;
char spec[3];
char code;
char compid[2];
char hdcode;
char pad1;
long hd; /* HD catalog number */
char name[3]; /* name[0]=alpha name[1]=2 name[3]=ori */
char nname; /* number of prose names */
/* 36 bytes to here */
};
typedef struct Mindexrec Mindexrec;
struct Mindexrec{ /* code knows the bit patterns in here; this is a long */
char m; /* M number */
char dummy;
short ngc;
};
typedef struct Bayerec Bayerec;
struct Bayerec{
long sao;
char name[3];
char pad;
};
/*
* Internal form
*/
typedef struct Namedrec Namedrec;
struct Namedrec{
char name[36];
};
typedef struct Namerec Namerec;
struct Namerec{
long sao;
long ngc;
long abell;
char name[36]; /* null terminated */
};
typedef struct Patchrec Patchrec;
struct Patchrec{
int nkey;
long key[60];
};
typedef struct Record Record;
struct Record{
Type type;
long index;
union{
SAOrec sao;
NGCrec ngc;
Abellrec abell;
Namedrec named;
Patchrec patch;
Planetrec planet;
/* PatchCrec is empty */
};
};
typedef struct Name Name;
struct Name{
char *name;
int type;
};
typedef struct Plate Plate;
struct Plate
{
char rgn[7];
char disk;
Angle ra;
Angle dec;
};
typedef struct Header Header;
struct Header
{
float param[Pnparam];
int amdflag;
float x;
float y;
float xi;
float eta;
};
typedef long Pix;
typedef struct Img Img;
struct Img
{
int nx;
int ny; /* ny is the fast-varying dimension */
Pix a[1];
};
#define RAD(x) ((x)*PI_180)
#define DEG(x) ((x)/PI_180)
#define ARCSECONDS_PER_RADIAN (DEG(1)*3600)
#define MILLIARCSEC (1000*60*60)
int nplate;
Plate plate[2000]; /* needs to go to 2000 when the north comes */
double PI_180;
double TWOPI;
double LN2;
int debug;
struct
{
float min;
float max;
float gamma;
float absgamma;
float mult1;
float mult2;
int neg;
} gam;
typedef struct Picture Picture;
struct Picture
{
int minx;
int miny;
int maxx;
int maxy;
char name[16];
uchar *data;
};
#ifndef _DRAW_H_
typedef struct Image Image;
#endif
extern double PI_180;
extern double TWOPI;
extern char *progname;
extern char *desctab[][2];
extern Name names[];
extern Record *rec;
extern long nrec;
extern Planetrec *planet;
/* for bbox: */
extern int folded;
extern DAngle ramin;
extern DAngle ramax;
extern DAngle decmin;
extern DAngle decmax;
extern Biobuf bout;
extern void saoopen(void);
extern void ngcopen(void);
extern void patchopen(void);
extern void mopen(void);
extern void constelopen(void);
extern void lowercase(char*);
extern void lookup(char*, int);
extern int typetab(int);
extern char*ngcstring(int);
extern char*skip(int, char*);
extern void prrec(Record*);
extern int equal(char*, char*);
extern int parsename(char*);
extern void radec(int, int*, int*, int*);
extern int btag(short);
extern long patcha(Angle, Angle);
extern long patch(int, int, int);
extern char*hms(Angle);
extern char*dms(Angle);
extern char*ms(Angle);
extern char*hm(Angle);
extern char*dm(Angle);
extern char*deg(Angle);
extern char*hm5(Angle);
extern long dangle(Angle);
extern Angle angle(DAngle);
extern void prdesc(char*, char*(*)[2], short*);
extern double xsqrt(double);
extern Angle dist(Angle, Angle, Angle, Angle);
extern Header* getheader(char*);
extern char* getword(char*, char*);
extern void amdinv(Header*, Angle, Angle, float, float);
extern void ppoinv(Header*, Angle, Angle);
extern void xypos(Header*, Angle, Angle, float, float);
extern void traneqstd(Header*, Angle, Angle);
extern Angle getra(char*);
extern Angle getdec(char*);
extern void getplates(void);
extern Img* dssread(char*);
extern void hinv(Pix*, int, int);
extern int input_bit(Biobuf*);
extern int input_nbits(Biobuf*, int);
extern int input_huffman(Biobuf*);
extern int input_nybble(Biobuf*);
extern void qtree_decode(Biobuf*, Pix*, int, int, int, int);
extern void start_inputing_bits(void);
extern Picture* image(Angle, Angle, Angle, Angle);
extern char* dssmount(int);
extern int dogamma(Pix);
extern void displaypic(Picture*);
extern void displayimage(Image*);
extern void plot(char*);
extern void astro(char*, int);
extern char* alpha(char*, char*);
extern char* skipbl(char*);
extern void flatten(void);
extern int bbox(long, long, int);
extern int inbbox(DAngle, DAngle);
extern char* nameof(Record*);
#define NINDEX 400

52
src/cmd/scat/strings.c Normal file
View file

@ -0,0 +1,52 @@
char *greek[]={ 0, /* 1-indexed */
"alpha", "beta", "gamma", "delta", "epsilon", "zeta", "eta", "theta",
"iota", "kappa", "lambda", "mu", "nu", "xsi", "omicron", "pi", "rho",
"sigma", "tau", "upsilon", "phi", "chi", "psi", "omega",
};
Rune greeklet[]={ 0,
0x3b1, 0x3b2, 0x3b3, 0x3b4, 0x3b5, 0x3b6, 0x3b7, 0x3b8, 0x3b9, 0x3ba, 0x3bb,
0x3bc, 0x3bd, 0x3be, 0x3bf, 0x3c0, 0x3c1, 0x3c3, 0x3c4, 0x3c5, 0x3c6, 0x3c7,
0x3c8, 0x3c9,
};
char *constel[]={ 0, /* 1-indexed */
"and", "ant", "aps", "aql", "aqr", "ara", "ari", "aur", "boo", "cae",
"cam", "cap", "car", "cas", "cen", "cep", "cet", "cha", "cir", "cma",
"cmi", "cnc", "col", "com", "cra", "crb", "crt", "cru", "crv", "cvn",
"cyg", "del", "dor", "dra", "equ", "eri", "for", "gem", "gru", "her",
"hor", "hya", "hyi", "ind", "lac", "leo", "lep", "lib", "lmi", "lup",
"lyn", "lyr", "men", "mic", "mon", "mus", "nor", "oct", "oph", "ori",
"pav", "peg", "per", "phe", "pic", "psa", "psc", "pup", "pyx", "ret",
"scl", "sco", "sct", "ser", "sex", "sge", "sgr", "tau", "tel", "tra",
"tri", "tuc", "uma", "umi", "vel", "vir", "vol", "vul",
};
Name names[]={
"gx", Galaxy,
"pl", PlanetaryN,
"oc", OpenCl,
"gb", GlobularCl,
"nb", DiffuseN,
"c+n",NebularCl,
"ast", Asterism,
"kt", Knot,
"***", Triple,
"d*", Double,
"*", Single,
"pd", PlateDefect,
"galaxy", Galaxy,
"planetary", PlanetaryN,
"opencluster", OpenCl,
"globularcluster", GlobularCl,
"nebula", DiffuseN,
"nebularcluster",NebularCl,
"asterism", Asterism,
"knot", Knot,
"triple", Triple,
"double", Double,
"single", Single,
"nonexistent", Nonexistent,
"unknown", Unknown,
"platedefect", PlateDefect,
0, 0,
};

368
src/cmd/scat/util.c Normal file
View file

@ -0,0 +1,368 @@
#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;
}