[comp.sources.amiga] Star Plotting program

doc@s.cc.purdue.edu.UUCP (06/09/87)

	Here is a star-plotting program submitted by Darrin West
(ihnp4!calgary!west).  He seems to be anxious to hear of any bugs
you have, or any fixes/enhancements, so be sure to send them to
him.  If you can't get them to him, I will gladly forward them.
This program requires the use of some datafiles which can be found
in comp.binaries.amiga along with the executable for this program.
    -Craig Norborg


#	This is a shell archive.
#	Remove everything above and includeing the cut line.
#	Then run the rest of the file through sh.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# Xshar: Extended Shell Archiver.
#	Run the following text with /bin/sh to create:
#	ReadMe
#	ReadMe.2
#	plot6.c
#	wind.c
# This archive created: Tue Jun  9 13:41:13 1987
# By: Craig Norborg (Purdue University Computing Center)
cat << \SHAR_EOF > ReadMe
Each .bin file is a list of "ent"'s, with the first four bytes being the
number of "ent"'s in the file.  Files may be plotted in any order, but should
be plotted from dimmest to brightest.  This is accomplished by passing the
file names to plot6 from brightest to dimmest (I read the files from the
end first).  To get any kind of speed out of this thing, the .bin files
should realy be in ram.

NB: If you modify or recompile this program, you MUST use motorola ffp format,
since the data files are written out that way.  If you do not use the -f -LF
flags (I run Lattice), then you will get very bizzare looking output.
I don't think it will crash and burn though.

I plan to add menu input and maybe a superbitmap to scroll through, and will
probably change the plot_star routine yet again (thus we have plotSIX!),
but feel free to make your own changes, and have fun star watching.  Oh yah,
REAL astronomers may be grossed out by the fact that you enter your view in
degrees (fractions too!), so I may change that to hours, minutes, seconds etc.
I just multiply hours by 15 in my head.

Try looking at ra 85 dc 0.  You should see Orion.  Does it look familiar?

I actually have not pushed this software very hard, so bugs may proliferate.
Let me know if you like/dislike anything, and send bug reports.

The close window gadget stops plotting immediately, in case something has
gone wrong.  Changing the window size causes scaling to occur, so might
mess things up.  (I should probably change this!)

I have not concentrated on making the code small, but have tried to make
it fast.  There are some more things I might do to speed it up (eg remove
stack checking), but it is quite acceptable now.  You should have seen it
when it read the original textual data files with scanf!!!!!

Make sure you remove my .signature before you run the files through sh.
I hope everyone has access to uudecode.

The stars colours are reminicent of the (b-v)? values, but are a guess,
and are a bit too startling.  I have some equations that I will try to
use for this next time.

Maybe after I finish the thesis. {:-|]  <-- stoic, determined look?

Darrin West.


SHAR_EOF
cat << \SHAR_EOF > ReadMe.2
You must type "plot6 to5.bin to5.5.bin ...", listing each file name by
decreasing order of magnitude (increasing number size) in order for
small stars to be plotted first, and allowing larger stars to be
plotted over the smaller ones.

I did find a bug in the input data.  If the width field is greater than
90 degrees (ie 180 degree field of view [I know, I should change
that]), then you look at the southern hemisphere instead of the
northern (given that that is where you wanted to look).  So, for a very
wide field of view, type 80 or 89 degrees.  Less field of view is less
distorted.  I usually use about 60 degrees, as this usually includes
the entire constellation I am looking at.  However, zooming in is fun
too.  (Did you know that one of the stars on the big dipper's handle is
a pair of stars [but not a double]?)

My favorite constellation is Orion.  He can be fould at 83 degrees ra
-3 degrees dc, view him with about 30 degrees of field of view.  Very
interesting!  To see the big dipper type in "0 90 60", it had better be
center stage above center.

	[ moderators note: The prompt for the program looks like it wants
	commas inbetween the data you type in.  This is not so, if you do
	not put spaces instead of commas, it gives you a blank screen ]
SHAR_EOF
cat << \SHAR_EOF > plot6.c
#include <stdio.h>
#include <math.h>

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 <file-list>\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;i<t;i++) {
	    check_input();
	    star = ar+i;
	    
	    c = sin(star->dc);
	    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);

}
SHAR_EOF
cat << \SHAR_EOF > wind.c
#include <intuition/intuition.h>
#include <stdio.h>

struct GfxBase *GfxBase;
struct RastPort *rp;
extern struct IntuitionBase *IntuitionBase;

struct Screen *myscreen,*OpenScreen();
struct NewScreen newscr = {
    0,0,640,200,4,
    0,1,
    HIRES,
    CUSTOMSCREEN,
    NULL,
    "Stars' Screen!",
    NULL,NULL
};

struct Window *mywindow,*OpenWindow();
struct NewWindow newwind = {
    0,0,640,200,
    1,2,
    CLOSEWINDOW | NEWSIZE /* |  REFRESHWINDOW */,
    WINDOWSIZING | WINDOWDEPTH | WINDOWCLOSE | WINDOWDRAG |
    ACTIVATE,
    NULL, NULL,
    "Stars! (c) Darrin West.",
    NULL, /* changed later */
    NULL,
    10,10,640,200,
    CUSTOMSCREEN
};

unsigned short colours[] = {
    0, 0x55f, 0x338, 0x224,	/* blue */
       0xfff, 0x888, 0x444,	/* white */
       0xff5, 0x883, 0x442,	/* yellow */
       0xf55, 0x833, 0x422,	/* red */
    0x0f0, 0x0f0, 0x0f0		/*junk */
};
    
int MAX_X,MAX_Y;
int rewrite = 0;

int window_init()
{
    GfxBase = (struct GFxBase *) OpenLibrary("graphics.library",0);
    if (GfxBase == NULL) {
        printf("Can't find graphics.library\n");
	return(0);
    }
    
    IntuitionBase = (struct IntuitionBase *) OpenLibrary("intuition.library",0);
    if (IntuitionBase == NULL){
        printf("Can't find intuition.library\n");
	return(0);
    }
    
    myscreen = OpenScreen(&newscr);
    if (myscreen == NULL){
        printf("Can't open screen\n");
	return(0);
    }
    
    myscreen->ViewPort.ColorMap = GetColorMap(16);
    LoadRGB4(&(myscreen->ViewPort),colours,16);
    
    newwind.Screen = myscreen;
    
    mywindow = OpenWindow(&newwind);
    if (mywindow == NULL){
        printf("Can't open window\n");
	return(0);
    }
    MAX_X = mywindow->Width;
    MAX_Y = mywindow->Height;
    rp = mywindow->RPort;
}

window_deinit()
{
    CloseWindow(mywindow);
    CloseScreen(myscreen);
    CloseLibrary(GfxBase);
    CloseLibrary(IntuitionBase);
}

check_input()
{
    int class,code;
    struct IntuiMessage *message;

	while (message = GetMsg(mywindow->UserPort)){
	    class = message->Class;
	    code = message->Code;
	    ReplyMsg(message);
	    
	    switch(class) {
		case CLOSEWINDOW:
		    printf("close\n");
		    fflush(stdout);
		    window_deinit();
		    exit(0);
		    break;
		case NEWSIZE:
		    MAX_X = mywindow->Width;
		    MAX_Y = mywindow->Height;
		    printf("size: %d %d\n",MAX_X,MAX_Y);
		    fflush(stdout);
		    rewrite = 1;
		    break;
		default:
		    printf("other\n");
		    fflush(stdout);
		    break;
	    }
	}
}

wt()
{
    Wait(1<<mywindow->UserPort->mp_SigBit);
}

point(x,y,cl)
int x,y,cl;
{
    SetAPen(rp,cl);
    WritePixel(rp,x,y);
}

wind_clear()
{
    SetRast(rp,0);
}
SHAR_EOF