Blob
1 #include <u.h>2 #include <libc.h>3 #include <bio.h>4 #include "sky.h"6 double PI_180 = 0.0174532925199432957692369;7 double TWOPI = 6.2831853071795864769252867665590057683943387987502;8 double LN2 = 0.69314718055994530941723212145817656807550013436025;9 static double angledangle=(180./PI)*MILLIARCSEC;11 #define rint scatrint13 int14 rint(char *p, int n)15 {16 int i=0;18 while(*p==' ' && n)19 p++, --n;20 while(n--)21 i=i*10+*p++-'0';22 return i;23 }25 DAngle26 dangle(Angle angle)27 {28 return angle*angledangle;29 }31 Angle32 angle(DAngle dangle)33 {34 return dangle/angledangle;35 }37 double38 rfloat(char *p, int n)39 {40 double i, d=0;42 while(*p==' ' && n)43 p++, --n;44 if(*p == '+')45 return rfloat(p+1, n-1);46 if(*p == '-')47 return -rfloat(p+1, n-1);48 while(*p == ' ' && n)49 p++, --n;50 if(n == 0)51 return 0.0;52 while(n-- && *p!='.')53 d = d*10+*p++-'0';54 if(n <= 0)55 return d;56 p++;57 i = 1;58 while(n--)59 d+=(*p++-'0')/(i*=10.);60 return d;61 }63 int64 sign(int c)65 {66 if(c=='-')67 return -1;68 return 1;69 }71 char*72 hms(Angle a)73 {74 static char buf[20];75 double x;76 int h, m, s, ts;78 x=DEG(a)/15;79 x += 0.5/36000.; /* round up half of 0.1 sec */80 h = floor(x);81 x -= h;82 x *= 60;83 m = floor(x);84 x -= m;85 x *= 60;86 s = floor(x);87 x -= s;88 ts = 10*x;89 sprint(buf, "%dh%.2dm%.2d.%ds", h, m, s, ts);90 return buf;91 }93 char*94 dms(Angle a)95 {96 static char buf[20];97 double x;98 int sign, d, m, s, ts;100 x = DEG(a);101 sign='+';102 if(a<0){103 sign='-';104 x=-x;105 }106 x += 0.5/36000.; /* round up half of 0.1 arcsecond */107 d = floor(x);108 x -= d;109 x *= 60;110 m = floor(x);111 x -= m;112 x *= 60;113 s = floor(x);114 x -= s;115 ts = floor(10*x);116 sprint(buf, "%c%d°%.2d'%.2d.%d\"", sign, d, m, s, ts);117 return buf;118 }120 char*121 ms(Angle a)122 {123 static char buf[20];124 double x;125 int d, m, s, ts;127 x = DEG(a);128 x += 0.5/36000.; /* round up half of 0.1 arcsecond */129 d = floor(x);130 x -= d;131 x *= 60;132 m = floor(x);133 x -= m;134 x *= 60;135 s = floor(x);136 x -= s;137 ts = floor(10*x);138 if(d != 0)139 sprint(buf, "%d°%.2d'%.2d.%d\"", d, m, s, ts);140 else141 sprint(buf, "%.2d'%.2d.%d\"", m, s, ts);142 return buf;143 }145 char*146 hm(Angle a)147 {148 static char buf[20];149 double x;150 int h, m, n;152 x = DEG(a)/15;153 x += 0.5/600.; /* round up half of tenth of minute */154 h = floor(x);155 x -= h;156 x *= 60;157 m = floor(x);158 x -= m;159 x *= 10;160 n = floor(x);161 sprint(buf, "%dh%.2d.%1dm", h, m, n);162 return buf;163 }165 char*166 hm5(Angle a)167 {168 static char buf[20];169 double x;170 int h, m;172 x = DEG(a)/15;173 x += 2.5/60.; /* round up 2.5m */174 h = floor(x);175 x -= h;176 x *= 60;177 m = floor(x);178 m -= m % 5;179 sprint(buf, "%dh%.2dm", h, m);180 return buf;181 }183 char*184 dm(Angle a)185 {186 static char buf[20];187 double x;188 int sign, d, m, n;190 x = DEG(a);191 sign='+';192 if(a<0){193 sign='-';194 x=-x;195 }196 x += 0.5/600.; /* round up half of tenth of arcminute */197 d = floor(x);198 x -= d;199 x *= 60;200 m = floor(x);201 x -= m;202 x *= 10;203 n = floor(x);204 sprint(buf, "%c%d°%.2d.%.1d'", sign, d, m, n);205 return buf;206 }208 char*209 deg(Angle a)210 {211 static char buf[20];212 double x;213 int sign, d;215 x = DEG(a);216 sign='+';217 if(a<0){218 sign='-';219 x=-x;220 }221 x += 0.5; /* round up half degree */222 d = floor(x);223 sprint(buf, "%c%d°", sign, d);224 return buf;225 }227 char*228 getword(char *ou, char *in)229 {230 int c;232 for(;;) {233 c = *in++;234 if(c == ' ' || c == '\t')235 continue;236 if(c == 0)237 return 0;238 break;239 }241 if(c == '\'')242 for(;;) {243 if(c >= 'A' && c <= 'Z')244 c += 'a' - 'A';245 *ou++ = c;246 c = *in++;247 if(c == 0)248 return 0;249 if(c == '\'') {250 *ou = 0;251 return in-1;252 }253 }254 for(;;) {255 if(c >= 'A' && c <= 'Z')256 c += 'a' - 'A';257 *ou++ = c;258 c = *in++;259 if(c == ' ' || c == '\t' || c == 0) {260 *ou = 0;261 return in-1;262 }263 }264 }266 /*267 * Read formatted angle. Must contain no embedded blanks268 */269 Angle270 getra(char *p)271 {272 Rune r;273 char *q;274 Angle f, d;275 int neg;277 neg = 0;278 d = 0;279 while(*p == ' ')280 p++;281 for(;;) {282 if(*p == ' ' || *p=='\0')283 goto Return;284 if(*p == '-') {285 neg = 1;286 p++;287 }288 if(*p == '+') {289 neg = 0;290 p++;291 }292 q = p;293 f = strtod(p, &q);294 if(q > p) {295 p = q;296 }297 p += chartorune(&r, p);298 switch(r) {299 default:300 Return:301 if(neg)302 d = -d;303 return RAD(d);304 case 'h':305 d += f*15;306 break;307 case 'm':308 d += f/4;309 break;310 case 's':311 d += f/240;312 break;313 case 0xB0: /* ° */314 d += f;315 break;316 case '\'':317 d += f/60;318 break;319 case '\"':320 d += f/3600;321 break;322 }323 }324 return 0;325 }327 double328 xsqrt(double a)329 {331 if(a < 0)332 return 0;333 return sqrt(a);334 }336 Angle337 dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2)338 {339 double a;341 a = sin(dec1) * sin(dec2) +342 cos(dec1) * cos(dec2) *343 cos(ra1 - ra2);344 a = atan2(xsqrt(1 - a*a), a);345 if(a < 0)346 a = -a;347 return a;348 }350 int351 dogamma(Pix c)352 {353 float f;355 f = c - gam.min;356 if(f < 1)357 f = 1;359 if(gam.absgamma == 1)360 c = f * gam.mult2;361 else362 c = exp(log(f*gam.mult1) * gam.absgamma) * 255;363 if(c > 255)364 c = 255;365 if(gam.neg)366 c = 255-c;367 return c;368 }