//Gravitational Wave Source Simulations

//Jeff Jauregui, April 2003

//code for implementing graphics borrowed from

//Andre LaMothe, "Windows game programming for dummies"

//***NOTE: This code may only compile correctly by creating

//a Win32 app in MS Visual C++

//header files needed to run windows program

#define WIN32_LEAN_AND_MEAN

#include <windows.h>

#include <windowsx.h>

#include <stdio.h>

#include <math.h>

#include <fstream.h>

//define a color and paint brush for drawing on the screen

COLORREF black = RGB(0,0,0);

HBRUSH bbrush = (HBRUSH__ *)GetStockObject(BLACK_BRUSH);

//define the number of radial and angular steps in the 2D polar grid

#define R_MAX 300

#define THETA_MAX 200

//define the number of time steps

#define T_MAX 100

#define pi 3.1416

//prototypes for various functions

double density_R_0(double, double, double); //returns a density as a function of r, theta, and time

double get_angle(double); //converts an angle index into an actual angle, in radians

double get_r(double); //converts an radius index into an actual value for the radius

double get_t(double); //converts a time index into a normalized time value

double calc_Q00(); //returns the first component of the mass quadrupole moment

void d2Q_dt2(); //calculates grav. wave amplitude using the second time derivative of Q00

//global variables

double grid[R_MAX][THETA_MAX]; //defines density(r,theta) polar grid

double Q00[T_MAX]; //an array to hold Q00 as a function of time

// defines for windows

#define WINDOW_CLASS_NAME "WINCLASS1"

// GLOBALS ////////////////////////////////////////////////

HWND main_window_handle = NULL; // save the window handle

//GRAPHICS FUNCTIONS:

LRESULT CALLBACK WindowProc(HWND hwnd,

UINT msg,

WPARAM wparam,

LPARAM lparam)

{

// this is the main message handler of the system

PAINTSTRUCT ps; // used in WM_PAINT

HDC hdc; // handle to a device context

// what is the message

int y = (int)HIWORD(lparam);

switch(msg)

{

case WM_CREATE:

{

// do initialization stuff here

return(0);

} break;

case WM_PAINT:

{

// simply validate the window

hdc = BeginPaint(hwnd,&ps);

EndPaint(hwnd,&ps);

return(0);

} break;

case WM_KEYDOWN:

PostMessage(hwnd,WM_DESTROY,NULL,NULL);

return 0;

break;

case WM_MOUSEMOVE:

break;

case WM_DESTROY:

{

// kill the application

DeleteObject(bbrush);

PostQuitMessage(0);

return(0);

} break;

default:break;

} // end switch

// process any messages that we didn't take care of

return (DefWindowProc(hwnd, msg, wparam, lparam));

} // end WinProc

// WINMAIN ////////////////////////////////////////////////

int WINAPI WinMain( HINSTANCE hinstance,

HINSTANCE hprevinstance,

LPSTR lpcmdline,

int ncmdshow)

{

WNDCLASS winclass; // this will hold the class we create

HWND hwnd; // generic window handle

MSG msg; // generic message

// first fill in the window class stucture

winclass.style = CS_DBLCLKS | CS_OWNDC |

CS_HREDRAW | CS_VREDRAW;

winclass.lpfnWndProc = WindowProc;

winclass.cbClsExtra = 0;

winclass.cbWndExtra = 0;

winclass.hInstance = hinstance;

winclass.hIcon = LoadIcon(NULL, IDI_APPLICATION);

winclass.hCursor = LoadCursor(NULL, IDC_ARROW);

winclass.hbrBackground = (HBRUSH__ *)GetStockObject(BLACK_BRUSH);

winclass.lpszMenuName = NULL;

winclass.lpszClassName = WINDOW_CLASS_NAME;

// register the window class

if (!RegisterClass(&winclass))

return(0);

// create the window

if (!(hwnd = CreateWindow(WINDOW_CLASS_NAME, // class

"Gravity", // title

WS_OVERLAPPEDWINDOW | WS_VISIBLE,

0,0, // x,y

320,200, // width, height

NULL, // handle to parent

NULL, // handle to menu

hinstance,// instance

NULL))) // creation parms

return(0);

ShowWindow(hwnd, SW_SHOWMAXIMIZED);

// get the window handle for graphics

main_window_handle = hwnd;

HDC hdc = GetDC(hwnd);

//***********Main body of program here***********

//create a rectangle to help erase the screen later

RECT rect;

rect.top = 0;

rect.bottom= 800;

rect.left = 0;

rect.right = 1200;

//loop through each time step, indexed by n

for(int n = 0; n < T_MAX; n++)

{

int r,o; //looping variables, r = radius, o = theta (angle)

//the min and max will be determined in order to perform color-mapping

double max = 0;

double min = 1000000;

//loop through grid for every point

for(r = 0; r < R_MAX; r++)

{

for(o = 0; o < THETA_MAX; o++)

{

//fill each discrete point on the grid with a density value

grid[r][o] = density_R_0(get_r(r),get_angle(o),(double)n);

if(grid[r][o] > max)

max = grid[r][o];

if(grid[r][o] < min)

min = grid[r][o];

}

}

//plot the density function in the window

bool draw_graphics = true; //turns on/off the displaying of graphics

if(draw_graphics)

{

//fill the screen with black in order to erase

FillRect(hdc,&rect,bbrush);

//loop through each point in the polar grid

for(r = 0; r < R_MAX; r++)

{

for(o = 0; o <= THETA_MAX; o++)

{

//draw a pixel whose intensity is determined by the density at that point

SetPixel(hdc, 1152/2+ 300*get_r(r)*cos(2*3.14*o/THETA_MAX),

864/2 + 300*get_r(r)*sin(2*3.14*o/THETA_MAX),

RGB(255/(max-min)*(grid[r][o]-min),0,0));

}

}

}

Q00[n] = calc_Q00();

}//end the main simulation loop

//calculate and output the gravitational wave amplitudes to a file wavedata.m

d2Q_dt2();

//***********End Main body of program***********

// enter windows event loop

while(1)

{

if (PeekMessage(&msg,NULL,0,0,PM_REMOVE))

{

// test if this is a quit

if (msg.message == WM_QUIT)

break;

// translate any accelerator keys

TranslateMessage(&msg);

// send the message to the window proc

DispatchMessage(&msg);

} // end if

} // end while

// return to Windows like this

return(msg.wParam);

} // end WinMain

//returns the density as a function of radius, angle o, and time

double density_R_0(double r, double o, double t)

{

double w = .12*get_t(t)*10; //angular frequency of rotation

//various possible density functions:

// return exp(-r*0.2);

// return exp(-r*(2+sin(2*o+w*t)));

//return exp(-r*0.2*t)*exp(-r*(1.1+sin(2*o+w*t)));

//let the mass be more concentrated toward the center over time,

//and let the distribution be asymmetric

return exp(-r*0.5*(get_t(t)+1)*(2+sin(2*o+w*get_t(t))));

}

//take in an angle "index", and return an angle on a scale of 0 to pi

double get_angle(double o)

{

return 2*pi*o/THETA_MAX;

}

//take in an radius "index", and return a radius on a scale of 0 to 1

double get_r(double r)

{

return r/R_MAX;

}

//take in an time "index",and return a time on a scale of 0 to 10

double get_t(double t)

{

return 10*t/T_MAX;

}

//calculate the Q00 component of the mass quadrupole moment

double calc_Q00()

{

//define the summing variable for numerical integration

double sum = 0;

//looping variables

int r, o;

//loop through the polargrid

for(r = 0; r < R_MAX; r++)

{

for(o = 0; o < THETA_MAX; o++)

{

//add to the sum the integrand (2x^2-y^2) times the density, times the area element

sum += pow(get_r(r),3)*(2*pow(cos(get_angle(o)),2)-pow(sin(get_angle(o)),2))*grid[r][o]*2*pi/THETA_MAX/R_MAX;//*10*2*pi/THETA_MAX;

}

}

//return the total sum

return sum;

}

//calculate the second time derivative the mass quadrupole moment, and output to a file

void d2Q_dt2()

{

//create and open the file for output

ofstream outfile;

outfile.open("wavedata.m");

//loop through the Q00 data and numerically evaluate the 2nd derivative

for(int i = 0; i < T_MAX-3; i++)

outfile << (Q00[i+2]-2*Q00[i+1]+Q00[i])*(double)(T_MAX*T_MAX) << '\n';

//close the file

outfile.close();

}