#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>  //for pow
//#include <vfw.h>  //for video
#include "avi.h"


#define VIDEO 0
#define NUMFRAMES 2800//60  // /30 for 1 second
#define FRAMESIZE 352*240*3
#define BMAPS 0
#define MAKEINDEX 0
#define CONSTVEL 1
#define GALAXY 0
#define TWOATOMS 0
#define PI 3.14159265

typedef struct _point3d
{
int a,b,d;
}point3d;

typedef struct _point3df
{
float a,b,d;
}point3df;

typedef struct _matrix3df
{
float a[3];
float b[3];
float d[3];
}matrix3df;

typedef struct tagAVIHEAD
{
unsigned long riff;
unsigned long fsize;
unsigned long avi;
unsigned long list;
unsigned long listlen; //?
unsigned long hdrl;
unsigned long avih;
//unsigned long more;  //?
}AVIHEAD;

typedef struct tagLISTSTR
{
unsigned long list;
unsigned long listlen;
unsigned long strl; //strl
unsigned long strh; //strh
unsigned long strhlen;
}LISTSTR;

typedef struct tagmovi
{
unsigned long list;  //LIST
unsigned long listsize;
unsigned long  movi;   //movi
}MOVI;

#if 0 
struct _aviindex_chunk {
FOURCC fcc;
DWORD cb;
WORD wLongsPerEntry; // size of each entry in aIndex array
BYTE bIndexSubType; // future use. must be 0
BYTE bIndexType; // one of AVI_INDEX_* codes
DWORD nEntriesInUse; // index of first unused member in aIndex array
DWORD dwChunkId; // fcc of what is indexed
DWORD dwReserved[3]; // meaning differs for each index
// type/subtype. 0 if unused
struct _aviindex_entry {
DWORD adw[wLongsPerEntry];
} aIndex[ ];
}AVISTDINDEX, * PAVISTDINDEX;
#endif

AVIINDEXENTRY avindex;


#if 0 
unsigned long parm[8];
unsigned long numframes;
		unsigned long caps;
		unsigned long streams;
		unsigned long framesize;
		unsigned long width;
		unsigned long height;
		unsigned long scale,rate;
		unsigned long dwlength;
		unsigned long editcount;
		unsigned long list2;
		unsigned long parm2[10];
		unsigned long fps;
		unsigned long fps0;
		unsigned long numframes2;
		unsigned long framesize2;
		unsigned long parm3[9];
		unsigned long width2;
		unsigned long height2;
		unsigned long weird1;
		unsigned long weird2;
		unsigned long framesize3;
		unsigned long weird3[4];
//		unsigned long junk;
//		unsigned long junklen;
//		unsigned long zero[964];
		unsigned long list3;
		unsigned long fsize1;
		unsigned long movi;
	}AVIHEAD;
#endif
AVIHEAD head;
MainAVIHeader head2;
LISTSTR head3;
AVIStreamHeader head4;
BITMAPINFOHEADER bmpihead;
unsigned long strf;
MOVI head5; 

#if 0 
typedef struct _BMPFRAME
{
unsigned int tag;
unsigned int framesize'
}BMPFRAME;
#endif

char head7[8]={ 0x30, 0x30, 0x64, 0x62,FRAMESIZE,FRAMESIZE>>8,FRAMESIZE>>16,FRAMESIZE>>24};
int curframe;

//prototypes
LRESULT CALLBACK WndProc (HWND, UINT, WPARAM, LPARAM);
void draw(void);
void init_points(void);
void transform_draw_points(void);
void update_screen(void);
void move_points(void);
void pmmult(point3df p1,point3df *p2,matrix3df m);
void mrotb(matrix3df *m, float a);
void mmult(matrix3df m1, matrix3df m2,matrix3df *m3);
void angle2matrix(point3df viewera,matrix3df *trot);


LPBYTE lpBuf,lpScreen,lpChr[4];
LPBITMAPINFO lpDIB;
LPRGBQUAD lpRGB;
HDC hdcWin,hdcMem;
HWND hwMain;
HBITMAP hBMP;
BYTE *screen;//[240*3][352];//[352][240*3];
DWORD dwTime,dwPreTime;

int DW=352;
int DH=240;
int DW2=176;
int DH2=120;
int numpoints;
int a,b,d;
point3df p[10000];
point3df v[10000];
point3df p2[10000];
point3df v2[10000];
point3df viewer;
float scale;
float gravity;

FILE *fptr,*fptr2;
matrix3df viewerm,ident;
point3df viewera;
float direk;  //amount of influence current direction has for constveloc
float cveloc;   //constant velocity of photons/pieces of matter
int dexp=2.0;  //distance exponent 2=distance squared
int bmnum;
char tstr[255],bmfile[255],bfn[20];
BITMAPFILEHEADER bmphead;
char lfile[255];
//unsigned char aviblock[515];
unsigned int aviblock[128];
int track;
int trackpoint;


//WINAPI
int PASCAL WinMain (HINSTANCE hInstance, HINSTANCE hPrevInstance,
                    PSTR szCmdLine, int iCmdShow){

	MSG  msg;
	WNDCLASS  wc ;

//	wc.cbSize        = sizeof(wndclass);
	wc.style         = CS_HREDRAW | CS_VREDRAW;
	wc.lpfnWndProc   = (WNDPROC)WndProc;
	wc.cbClsExtra    = 0;
	wc.cbWndExtra    = 0;
	wc.hInstance     = hInstance;
	wc.hIcon         = NULL;
	wc.hCursor       = LoadCursor(NULL,IDC_ARROW);
	wc.hbrBackground = (HBRUSH)GetStockObject(BLACK_BRUSH);
	wc.lpszMenuName  = NULL;
	wc.lpszClassName = "Draw3D";  //"CWindow";
	wc.hIcon       = NULL;
	RegisterClass (&wc);


    hwMain = CreateWindowEx(
        WS_EX_APPWINDOW,
        "Draw3D",
        "Draw3D",
        WS_OVERLAPPED | WS_CAPTION | WS_SYSMENU | 
		WS_THICKFRAME | WS_MINIMIZEBOX | WS_MAXIMIZEBOX,
		CW_USEDEFAULT, CW_USEDEFAULT,
        DW,//WinX,640,//GetSystemMetrics(SM_CXSCREEN), //100
        DH,//WinY,480,//GetSystemMetrics(SM_CYSCREEN),  //100
        NULL,
        NULL,
        hInstance,
        NULL );

	screen=(unsigned char *)malloc(DW*DH*3);

	init_points();


#if VIDEO
//open the video file
	fptr=fopen("c:\\temparea\\test.avi","wb");
//add avi header
	head.riff=0x46464952;
	head.fsize=0x0;
  head.fsize=sizeof(AVIHEAD)+
		sizeof(MainAVIHeader)+
		sizeof(LISTSTR)+sizeof(AVIStreamHeader)+4+
		sizeof(BITMAPINFOHEADER)+sizeof(MOVI)+4+FRAMESIZE*NUMFRAMES;
	//this number is 0x738 (1848) too large
  //head.fsize=0x7ecf0;


	head.avi=0x20495641;
	head.list=0x5453494c;
  head.listlen=0xc0;//0xc0; //0x16e;  //0xc0;  //what is this? got to be c0!
	head.hdrl = 0x6c726468;  //hdrl
	head.avih = 0x68697661;  //avih
  //head.more=0;

	fwrite(&head,sizeof(AVIHEAD),1,fptr);

  head2.dwMicroSecPerFrame	=0x38;
	head2.dwMaxBytesPerSec =0x1046b;//0x8235;//0x30D40;//0x1046b;
	head2.dwPaddingGranularity =0;//0xD160;//0;//dwReserved1
	head2.dwFlags =0;//0x810; //0;  //?
  head2.dwTotalFrames =0x810;//NUMFRAMES;
  head2.dwInitialFrames =NUMFRAMES;//0;//0x5f;//0;
	
	head2.noidea = 0;
	
	head2.dwStreams = 1;
  head2.dwSuggestedBufferSize = FRAMESIZE;
	head2.dwWidth =DW;
	head2.dwHeight =DH;
	
	head2.dwScale =0;//0xbff6241a;
  head2.dwRate = 0;//0xbff6241a;
  head2.dwStart = 0;
  head2.dwLength = 0; //?
  //head2.dwReserved2 =0;
	//memset(head2.dwReserved,0,32);

	fwrite(&head2,sizeof(MainAVIHeader),1,fptr);

  head3.list = 	0x5453494c;
  head3.listlen = sizeof(LISTSTR)+sizeof(AVIStreamHeader)+sizeof(bmpihead);
									//this LIST minus LIST and len plus streamheader + strf + 0x28 + bitmapinfoheader 
	
  head3.strl = 0x6c727473;  //strl
	head3.strh= 0x68727473; //strh
  head3.strhlen = 0x38;//0x40;

	fwrite(&head3,sizeof(LISTSTR),1,fptr);

#if 0 
AVISTREAMINFO headt;

	//stream header
  headt.fccType = 0x73646976;  //vids
  headt.fccHandler =  0x6d737663; //DIB 0x20424944; //msvc 
  headt.dwFlags = 0;
	headt.dwCaps =0;
	headt.wPriority =0;
	headt.wLanguage=0;
	headt.dwScale =1;
	headt.dwRate = 30;  //fps
  headt.dwStart = 0;
	headt.dwLength=NUMFRAMES;
	headt.dwInitialFrames=0;
	headt.dwSuggestedBufferSize = FRAMESIZE;
  headt.dwQuality =0;
	headt.dwSampleSize =0;
	
RECT rctB = {0,0,352,240};

	headt.rcFrame=rctB;
  //headt.EditCount=0;
	//headt.dwFormatChangeCount=0;
	//memcpy(headt.szName,0,32);


	//fwrite(&headt,sizeof(AVISTREAMINFO),1,fptr);

#endif

//
//  head4.dwReserved2 =DW;//0;
//  head4.dwReserved3 =DH;//0;
  //head4.dwReserved4=0;

//#if 0 
	//stream header
  head4.fccType = 0x73646976;  //vids
  //head4.fccHandler = 0x20424944; //DIB 0x6376736d;//msvc 
	head4.fccHandler = 0;//tph - 0x20424944; //DIB 0x6376736d;//msvc 
  head4.dwFlags = 0;
	head4.dwReserved1 =0;
	head4.dwInitialFrames =0;
	head4.dwScale =1;
	head4.dwRate = 30;  //fps
  head4.dwStart = 0;
	head4.dwLength=NUMFRAMES;
	head4.dwSuggestedBufferSize = FRAMESIZE;
  head4.dwQuality =0;
	head4.dwSampleSize =0;
  //head4.dwReserved2 =DW;//0;
  //head4.dwReserved3 =DH;//0;
	// RECT rctB = {0,0,DW,DH};
//	head4.rcFrame=rctB;
  head4.sx =0;
	head4.sy =0;
	head4.ex =DW;
	head4.ey =DH;

  //head4.dwReserved4=0;
	
	fwrite(&head4,sizeof(AVIStreamHeader),1,fptr);
//#endif
	
	//stream format
  strf=0x66727473;  //strf
	fwrite(&strf,4,1,fptr);

	memset(&bmpihead,0,sizeof(bmpihead));
bmpihead.biSize=0x28;
bmpihead.biWidth=DW;
bmpihead.biHeight=DH;
bmpihead.biPlanes=1;
bmpihead.biBitCount=24;
bmpihead.biCompression=0;
bmpihead.biSizeImage=DW*DH*3;//0;
bmpihead.biXPelsPerMeter=0;//0xb12;
bmpihead.biYPelsPerMeter=0;//0xb12;
bmpihead.biClrUsed=0;//0x100;
bmpihead.biClrImportant=0;//0x100;

fwrite(&bmpihead.biSize,4,1,fptr);  //0x28 = size of header?
fwrite(&bmpihead,sizeof(BITMAPINFOHEADER),1,fptr);
	

//add JUNK
unsigned int junk;
unsigned char dt;
int junksize;
junk=0x4b4e554a;
dt=0;
junksize=1816;  //0x718

fwrite(&junk,4,1,fptr);
fwrite(&junksize,4,1,fptr);
for (a=0;a<1816;a++)  //should be a 0x7f4  (2048)
	{
	fwrite(&dt,1,1,fptr);
	}

//END AVI HEADER


	head5.list = 0x5453494c;  //LIST
	head5.listsize = 4+sizeof(head7)+FRAMESIZE*NUMFRAMES;//0;//head7=8  //?was 12+
	head5.movi = 0x69766f6d;   //movi

		
		//NUMFRAMES*FRAMESIZE+
		//NUMFRAMES*sizeof(head7)-8;

		fwrite(&head5,sizeof(MOVI),1,fptr);



		curframe=0;
#endif


	ShowWindow (hwMain,iCmdShow); 
	UpdateWindow (hwMain);
	SetFocus( hwMain);


memset(screen,0,DH*DW*3);


		hdcWin=GetDC(hwMain);
	  hdcMem = CreateCompatibleDC(hdcWin);  //NULL);



	while (1) { 

		if (PeekMessage (&msg,NULL,0,0,PM_NOREMOVE)) {

			if (!GetMessage (&msg,NULL,0,0)) 
				return msg.wParam ;

			TranslateMessage(&msg);
			DispatchMessage(&msg);

		} 
			transform_draw_points();  
			move_points();
		  update_screen();

	}

#if BMAPS
  //fclose(fptr2);
#endif

}

void move_points(void)
{
float dist,dist2;
float da,db,dd;
float denom,denom2;

  for (a=0;a<numpoints;a++)
		{
  
		v2[a].a = 0;
		v2[a].b = 0;
		v2[a].d = 0;
		for (b=0;b<numpoints;b++)
			{
			if (b != a)
			{
      da = p[b].a - p[a].a;
      db = p[b].b - p[a].b;
      dd = p[b].d - p[a].d;


			dist = da*da+db*db+dd*dd;
			dist = pow(dist,0.5);

      dist = pow(dist,dexp);

      if (dist < 1.0) {dist = 1.0;}
      //If (dist = 0) Then
       //this is bounce off each other
      // momen2(b).a = -1 * momen2(b).a
      //  momen2(b).b = -1 * momen2(b).b
      //  momen2(b).d = -1 * momen2(b).d
        
      //this is stuck to each other
      //momen(b).a = 0
       //momen(b).b = 0
      //Else
//        momen2(b).a = (points(a).a - points(b).a) / dist
//        momen2(b).b = (points(a).b - points(b).b) / dist
//        momen2(b).d = (points(a).d - points(b).d) / dist

#if CONSTVEL
			  //sum each normalized influence (for example .1 .2 .7 = 1.0)
			  //of all pieces of matter on 1 piece
			//this final velocity will be > 1.0 because the influences are added
//			  dist=dist*dist;  //for inverse distance squared
			  denom=fabs(da)+fabs(db)+fabs(dd);
				if (denom==0) {denom=1.0;}
        v2[a].a = v2[a].a + (da / denom)/dist;
        v2[a].b = v2[a].b + (db / denom)/dist;
        v2[a].d = v2[a].d + (dd / denom)/dist;

				//more distance between 2 photons = less influence

#else

        //dist2 = dist * dist;
      
        v2[a].a = v2[a].a + (da / dist);
        v2[a].b = v2[a].b + (db / dist);
        v2[a].d = v2[a].d + (dd / dist);
#endif
      
      //'End If
			}  //end if different photons
		}  //end for b
	} //end for a

for (a=0;a<numpoints;a++)
	{  
#if CONSTVEL
//the influence of photons compared to 
//current distribution of velocity should not change 
	//depending on number of photons
	v2[a].a/=numpoints;  
  v2[a].b/=numpoints;
  v2[a].d/=numpoints;


		//normalize normalized sum so final overall influence = 
	// (a + b + d = 1.0)  1.0 is the velocity the photon must move
	denom=fabs(v2[a].a) + fabs(v2[a].b) + fabs(v2[a].d);

  

#if 1
	//current direction has more effect than new direction

	denom2=(1-direk)*denom+direk*(fabs(v[a].a) +fabs(v[a].b) +fabs(v[a].d)); 
	if (denom2==0) {denom2=1.0;}
	denom2/=cveloc;  //will make base velocity = cveloc
//  v[a].a = direk*(v[a].a/denom2)+(1-direk)*v2[a].a/denom;
//  v[a].b = direk*(v[a].b/denom2)+(1-direk)*v2[a].b/denom;
//  v[a].d = direk*(v[a].d/denom2)+(1-direk)*v2[a].d/denom;
	//take some part of the current velocity and the new velocity
	//normalized for a third time to = 1.0
  v[a].a = (direk*v[a].a+(1-direk)*v2[a].a)/denom2;
  v[a].b = (direk*v[a].b+(1-direk)*v2[a].b)/denom2;
  v[a].d = (direk*v[a].d+(1-direk)*v2[a].d)/denom2;

#else
  //only other matter effects direction
	//if (denom==0) {denom=1.0;}

  v[a].a = v2[a].a/denom;
  v[a].b = v2[a].b/denom;
  v[a].d = v2[a].d/denom;
#endif


#else
  v[a].a = v[a].a + v2[a].a/gravity;
  v[a].b = v[a].b + v2[a].b/gravity;
  v[a].d = v[a].d + v2[a].d/gravity;
#endif
  
  p[a].a = p[a].a + v[a].a;
  p[a].b = p[a].b + v[a].b;
  p[a].d = p[a].d + v[a].d;
	}  //end for a  
   

//move viewer if tracking a point
  if (track)
		{
  viewer.a = p[trackpoint].a;
	viewer.b = p[trackpoint].b;
	//viewer.d = p[0].d;
	}


}

void update_screen(void)
{
int a,b;
	
	char byt;
	byt=0xff;

	
#if VIDEO

  if (curframe<NUMFRAMES)
	{
#if 0 
		//add the bitmap to the video file
			fwrite(&head7,8,1,fptr);   

			//also need to put bitmap header?
fwrite(&bmpihead,sizeof(BITMAPINFOHEADER),1,fptr);
#endif

//added tph
//add the bitmap to the video file
		  fwrite(&head7,8,1,fptr);   

			//also need to put bitmap header?
      fwrite(&bmpihead,sizeof(BITMAPINFOHEADER),1,fptr);

			fwrite(screen,FRAMESIZE,1,fptr);

#if 0 
			for(b=DH-1;b>-1;b--)
			  {
				for(a=0;a<DW*3;a++)
				{
					byt=*(screen+b*DW*3+a);
				fwrite(&byt,1,1,fptr);  
//screen+a+b*DW*3;        
	//			fwrite(&byt,1,1,fptr);  
				}
			  }
#endif
			curframe++;
	}
	else
		{
		if (curframe==NUMFRAMES)
			{	
#if MAKEINDEX
unsigned int temp;

//write index

memset(aviblock,0,128);

avindex.ckid = 0x31786469;  //idx1

//if NUMFRAMES>1 add Chunklength=0x20
//for each frame
for (a=1;a<NUMFRAMES;a++)
{
avindex.dwChunkLength = 0x10;  //
avindex.dwChunkOffset =0x62643030; //00db 
avindex.dwFlags =0x20;//0x10;  //



	fwrite(&avindex,sizeof(AVIINDEXENTRY),1,fptr);
	temp=4;
	fwrite(&temp,4,1,fptr);
	temp=FRAMESIZE;
	fwrite(&temp,4,1,fptr);
}
	fwrite(&aviblock,512,1,fptr);  //128 32 bits=512 bytes
#endif //if MAKEINDEX

	//close the video file
  fclose(fptr);	


			MessageBox(hwMain,"Done making movie","Movie Done",0);
			curframe++;
			}
		}
#endif




	  hBMP = CreateBitmap (DW, DH, 1, 24, screen) ;
    SelectObject(hdcMem, hBMP);

    StretchBlt(hdcWin,0,0,DW,DH,hdcMem,0,0,DW,DH,SRCCOPY);


#if BMAPS
//strcpy(bmfile,"images\\bm");
strcpy(bmfile,"images\\");
//itoa(bmnum,tstr,10);
//strcat(bmfile,tstr);
strcat(bmfile,bfn);
strcat(bmfile,".bmp");

//97 to 122
//add filname
if (bfn[7]==122)
	{
	if (bfn[6]==122)
		{
		if (bfn[5]==122)
			{
			bfn[7]=97;
			bfn[6]=97;
			bfn[5]=97;
			bfn[4]++;
			}
		bfn[7]=97;
		bfn[6]=97;
		bfn[5]++;
		}
	bfn[7]=97;
	bfn[6]++;
	}
else
	{
	bfn[7]++;
	}


//fwrite(&bmfile,1,sizeof(bmfile),fptr2);  


fptr=fopen(bmfile,"wb");
fwrite(&bmphead,sizeof(BITMAPFILEHEADER),1,fptr);
fwrite(&bmpihead,sizeof(BITMAPINFOHEADER),1,fptr);
//fwrite(&screen,0,DW*DH*3,fptr);

			for(b=DH-1;b>-1;b--)
			  {
				for(a=0;a<DW*3;a++)
				{
					byt=*(screen+b*DW*3+a);
				fwrite(&byt,1,1,fptr);  
//screen+a+b*DW*3;        
	//			fwrite(&byt,1,1,fptr);  
				}
			  }

  fclose(fptr);

		
bmnum++;

#endif

//		InvalidateRect(hwMain,NULL,FALSE);
		UpdateWindow (hwMain);  

	  DeleteObject(hBMP);


}//end update_screen


void init_points(void)
{
float x,y,z,xi,yi,zi;
int np2,np3;
int xs,ys,xm,ym,zm,zs;
float denom;
float xwidth,ywidth,zwidth;
float xmin,ymin,zmin,xmax,ymax,zmax;


track=0;

ident.a[0]=1.0;  ident.a[1]=0.0;  ident.a[2]=0.0;
ident.b[0]=0.0;  ident.b[1]=1.0;  ident.b[2]=0.0;
ident.d[0]=0.0;  ident.d[1]=0.0;  ident.d[2]=1.0;

//viewerm=ident;

viewer.a = 0.0;  viewer.b = 0.0; viewer.d = -8000.0;
//viewera.a =0.0;  viewera.b =PI/2.0;  viewera.d =0;//viewer angle
viewera.a =0.0;  viewera.b =0.0;  viewera.d =0;//viewer angle
angle2matrix(viewera,&viewerm);


#if CONSTVEL
//viewer.a = -8000.0;  viewer.b = 0.0; viewer.d = 0.0;
//viewera.a =0.0;  viewera.b =PI/2.0;  viewera.d =0.0;//viewer angle
//.999 sends photons out
//.99 keeps photons in (!)  .99*old direction  .01*new direction
//direk=.5;  //for cveloc=1.0 
//.0005 10.0 400=interesting (nucleus)
//.0005  1.0 400=interesting (galaxies)
//.01 1.0 30=fall to a certain distance and stop
//.98 1.0 30 1.0 .5  = looks like gravity  square root inv dist
//.38 1.0 10 1.0 1.0 = looks kind of like gravity inv dist
//.12 1.0 50 1.0 2.0 = eventually may look like gravity
//above is used for movie
//interesting watch little cluster - is slower than other single photons 
//because the more photons tangled the slower the velocity of the entire group
//then one photon is thrown out- gets free
direk=0.12;//.000298;//.1;  //more small more attraction
//0= no momentum of direction
cveloc=1.0;  //constant velocity of photons
numpoints=200;//3;//400;
gravity=1.0;  //more gravity - more clumping new vel:old vel
dexp=2.0;
//np2=(int)pow(numpoints,0.5);
np2=(int)pow(numpoints,0.333);

xwidth=200.0;
ywidth=200.0;
zwidth=200.0;


xmin=-xwidth/2;
ymin=-ywidth/2;
zmin=-zwidth/2;
xmax=xwidth/2;
ymax=ywidth/2;
zmax=zwidth/2;

x=xmin;
y=ymin;
z=zmin;
//x=-100;  
//y=-100;
//z=-100;

//xi=(float)200/np2;
//yi=(float)200/np2;
//zi=(float)200/np2;

xi=xwidth/np2;
yi=ywidth/np2;
zi=zwidth/np2;


for (a=0;a<numpoints;a++)
	{
	x+=xi;
	//z+=zi;
	//if (x>100) {x=-100; y+=yi;}
	if (x>xmax) {x=xmin; y+=yi;}

	//if (z>100) {z=-100; y+=yi;}
	//if (y>100) 
		//{x=-100; y=-100; z+=zi;}
	if (y>ymax) 
		{x=xmin; y=ymin; z+=zi;}
	//if (y>100) 
		//{z=-100; y=-100; x+=xi;}
  p[a].a = -x;
	p[a].b = y;
  p[a].d = z;

	denom=abs(x)+abs(y)+abs(z);
 	v[a].a=0;//(x*2)/denom;
	v[a].b=0;//(y*2)/denom;
	v[a].d=0;//(z*2)/denom;

	}

p[numpoints-1].a = -50000.0;

#if 0 //line of points to left
//numpoints=500;
for (a=numpoints;a<numpoints+100;a++)
	{
  p[a].a = -800;
	p[a].b = 0;
  p[a].d = a*2;

 	v[a].a=0.0;
	v[a].b=0.0;
	v[a].d=0.0;

}
numpoints=numpoints+100;
#endif

#endif  //CONSTVEL


#if GALAXY
//viewer.a = -8000.0;  viewer.b = 0.0; viewer.d = 0.0;
//viewera.a =0.0;  viewera.b =PI/2.0;  viewera.d =0.0;//viewer angle

numpoints=400;
gravity=1;  //more gravity - more clumping
//np2=(int)pow(numpoints,0.5);
np2=(int)pow(numpoints,0.333);
x=-100;  
y=-100;
z=-100;

xi=(float)200/np2;
yi=(float)200/np2;
zi=(float)200/np2;

for (a=0;a<numpoints;a++)
	{
	x+=xi;
	if (x>100) {x=-100; y+=yi;}
	if (y>100) 
		{x=-100; y=-100; z+=zi;}
  p[a].a = -x;
	p[a].b = y;
  p[a].d = z;

	v[a].a=0.0;
	v[a].b=0.0;
	v[a].d=0.0;

	v[a].a=-y/8;//(ym-y)/20;
	v[a].b=-x/8;  //-(xm-x)/20;


	}

a=numpoints;
numpoints+=100;  //200 more points
for(b=a;b<numpoints;b++)
{
p[b].a=-20.0;//-2000.0;//(float)(numpoints-a)/b;
p[b].b=0.0;//(float)(numpoints-a)/b;
p[b].d=-5000.0;//0.0;

v[b].a=0.0;//10.0;
v[b].b=0.0;
v[b].d=30.0;//-200.0;
}

//p[numpoints-1].a = -50000.0;
#endif  //GALAXY

#if TWOATOMS
//scale=1.0/4.0;
viewer.a=1000;
viewer.d = -5000.0;

numpoints=300;
gravity=1;//numpoints;//*1000;


np3=(int)(numpoints/2);
np2=(int)pow(np3,0.5);

xi=15.0;
yi=15.0;

xm=-2000;   //-DW2*1000;  
ym=0;    //-DH2*1000;

xs=xm-(np2*xi)/2;
ys=ym-(np2*yi)/2;

x=xs;  
y=ys;

for (a=0;a<np3;a++)
	{
	x+=xi;
	if (x>(xs+np2*xi)) {x=xs; y+=yi;}
  p[a].a = x;
	p[a].b = y;
  p[a].d = 0;

	v[a].a=(ym-y)/20;
	v[a].b=-(xm-x)/20;
	}


np3=(int)(numpoints/2);
np2=(int)pow(np3,0.5);

xi=15.0;
yi=15.0;

xm=2000;   //-DW2*1000;  
ym=0;    //-DH2*1000;

xs=xm-(np2*xi)/2;
ys=ym-(np2*yi)/2;

x=xs;  
y=ys;

for (a=np3;a<numpoints;a++)
	{
	x+=xi;
	if (x>(xs+np2*xi)) {x=xs; y+=yi;}
  p[a].a = x;
	p[a].b = y;
  p[a].d = 0;

	v[a].a=(ym-y)/20;
	v[a].b=-(xm-x)/20;
	}



#endif  //TWOATOMS


#if BMAPS
	bmnum=0;

bmphead.bfOffBits =sizeof(BITMAPFILEHEADER);
bmphead.bfReserved1 =0;
bmphead.bfReserved2 =0;
bmphead.bfSize = DW*DH*3+sizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER);
bmphead.bfType =19778;  // "BM"  19778



	memset(&bmpihead,0,sizeof(bmpihead));
bmpihead.biSize=0x28;
bmpihead.biWidth=DW;
bmpihead.biHeight=DH;
bmpihead.biPlanes=1;
bmpihead.biBitCount=24;
bmpihead.biCompression=0;
bmpihead.biSizeImage=DW*DH*3;//0;
bmpihead.biXPelsPerMeter=0;//0xb12;
bmpihead.biYPelsPerMeter=0;//0xb12;
bmpihead.biClrUsed=0;//0x100;
bmpihead.biClrImportant=0;//0x100;



//write out to bitmap files
strcpy(bfn,"aaaaaaaa");


//fptr2=fopen("bmplist.txt","wb");



#endif



}


void transform_draw_points(void)
{
int ox,oy,oz;
unsigned char * offst;
float tz;
point3df p2,p3;


//clear screen
memset(screen,0,FRAMESIZE);//DH*DW*3);

for(a=0;a<numpoints;a++)
	{
  //ox=p[a].a*scale+DW2;
	//oy=p[a].b*scale+DH2;

	//viewer angle - both above and this could be done with 4x4
	pmmult(p[a],&p2,viewerm);

//viewer location
	p3.a=p2.a - viewer.a;
	p3.b=p2.b - viewer.b;
	p3.d=p2.d - viewer.d;


	//tz=p2.d-viewer.d;
	//perspective and screen offset
	ox=p3.a*2000/p3.d  + DW2;
	oy=p3.b*2000/p3.d  + DH2;

	//ox=p3.a  + DW2;
	//oy=p3.b  + DH2;

	


	if (ox>-1 && ox < DW && oy>-1 && oy<DH)  //clip
		{
		offst=screen+ox*3+oy*DW*3;
		*(offst)=0xff;
		*(offst+1)=0xff;
		*(offst+2)=0xff;
		}
	}




}


void pmmult(point3df p1,point3df *p2,matrix3df m)
{
p2->a=p1.a*m.a[0] + p1.b*m.b[0] + p1.d*m.d[0]; 
p2->b=p1.a*m.a[1] + p1.b*m.b[1] + p1.d*m.d[1];
p2->d=p1.a*m.a[2] + p1.b*m.b[2] + p1.d*m.d[2];
}


void mrota(matrix3df *m, float a)
{
	//rotation matrix
m->a[0]=1.0;  m->a[1]=0.0;  m->a[2]=sin(a);
m->b[0]=0.0;  m->b[1]=cos(a);  m->b[2]=-sin(a);
m->d[0]=0.0;  m->d[1]=sin(a);  m->d[2]=cos(a);
}


void mrotb(matrix3df *m, float a)
{
	//rotation matrix
m->a[0]=cos(a);  m->a[1]=0.0;  m->a[2]=sin(a);
m->b[0]=0.0;  m->b[1]=1.0;  m->b[2]=0.0;
m->d[0]=-sin(a);  m->d[1]=0.0;  m->d[2]=cos(a);

}

void mrotd(matrix3df *m, float a)
{
	//rotation matrix
m->a[0]=cos(a);  m->a[1]=-sin(a);  m->a[2]=0.0;
m->b[0]=sin(a);  m->b[1]=cos(a);  m->b[2]=0.0;
m->d[0]=0.0;  m->d[1]=0.0;  m->d[2]=1.0;
}



void mmult(matrix3df m1, matrix3df m2,matrix3df *m3)
{
m3->a[0]=m1.a[0]*m2.a[0]+m1.a[1]*m2.b[0]+m1.a[2]*m2.d[0];
m3->a[1]=m1.a[0]*m2.a[1]+m1.a[1]*m2.b[1]+m1.a[2]*m2.d[1];
m3->a[2]=m1.a[0]*m2.a[2]+m1.a[1]*m2.b[2]+m1.a[2]*m2.d[2];

m3->b[0]=m1.b[0]*m2.a[0]+m1.b[1]*m2.b[0]+m1.b[2]*m2.d[0];
m3->b[1]=m1.b[0]*m2.a[1]+m1.b[1]*m2.b[1]+m1.b[2]*m2.d[1];
m3->b[2]=m1.b[0]*m2.a[2]+m1.b[1]*m2.b[2]+m1.b[2]*m2.d[2];

m3->d[0]=m1.d[0]*m2.a[0]+m1.a[1]*m2.d[0]+m1.a[2]*m2.d[0];
m3->d[1]=m1.d[0]*m2.a[1]+m1.a[1]*m2.d[1]+m1.a[2]*m2.d[1];
m3->d[2]=m1.d[0]*m2.a[2]+m1.a[1]*m2.d[2]+m1.a[2]*m2.d[2];
}

void angle2matrix(point3df viewera,matrix3df *trot)
{
matrix3df tm,tm2;

mrota(&tm,viewera.a);
*trot=tm;
mrotb(&tm,viewera.b);
mmult(*trot,tm,&tm2);
mrotd(&tm,viewera.d);
mmult(tm2,tm,trot);

}


LRESULT CALLBACK WndProc (HWND hwnd, UINT iMsg, WPARAM wParam, LPARAM lParam) 
{

matrix3df trot,tview;
float difdist0,difdist1;


	switch (iMsg) {

    case WM_KEYDOWN://0x0100:  //WM_KEYDOWN
      switch( wParam )  //0x84
        {
        case VK_ESCAPE:
        case VK_F12:
			  //QUIT=1;
			  PostMessage(hwnd, WM_CLOSE, 0, 0);
			  //free(MemoryStart);
            break;
		    case VK_NUMPAD5:   //Z+
					//scale*=2.0;
					viewer.d/=2.0;
			     break;
		    case VK_NUMPAD0:   //Z-
		    case VK_INSERT:
					//scale/=2.0;
					  viewer.d*=2.0;
			    break;
				case VK_DOWN:
				case VK_NUMPAD2:
					  viewer.b-=(viewer.d/100.0);
					break;
				case VK_UP:
        case VK_NUMPAD8:
					  viewer.b+=(viewer.d/100.0);
					break;
				case VK_LEFT:
        case VK_NUMPAD4:
					  viewer.a+=(viewer.d/100.0);
					break;
				case VK_RIGHT:
        case VK_NUMPAD6:
					  viewer.a-=(viewer.d/100.0);
					break;
				case 99:  //rt dn
					  viewer.a-=(viewer.d/100.0);
					  viewer.b-=(viewer.d/100.0);
					break;
				case 105:  //rt up
					  viewer.a-=(viewer.d/100.0);
					  viewer.b+=(viewer.d/100.0);
					break;
				case 103:  //lt up
					  viewer.a+=(viewer.d/100.0);
					  viewer.b+=(viewer.d/100.0);
					break;
				case 97:  //lt dn
					  viewer.a+=(viewer.d/100.0);
					  viewer.b-=(viewer.d/100.0);
					break;

				case 65:  //a
					  viewera.b+=0.1;
						if (viewera.b>2*PI) 
							{viewera.b=0.0;}
						angle2matrix(viewera,&trot);
						mmult(ident,trot,&viewerm);
						//rotate on y axis
					break;
				case 68:  //d
					  viewera.b-=0.1;
						if (viewera.b<-2*PI) 
							{viewera.b=0.0;}
						angle2matrix(viewera,&trot);
						mmult(ident,trot,&viewerm);
						//rotate on y axis
					break;
        case 84:  //t  track a point
					  track=!track;
						if (track==1)
						{
							difdist0=0;
							
						for(a=0;a<numpoints;a++)
						{
							difdist1=pow((viewer.a - p[a].a),2.0);
							difdist1+=pow((viewer.b - p[a].b),2.0);
							difdist1=pow(difdist1,0.5);
						if (difdist1 < difdist0 || a==0)
						{
						difdist0=difdist1;
						trackpoint=a;
						}
						}  //end for a
						}  //end track==1
					  break;
				case 87:  //w
					  viewera.a+=0.1;
						if (viewera.a>2*PI) 
							{viewera.a=0.0;}
						angle2matrix(viewera,&trot);
						mmult(ident,trot,&viewerm);
						//rotate on x axis
					break;
				case 88:  //x
					  viewera.a-=0.1;
						if (viewera.a<-2*PI) 
							{viewera.a=0.0;}
						angle2matrix(viewera,&trot);
						mmult(ident,trot,&viewerm);
						//rotate on x axis
					break;

				}
       break;
		case WM_DESTROY : 

			//test
	  ReleaseDC(hwMain,hdcWin);
    DeleteDC(hdcMem);
//	  DeleteObject(hBMP);


			PostQuitMessage(0);
			return 0;
		break;
	}

	return DefWindowProc (hwnd, iMsg, wParam, lParam) ;

}