dots

Personal dotfiles
git clone git://git.gormless.xyz/dots.git
Log | Files | Refs

occ.c (3236B)


      1 #include "astro.h"
      2 
      3 Occ	 o1, o2;
      4 Obj2	 xo1, xo2;
      5 
      6 void
      7 occult(Obj2 *p1, Obj2 *p2, double d)
      8 {
      9 	int i, i1, N;
     10 	double d1, d2, d3, d4;
     11 	double x, di, dx, x1;
     12 
     13 	USED(d);
     14 
     15 	d3 = 0;
     16 	d2 = 0;
     17 	occ.t1 = -100;
     18 	occ.t2 = -100;
     19 	occ.t3 = -100;
     20 	occ.t4 = -100;
     21 	occ.t5 = -100;
     22 	for(i=0; i<=NPTS+1; i++) {
     23 		d1 = d2;
     24 		d2 = d3;
     25 		d3 = dist(&p1->point[i], &p2->point[i]);
     26 		if(i >= 2 && d2 <= d1 && d2 <= d3)
     27 			goto found;
     28 	}
     29 	return;
     30 
     31 found:
     32 	N = 2880*PER/NPTS; /* 1 min steps */
     33 	i -= 2;
     34 	set3pt(p1, i, &o1);
     35 	set3pt(p2, i, &o2);
     36 
     37 	di = i;
     38 	x = 0;
     39 	dx = 2./N;
     40 	for(i=0; i<=N; i++) {
     41 		setpt(&o1, x);
     42 		setpt(&o2, x);
     43 		d1 = d2;
     44 		d2 = d3;
     45 		d3 = dist(&o1.act, &o2.act);
     46 		if(i >= 2 && d2 <= d1 && d2 <= d3)
     47 			goto found1;
     48 		x += dx;
     49 	}
     50 	fprint(2, "bad 1 \n");
     51 	return;
     52 
     53 found1:
     54 	if(d2 > o1.act.semi2+o2.act.semi2+50)
     55 		return;
     56 	di += x - 3*dx;
     57 	x = 0;
     58 	for(i=0; i<3; i++) {
     59 		setime(day + deld*(di + x));
     60 		(*p1->obj)();
     61 		setobj(&xo1.point[i]);
     62 		(*p2->obj)();
     63 		setobj(&xo2.point[i]);
     64 		x += 2*dx;
     65 	}
     66 	dx /= 60;
     67 	x = 0;
     68 	set3pt(&xo1, 0, &o1);
     69 	set3pt(&xo2, 0, &o2);
     70 	for(i=0; i<=240; i++) {
     71 		setpt(&o1, x);
     72 		setpt(&o2, x);
     73 		d1 = d2;
     74 		d2 = d3;
     75 		d3 = dist(&o1.act, &o2.act);
     76 		if(i >= 2 && d2 <= d1 && d2 <= d3)
     77 			goto found2;
     78 		x += 1./120.;
     79 	}
     80 	fprint(2, "bad 2 \n");
     81 	return;
     82 
     83 found2:
     84 	if(d2 > o1.act.semi2 + o2.act.semi2)
     85 		return;
     86 	i1 = i-1;
     87 	x1 = x - 1./120.;
     88 	occ.t3 = di + i1 * dx;
     89 	occ.e3 = o1.act.el;
     90 	d3 = o1.act.semi2 - o2.act.semi2;
     91 	if(d3 < 0)
     92 		d3 = -d3;
     93 	d4 = o1.act.semi2 + o2.act.semi2;
     94 	for(i=i1,x=x1;; i++) {
     95 		setpt(&o1, x);
     96 		setpt(&o2, x);
     97 		d1 = d2;
     98 		d2 = dist(&o1.act, &o2.act);
     99 		if(i != i1) {
    100 			if(d1 <= d3 && d2 > d3) {
    101 				occ.t4 = di + (i-.5) * dx;
    102 				occ.e4 = o1.act.el;
    103 			}
    104 			if(d2 > d4) {
    105 				if(d1 <= d4) {
    106 					occ.t5 = di + (i-.5) * dx;
    107 					occ.e5 = o1.act.el;
    108 				}
    109 				break;
    110 			}
    111 		}
    112 		x += 1./120.;
    113 	}
    114 	for(i=i1,x=x1;; i--) {
    115 		setpt(&o1, x);
    116 		setpt(&o2, x);
    117 		d1 = d2;
    118 		d2 = dist(&o1.act, &o2.act);
    119 		if(i != i1) {
    120 			if(d1 <= d3 && d2 > d3) {
    121 				occ.t2 = di + (i+.5) * dx;
    122 				occ.e2 = o1.act.el;
    123 			}
    124 			if(d2 > d4) {
    125 				if(d1 <= d4) {
    126 					occ.t1 = di + (i+.5) * dx;
    127 					occ.e1 = o1.act.el;
    128 				}
    129 				break;
    130 			}
    131 		}
    132 		x -= 1./120.;
    133 	}
    134 }
    135 
    136 void
    137 setpt(Occ *o, double x)
    138 {
    139 	double y;
    140 
    141 	y = x * (x-1);
    142 	o->act.ra = o->del0.ra +
    143 		x*o->del1.ra + y*o->del2.ra;
    144 	o->act.decl2 = o->del0.decl2 +
    145 		x*o->del1.decl2 + y*o->del2.decl2;
    146 	o->act.semi2 = o->del0.semi2 +
    147 		x*o->del1.semi2 + y*o->del2.semi2;
    148 	o->act.el = o->del0.el +
    149 		x*o->del1.el + y*o->del2.el;
    150 }
    151 
    152 
    153 double
    154 pinorm(double a)
    155 {
    156 
    157 	while(a < -pi)
    158 		a += pipi;
    159 	while(a > pi)
    160 		a -= pipi;
    161 	return a;
    162 }
    163 
    164 void
    165 set3pt(Obj2 *p, int i, Occ *o)
    166 {
    167 	Obj1 *p1, *p2, *p3;
    168 	double a;
    169 
    170 	p1 = p->point+i;
    171 	p2 = p1+1;
    172 	p3 = p2+1;
    173 
    174 	o->del0.ra = p1->ra;
    175 	o->del0.decl2 = p1->decl2;
    176 	o->del0.semi2 = p1->semi2;
    177 	o->del0.el = p1->el;
    178 
    179 	a = p2->ra - p1->ra;
    180 	o->del1.ra = pinorm(a);
    181 	a = p2->decl2 - p1->decl2;
    182 	o->del1.decl2 = pinorm(a);
    183 	o->del1.semi2 = p2->semi2 - p1->semi2;
    184 	o->del1.el = p2->el - p1->el;
    185 
    186 	a = p1->ra + p3->ra - 2*p2->ra;
    187 	o->del2.ra = pinorm(a)/2;
    188 	a = p1->decl2 + p3->decl2 - 2*p2->decl2;
    189 	o->del2.decl2 = pinorm(a)/2;
    190 	o->del2.semi2 = (p1->semi2 + p3->semi2 - 2*p2->semi2) / 2;
    191 	o->del2.el = (p1->el + p3->el - 2*p2->el) / 2;
    192 }