#include #include struct ent{ float ra,dc; short mg,cl; }; extern int MAX_X,MAX_Y; /* 18.5cm/screen Y 26.2cm/screen X -> aspect ratio 1.4162 */ #define ASPECT 1.4162 /* MATH! spherical to rectangular coordinates. earth's axis is z axis dc measured from positive (north) z axis ra measured from x axis (we'll fix that later) All stars are on a unit sphere. x = sin(pi/2-dc)cos(ra) y = sin(pi/2-dc)sin(ra) z = cos(pi/2-dc) rotation rule (about origin [ie about z axis]) a' = a cos(th) - b sin(th) b' = a sin(th) + b cos(th) rotate to see desired line of longitude(ra0). ie line up y axis with chosen longitude (remember we were measuring from x) rotate counter clockwise pi/2, rotate clockwise ra0 about z (ie pi/2-ra0) x' = sin(pi/2-dc)cos(ra)cos(pi/2-ra0) - sin(pi/2-dc)sin(ra)sin(pi/2-ra0) = sin(pi/2-dc)cos(ra+pi/2-ra0) [believe me!] y' = sin(pi/2-dc)cos(ra)sin(pi/2-ra0) + sin(pi/2-dc)sin(ra)cos(pi/2-ra0) = sin(pi/2-dc)sin(ra+pi/2-ra0) z' = cos(pi/2-dc) now rotate down from the north pole along your longitude to desired latitude ie rotate down pi/2 then up by dc0 we are rotating around the x axis x" = sin(pi/2-dc)cos(ra+pi/2-ra0) y" = sin(pi/2-dc)sin(ra+pi/2-ra0)cos(pi/2-dc0)-cos(pi/2-dc)sin(pi/2-dc0) = cos(dc)sin(ra+pi/2-ra0)cos(pi/2-dc0)-sin(dc)sin(pi/2-dc0) z" = sin(pi/2-dc)sin(ra+pi/2-ra0)sin(pi/2-dc0)+cos(pi/2-dc)cos(pi/2-dc0) = cos(dc)sin(ra+pi/2-ra0)sin(pi/2-dc0) + sin(dc)cos(pi/2-dc0) I set dc0 to pi/2-dc0 and ra0 ot pi/2-ra0. a = sin(dc) b = cos(dc) c = sin(dc0) d = cos(dc0) e = sin(ra + ra0) f = cos(ra + ra0) x" = df; y" = deb-ca; z" = dea+cb; Now! viewing transforms. I assume we are at (0,0,0) looking at ra0 x dc0. We are looking down the z axis with y up. our famous Ys = d/s Ye/Ze view box transform gives us: (ie Y in screen coordinates from eye coordinates) Xs = df/(dea+cb); Ys = (deb-ca)/(dea+cb); first we discard any stars behind us (ie z < 0) if the star is more than half a screen to the side, discard it. NB: we need not calculate f until the last moment. divide x by the aspect ratio, then map to physical coordinates. the rest is a subjective star plotting mechanism */ #define mabs(x) (((x)<0.0)?-(x):(x)) #define mfloor(x) (((x)<0.0)?-(int)-(x):(int)(x)) main(argc,argv) int argc; char **argv; { FILE *fp; struct ent *ar,*star; int t,i,num; double a,b,c,d,e,f,de,ra0,de0,wi,xt,yt; double HLF; double z; double TO_RADS; if (argc <= 1) { printf("usage plot \n"); exit(0); } /* beware the ffp constant expression bug! */ TO_RADS = (TO_RADS=2.0) * 3.1415926535 / (TO_RADS=360.0); HLF = 90.0 * TO_RADS; printf("ra,de,wi [degrees]:");fflush(stdout); scanf("%lf %lf %lf",&ra0,&de0,&wi); ra0 = HLF-ra0*TO_RADS; de0 = HLF-de0*TO_RADS; wi = tan(mabs(wi)*TO_RADS); window_init(); a = sin(de0); b = cos(de0); while (argc-- > 1){ fp = fopen(argv[argc],"r"); if (fp == NULL){ printf("cant find %s\n",argv[argc]); exit(0); } t = fread((char*)&num,4,1,fp); if (t != 1) { printf("didn't read file size!\n"); exit(0); } printf("reading %d records\n",num); fflush(stdout); ar = (struct ent*)malloc(num*sizeof(struct ent)); if (ar == NULL){ printf("couldn't malloc %d bytes\n",num*sizeof(struct ent)); exit(0); } t = fread((char*)ar,sizeof(struct ent),num,fp); fclose(fp); printf("read %d\n",t); for (i = 0;idc); d = cos(star->dc); e = sin(star->ra + ra0); de = d * e; z = wi*(de * a + c * b); if (z <= 1.0e-16) continue; /* behind us, so skip it */ yt = (de * b - c * a) / z; if (yt < -0.5 || yt > 0.5) continue; f = cos(star->ra + ra0); xt = d * f / ASPECT / z; if (xt < -0.5 || xt > 0.5) continue; plot_star((int)((xt + 0.5)*MAX_X), (int)((yt + 0.5)*MAX_Y),star->mg,star->cl); } printf("done\n"); fflush(stdout); free (ar); } while (1){ wt(); check_input(); } } plot_star(x,y,mg,cl) int x,y; short mg,cl; { int c,colour,m; m = (int)mfloor(((double) mg)/100.0+0.49); c = (int)mfloor(cl/60.0+0.49) + 1; if (c < 0) c = 0; else if (c > 3) c = 3; c = 1 + 3 * c + 2;/*dimmest of range*/ colour = c-1;/* medium */ switch(m){ case 0: point(x-3,y,colour); point(x-4,y,colour); point(x+3,y,colour); point(x+4,y,colour); point(x+5,y,colour); point(x-5,y,colour); colour = c-2; point(x+2,y+1,colour); point(x+2,y-1,colour); point(x-2,y+1,colour); point(x-2,y-1,colour); case 1: point(x,y+2,colour); point(x,y-2,colour); colour = c-2; case 2: /* in bright colour */ point(x+3,y,colour); point(x-3,y,colour); point(x+2,y,colour); point(x-2,y,colour); colour = c-2; point(x+1,y+1,colour); point(x+1,y-1,colour); point(x-1,y+1,colour); point(x-1,y-1,colour); /* in brightest colour */ case 3: /* in dimmest colour */ colour = c-2; /* brighten colour */ case 4: /* in dimmest colour*/ point(x,y+1,colour); point(x,y-1,colour); point(x+1,y,colour); point(x-1,y,colour); point(x+2,y,colour); point(x-2,y,colour); case 5: colour = c-2;/* set to brightest */ break; case 6: colour = c-1; break; default: colour = c;/* in dimmest colour */ } point(x,y,colour); }