Blob
1 #include "astro.h"3 char* herefile;5 int6 main(int argc, char *argv[])7 {8 int i, j;9 double d;11 pi = atan(1.0)*4;12 pipi = pi*2;13 radian = pi/180;14 radsec = radian/3600;15 converge = 1.0e-14;17 startab = unsharp("#9/sky/estartab");18 herefile = unsharp("#9/sky/here");20 fmtinstall('R', Rconv);21 fmtinstall('D', Dconv);23 per = PER;24 deld = PER/NPTS;25 init();26 args(argc, argv);27 init();29 loop:30 d = day;31 pdate(d);32 if(flags['p'] || flags['e']) {33 print(" ");34 ptime(d);35 pstime(d);36 }37 print("\n");38 for(i=0; i<=NPTS+1; i++) {39 setime(d);41 for(j=0; objlst[j]; j++) {42 (*objlst[j]->obj)();43 setobj(&objlst[j]->point[i]);44 if(flags['p']) {45 if(flags['m'])46 if(strcmp(objlst[j]->name, "Comet"))47 continue;48 output(objlst[j]->name, &objlst[j]->point[i]);49 }50 }51 if(flags['e']) {52 d = dist(&eobj1->point[i], &eobj2->point[i]);53 print("dist %s to %s = %.4f\n", eobj1->name, eobj2->name, d);54 }55 /* if(flags['p']) { */56 /* pdate(d); */57 /* print(" "); */58 /* ptime(d); */59 /* print("\n"); */60 /* } */61 if(flags['p'] || flags['e'])62 break;63 d += deld;64 }65 if(!(flags['p'] || flags['e']))66 search();67 day += per;68 nperiods -= 1;69 if(nperiods > 0)70 goto loop;71 exits(0);72 return 0; /* gcc */73 }75 void76 args(int argc, char *argv[])77 {78 char *p;79 long t;80 int f, i;81 Obj2 *q;83 memset(flags, 0, sizeof(flags));84 ARGBEGIN {85 default:86 fprint(2, "astro [-adeklmopst] [-c nperiod] [-C tperiod]\n");87 exits(0);89 case 'c':90 nperiods = 1;91 p = ARGF();92 if(p)93 nperiods = atol(p);94 flags['c']++;95 break;96 case 'C':97 p = ARGF();98 if(p)99 per = atof(p);100 break;101 case 'e':102 eobj1 = nil;103 eobj2 = nil;104 p = ARGF();105 if(p) {106 for(i=0; q=objlst[i]; i++) {107 if(strcmp(q->name, p) == 0)108 eobj1 = q;109 if(strcmp(q->name1, p) == 0)110 eobj1 = q;111 }112 p = ARGF();113 if(p) {114 for(i=0; q=objlst[i]; i++) {115 if(strcmp(q->name, p) == 0)116 eobj2 = q;117 if(strcmp(q->name1, p) == 0)118 eobj2 = q;119 }120 }121 }122 if(eobj1 && eobj2) {123 flags['e']++;124 break;125 }126 fprint(2, "cant recognize eclipse objects\n");127 exits("eflag");129 case 'a':130 case 'd':131 case 'j':132 case 'k':133 case 'l':134 case 'm':135 case 'o':136 case 'p':137 case 's':138 case 't':139 flags[ARGC()]++;140 break;141 } ARGEND142 if(*argv){143 fprint(2, "usage: astro [-dlpsatokm] [-c nday] [-e obj1 obj2]\n");144 exits("usage");145 }147 t = time(0);148 day = t/86400. + 25567.5;149 if(flags['d'])150 day = readate();151 if(flags['j'])152 print("jday = %.4f\n", day);153 deltat = day * .001704;154 if(deltat > 32.184) /* assume date is utc1 */155 deltat = 32.184; /* correct by leap sec */156 if(flags['t'])157 deltat = readdt();159 if(flags['l']) {160 fprint(2, "nlat wlong elev\n");161 readlat(0);162 } else {163 f = open(herefile, OREAD);164 if(f < 0) {165 fprint(2, "%s?\n", herefile);166 /* btl mh */167 nlat = (40 + 41.06/60)*radian;168 awlong = (74 + 23.98/60)*radian;169 elev = 150 * 3.28084;170 } else {171 readlat(f);172 close(f);173 }174 }175 }177 double178 readate(void)179 {180 int i;181 Tim t;183 fprint(2, "year mo da hr min\n");184 rline(0);185 for(i=0; i<5; i++)186 t.ifa[i] = atof(skip(i));187 return convdate(&t);188 }190 double191 readdt(void)192 {194 fprint(2, "ΔT (sec) (%.3f)\n", deltat);195 rline(0);196 return atof(skip(0));197 }199 double200 etdate(long year, int mo, double day)201 {202 Tim t;204 t.ifa[0] = year;205 t.ifa[1] = mo;206 t.ifa[2] = day;207 t.ifa[3] = 0;208 t.ifa[4] = 0;209 return convdate(&t) + 2415020;210 }212 void213 readlat(int f)214 {216 rline(f);217 nlat = atof(skip(0)) * radian;218 awlong = atof(skip(1)) * radian;219 elev = atof(skip(2)) * 3.28084;220 }222 double223 fmod(double a, double b)224 {225 return a - floor(a/b)*b;226 }