metropolis code
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <SDL3/SDL.h>
#include <SDL3/SDL_main.h>
#include "sdl_utils.c"
void draw_dipole(SDL_Object* sdl, int** s, int size, int x, int y){
// Pick "appropriate" colour for up/down dipole
if(s[x][y]==1) SDL_SetRenderDrawColor(sdl->renderer, 0, 0, 0, SDL_ALPHA_OPAQUE);
else SDL_SetRenderDrawColor(sdl->renderer, 255, 255, 255, SDL_ALPHA_OPAQUE);
// Draw the dipole
float rim = (float)100/size;
float width = (float)(sdl->window_size)/size;
SDL_FRect dipole = {x * width + rim/2, y * width + rim/2, width - rim, width - rim};
SDL_RenderFillRect(sdl->renderer, &dipole);
}
bool init_grid(SDL_Object* sdl, int** s, int size, bool no_graphics){ // Fill grid randomly and draw to texture and render once for each new row
for(int x=0; x<size; x++){
for(int y=0; y<size; y++){
s[x][y] = 2 * (rand()%2) - 1;
draw_dipole(sdl, s, size, x, y);
if(!y && !no_graphics) render(sdl);
if(!event_handling(sdl)) return false;
}
}
return true;
}
int delta_U(int** s, int size, int x, int y){ // Calculate energy difference for flipping dipole at (x,y)
int dU = 0;
for(int n=0; n<2; n++){
dU += 2*s[x][y]*(s[(x+2*n-1+size)%size][y] + s[x][(y+2*n-1+size)%size]); // Implementing the periodic boundary conditions in a mathematically elegant way
}
return dU;
}
void calc_U(int* U, int** s, int size){ // Calculate the inner energy for the entire grid.
for(int x=0; x<size; x++){
for(int y=0; y<size; y++){
*U += -s[x][y] * (s[(x+1+size)%size][y] + s[x][(y+1+size)%size]);
}
}
}
void metropolis(int** s, int size, int iterations, float T, int* U, int* M, SDL_Object* sdl, bool no_graphics){
int x; int y;
for(int i=0; i<iterations; i++){ // main loop
// Render and present changes after every some amount of flips
if(!(i%(100*(size^2))) && !no_graphics) render(sdl);
// Pick random dipole
x = (rand() % size+1)-1;
y = (rand() % size+1)-1;
int dU = delta_U(s, size, x, y); // Calculate energy difference
float random_number = ((double)rand()/(double)RAND_MAX);
// Flip dipole at (x,y) if energy difference is negative,
// zero or if positive gets flipped with a probability matching
// the boltzmann weight ratio and update texture
if((dU<=0) || random_number<(exp(-dU/T))){
U[i+1] = U[i] + dU; // Update inner energy
s[x][y] *= -1;
M[i+1] = M[i] + 2 * s[x][y]; // Update magnetization
draw_dipole(sdl, s, size, x, y);
}
else{
U[i+1]=U[i];
M[i+1]=M[i];
}
if(!event_handling(sdl)) break; // Check if user wants to quit
}
}
void temperature_axis_mode(int** s, int size, int iterations, float step, float T, int* U, int* U_average, int* M, SDL_Object* sdl){
printf("Alt Mode live ticker:\n");
int num_values = (int)(3 * (1/step) + 1);
float* T_axis = malloc(num_values * sizeof(float));
int* U_avg_axis = malloc(num_values * sizeof(int));
for(T=4; T>1; T-=step) {
printf("Debug: %d\n", (int)(1/step*T-1/step));
printf("We have a temperature of %1.2f.\n", T);
for(int i=0; i<iterations+1; i++){
U[i] = 0;
M[i] = 0;
}
metropolis(s, size, iterations, T, U, M, sdl, true);
*U_average = 0;
for(int i=0; i<iterations+1; i++) *U_average += U[i];
*U_average = *U_average/(iterations+1);
T_axis[(int)(1/step*T-1/step)] = T;
U_avg_axis[(int)(1/step*T-1/step)] = *U_average;
printf("Our Ising run decided on an average energy of %d.\n", *U_average);
}
FILE* fp = fopen("./alt_mode_data.txt", "w");
fprintf(fp, "T, U\n");
for(int i=0; i<num_values; i++){ fprintf(fp,"%1.2f, %d\n", T_axis[i], U_avg_axis[i]);
fclose(fp);
}
free(T_axis);
free(U_avg_axis);
}
int main(int argc, char** argv){
// program init variables(size of grid and temperature in units of epsilon/k)
// TODO: Let SDL handle fetching these from user
float T = atof(argv[1]);
int size = atoi(argv[2]);
int scale = atoi(argv[3]);
bool no_graphics = false;
if(argv[4]) no_graphics = true;
srand(time(NULL)); // Set up rand
SDL_Object sdl; // Declare object containing window-, renderer-reference and event
if(!no_graphics) SDL_Object_Init(&sdl, 900, "Metropolis Simulation"); // Initialize SDL
// Init grid
int** s = malloc(size*sizeof(int*));
for(int i=0; i<size; i++){
s[i] = malloc(size*sizeof(int));
}
// Stop if user hits "esc" in the process of initializing the grid
if(!init_grid(&sdl, s, size, no_graphics)) goto jumppoint;
if(!no_graphics) render(&sdl);
int iterations = pow(10, scale)*pow(size, 2); // Set number of iterations
// Declare and init inner energy
int* U = malloc((iterations+1) * sizeof(int));
for(int i=0; i<iterations+1; i++) U[i] = 0;
calc_U(&U[0], s, size); // Calc initial inner energy
// Declare and init magnetization
int* M = malloc((iterations+1) * sizeof(int));
for(int i=0; i<iterations; i++) M[i] = 0;
for(int x=0; x<size; x++){
for(int y=0; y<size; y++){ // Calc initial magnetization
M[0] += s[x][y];
}
}
int U_average = 0;
int M_average = 0;
printf("T: %2.1f, size: %d\n", T, size);
float step;
if(no_graphics && argv[5]){
step = atof(argv[5]);
temperature_axis_mode(s, size, iterations, step, T, U, &U_average, M, &sdl);
goto jumppoint;
}
// Call subroutine containing the Metropolis Algorithm
metropolis(s, size, iterations, T, U, M, &sdl, no_graphics);
for(int i=0; i<iterations+1; i++) U_average += U[i];
U_average = U_average/(iterations+1);
for(int i=0; i<iterations+1; i++) M_average += M[i];
M_average = M_average/(iterations+1);
if(!no_graphics) render(&sdl); // Render the final state
printf("Initial energy: %d\n", U[0]); // Print initial energy
printf("Average energy: %d\n", U_average); // Print average energy
printf("Final energy: %d\n", U[iterations]); // Print final energy
// Print results for the magnetization
printf("Initial magnetization: %d\n", M[0]);
printf("Average magnetization: %d\n", M_average);
printf("Final magnetization: %d\n", M[iterations]);
// Keep window open at the end until user wants to quit
if(!no_graphics) while(event_handling(&sdl));
jumppoint:
// Free memory
SDL_DestroyWindow(sdl.window);
SDL_DestroyRenderer(sdl.renderer);
free(U);
free(M);
for(int i=0; i<size; i++){
free(s[i]);
}
free(s);
// End
SDL_Quit();
return 0;
}
some sdl util code
#include <SDL3/SDL.h>
#include <SDL3/SDL_main.h>
typedef struct{
SDL_Window* window;
SDL_Renderer* renderer;
SDL_Texture* texture;
SDL_Event event;
int window_size;
char* window_title;
} SDL_Object;
void SDL_Object_Init(SDL_Object* sdl, size_t window_size, char* window_title){
SDL_Init(SDL_INIT_VIDEO);
sdl->window_size = window_size;
sdl->window_title = window_title;
sdl->window = SDL_CreateWindow(sdl->window_title, sdl->window_size, sdl->window_size, 0);
sdl->renderer = SDL_CreateRenderer(sdl->window, NULL);
sdl->texture = SDL_CreateTexture(sdl->renderer, SDL_PIXELFORMAT_RGBA8888, SDL_TEXTUREACCESS_TARGET, sdl->window_size, sdl->window_size);
SDL_SetRenderTarget(sdl->renderer, sdl->texture);
}
bool event_handling(SDL_Object* sdl){
while(SDL_PollEvent(&sdl->event)){
switch(sdl->event.type){
case SDL_EVENT_WINDOW_CLOSE_REQUESTED:{
if(sdl->window){
SDL_DestroyWindow(sdl->window);
sdl->window=NULL;
}
}
break;
case SDL_EVENT_KEY_DOWN:{
switch(sdl->event.key.key){
case SDLK_ESCAPE:
return false;
break;
}
}
break;
case SDL_EVENT_QUIT:
break;
}
}
return true;
}
void render(SDL_Object* sdl){ // Set renderer onto window, set bg-colour, draw texture on background, present and set renderer back onto texture
SDL_SetRenderTarget(sdl->renderer, NULL);
SDL_SetRenderDrawColor(sdl->renderer, 127, 127, 127, SDL_ALPHA_OPAQUE);
SDL_RenderClear(sdl->renderer);
SDL_FRect dstrect = {0, 0, sdl->window_size, sdl->window_size};
SDL_RenderTexture(sdl->renderer, sdl->texture, NULL, &dstrect);
SDL_RenderPresent(sdl->renderer);
SDL_Delay(50);
SDL_SetRenderTarget(sdl->renderer, sdl->texture);
}
Thank you for having looked at my code.