# HG changeset patch # User pouya@nohup.io # Date 1612133119 -3600 # Sun Jan 31 23:45:19 2021 +0100 # Node ID be34223b93a0ba5ff53525c951ec06dba7d846da # Parent 23ecfc48ea840e6bb12bd1e48afae63ea641a770 first (untested) commit of the 'sky' mini-collection of astronomical tools diff --git a/sky/LICENCE b/sky/LICENCE new file mode 100644 --- /dev/null +++ b/sky/LICENCE @@ -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. + diff --git a/sky/cmd/skywatch.c b/sky/cmd/skywatch.c new file mode 100644 diff --git a/sky/include/lex.h b/sky/include/lex.h new file mode 100644 --- /dev/null +++ b/sky/include/lex.h @@ -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); diff --git a/sky/include/sky.h b/sky/include/sky.h new file mode 100644 --- /dev/null +++ b/sky/include/sky.h @@ -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 ) + diff --git a/sky/liblex/README b/sky/liblex/README new file mode 100644 --- /dev/null +++ b/sky/liblex/README @@ -0,0 +1,1 @@ +Basic lexer after Rob Pike. Implemented but not yet tested! diff --git a/sky/liblex/lex.c b/sky/liblex/lex.c new file mode 100644 --- /dev/null +++ b/sky/liblex/lex.c @@ -0,0 +1,112 @@ +#include +#include +#include +#include + +#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; +} diff --git a/sky/liblex/mkfile b/sky/liblex/mkfile new file mode 100644 --- /dev/null +++ b/sky/liblex/mkfile @@ -0,0 +1,20 @@ + +#include + +#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; +}