Blob
1 #include "astro.h"3 double4 dist(Obj1 *p, Obj1 *q)5 {6 double a;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 }14 int15 rline(int f)16 {17 char *p;18 int c;19 static char buf[1024];20 static int bc, bn, bf;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 }41 double42 sunel(double t)43 {44 int i;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 }54 double55 rise(Obj2 *op, double el)56 {57 Obj2 *p;58 int i;59 double e1, e2;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;71 found:72 return i - 1 + (el-e1)/(e2-e1);73 }75 double76 set(Obj2 *op, double el)77 {78 Obj2 *p;79 int i;80 double e1, e2;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;92 found:93 return i - 1 + (el-e1)/(e2-e1);94 }96 double97 solstice(int n)98 {99 int i;100 double d1, d2, d3;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;119 found:120 return i - (d3-d2)/(d1-d2);121 }123 double124 betcross(double b)125 {126 int i;127 double d1, d2;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;138 found:139 return i - (b-d2)/(d1-d2);140 }142 double143 melong(Obj2 *op)144 {145 Obj2 *p;146 int i;147 double d1, d2, d3;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;161 found:162 return i - 2;163 }165 #define NEVENT 100166 Event events[NEVENT];167 Event* eventp = 0;169 void170 event(char *format, char *arg1, char *arg2, double tim, int flag)171 {172 Event *p;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 }195 int evcomp(const void*, const void*);197 void198 evflush(void)199 {200 Event *p;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 }216 int217 evcomp(const void *a1, const void *a2)218 {219 double t1, t2;220 Event *p1, *p2;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 }237 double238 pyth(double x)239 {241 x *= x;242 if(x > 1)243 x = 1;244 return sqrt(1-x);245 }247 char*248 skip(int n)249 {250 int i;251 char *cp;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 }