Zadanie 1 polepszenie pracy algorytmu
This commit is contained in:
parent
0c5a465906
commit
23be099698
@ -1,225 +1,167 @@
|
|||||||
#include <stdio.h>
|
|
||||||
#include <iostream>
|
|
||||||
#include "shader.h"
|
|
||||||
#include <cmath>
|
|
||||||
#include <GL/glew.h>
|
#include <GL/glew.h>
|
||||||
#include <GLFW/glfw3.h>
|
#include <GLFW/glfw3.h>
|
||||||
#include <glm/glm.hpp>
|
#include <cmath>
|
||||||
#include <glm/gtc/matrix_transform.hpp>
|
#include <iostream>
|
||||||
#include <glm/gtc/type_ptr.hpp>
|
#include <vector>
|
||||||
GLFWwindow* window;
|
|
||||||
|
|
||||||
#define PI 3.141592 // PI approximate value
|
// Konstante
|
||||||
#define HEIGHT 0.0025 // step length for numerical integration
|
const float g = 9.81f; // przyspieszenie ziemskie
|
||||||
#define ROD_LENGHT 0.5 // length of rod
|
const float r = 1.0f; // długość nici
|
||||||
#define GRAVITY 9.81 // gravitational constant
|
const float dt = 0.01f; // krok czasowy
|
||||||
#define theta_0 PI / 2 // Initial angle
|
|
||||||
#define omega_0 0 // Initial angular velocity
|
|
||||||
#define time_0 0 // Initial time
|
|
||||||
#define RADIUS 0.15 // Radius of pendulum circle
|
|
||||||
|
|
||||||
#define A 1.4 // Amplitude of the driving force
|
// Zmienne globalne
|
||||||
#define k 0.67 // Frequency of the driving force
|
float theta = 0.5f; // początkowy kąt
|
||||||
|
float omega = 0.0f; // początkowa prędkość kątowa
|
||||||
|
|
||||||
// Function for angular velocity (f function for numerical integration)
|
GLuint VBO, VAO, EBO;
|
||||||
float f(float time, float theta, float omega) {
|
|
||||||
return omega;
|
void initOpenGL()
|
||||||
|
{
|
||||||
|
// Setup Vertex Array Object and Vertex Buffer Object
|
||||||
|
glGenVertexArrays(1, &VAO);
|
||||||
|
glGenBuffers(1, &VBO);
|
||||||
|
|
||||||
|
glBindVertexArray(VAO);
|
||||||
|
|
||||||
|
glBindBuffer(GL_ARRAY_BUFFER, VBO);
|
||||||
|
glBufferData(GL_ARRAY_BUFFER, sizeof(float) * 4, nullptr, GL_DYNAMIC_DRAW);
|
||||||
|
|
||||||
|
// Pozycja
|
||||||
|
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 2 * sizeof(float), (void*)0);
|
||||||
|
glEnableVertexAttribArray(0);
|
||||||
|
|
||||||
|
glBindBuffer(GL_ARRAY_BUFFER, 0);
|
||||||
|
glBindVertexArray(0);
|
||||||
}
|
}
|
||||||
|
|
||||||
// Function for angular acceleration (g function for numerical integration)
|
// Metoda Eulera
|
||||||
float g(float time, float theta, float omega) {
|
void updateEuler()
|
||||||
return -(GRAVITY / ROD_LENGHT) * sin(theta) + (A * cos(k * time)) / ROD_LENGHT;
|
{
|
||||||
|
float alpha = -(g / r) * sin(theta); // przyspieszenie kątowe
|
||||||
|
omega += alpha * dt; // aktualizacja prędkości kątowej
|
||||||
|
theta += omega * dt; // aktualizacja kąta
|
||||||
}
|
}
|
||||||
|
|
||||||
// Function to draw the circle (pendulum ball)
|
// Metoda Verleta
|
||||||
void drawCircle(float array[]) {
|
void updateVerlet()
|
||||||
int corner_one, corner_two, corner_three; // corners of triangles. GL_TRIANGLES starts to draw counterclockwise.
|
{
|
||||||
corner_one = -6;
|
static float prev_theta = theta;
|
||||||
corner_two = -4;
|
float alpha = -(g / r) * sin(theta); // przyspieszenie kątowe
|
||||||
corner_three = -2;
|
float new_theta = 2 * theta - prev_theta + alpha * dt * dt; // aktualizacja kąta
|
||||||
|
prev_theta = theta;
|
||||||
for (int angle = 1; angle <= 360; angle++) {
|
theta = new_theta;
|
||||||
corner_one = corner_one + 6;
|
|
||||||
corner_two = corner_two + 6;
|
|
||||||
corner_three = corner_three + 6;
|
|
||||||
|
|
||||||
array[corner_one] = 0.0f;
|
|
||||||
array[corner_one + 1] = 0.0f;
|
|
||||||
|
|
||||||
array[corner_two] = RADIUS * cos((angle - 1) * PI / 180);
|
|
||||||
array[corner_two + 1] = RADIUS * sin((angle - 1) * PI / 180);
|
|
||||||
|
|
||||||
array[corner_three] = RADIUS * cos(angle * PI / 180);
|
|
||||||
array[corner_three + 1] = RADIUS * sin(angle * PI / 180);
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
void RungeKuttaIntegration(float& theta, float& omega, float& time) {
|
// Metoda Rungego-Kutty rzędu 4
|
||||||
float h = HEIGHT; // Step size
|
void updateRungeKutta()
|
||||||
float k1_theta = h * f(time, theta, omega);
|
{
|
||||||
float k1_omega = h * g(time, theta, omega);
|
auto f = [](float theta, float omega) { return -(g / r) * sin(theta); };
|
||||||
float k2_theta = h * f(time + h / 2, theta + k1_theta / 2, omega + k1_omega / 2);
|
|
||||||
float k2_omega = h * g(time + h / 2, theta + k1_theta / 2, omega + k1_omega / 2);
|
|
||||||
float k3_theta = h * f(time + h / 2, theta + k2_theta / 2, omega + k2_omega / 2);
|
|
||||||
float k3_omega = h * g(time + h / 2, theta + k2_theta / 2, omega + k2_omega / 2);
|
|
||||||
float k4_theta = h * f(time + h, theta + k3_theta, omega + k3_omega);
|
|
||||||
float k4_omega = h * g(time + h, theta + k3_theta, omega + k3_omega);
|
|
||||||
|
|
||||||
// Update theta and omega
|
float k1_theta = omega;
|
||||||
theta = theta + h/6 * (k1_theta + 2*k2_theta + 2*k3_theta + k4_theta);
|
float k1_omega = f(theta, omega);
|
||||||
omega = omega + h/6 * (k1_omega + 2*k2_omega + 2*k3_omega + k4_omega);
|
|
||||||
|
|
||||||
// Keep theta in the range of -2PI to 2PI
|
float k2_theta = omega + 0.5f * dt * k1_omega;
|
||||||
if (theta > 2 * PI) theta -= 2 * PI;
|
float k2_omega = f(theta + 0.5f * dt * k1_theta, omega + 0.5f * dt * k1_omega);
|
||||||
if (theta < -2 * PI) theta += 2 * PI;
|
|
||||||
|
|
||||||
time += h; // Increment time
|
float k3_theta = omega + 0.5f * dt * k2_omega;
|
||||||
|
float k3_omega = f(theta + 0.5f * dt * k2_theta, omega + 0.5f * dt * k2_omega);
|
||||||
|
|
||||||
|
float k4_theta = omega + dt * k3_omega;
|
||||||
|
float k4_omega = f(theta + dt * k3_theta, omega + dt * k3_omega);
|
||||||
|
|
||||||
|
theta += (dt / 6.0f) * (k1_theta + 2.0f * k2_theta + 2.0f * k3_theta + k4_theta);
|
||||||
|
omega += (dt / 6.0f) * (k1_omega + 2.0f * k2_omega + 2.0f * k3_omega + k4_omega);
|
||||||
}
|
}
|
||||||
|
|
||||||
void EulerIntegration(float& theta, float& omega, float& time) {
|
void drawPendulum()
|
||||||
float h = HEIGHT; // Step size
|
{
|
||||||
float theta_new = theta + h * f(time, theta, omega);
|
float x = r * sin(theta);
|
||||||
float omega_new = omega + h * g(time, theta, omega);
|
float y = -r * cos(theta);
|
||||||
|
|
||||||
// Update theta and omega
|
float vertices[] = {
|
||||||
theta = theta_new;
|
0.0f, 0.0f, // Punkt zaczepienia
|
||||||
omega = omega_new;
|
x, y // Punkt masy
|
||||||
|
};
|
||||||
|
|
||||||
// Keep theta in the range of -2PI to 2PI
|
glBindBuffer(GL_ARRAY_BUFFER, VBO);
|
||||||
if (theta > 2 * PI) theta -= 2 * PI;
|
glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(vertices), vertices);
|
||||||
if (theta < -2 * PI) theta += 2 * PI;
|
glBindBuffer(GL_ARRAY_BUFFER, 0);
|
||||||
|
|
||||||
time += h; // Increment time
|
glClear(GL_COLOR_BUFFER_BIT);
|
||||||
}
|
glBindVertexArray(VAO);
|
||||||
|
|
||||||
void VerletIntegration(float& theta, float& omega, float& time) {
|
glDrawArrays(GL_LINES, 0, 2);
|
||||||
float h = HEIGHT; // Step size
|
glDrawArrays(GL_POINTS, 1, 1);
|
||||||
float theta_new = theta + h * omega + 0.5 * h * h * g(time, theta, omega);
|
|
||||||
float omega_new = (theta_new - theta) / h;
|
|
||||||
|
|
||||||
// Update theta and omega
|
|
||||||
theta = theta_new;
|
|
||||||
omega = omega_new;
|
|
||||||
|
|
||||||
// Keep theta in the range of -2PI to 2PI
|
|
||||||
if (theta > 2 * PI) theta -= 2 * PI;
|
|
||||||
if (theta < -2 * PI) theta += 2 * PI;
|
|
||||||
|
|
||||||
time += h; // Increment time
|
|
||||||
|
|
||||||
|
glBindVertexArray(0);
|
||||||
}
|
}
|
||||||
|
|
||||||
int main()
|
int main()
|
||||||
{
|
{
|
||||||
if (!glfwInit())
|
if (!glfwInit())
|
||||||
{
|
{
|
||||||
fprintf( stderr, "Failed to initialize GLFW\n" );
|
std::cerr << "Nie można zainicjalizować GLFW" << std::endl;
|
||||||
getchar();
|
|
||||||
return -1;
|
return -1;
|
||||||
}
|
}
|
||||||
|
|
||||||
glfwWindowHint(GLFW_SAMPLES, 4);
|
GLFWwindow* window = glfwCreateWindow(800, 600, "Wahadło Matematyczne", NULL, NULL);
|
||||||
glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
|
if (!window)
|
||||||
glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
|
{
|
||||||
glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
|
std::cerr << "Nie można utworzyć okna GLFW" << std::endl;
|
||||||
|
|
||||||
window = glfwCreateWindow( 630, 600, "ZADANIE 1", NULL, NULL);
|
|
||||||
if( window == NULL ){
|
|
||||||
fprintf( stderr, "Failed to open GLFW window. If you have an Intel GPU, they are not 3.3 compatible. Try the 2.1 version of the tutorials.\n" );
|
|
||||||
getchar();
|
|
||||||
glfwTerminate();
|
glfwTerminate();
|
||||||
return -1;
|
return -1;
|
||||||
}
|
}
|
||||||
|
|
||||||
glfwMakeContextCurrent(window);
|
glfwMakeContextCurrent(window);
|
||||||
glewExperimental = true;
|
glewExperimental = GL_TRUE;
|
||||||
if (glewInit() != GLEW_OK) {
|
if (glewInit() != GLEW_OK)
|
||||||
fprintf(stderr, "Failed to initialize GLEW\n");
|
{
|
||||||
getchar();
|
std::cerr << "Nie można zainicjalizować GLEW" << std::endl;
|
||||||
glfwTerminate();
|
|
||||||
return -1;
|
return -1;
|
||||||
}
|
}
|
||||||
|
|
||||||
glfwSetInputMode(window, GLFW_STICKY_KEYS, GL_TRUE);
|
glMatrixMode(GL_PROJECTION);
|
||||||
glClearColor(1.0f, 0.8f, 0.0f, 0.0f);
|
glLoadIdentity();
|
||||||
|
glOrtho(-2, 2, -2, 2, -1, 1);
|
||||||
|
glMatrixMode(GL_MODELVIEW);
|
||||||
|
|
||||||
Shader myshader("pendulum_vs.glsl" , "pendulum_fs.glsl");
|
glPointSize(10.0f);
|
||||||
unsigned int shaderProgram = myshader.programID();
|
|
||||||
float vertices[2160];
|
|
||||||
float vertices2[] = { //vertices2 gives us the rod of the pendulum.
|
|
||||||
-0.01f, 0.0f,
|
|
||||||
0.01f, 0.0f,
|
|
||||||
0.01f, 0.8f,
|
|
||||||
-0.01f, 0.8f,
|
|
||||||
-0.01, 0.0f
|
|
||||||
};
|
|
||||||
|
|
||||||
drawCircle(vertices); //draws the pendulum ball.
|
initOpenGL();
|
||||||
|
|
||||||
unsigned int VBO, VAO, VAO2;
|
int method = 1; // Domyślnie metoda Eulera
|
||||||
glGenVertexArrays(1, &VAO);
|
|
||||||
glGenBuffers(1, &VBO);
|
|
||||||
glBindVertexArray(VAO);
|
|
||||||
glBindBuffer(GL_ARRAY_BUFFER, VBO);
|
|
||||||
glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
|
|
||||||
|
|
||||||
//position attribute for 'vertices'(ball of the pendulum)
|
std::cout << "Wybierz metodę: (1) Euler, (2) Verlet, (3) Runge-Kutta: ";
|
||||||
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 2 * sizeof(float), (void*)0);
|
std::cin >> method;
|
||||||
glEnableVertexAttribArray(0);
|
|
||||||
|
|
||||||
glGenVertexArrays(1, &VAO2);
|
while (!glfwWindowShouldClose(window))
|
||||||
glGenBuffers(1, &VBO);
|
{
|
||||||
glBindVertexArray(VAO2);
|
switch (method)
|
||||||
glBindBuffer(GL_ARRAY_BUFFER, VBO);
|
{
|
||||||
glBufferData(GL_ARRAY_BUFFER, sizeof(vertices2), vertices2, GL_STATIC_DRAW);
|
case 1:
|
||||||
|
updateEuler();
|
||||||
|
break;
|
||||||
|
case 2:
|
||||||
|
updateVerlet();
|
||||||
|
break;
|
||||||
|
case 3:
|
||||||
|
updateRungeKutta();
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
updateEuler();
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
//position attribute for 'vertices2'(rod of the pendulum)
|
drawPendulum();
|
||||||
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 2 * sizeof(float), (void*)0);
|
|
||||||
glEnableVertexAttribArray(0);
|
|
||||||
|
|
||||||
float theta= theta_0;
|
|
||||||
float omega= omega_0;
|
|
||||||
float time = time_0;
|
|
||||||
|
|
||||||
float k1,k2,k3,k4,l1,l2,l3,l4;
|
|
||||||
float driving_force;
|
|
||||||
|
|
||||||
//Initialize f and g functions.
|
|
||||||
f(time,theta,omega);
|
|
||||||
g(time,theta,omega);
|
|
||||||
float current_angle;
|
|
||||||
do{
|
|
||||||
driving_force = A * cos(k * time);
|
|
||||||
current_angle = theta * 180 / PI; //converts theta(radian) to degree
|
|
||||||
|
|
||||||
glClear( GL_COLOR_BUFFER_BIT );
|
|
||||||
glUseProgram(shaderProgram);
|
|
||||||
|
|
||||||
glm::mat4 model = glm::mat4(1.0f);
|
|
||||||
glm::mat4 projection = glm::mat4(1.0f);
|
|
||||||
glm::mat4 view = glm::mat4(1.0f);
|
|
||||||
view = glm::rotate(view, glm::radians(current_angle), glm::vec3(0.0f, 0.0f, 1.0f));
|
|
||||||
view = glm::translate(view, glm::vec3(0.0f, -0.8f, 0.0f));
|
|
||||||
|
|
||||||
glUniformMatrix4fv(glGetUniformLocation(shaderProgram, "projection"), 1, GL_FALSE, glm::value_ptr(projection));
|
|
||||||
glUniformMatrix4fv(glGetUniformLocation(shaderProgram, "view"), 1, GL_FALSE, glm::value_ptr(view));
|
|
||||||
glUniformMatrix4fv(glGetUniformLocation(shaderProgram, "model"), 1, GL_FALSE, glm::value_ptr(model));
|
|
||||||
|
|
||||||
// RungeKuttaIntegration(theta, omega, time);
|
|
||||||
// VerletIntegration(theta, omega, time);
|
|
||||||
EulerIntegration(theta, omega, time);
|
|
||||||
|
|
||||||
glBindVertexArray(VAO);
|
|
||||||
glDrawArrays(GL_TRIANGLES, 0, 1080);
|
|
||||||
glBindVertexArray(VAO2);
|
|
||||||
glDrawArrays(GL_TRIANGLE_STRIP, 0, 5);
|
|
||||||
|
|
||||||
glfwSwapBuffers(window);
|
glfwSwapBuffers(window);
|
||||||
glfwPollEvents();
|
glfwPollEvents();
|
||||||
}
|
}
|
||||||
while( glfwGetKey(window, GLFW_KEY_ESCAPE ) != GLFW_PRESS &&
|
|
||||||
glfwWindowShouldClose(window) == 0);
|
|
||||||
|
|
||||||
|
glDeleteVertexArrays(1, &VAO);
|
||||||
|
glDeleteBuffers(1, &VBO);
|
||||||
|
|
||||||
|
glfwDestroyWindow(window);
|
||||||
glfwTerminate();
|
glfwTerminate();
|
||||||
|
|
||||||
return 0;
|
return 0;
|
||||||
|
Loading…
Reference in New Issue
Block a user