@@ 0,0 1,27 @@ 
+Copyright (c) 2021 Pouya Tafti.  All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+1.  Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2.  Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the
+   distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC.  AND
+CONTRIBUTORS ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
+INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
+IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+

          

@@ 0,0 1,42 @@ 
+typedef struct State State;
+typedef struct Tok Tok;
+typedef struct Lexer Lexer;
+
+enum {
+	Ltoksize = 256,
+
+	Leof = -1,
+	Loverflow = -10,
+};
+
+/* wrapper to simulate recursive function type */
+struct State {
+	State (*next)(Lexer *l);
+};
+
+struct Tok {
+	int typ;
+	Rune str[Ltoksize+1];
+};
+	
+struct Lexer {
+	Biobuf *bp;
+	Rune b[Ltoksize+1], *cur;
+	int start, pos;
+	Channel *c;
+	State state;
+};
+
+int linit(Lexer *l, int fd, State state);
+void lterm(Lexer *l);
+
+long lnext(Lexer *l);
+void lbackup(Lexer *l);
+long lpeek(Lexer *l);
+
+void lignore(Lexer *l);
+int laccept(Lexer *l, Rune *s);
+void lacceptrun(Lexer *l, Rune *s);
+
+void lemit(Lexer *l, int typ);
+Tok lnexttok(Lexer *l);

          

@@ 0,0 1,103 @@ 
+/** dat **/
+
+typedef struct Gdloc Gdloc;
+typedef struct Eqloc Eqloc;
+typedef struct Hzloc Hzloc;
+typedef struct Gxloc Gxloc; 
+typedef struct Planar Planar;
+typedef struct Obs Obs;
+typedef struct Viewangle Viewangle;
+typedef struct Camera Camera;
+
+/* geodetic location */
+struct Gdloc {
+	double lat;	/* in rad */
+	double lng;	/* in rad */
+	double alt;	/* in meters from mean sea level */
+};
+
+/* equatorial location */
+struct Eqloc {
+	double ra;		/* right ascension in rad */
+	double dec;	/* declination in rad */
+};
+
+/* horizontal location */
+struct Hzloc {
+	double az;	/* azimuth in rad */
+	double z;		/* zenith in rad (complement of elevation) */
+};
+
+/* galactic location */
+struct Gxloc {
+	double b;		/* galactic latitude in rad */
+	double l;		/* galactic longitude in rad */
+};
+
+/* 2D projection */
+struct Planar {
+	double x, y;
+};
+
+struct Obs {
+	Gdloc g;
+	double P;		/* pressure in kPa */
+	double T;		/* temperature in K */
+	double rh;		/* relative humidity */
+};
+
+struct Viewangle {
+	double az;	/* azimuth in rad */
+	double z;		/* zenith in rad (complement of elevation) */
+	double roll;	/* roll in rad */
+};
+
+struct Camera {
+	double f;		/* focal length in mm */
+	double w, h;	/* sensor dimensions in mm */
+	int wpx, hpx;	/* sensor dimensions in pixels */		
+	Viewangle angle;
+};
+
+#define TWOPI (PI+PI)
+#define DEG 1.745329251994329576923e-2
+#define DMIN 2.908882086657215961539e-4
+#define DSEC 4.848136811095359935899e-6
+#define HOUR 2.617993877991494365386e-1
+#define MIN 4.363323129985823942309e-3
+#define SEC 7.272205216643039903849e-5
+
+#define JD00 2451545.0
+#define JDX100 36525.0
+
+/* atmospheric conditions for refraction */
+#define STP_kPa 100.0d
+#define STP_K 273.15d
+#define NTP_kPa 101.325d
+#define NTP_K 293.15d
+
+
+/** fns **/
+
+Eqloc gx2eq00(Gxloc gx);
+Gxloc eq2gx00(Eqloc eq);
+
+Hzloc eq2hz(Eqloc eq, Gdloc geod, double jdutc, double dut1, double dt);
+Eqloc hz2eq(Hzloc hz, Gdloc geod, double jdutc, double dut1, double dt);
+
+Planar hz2pxl(Camera *c, Hzloc hz);
+Hzloc pxl2hz(Camera *c, Planar px);
+
+double gast(double jd, double dut1, double dt);
+double eqeqx(double jd, double dut1, double dt);
+
+double saemundsson(Obs *o, double z);
+
+#define deg2rad(d) ( (d) * DEG )
+#define rad2deg(a) ( (a) / DEG )
+#define h2rad(h) ( (h) * HOUR )
+#define rad2h(a) ( (a) / HOUR )
+#define hms2h(h, m, s) ( (h) + (m)/60.0d + (s)/3600.0d )
+#define hms2rad(h, m, s) ( (h)*HOUR+(m)*MIN+(s)*SEC )
+#define degms2rad(d, m, s) ( (d)*DEG+(m)*DMIN+(s)*DSEC )
+

          

@@ 0,0 1,1 @@ 
+Basic lexer after Rob Pike.  Implemented but not yet tested!

          

@@ 0,0 1,112 @@ 
+#include <u.h>
+#include <libc.h>
+#include <thread.h>
+#include <bio.h>
+
+#include "lex.h"
+
+int
+linit(Lexer *l, int fd, State state)
+{
+	if((l->c=chancreate(sizeof(Tok), 2)) == nil)
+		return Leof;
+
+	if(Binit(l->bp, fd, OREAD) == Beof) {
+		chanfree(l->c);
+		return Leof;
+	}
+
+	l->start = l->pos = 0;
+	l->cur = l->b;
+	l->state = state;
+
+	return 0;
+}
+
+void
+lterm(Lexer *l)
+{
+	chanfree(l->c);
+	Bterm(l->bp);
+}
+
+long
+lnext(Lexer *l)
+{
+	long c;
+
+	if(l->pos-l->start > Ltoksize) return Loverflow;
+
+	if((c = Bgetrune(l->bp)) == Beof) return Leof;
+
+	l->pos++;
+	*l->cur++ = c;
+	return c; 
+}
+
+void
+lbackup(Lexer *l)
+{
+	if(l->pos <= l->start) return;
+	Bungetrune(l->bp);
+	l->pos--;
+	l->cur--;
+}
+
+long
+lpeek(Lexer *l)
+{
+	long c;
+
+	c = lnext(l);
+	lbackup(l);
+
+	return c;
+}
+
+void
+lignore(Lexer *l)
+{
+	l->start = l->pos;
+	l->cur = l->b;
+}
+
+int
+laccept(Lexer *l, Rune *s)
+{
+	long c;
+
+	if((c = lnext(l)) < 0) return 0;
+
+	if(runestrchr(s, c) != 0) return 1;
+
+	lbackup(l);
+	return 0;
+}
+
+void
+lacceptrun(Lexer *l, Rune *s)
+{
+	while(laccept(l, s));
+}
+
+void
+lemit(Lexer *l, int typ)
+{
+	Tok t;
+
+	t.typ = typ;
+	runestrncpy(t.str, l->b, l->pos-l->start);
+	t.str[Ltoksize] = 0;
+
+	send(l->c, &t);
+}
+
+Tok
+lnexttok(Lexer *l)
+{
+	Tok t;
+
+	while(!nbrecv(l->c, &t)) l->state = l->state.next(l);
+	return t;	
+}

          

@@ 0,0 1,20 @@ 
+</$objtype/mkfile
+
+LIB=liblex.a
+
+INCDIR=../include
+
+CFLAGS=-I$INCDIR
+
+OFILES=\
+	${LIB:lib%.a=%.$O}\
+
+HFILES=\
+	$INCDIR/${LIB:lib%.a=%.h}\
+
+UPDATE=\
+	mkfile\
+	$HFILES\
+	${OFILES:%.$O=%.c}
+
+</sys/src/cmd/mklib
  No newline at end of file

          

@@ 0,0 1,19 @@ 
+This library provides a minimally complete set of functions to convert
+between various celestial coordinates, often favouring simplicity over
+precision.
+
+Starting with the B2000 galactic coordinates of an astronomical
+object, you can obtain its coordinates in the equatorial and
+horizontal systems, and its pixel coordinates through a simple
+telescope, using the sequence of functions gx2eq, eq2hz, and hz2pxl.
+You can further use Saemundsson's formula, provided by the function
+saemundsson, to correct for atmospheric refraction (before the optical
+projection step).
+
+A few other lower level functions and macros are also available.  For
+details see the source code.
+
+The code is in the Plan 9 dialect of C, but it should be usable on
+most systems with a standard C compiler and support for floating point
+arithmetic with almost no change other than replacing the standard
+header files (and potentially the definition of NaN).

          

@@ 0,0 1,20 @@ 
+</$objtype/mkfile
+
+LIB=libsky.a
+
+INCDIR=../include
+
+CFLAGS=-I$INCDIR
+
+OFILES=\
+	${LIB:lib%.a=%.$O}\
+
+HFILES=\
+	$INCDIR/${LIB:lib%.a=%.h}\
+
+UPDATE=\
+	mkfile\
+	$HFILES\
+	${OFILES:%.$O=%.c}
+
+</sys/src/cmd/mklib
  No newline at end of file

          

@@ 0,0 1,345 @@ 
+#include <u.h>
+#include <libc.h>
+
+#include "sky.h"
+
+double gammab_fw(double jd, double dut1, double dt);
+double phib_fw(double jd, double dut1, double dt);
+double psib_fw(double jd, double dut1, double dt);
+double epsA_fw(double jd, double dut1, double dt);
+void sXY(double jd, double dut1, double dt, double *s, double *X, double *Y);
+
+double _sdG = 0;
+double _cdG = 0;
+
+Eqloc
+gx2eq00(Gxloc gx)
+{
+	double sb, cb, slCPl, clCPl, sd, saaGcd, caaGcd;
+	Eqloc eq;
+
+	if(_sdG==0) _sdG = sin(deg2rad(27.12825));
+	if(_cdG==0) _cdG = cos(deg2rad(27.12825));
+
+	sb = sin(gx.b);
+	cb = cos(gx.b);
+	slCPl = sin(deg2rad(122.93192)-gx.l);
+	clCPl = cos(deg2rad(122.93192)-gx.l);
+
+	sd = sb*_sdG + cb*_cdG*clCPl;
+
+	saaGcd = cb*slCPl;
+	caaGcd = sb*_cdG - cb*_sdG*clCPl;
+
+	eq.dec = asin(sd);
+	eq.ra = atan2(saaGcd, caaGcd) + deg2rad(192.85948);
+	while(eq.ra<-PIO2) eq.ra+=PI;
+	while(eq.ra>PIO2) eq.ra -= PI;
+
+	return eq;
+}
+
+Gxloc
+eq2gx00(Eqloc eq)
+{
+	double sd, cd, saaG, caaG, sb, slCPlcb, clCPlcb;
+	Gxloc gx;
+
+	if(_sdG==0) _sdG = sin(deg2rad(27.12825));
+	if(_cdG==0) _cdG = cos(deg2rad(27.12825));
+
+	sd = sin(eq.dec);
+	cd = cos(eq.dec);
+	saaG = sin(eq.ra - deg2rad(192.85948));
+	caaG = cos(eq.ra - deg2rad(192.85948));
+
+	sb = sd*_sdG + cd*_cdG*caaG;
+
+	slCPlcb = cd*saaG;
+	clCPlcb = sd*_cdG - cd*_sdG*caaG;
+
+	gx.b = asin(sb);
+	gx.l = deg2rad(122.93192) - atan2(slCPlcb, clCPlcb);
+	while(gx.l<-PIO2) gx.l+=PI;
+	while(gx.l>PIO2) gx.l -= PI;
+
+	return gx;
+}
+
+Hzloc
+eq2hz(Eqloc eq, Gdloc geod, double jdutc, double dut1, double dt)
+{
+	double gst, /* tU, era, */ h, sph, cph, sd, cd, sh, ch, u, v, w;
+	Hzloc hz;
+
+	gst = gast(jdutc, dut1, dt);
+	h = gst + geod.lng - eq.ra;
+/*
+	tU = jd+dut1-JD00;
+	era = TWOPI * (tU -(long)tU + 0.7790572732640 + 0.00273781191135448 * tU);
+	h = era + geod.lng - eq.ra;
+*/
+	sph = sin(geod.lat);
+	cph = cos(geod.lat);
+	sd = sin(eq.dec);
+	cd = cos(eq.dec);
+	sh = sin(h);
+	ch = cos(h);
+
+	u = sd*cph - cd*ch*sph;
+	v = -cd*sh;
+	w = sd*sph + cd*ch*cph;
+
+	hz.az = atan2(v, u);
+	hz.az = hz.az < 0.0 ? hz.az+TWOPI : hz.az;
+
+	hz.z = acos(w);
+	//hz.z = PIO2 - atan2(w, sqrt(u*u+v*v));
+
+	return hz;
+}
+
+Eqloc
+hz2eq(Hzloc hz, Gdloc geod, double jdutc, double dut1, double dt)
+{
+	double sph, cph, saz, caz, sz, cz, u, v, w, h, gst;
+	Eqloc eq;
+
+	sph = sin(geod.lat);
+	cph = cos(geod.lat);
+	saz = sin(hz.az);
+	caz = cos(hz.az);
+	sz = sin(hz.z);
+	cz = cos(hz.z);
+
+	u = cz*cph - caz*sz*sph;
+	v = -saz*sz;
+	w = cz*sph + caz*sz*cph;
+
+	h = atan2(v, u);
+	gst = gast(jdutc, dut1, dt);
+	
+	eq.ra = gst + geod.lng - h;
+	eq.dec = asin(w);
+	//eq.dec = atan2(w, sqrt(u*u+v*v));
+
+	return eq;
+}
+
+Planar
+hz2pxl(Camera *c, Hzloc hz)
+{
+	double dz, daz, tz, taz, sr, cr;
+	Planar pt, ptr, px;
+
+	dz = hz.z - (c->angle).z;
+	daz = hz.az - (c->angle).az;
+
+	tz = tan(dz);
+	taz = tan(daz);
+
+	sr = sin((c->angle).roll);
+	cr = cos((c->angle).roll);
+
+	pt.x = -c->f*tz;
+	pt.y = c->f*taz;
+	ptr.x = pt.x*cr-pt.y*sr;
+	ptr.y = pt.x*sr+pt.y*cr;
+
+	while(dz < -PI) dz += TWOPI;
+	while(dz > PI) dz -= TWOPI;
+	while(daz < -PI) daz += TWOPI;
+	while(daz > PI) daz -= TWOPI;
+
+	/* behind the camera */
+	if(dz < -PIO2 || dz > PIO2 || daz < -PIO2 || daz > PIO2){
+		px.x = NaN();
+		px.y = NaN();
+		return px;
+	}
+
+	px.x = ptr.x*c->wpx/c->w;
+	px.y = ptr.y*c->hpx/c->h;
+	return px;
+}
+
+Hzloc
+pxl2hz(Camera *c, Planar px)
+{
+	double sr, cr, tz, taz;
+	Planar pt, ptr;
+	Hzloc hz;
+
+	sr = sin(c->angle.roll);
+	cr = cos(c->angle.roll);
+
+	ptr.x = px.x*c->w/c->wpx;
+	ptr.y = px.y*c->h/c->hpx;
+	pt.x = ptr.x*cr+ptr.y*sr;
+	pt.y = -ptr.x*sr+ptr.y*cr;
+
+	tz = -pt.x / c->f;
+	taz = -pt.y / c->f;
+
+	hz.z = c->angle.z + atan(tz);
+	hz.az = c->angle.az + atan(taz);
+
+	return hz;
+}
+
+double
+gast(double jd, double dut1, double dt)
+{
+	double tU, t, era, gmst;
+
+	tU = jd+dut1-JD00;
+	t = (tU+dt)/JDX100;
+
+	/* earth rotation angle is used for the CIO-based hour angle */
+	era = TWOPI * (tU -(long)tU + 0.7790572732640 + 0.00273781191135448 * tU);
+	
+	/* from Fukushima 2003 */
+	gmst = era +
+		(0.012911569 + t*4612.160517397) * DSEC +
+		t*t*(1.391542507 + t*(-0.000124849 + t*(-0.000004991 + t*(-0.000000479)))) * DSEC;
+
+//	/* GMST approx from Capitaine et al. 2005 */
+//	gmst = era + DSEC * (((((
+//		-0.0000000368)*t +
+//		-0.000029956)*t +
+//		-0.00000044)*t +
+//		1.3915817)*t +
+//		4612.156534)*t +
+//		0.014506;
+
+	return gmst+eqeqx(jd, dut1, dt);
+}
+
+double
+eqeqx(double jd, double dut1, double dt)
+{
+	double tU, t, om, lm, ls, gs, gls, dpsi, epsbar;
+
+	tU = jd+dut1-JD00;
+	t = (tU+dt)/JDX100;
+
+	/* longitude of the scending node of moon's orbit */
+	om = (125.04452 + t*(-134.136261 + t*(0.0020708 + t/450000.0))) * DEG;
+
+	/* mean longitude of the moon */
+	lm = (218.31654591 + t*(481267.88134236 + t*(-0.00163 + t*(1./538841.0 - t/65194000.0)))) * DEG;
+
+	/* mean longitude of the sun */
+	ls = (280.46645 + t*(36000.76983 + t*0.0003032)) * DEG;
+
+	/* mean anomaly of the sun */
+	gs = (357.52910 + t*(35999.05030 + t*( -0.0001559 + t*(-0.00000048)))) * DEG;
+
+	/* geometric longitude of the sun */
+	gls = ls + ( (1.9146000 + t*(-0.004817 + t*(- 0.000014)))*sin(gs) + (0.019993 - 0.000101*t)*sin(2.0*gs) + 0.000290*sin(3.0*gs) ) * DEG;
+	
+	dpsi =  ( -17.20*sin(om) - 1.32*sin(gls+gls) - 0.23*sin(lm+lm) + 0.21*sin(om+om) ) * DSEC;
+	epsbar = (84381.448 + t*(-46.8150 + t*(-0.00059 + t*0.001813))) * DSEC;
+
+	return dpsi*cos(epsbar) + (0.00264*sin(om) + 0.000063*sin(om+om)) * DSEC;	
+}
+
+double
+saemundsson(Obs *o, double z)
+{
+	return -h2rad(1.02/60.0) *
+		o->P/101.0 *
+		283.0/o->T /
+		tan(PIO2-z+deg2rad(10.3)/
+			(PIO2-z+deg2rad(5.11)));
+}
+
+
+/* functions below this line are not used at the moment but kept for future use */
+
+double
+gammab_fw(double jd, double dut1, double dt)
+{
+	double t;
+
+	t = (jd+dut1+dt-JD00) / JDX100;
+
+	/* from Wallace & Capitaine, 2006 (Fukushima-Williams equations) */
+	return ((((((
+		0.0000000260)*t - 
+		0.000002788)*t -
+		0.00031238)*t +
+		0.4932044)*t +
+		10.556378)*t -
+		0.052928) * DSEC;
+}
+
+double
+phib_fw(double jd, double dut1, double dt)
+{
+	double t;
+
+	t = (jd+dut1+dt-JD00) / JDX100;
+
+	/* from Wallace & Capitaine, 2006 (Fukushima-Williams equations) */
+	return ((((((
+		-0.0000000176)*t -
+		0.000000440)*t +
+		0.00053289)*t +
+		0.0511268)*t -
+		46.811016)*t +
+		84381.412819) * DSEC;
+}
+
+double
+psib_fw(double jd, double dut1, double dt)
+{
+	double t;
+
+	t = (jd+dut1+dt-JD00) / JDX100;
+
+	/* from Wallace & Capitaine, 2006 (Fukushima-Williams equations) */
+	return ((((((
+		-0.0000000148)*t +
+		0.000026452)*t -
+		0.00018522)*t -
+		1.5584175)*t +
+		5038.481484)*t -
+		0.041775) * DSEC;
+}
+
+double
+epsA_fw(double jd, double dut1, double dt)
+{
+	double t;
+
+	t = (jd+dut1+dt-JD00) / JDX100;
+
+	/* from Wallace & Capitaine, 2006 (Fukushima-Williams equations) */
+	return ((((((
+		-0.0000000434)*t -
+		0.000000576)*t +
+		0.00200340)*t - 
+		0.0001831)*t -
+		46.836769)*t +
+		84381.406) * DSEC;
+}
+
+void
+sXY(double jd, double dut1, double dt, double *s, double *X, double *Y)
+{
+	double tau, t, sXY2, Omega;
+
+	tau = jd+dut1+dt-JD00;
+	t = tau / JDX100;
+
+	sXY2 = (((((15.62*t + 27.98)*t - 72574.11)*t - 122.68)*t + 3808.65)*t + 94.0);
+	Omega = 2.182 - 9.242e-4*t;
+
+	*X = 2.6603e-7*t - 33.2e-6*sin(Omega);
+	*Y = -8.14e-14*t*t+44.6e-6*cos(Omega);
+	
+	*s = sXY2 - (*X)*(*Y)/2;
+
+	return;	
+}