13e12c5d1SDavid du Colombier #include <u.h> 23e12c5d1SDavid du Colombier #include <libc.h> 3*219b2ee8SDavid du Colombier #include <bio.h> 43e12c5d1SDavid du Colombier #include "sky.h" 53e12c5d1SDavid du Colombier 63e12c5d1SDavid du Colombier double PI_180 = 0.0174532925199432957692369; 73e12c5d1SDavid du Colombier double TWOPI = 6.2831853071795864769252867665590057683943387987502; 8*219b2ee8SDavid du Colombier double LN2 = 0.69314718055994530941723212145817656807550013436025; 9*219b2ee8SDavid du Colombier static double angledangle=(180./PI)*60.*60.*1000.; 103e12c5d1SDavid du Colombier 113e12c5d1SDavid du Colombier int 123e12c5d1SDavid du Colombier rint(char *p, int n) 133e12c5d1SDavid du Colombier { 143e12c5d1SDavid du Colombier int i=0; 153e12c5d1SDavid du Colombier 163e12c5d1SDavid du Colombier while(*p==' ' && n) 173e12c5d1SDavid du Colombier p++, --n; 183e12c5d1SDavid du Colombier while(n--) 193e12c5d1SDavid du Colombier i=i*10+*p++-'0'; 203e12c5d1SDavid du Colombier return i; 213e12c5d1SDavid du Colombier } 223e12c5d1SDavid du Colombier 233e12c5d1SDavid du Colombier DAngle 243e12c5d1SDavid du Colombier dangle(Angle angle) 253e12c5d1SDavid du Colombier { 263e12c5d1SDavid du Colombier return angle*angledangle; 273e12c5d1SDavid du Colombier } 283e12c5d1SDavid du Colombier 293e12c5d1SDavid du Colombier Angle 303e12c5d1SDavid du Colombier angle(DAngle dangle) 313e12c5d1SDavid du Colombier { 323e12c5d1SDavid du Colombier return dangle/angledangle; 333e12c5d1SDavid du Colombier } 343e12c5d1SDavid du Colombier 353e12c5d1SDavid du Colombier double 363e12c5d1SDavid du Colombier rfloat(char *p, int n) 373e12c5d1SDavid du Colombier { 383e12c5d1SDavid du Colombier double i, d=0; 393e12c5d1SDavid du Colombier 403e12c5d1SDavid du Colombier while(*p==' ' && n) 413e12c5d1SDavid du Colombier p++, --n; 423e12c5d1SDavid du Colombier if(*p == '+') 433e12c5d1SDavid du Colombier return rfloat(p+1, n-1); 443e12c5d1SDavid du Colombier if(*p == '-') 453e12c5d1SDavid du Colombier return -rfloat(p+1, n-1); 463e12c5d1SDavid du Colombier while(*p == ' ' && n) 473e12c5d1SDavid du Colombier p++, --n; 483e12c5d1SDavid du Colombier if(n == 0) 493e12c5d1SDavid du Colombier return 0.0; 503e12c5d1SDavid du Colombier while(n-- && *p!='.') 513e12c5d1SDavid du Colombier d = d*10+*p++-'0'; 523e12c5d1SDavid du Colombier if(n <= 0) 533e12c5d1SDavid du Colombier return d; 543e12c5d1SDavid du Colombier p++; 553e12c5d1SDavid du Colombier i = 1; 563e12c5d1SDavid du Colombier while(n--) 573e12c5d1SDavid du Colombier d+=(*p++-'0')/(i*=10.); 583e12c5d1SDavid du Colombier return d; 593e12c5d1SDavid du Colombier } 603e12c5d1SDavid du Colombier 613e12c5d1SDavid du Colombier int 623e12c5d1SDavid du Colombier sign(int c) 633e12c5d1SDavid du Colombier { 643e12c5d1SDavid du Colombier if(c=='-') 653e12c5d1SDavid du Colombier return -1; 663e12c5d1SDavid du Colombier return 1; 673e12c5d1SDavid du Colombier } 683e12c5d1SDavid du Colombier 693e12c5d1SDavid du Colombier char* 703e12c5d1SDavid du Colombier hms(Angle a) 713e12c5d1SDavid du Colombier { 723e12c5d1SDavid du Colombier static char buf[20]; 73*219b2ee8SDavid du Colombier double x; 743e12c5d1SDavid du Colombier int h, m, s, ts; 753e12c5d1SDavid du Colombier 76*219b2ee8SDavid du Colombier x=DEG(a)/15; 77*219b2ee8SDavid du Colombier x += 0.5/36000.; /* round up half of 0.1 sec */ 783e12c5d1SDavid du Colombier h = floor(x); 793e12c5d1SDavid du Colombier x -= h; 803e12c5d1SDavid du Colombier x *= 60; 813e12c5d1SDavid du Colombier m = floor(x); 823e12c5d1SDavid du Colombier x -= m; 833e12c5d1SDavid du Colombier x *= 60; 843e12c5d1SDavid du Colombier s = floor(x); 853e12c5d1SDavid du Colombier x -= s; 863e12c5d1SDavid du Colombier ts = 10*x; 873e12c5d1SDavid du Colombier sprint(buf, "%dh%.2dm%.2d.%ds", h, m, s, ts); 883e12c5d1SDavid du Colombier return buf; 893e12c5d1SDavid du Colombier } 903e12c5d1SDavid du Colombier 913e12c5d1SDavid du Colombier char* 923e12c5d1SDavid du Colombier dms(Angle a) 933e12c5d1SDavid du Colombier { 943e12c5d1SDavid du Colombier static char buf[20]; 95*219b2ee8SDavid du Colombier double x; 963e12c5d1SDavid du Colombier int sign, d, m, s, ts; 973e12c5d1SDavid du Colombier 98*219b2ee8SDavid du Colombier x = DEG(a); 993e12c5d1SDavid du Colombier sign='+'; 1003e12c5d1SDavid du Colombier if(a<0){ 1013e12c5d1SDavid du Colombier sign='-'; 1023e12c5d1SDavid du Colombier x=-x; 1033e12c5d1SDavid du Colombier } 104*219b2ee8SDavid du Colombier x += 0.5/36000.; /* round up half of 0.1 arcsecond */ 1053e12c5d1SDavid du Colombier d = floor(x); 1063e12c5d1SDavid du Colombier x -= d; 1073e12c5d1SDavid du Colombier x *= 60; 1083e12c5d1SDavid du Colombier m = floor(x); 1093e12c5d1SDavid du Colombier x -= m; 1103e12c5d1SDavid du Colombier x *= 60; 1113e12c5d1SDavid du Colombier s = floor(x); 1123e12c5d1SDavid du Colombier x -= s; 1133e12c5d1SDavid du Colombier ts = floor(10*x); 1143e12c5d1SDavid du Colombier sprint(buf, "%c%d°%.2d'%.2d.%d\"", sign, d, m, s, ts); 1153e12c5d1SDavid du Colombier return buf; 1163e12c5d1SDavid du Colombier } 1173e12c5d1SDavid du Colombier 1183e12c5d1SDavid du Colombier char* 119*219b2ee8SDavid du Colombier ms(Angle a) 120*219b2ee8SDavid du Colombier { 121*219b2ee8SDavid du Colombier static char buf[20]; 122*219b2ee8SDavid du Colombier double x; 123*219b2ee8SDavid du Colombier int d, m, s, ts; 124*219b2ee8SDavid du Colombier 125*219b2ee8SDavid du Colombier x = DEG(a); 126*219b2ee8SDavid du Colombier x += 0.5/36000.; /* round up half of 0.1 arcsecond */ 127*219b2ee8SDavid du Colombier d = floor(x); 128*219b2ee8SDavid du Colombier x -= d; 129*219b2ee8SDavid du Colombier x *= 60; 130*219b2ee8SDavid du Colombier m = floor(x); 131*219b2ee8SDavid du Colombier x -= m; 132*219b2ee8SDavid du Colombier x *= 60; 133*219b2ee8SDavid du Colombier s = floor(x); 134*219b2ee8SDavid du Colombier x -= s; 135*219b2ee8SDavid du Colombier ts = floor(10*x); 136*219b2ee8SDavid du Colombier if(d != 0) 137*219b2ee8SDavid du Colombier sprint(buf, "%d°%.2d'%.2d.%d\"", d, m, s, ts); 138*219b2ee8SDavid du Colombier else 139*219b2ee8SDavid du Colombier sprint(buf, "%.2d'%.2d.%d\"", m, s, ts); 140*219b2ee8SDavid du Colombier return buf; 141*219b2ee8SDavid du Colombier } 142*219b2ee8SDavid du Colombier 143*219b2ee8SDavid du Colombier char* 1443e12c5d1SDavid du Colombier hm(Angle a) 1453e12c5d1SDavid du Colombier { 1463e12c5d1SDavid du Colombier static char buf[20]; 147*219b2ee8SDavid du Colombier double x; 148*219b2ee8SDavid du Colombier int h, m, n; 1493e12c5d1SDavid du Colombier 150*219b2ee8SDavid du Colombier x = DEG(a)/15; 151*219b2ee8SDavid du Colombier x += 0.5/600.; /* round up half of tenth of minute */ 1523e12c5d1SDavid du Colombier h = floor(x); 1533e12c5d1SDavid du Colombier x -= h; 1543e12c5d1SDavid du Colombier x *= 60; 1553e12c5d1SDavid du Colombier m = floor(x); 156*219b2ee8SDavid du Colombier x -= m; 157*219b2ee8SDavid du Colombier x *= 10; 158*219b2ee8SDavid du Colombier n = floor(x); 159*219b2ee8SDavid du Colombier sprint(buf, "%dh%.2d.%1dm", h, m, n); 1603e12c5d1SDavid du Colombier return buf; 1613e12c5d1SDavid du Colombier } 1623e12c5d1SDavid du Colombier 1633e12c5d1SDavid du Colombier char* 1643e12c5d1SDavid du Colombier dm(Angle a) 1653e12c5d1SDavid du Colombier { 1663e12c5d1SDavid du Colombier static char buf[20]; 167*219b2ee8SDavid du Colombier double x; 1683e12c5d1SDavid du Colombier int sign, d, m, n; 1693e12c5d1SDavid du Colombier 170*219b2ee8SDavid du Colombier x = DEG(a); 1713e12c5d1SDavid du Colombier sign='+'; 1723e12c5d1SDavid du Colombier if(a<0){ 1733e12c5d1SDavid du Colombier sign='-'; 1743e12c5d1SDavid du Colombier x=-x; 1753e12c5d1SDavid du Colombier } 176*219b2ee8SDavid du Colombier x += 0.5/600.; /* round up half of tenth of arcminute */ 1773e12c5d1SDavid du Colombier d = floor(x); 1783e12c5d1SDavid du Colombier x -= d; 1793e12c5d1SDavid du Colombier x *= 60; 1803e12c5d1SDavid du Colombier m = floor(x); 1813e12c5d1SDavid du Colombier x -= m; 1823e12c5d1SDavid du Colombier x *= 10; 1833e12c5d1SDavid du Colombier n = floor(x); 1843e12c5d1SDavid du Colombier sprint(buf, "%c%d°%.2d.%.1d'", sign, d, m, n); 1853e12c5d1SDavid du Colombier return buf; 1863e12c5d1SDavid du Colombier } 187*219b2ee8SDavid du Colombier 188*219b2ee8SDavid du Colombier char* 189*219b2ee8SDavid du Colombier getword(char *ou, char *in) 190*219b2ee8SDavid du Colombier { 191*219b2ee8SDavid du Colombier int c; 192*219b2ee8SDavid du Colombier 193*219b2ee8SDavid du Colombier for(;;) { 194*219b2ee8SDavid du Colombier c = *in++; 195*219b2ee8SDavid du Colombier if(c == ' ' || c == '\t') 196*219b2ee8SDavid du Colombier continue; 197*219b2ee8SDavid du Colombier if(c == 0) 198*219b2ee8SDavid du Colombier return 0; 199*219b2ee8SDavid du Colombier break; 200*219b2ee8SDavid du Colombier } 201*219b2ee8SDavid du Colombier 202*219b2ee8SDavid du Colombier if(c == '\'') 203*219b2ee8SDavid du Colombier for(;;) { 204*219b2ee8SDavid du Colombier if(c >= 'A' && c <= 'Z') 205*219b2ee8SDavid du Colombier c += 'a' - 'A'; 206*219b2ee8SDavid du Colombier *ou++ = c; 207*219b2ee8SDavid du Colombier c = *in++; 208*219b2ee8SDavid du Colombier if(c == 0) 209*219b2ee8SDavid du Colombier return 0; 210*219b2ee8SDavid du Colombier if(c == '\'') { 211*219b2ee8SDavid du Colombier *ou = 0; 212*219b2ee8SDavid du Colombier return in-1; 213*219b2ee8SDavid du Colombier } 214*219b2ee8SDavid du Colombier } 215*219b2ee8SDavid du Colombier for(;;) { 216*219b2ee8SDavid du Colombier if(c >= 'A' && c <= 'Z') 217*219b2ee8SDavid du Colombier c += 'a' - 'A'; 218*219b2ee8SDavid du Colombier *ou++ = c; 219*219b2ee8SDavid du Colombier c = *in++; 220*219b2ee8SDavid du Colombier if(c == ' ' || c == '\t' || c == 0) { 221*219b2ee8SDavid du Colombier *ou = 0; 222*219b2ee8SDavid du Colombier return in-1; 223*219b2ee8SDavid du Colombier } 224*219b2ee8SDavid du Colombier } 225*219b2ee8SDavid du Colombier } 226*219b2ee8SDavid du Colombier 227*219b2ee8SDavid du Colombier /* 228*219b2ee8SDavid du Colombier * Read formatted angle. Must contain no embedded blanks 229*219b2ee8SDavid du Colombier */ 230*219b2ee8SDavid du Colombier Angle 231*219b2ee8SDavid du Colombier getra(char *p) 232*219b2ee8SDavid du Colombier { 233*219b2ee8SDavid du Colombier Rune r; 234*219b2ee8SDavid du Colombier char *q; 235*219b2ee8SDavid du Colombier Angle f, d; 236*219b2ee8SDavid du Colombier int neg; 237*219b2ee8SDavid du Colombier 238*219b2ee8SDavid du Colombier neg = 0; 239*219b2ee8SDavid du Colombier d = 0; 240*219b2ee8SDavid du Colombier while(*p == ' ') 241*219b2ee8SDavid du Colombier p++; 242*219b2ee8SDavid du Colombier for(;;) { 243*219b2ee8SDavid du Colombier if(*p == ' ') 244*219b2ee8SDavid du Colombier goto Return; 245*219b2ee8SDavid du Colombier if(*p == '-') { 246*219b2ee8SDavid du Colombier neg = 1; 247*219b2ee8SDavid du Colombier p++; 248*219b2ee8SDavid du Colombier } 249*219b2ee8SDavid du Colombier q = p; 250*219b2ee8SDavid du Colombier f = strtod(p, &q); 251*219b2ee8SDavid du Colombier if(q > p) { 252*219b2ee8SDavid du Colombier p = q; 253*219b2ee8SDavid du Colombier } 254*219b2ee8SDavid du Colombier p += chartorune(&r, p); 255*219b2ee8SDavid du Colombier switch(r) { 256*219b2ee8SDavid du Colombier default: 257*219b2ee8SDavid du Colombier Return: 258*219b2ee8SDavid du Colombier if(neg) 259*219b2ee8SDavid du Colombier d = -d; 260*219b2ee8SDavid du Colombier return RAD(d); 261*219b2ee8SDavid du Colombier case 'h': 262*219b2ee8SDavid du Colombier d += f*15; 263*219b2ee8SDavid du Colombier break; 264*219b2ee8SDavid du Colombier case 'm': 265*219b2ee8SDavid du Colombier d += f/4; 266*219b2ee8SDavid du Colombier break; 267*219b2ee8SDavid du Colombier case 's': 268*219b2ee8SDavid du Colombier d += f/240; 269*219b2ee8SDavid du Colombier break; 270*219b2ee8SDavid du Colombier case L'°': 271*219b2ee8SDavid du Colombier d += f; 272*219b2ee8SDavid du Colombier break; 273*219b2ee8SDavid du Colombier case '\'': 274*219b2ee8SDavid du Colombier d += f/60; 275*219b2ee8SDavid du Colombier break; 276*219b2ee8SDavid du Colombier case '\"': 277*219b2ee8SDavid du Colombier d += f/3600; 278*219b2ee8SDavid du Colombier break; 279*219b2ee8SDavid du Colombier } 280*219b2ee8SDavid du Colombier } 281*219b2ee8SDavid du Colombier return 0; 282*219b2ee8SDavid du Colombier } 283*219b2ee8SDavid du Colombier 284*219b2ee8SDavid du Colombier double 285*219b2ee8SDavid du Colombier xsqrt(double a) 286*219b2ee8SDavid du Colombier { 287*219b2ee8SDavid du Colombier 288*219b2ee8SDavid du Colombier if(a < 0) 289*219b2ee8SDavid du Colombier return 0; 290*219b2ee8SDavid du Colombier return sqrt(a); 291*219b2ee8SDavid du Colombier } 292*219b2ee8SDavid du Colombier 293*219b2ee8SDavid du Colombier Angle 294*219b2ee8SDavid du Colombier dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2) 295*219b2ee8SDavid du Colombier { 296*219b2ee8SDavid du Colombier double a; 297*219b2ee8SDavid du Colombier 298*219b2ee8SDavid du Colombier a = sin(dec1) * sin(dec2) + 299*219b2ee8SDavid du Colombier cos(dec1) * cos(dec2) * 300*219b2ee8SDavid du Colombier cos(ra1 - ra2); 301*219b2ee8SDavid du Colombier a = atan2(xsqrt(1 - a*a), a); 302*219b2ee8SDavid du Colombier if(a < 0) 303*219b2ee8SDavid du Colombier a = -a; 304*219b2ee8SDavid du Colombier return a; 305*219b2ee8SDavid du Colombier } 306*219b2ee8SDavid du Colombier 307*219b2ee8SDavid du Colombier int 308*219b2ee8SDavid du Colombier dogamma(Pix c) 309*219b2ee8SDavid du Colombier { 310*219b2ee8SDavid du Colombier float f; 311*219b2ee8SDavid du Colombier 312*219b2ee8SDavid du Colombier f = c - gam.min; 313*219b2ee8SDavid du Colombier if(f < 1) 314*219b2ee8SDavid du Colombier f = 1; 315*219b2ee8SDavid du Colombier 316*219b2ee8SDavid du Colombier if(gam.absgamma == 1) 317*219b2ee8SDavid du Colombier c = f * gam.mult2; 318*219b2ee8SDavid du Colombier else 319*219b2ee8SDavid du Colombier c = exp(log(f*gam.mult1) * gam.absgamma) * 255; 320*219b2ee8SDavid du Colombier if(c > 255) 321*219b2ee8SDavid du Colombier c = 255; 322*219b2ee8SDavid du Colombier if(gam.neg) 323*219b2ee8SDavid du Colombier c = 255-c; 324*219b2ee8SDavid du Colombier return c; 325*219b2ee8SDavid du Colombier } 326