URI:
       tdist.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tdist.c (3488B)
       ---
            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 int        evcomp(const void*, const void*);
          196 
          197 void
          198 evflush(void)
          199 {
          200         Event *p;
          201 
          202         if(eventp == 0)
          203                 return;
          204         qsort(events, eventp-events, sizeof *p, evcomp);
          205         for(p = events; p<eventp; p++) {
          206                 if((p->flag&SIGNIF) && flags['s'])
          207                         print("ding ding ding ");
          208                 print(p->format, p->arg1, p->arg2);
          209                 if(p->flag & PTIME)
          210                         ptime(day + p->tim*deld);
          211                 print("\n");
          212         }
          213         eventp = 0;
          214 }
          215 
          216 int
          217 evcomp(const void *a1, const void *a2)
          218 {
          219         double t1, t2;
          220         Event *p1, *p2;
          221 
          222         p1 = (Event*)a1;
          223         p2 = (Event*)a2;
          224         t1 = p1->tim;
          225         t2 = p2->tim;
          226         if(p1->flag & SIGNIF)
          227                 t1 -= 1000.;
          228         if(p2->flag & SIGNIF)
          229                 t2 -= 1000.;
          230         if(t1 > t2)
          231                 return 1;
          232         if(t2 > t1)
          233                 return -1;
          234         return 0;
          235 }
          236 
          237 double
          238 pyth(double x)
          239 {
          240 
          241         x *= x;
          242         if(x > 1)
          243                 x = 1;
          244         return sqrt(1-x);
          245 }
          246 
          247 char*
          248 skip(int n)
          249 {
          250         int i;
          251         char *cp;
          252 
          253         cp = line;
          254         for(i=0; i<n; i++) {
          255                 while(*cp == ' ' || *cp == '\t')
          256                         cp++;
          257                 while(*cp != '\n' && *cp != ' ' && *cp != '\t')
          258                         cp++;
          259         }
          260         while(*cp == ' ' || *cp == '\t')
          261                 cp++;
          262         return cp;
          263 }