sys/src/cmd/astro/dist.c

   1 #include "astro.h"
   2
   3 double
   4 dist(Obj1 *p, Obj1 *q)
   5 {
   6         double a;
   7
   8         a = sin(p->decl2)*sin(q->decl2) +
   9                 cos(p->decl2)*cos(q->decl2)*cos(p->ra-q->ra);
  10         a = fabs(atan2(pyth(a), a)) / radsec;
  11         return a;
  12 }
  13
  14 int
  15 rline(int f)
  16 {
  17         char *p;
  18         int c;
  19         static char buf[1024];
  20         static int bc, bn, bf;
  21
  22         if(bf != f) {
  23                 bf = f;
  24                 bn = 0;
  25         }
  26         p = line;
  27         do {
  28                 if(bn <= 0) {
  29                         bn = read(bf, buf, sizeof(buf));
  30                         if(bn <= 0)
  31                                 return 1;
  32                         bc = 0;
  33                 }
  34                 c = buf[bc];
  35                 bn--; bc++;
  36                 *p++ = c;
  37         } while(c != '\n');
  38         return 0;
  39 }
  40
  41 double
  42 sunel(double t)
  43 {
  44         int i;
  45
  46         i = floor(t);
  47         if(i < 0 || i > NPTS+1)
  48                 return -90;
  49         t = osun.point[i].el +
  50                 (t-i)*(osun.point[i+1].el - osun.point[i].el);
  51         return t;
  52 }
  53
  54 double
  55 rise(Obj2 *op, double el)
  56 {
  57         Obj2 *p;
  58         int i;
  59         double e1, e2;
  60
  61         e2 = 0;
  62         p = op;
  63         for(i=0; i<=NPTS; i++) {
  64                 e1 = e2;
  65                 e2 = p->point[i].el;
  66                 if(i >= 1 && e1 <= el && e2 > el)
  67                         goto found;
  68         }
  69         return -1;
  70
  71 found:
  72         return i - 1 + (el-e1)/(e2-e1);
  73 }
  74
  75 double
  76 set(Obj2 *op, double el)
  77 {
  78         Obj2 *p;
  79         int i;
  80         double e1, e2;
  81
  82         e2 = 0;
  83         p = op;
  84         for(i=0; i<=NPTS; i++) {
  85                 e1 = e2;
  86                 e2 = p->point[i].el;
  87                 if(i >= 1 && e1 > el && e2 <= el)
  88                         goto found;
  89         }
  90         return -1;
  91
  92 found:
  93         return i - 1 + (el-e1)/(e2-e1);
  94 }
  95
  96 double
  97 solstice(int n)
  98 {
  99         int i;
 100         double d1, d2, d3;
 101
 102         d3 = (n*pi)/2 - pi;
 103         if(n == 0)
 104                 d3 += pi;
 105         d2 = 0.;
 106         for(i=0; i<=NPTS; i++) {
 107                 d1 = d2;
 108                 d2 = osun.point[i].ra;
 109                 if(n == 0) {
 110                         d2 -= pi;
 111                         if(d2 < -pi)
 112                                 d2 += pipi;
 113                 }
 114                 if(i >= 1 && d3 >= d1 && d3 < d2)
 115                         goto found;
 116         }
 117         return -1;
 118
 119 found:
 120         return i - (d3-d2)/(d1-d2);
 121 }
 122
 123 double
 124 betcross(double b)
 125 {
 126         int i;
 127         double d1, d2;
 128
 129         d2 = 0;
 130         for(i=0; i<=NPTS; i++) {
 131                 d1 = d2;
 132                 d2 = osun.point[i].mag;
 133                 if(i >= 1 && b >= d1 && b < d2)
 134                         goto found;
 135         }
 136         return -1;
 137
 138 found:
 139         return i - (b-d2)/(d1-d2);
 140 }
 141
 142 double
 143 melong(Obj2 *op)
 144 {
 145         Obj2 *p;
 146         int i;
 147         double d1, d2, d3;
 148
 149         d2 = 0;
 150         d3 = 0;
 151         p = op;
 152         for(i=0; i<=NPTS; i++) {
 153                 d1 = d2;
 154                 d2 = d3;
 155                 d3 = dist(&p->point[i], &osun.point[i]);
 156                 if(i >= 2 && d2 >= d1 && d2 >= d3)
 157                         goto found;
 158         }
 159         return -1;
 160
 161 found:
 162         return i - 2;
 163 }
 164
 165 #define NEVENT  100
 166 Event   events[NEVENT];
 167 Event*  eventp = 0;
 168
 169 void
 170 event(char *format, char *arg1, char *arg2, double tim, int flag)
 171 {
 172         Event *p;
 173
 174         if(flag & DARK)
 175                 if(sunel(tim) > -12)
 176                         return;
 177         if(flag & LIGHT)
 178                 if(sunel(tim) < 0)
 179                         return;
 180         if(eventp == 0)
 181                 eventp = events;
 182         p = eventp;
 183         if(p >= events+NEVENT) {
 184                 fprint(2, "too many events\n");
 185                 return;
 186         }
 187         eventp++;
 188         p->format = format;
 189         p->arg1 = arg1;
 190         p->arg2 = arg2;
 191         p->tim = tim;
 192         p->flag = flag;
 193 }
 194
 195 void
 196 evflush(void)
 197 {
 198         Event *p;
 199
 200         if(eventp == 0)
 201                 return;
 202         qsort(events, eventp-events, sizeof *p, evcomp);
 203         for(p = events; p<eventp; p++) {
 204                 if((p->flag&SIGNIF) && flags['s'])
 205                         print("ding ding ding ");
 206                 print(p->format, p->arg1, p->arg2);
 207                 if(p->flag & PTIME)
 208                         ptime(day + p->tim*deld);
 209                 print("\n");
 210         }
 211         eventp = 0;
 212 }
 213
 214 int
 215 evcomp(void *a1, void *a2)
 216 {
 217         double t1, t2;
 218         Event *p1, *p2;
 219
 220         p1 = a1;
 221         p2 = a2;
 222         t1 = p1->tim;
 223         t2 = p2->tim;
 224         if(p1->flag & SIGNIF)
 225                 t1 -= 1000.;
 226         if(p2->flag & SIGNIF)
 227                 t2 -= 1000.;
 228         if(t1 > t2)
 229                 return 1;
 230         if(t2 > t1)
 231                 return -1;
 232         return 0;
 233 }
 234
 235 double
 236 pyth(double x)
 237 {
 238
 239         x *= x;
 240         if(x > 1)
 241                 x = 1;
 242         return sqrt(1-x);
 243 }
 244
 245 char*
 246 skip(int n)
 247 {
 248         int i;
 249         char *cp;
 250
 251         cp = line;
 252         for(i=0; i<n; i++) {
 253                 while(*cp == ' ' || *cp == '\t')
 254                         cp++;
 255                 while(*cp != '\n' && *cp != ' ' && *cp != '\t')
 256                         cp++;
 257         }
 258         while(*cp == ' ' || *cp == '\t')
 259                 cp++;
 260         return cp;
 261 }