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