diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000..12b45a4 --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,56 @@ +name: Build and Test + +on: + push: + branches: [ '**' ] + pull_request: + branches: [ '**' ] + +jobs: + build-linux: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Install Dependencies + run: | + sudo apt-get update + sudo apt-get install -y libsdl2-dev build-essential + + - name: Compile + run: make clean && make + + - name: Check binary + run: | + if [ ! -f bin/singularity ]; then + echo "Binary not found!" + exit 1 + fi + echo "Linux build successful." + + build-windows: + runs-on: windows-latest + steps: + - uses: actions/checkout@v4 + + - name: Set up MSYS2 + uses: msys2/setup-msys2@v2 + with: + msystem: MINGW64 + install: >- + mingw-w64-x86_64-gcc + mingw-w64-x86_64-SDL2 + make + + - name: Compile + shell: msys2 {0} + run: make clean && make + + - name: Check binary + shell: msys2 {0} + run: | + if [ ! -f bin/singularity.exe ]; then + echo "Binary not found!" + exit 1 + fi + echo "Windows build successful." diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml new file mode 100644 index 0000000..78a1922 --- /dev/null +++ b/.github/workflows/release.yml @@ -0,0 +1,91 @@ +name: Release + +on: + push: + tags: + - 'v*' + +jobs: + build-linux: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Install Dependencies + run: | + sudo apt-get update + sudo apt-get install -y libsdl2-dev build-essential + + - name: Compile + run: make clean && make + + - name: Package + run: | + mkdir -p singularity-linux + cp bin/singularity singularity-linux/ + tar -czvf singularity-linux.tar.gz singularity-linux/ + + - name: Upload Artifact + uses: actions/upload-artifact@v4 + with: + name: linux-artifact + path: singularity-linux.tar.gz + + build-windows: + runs-on: windows-latest + steps: + - uses: actions/checkout@v4 + + - name: Set up MSYS2 + uses: msys2/setup-msys2@v2 + with: + msystem: MINGW64 + install: >- + mingw-w64-x86_64-gcc + mingw-w64-x86_64-SDL2 + make + + - name: Compile & Bundle + shell: msys2 {0} + run: make clean && make bundle-sdl2 + + - name: Package + shell: bash + run: | + mkdir -p singularity-windows + cp bin/singularity.exe singularity-windows/ + cp bin/SDL2.dll singularity-windows/ + 7z a singularity-windows.zip singularity-windows/ + + - name: Upload Artifact + uses: actions/upload-artifact@v4 + with: + name: windows-artifact + path: singularity-windows.zip + + create-release: + needs: [build-linux, build-windows] + runs-on: ubuntu-latest + permissions: + contents: write + steps: + - name: Download Linux Artifact + uses: actions/download-artifact@v4 + with: + name: linux-artifact + + - name: Download Windows Artifact + uses: actions/download-artifact@v4 + with: + name: windows-artifact + + - name: Create GitHub Release + uses: softprops/action-gh-release@v2 + with: + files: | + singularity-linux.tar.gz + singularity-windows.zip + draft: false + prerelease: false + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/.gitignore b/.gitignore index dbe9c82..76dd37e 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,7 @@ -.vscode/ \ No newline at end of file +.vscode/ +bin/ +build/ +output.txt +*.txt +*.bmp +*.log \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..859ed22 --- /dev/null +++ b/README.md @@ -0,0 +1,77 @@ +# Singularity + +Singularity is a high-performance **S-plane visualizer** and experimental music instrument that translates complex transfer functions, formed by the interaction of **poles and zeros**, into a symphony of color and sound. + +## Purpose & Inspiration + +This project is a tribute to the logical elegance of **Laplace Transformations** and control theory. It provides a real-time, interactive window into the **S-Domain**, where the stability and behavior of linear systems are dictated by the landscape of singularities. + +Influenced by the work and spirit of the late developer **Terry Davis**, Singularity was built with a dedication to first-principles development and the creation of unique, personal digital artifacts. + +The goal is to move beyond the textbook definitions of transfer functions. It allows users to **see the colors of math** and explore the infinite intricacy of the complex plane. By mapping mathematical magnitude and phase to vibrant visual gradients and synthetic soundscapes, we turn the screen into an auditory and visual representation of the **Laplace landscape**. + +## Features + +- **High-Performance Math Engine**: Built in C with a Structure of Arrays (SoA) layout for maximum SIMD efficiency. +- **Dynamic Sonification**: The field magnitude and phase are mapped to a procedurally generated synthesizer, turning the visualizer into a playable instrument. +- **Advanced Visual Overlays**: + - **Phase Trails**: Real-time vector field streamlines tracing the flow of the complex plane. + - **Domain Coloring**: A warped checkerboard pattern illustrating the transformation of space. +- **Cinematic Experience**: Supports auto-orbit, pulsing animations, and high-resolution screenshot export. +- **Cross-Platform**: Compiles natively on Windows (MinGW) and Linux (Ubuntu). + +## Usage & Controls + +### General Navigation +| Key | Action | +|---|---| +| `Alt + Enter` | Toggle Fullscreen | +| `Escape` | Exit Application | +| `H` | Toggle HUD Overlays | +| `S` / `F1` | Open Settings Panel | +| `P` | Export Screenshot (saved as `.bmp`) | +| `Ctrl + R` | Randomize all singularities | +| `Ctrl + Z / Y` | Undo / Redo changes | + +### Interaction & Manipulation +| Input | Action | +|---|---| +| `G` | Toggle Select ↔ Move Mode | +| `Z` / `X` | Add Zero / Pole at cursor (Move Mode) | +| `Delete` | Remove selected singularity | +| `Mouse Wheel` | Zoom in/out at cursor position | +| `Middle Click Drag` | Pan around the complex plane | +| `Right Click Drag` | Drag to select a group of singularities | +| `1` - `9` | Load mathematical presets (Ring, Elliptic, Cross, etc.) | + +### Visual & Audio Effects +| Key | Action | +|---|---| +| `Tab` | Cycle through 6 distinct Color Maps | +| `Backspace`| Cycle Normalization (Linear, Log, Steps) | +| `C` | Cycle Entity Rendering (Solid, Transparent, Dashed, etc.) | +| `O` | Toggle Auto-Orbit Animation | +| `U` | Toggle Magnitude Pulsation | +| `F` | Toggle Vector Field / Phase Trails | +| `D` | Toggle Domain Coloring Overlay | +| `A` | **Toggle Audio Sonification** (Theremin Mode) | + +## Building + +### Windows (MinGW/MSYS2) +Ensure you have `gcc`, `make`, and `SDL2` installed via MSYS2. +```bash +make +./bin/singularity.exe +``` + +### Linux (Ubuntu) +```bash +sudo apt install libsdl2-dev build-essential +make +./bin/singularity +``` + +--- + +*Visualization is not just about seeing; it's about understanding the rhythmic beauty of thought.* diff --git a/images/Screenshot 2025-08-15 200518.png b/images/Screenshot 2025-08-15 200518.png new file mode 100644 index 0000000..2b965ea Binary files /dev/null and b/images/Screenshot 2025-08-15 200518.png differ diff --git a/images/Screenshot 2025-08-15 200537.png b/images/Screenshot 2025-08-15 200537.png new file mode 100644 index 0000000..325fe85 Binary files /dev/null and b/images/Screenshot 2025-08-15 200537.png differ diff --git a/images/Screenshot 2025-08-15 200545.png b/images/Screenshot 2025-08-15 200545.png new file mode 100644 index 0000000..6423cbf Binary files /dev/null and b/images/Screenshot 2025-08-15 200545.png differ diff --git a/images/Screenshot 2025-08-15 200554.png b/images/Screenshot 2025-08-15 200554.png new file mode 100644 index 0000000..57c6987 Binary files /dev/null and b/images/Screenshot 2025-08-15 200554.png differ diff --git a/images/Screenshot 2025-08-15 200601.png b/images/Screenshot 2025-08-15 200601.png new file mode 100644 index 0000000..0ac1594 Binary files /dev/null and b/images/Screenshot 2025-08-15 200601.png differ diff --git a/images/Screenshot 2025-08-15 200607.png b/images/Screenshot 2025-08-15 200607.png new file mode 100644 index 0000000..e77196f Binary files /dev/null and b/images/Screenshot 2025-08-15 200607.png differ diff --git a/include/app.h b/include/app.h new file mode 100644 index 0000000..26a2ee4 --- /dev/null +++ b/include/app.h @@ -0,0 +1,38 @@ +#ifndef APP_H +#define APP_H + +#include +#include "img_utils.h" + +typedef struct { + SDL_Window *window; + SDL_Renderer *renderer; + SDL_Texture *H_tex; + SDL_Texture *UI_tex; + img_t H_img; + img_t UI_img; + + int screen_w, screen_h; /* current window pixel dimensions */ + int grid_w, grid_h; /* low-res math grid (screen / GRID_DIVISOR) */ + int is_fullscreen; + int windowed_w, windowed_h; /* saved windowed size for toggle */ + + double x_range[2]; /* [x_min, x_max] */ + double y_range[2]; /* [y_min, y_max] */ + + /* Camera for zoom/pan */ + double cam_cx, cam_cy; /* center in complex space */ + double cam_zoom; /* multiplier: >1 = zoomed in */ + + SDL_Cursor *cursor_arrow; + SDL_Cursor *cursor_hand; + SDL_Cursor *cursor_move; +} App; + +/* Returns 0 on success, non-zero on failure */ +int app_init(App *a); +void app_toggle_fullscreen(App *a); +void app_handle_resize(App *a, int w, int h); +void app_destroy(App *a); + +#endif /* APP_H */ diff --git a/include/audio.h b/include/audio.h new file mode 100644 index 0000000..666477e --- /dev/null +++ b/include/audio.h @@ -0,0 +1,21 @@ +#ifndef AUDIO_H +#define AUDIO_H + +#include + +/* Initialize SDL audio and start the callback thread */ +void audio_init(void); + +/* Shutdown the audio device */ +void audio_shutdown(void); + +/* Smoothly update the target frequency and amplitude. + freq: in Hz. amp: 0.0 to 1.0. */ +void audio_set_params(double freq, double amp); + +/* Toggle synthesis on/off */ +void audio_set_enabled(int enabled); + +void audio_get_params(double *freq, double *amp); + +#endif /* AUDIO_H */ diff --git a/include/cpu_detect.h b/include/cpu_detect.h new file mode 100644 index 0000000..b65fcb0 --- /dev/null +++ b/include/cpu_detect.h @@ -0,0 +1,17 @@ +#ifndef _CPU_DETECT_HH +#define _CPU_DETECT_HH + +#include + +typedef enum { + CPU_FEATURE_NONE = 0, + CPU_FEATURE_SSE2 = 1 << 0, + CPU_FEATURE_AVX = 1 << 1, + CPU_FEATURE_AVX2 = 1 << 2, + CPU_FEATURE_AVX512F = 1 << 3, +} cpu_features_t; + +cpu_features_t detect_cpu_features(void); +const char* get_simd_level_name(cpu_features_t features); + +#endif diff --git a/include/font_utils.h b/include/font_utils.h new file mode 100644 index 0000000..5e96f80 --- /dev/null +++ b/include/font_utils.h @@ -0,0 +1,15 @@ +#ifndef _FONT_UTILS_HH +#define _FONT_UTILS_HH + +#include +#include "img_utils.h" + +void draw_char(img_t *img, char c, uint32_t x, uint32_t y); +void draw_char_colored(img_t *img, char c, uint32_t x, uint32_t y, uint32_t color); +void draw_string_colored(img_t *img, const char *str, uint32_t x, uint32_t y, uint32_t color); +void draw_string(img_t *img, const char *str, uint32_t x, uint32_t y); + +void draw_char_dynamic(img_t *img, char c, uint32_t x, uint32_t y, int primary_channel); +void draw_string_dynamic(img_t *img, const char *str, uint32_t x, uint32_t y, int primary_channel); + +#endif /* _FONT_UTILS_HH */ diff --git a/include/img_utils.h b/include/img_utils.h new file mode 100644 index 0000000..378a325 --- /dev/null +++ b/include/img_utils.h @@ -0,0 +1,59 @@ +#ifndef _IMG_UTILS_HH +#define _IMG_UTILS_HH + +#include +#include +#include + +struct img_t; + +typedef struct img_t { + uint32_t *data; + size_t height; + size_t width; + struct img_t *bg_ref; +} img_t; + +img_t *create_img(size_t width, size_t height, img_t *img); + +static inline void channel_color(img_t *img, int32_t x, int32_t y, int ch, uint8_t *tr, uint8_t *tg, uint8_t *tb) { + uint32_t bg_pixel = 0; + if (img->bg_ref) { + /* Sample from background reference (scaled) */ + int32_t bx = x * (int32_t)img->bg_ref->width / (int32_t)img->width; + int32_t by = y * (int32_t)img->bg_ref->height / (int32_t)img->height; + if ((uint32_t)bx < img->bg_ref->width && (uint32_t)by < img->bg_ref->height) { + bg_pixel = img->bg_ref->data[by * img->bg_ref->width + bx]; + } + } else { + /* Sample from current image */ + uint32_t idx = (uint32_t)y * img->width + (uint32_t)x; + bg_pixel = img->data[idx]; + } + + uint8_t r = (bg_pixel >> 16) & 0xFF, g = (bg_pixel >> 8) & 0xFF, b = bg_pixel & 0xFF; + if (ch == 1) { *tr = 255; *tg = (255-g)/2; *tb = (255-b)/2; } + else if (ch == 3) { *tr = (255-r)/2; *tg = (255-g)/2; *tb = 255; } + else { *tr = (255-r); *tg = 255; *tb = (255-b); } +} + +/* Basic shapes */ +void draw_circle_aa(img_t *img, int32_t cx, int32_t cy, int32_t radius, uint32_t color); +void draw_circle_aa_dynamic(img_t *img, int32_t cx, int32_t cy, int32_t radius, int primary_channel); +void draw_circle_aa_alpha(img_t *img, int32_t cx, int32_t cy, int32_t radius, int primary_channel, uint8_t alpha_scale); +void draw_circle_dashed_aa_dynamic(img_t *img, int32_t cx, int32_t cy, int32_t radius, int primary_channel); +void draw_circle_glow(img_t *img, int32_t cx, int32_t cy, int32_t radius, int primary_channel); +void draw_circle_glow_color(img_t *img, int32_t cx, int32_t cy, int32_t radius, uint32_t color); + +void fill_circle(img_t *img, int32_t cx, int32_t cy, int32_t radius, uint32_t color); +void fill_circle_dynamic(img_t *img, int32_t cx, int32_t cy, int32_t radius, int primary_channel); +void fill_circle_glow(img_t *img, int32_t cx, int32_t cy, int32_t radius, int primary_channel); + +void draw_rect_fill(img_t *img, int32_t x, int32_t y, int32_t w, int32_t h, uint32_t color); +void draw_rect_blend(img_t *img, int32_t x, int32_t y, int32_t w, int32_t h, uint32_t color); +void draw_rect_stroke(img_t *img, int32_t x, int32_t y, int32_t w, int32_t h, uint32_t color); + +void draw_rounded_rect_blend(img_t *img, int32_t x, int32_t y, int32_t w, int32_t h, int32_t r, uint32_t color); +void draw_rounded_rect_fill(img_t *img, int32_t x, int32_t y, int32_t w, int32_t h, int32_t r, uint32_t color); + +#endif /* _IMG_UTILS_HH */ \ No newline at end of file diff --git a/include/presets.h b/include/presets.h new file mode 100644 index 0000000..12470da --- /dev/null +++ b/include/presets.h @@ -0,0 +1,12 @@ +#ifndef PRESETS_H +#define PRESETS_H + +#include "simulation.h" + +/* + * Load a named preset configuration into the simulation. + * n should be between 1 and 9. + */ +void load_preset(Simulation *sim, int n); + +#endif /* PRESETS_H */ diff --git a/include/simd_avx.h b/include/simd_avx.h new file mode 100644 index 0000000..b6edacf --- /dev/null +++ b/include/simd_avx.h @@ -0,0 +1,13 @@ +#ifndef _SIMD_AVX_HH +#define _SIMD_AVX_HH + +#include "transfer_function.h" +#include + +float **avx_compute_H( + float x_range[2], float y_range[2], + new_singularity_array_t zeros_arr, new_singularity_array_t poles_arr, + uint16_t height, uint16_t width, + float **H); + +#endif \ No newline at end of file diff --git a/include/simd_avx512.h b/include/simd_avx512.h new file mode 100644 index 0000000..474d549 --- /dev/null +++ b/include/simd_avx512.h @@ -0,0 +1,13 @@ +#ifndef _SIMD_AVX512_HH +#define _SIMD_AVX512_HH + +#include "transfer_function.h" +#include + +float **avx512_compute_H( + float x_range[2], float y_range[2], + new_singularity_array_t zeros_arr, new_singularity_array_t poles_arr, + uint16_t height, uint16_t width, + float **H); + +#endif diff --git a/include/simd_sse.h b/include/simd_sse.h new file mode 100644 index 0000000..59d4f5d --- /dev/null +++ b/include/simd_sse.h @@ -0,0 +1,19 @@ +#ifndef _SIMD_SSE_HH +#define _SIMD_SSE_HH + +#include "transfer_function.h" +#include "img_utils.h" +#include +#include + +float **sse2_compute_H( + float x_range[2], float y_range[2], + new_singularity_array_t zeros_arr, new_singularity_array_t poles_arr, + uint16_t height, uint16_t width, + float **H); + +float **sse_normalize_H( + const float **H, size_t size, + float **normalized); + +#endif \ No newline at end of file diff --git a/include/simulation.h b/include/simulation.h new file mode 100644 index 0000000..a1ea1ee --- /dev/null +++ b/include/simulation.h @@ -0,0 +1,44 @@ +#ifndef SIMULATION_H +#define SIMULATION_H + +#include "transfer_function.h" +#include + +typedef struct { + singularity_array_t zeros; + singularity_array_t poles; + + /* Persistent SOA buffers — grown with realloc, never freed until shutdown */ + singularity_soa_t zeros_soa; + singularity_soa_t poles_soa; + size_t soa_cap_z; + size_t soa_cap_p; + + /* Output buffers for normalised H (persistent) */ + float **H_simd; + double *nH_re; + double *nH_im; + size_t nH_cap; + + /* s-grid (rebuilt on window resize) */ + s_grid_soa_t s_grid; +} Simulation; + +void sim_init(Simulation *s, double x_range[2], double y_range[2], + int grid_w, int grid_h); +void sim_rebuild_grid(Simulation *s, double x_range[2], double y_range[2], + int grid_w, int grid_h); +void sim_randomise(Simulation *s, double x_range[2], double y_range[2]); + +singularity_t *sim_add_zero(Simulation *s, double complex pos); +singularity_t *sim_add_pole(Simulation *s, double complex pos); +int sim_delete_entity(Simulation *s, singularity_t *entity); +singularity_t *sim_find_entity(Simulation *s, double complex pos, + double threshold, int *out_type); + +void sim_compute(Simulation *s, int grid_w, int grid_h); +void sim_destroy(Simulation *s); + +double complex sim_eval_H(Simulation *s, double complex pos); + +#endif /* SIMULATION_H */ diff --git a/include/transfer_function.h b/include/transfer_function.h index c6dee4b..838c535 100644 --- a/include/transfer_function.h +++ b/include/transfer_function.h @@ -1,33 +1,103 @@ -#ifndef _TF_HH -#define _TF_HH - -#include -#include -#include -#include - -// for s -double complex *generate_s_grid(uint16_t height, uint16_t width, double y_range[2], double x_range[2]); - -// for poles and zeros -typedef struct -{ - double complex val; - uint8_t e; - double m; - double complex c; -} singularity_t; - -typedef struct -{ - singularity_t *data; - size_t size; - size_t capacity; -} singularity_array_t; - -void create_singularities(singularity_array_t *arr, size_t initial_capacity); -void add_singularity(singularity_array_t *arr, singularity_t s); -void remove_singularity(singularity_array_t *arr, size_t index); -void free_singularities(singularity_array_t *arr); - +#ifndef _TF_HH +#define _TF_HH + +#include +#include +#include +#include +#include +#include +#include +#include "img_utils.h" + +#define MAX(a, b) ((a) > (b) ? (a) : (b)) + +// for s +double complex *generate_s_grid(uint16_t height, uint16_t width, double y_range[2], double x_range[2]); + +// SoA (Structure of Arrays) version for better SIMD performance +typedef struct { + float *real; + float *imag; + uint32_t length; +} s_grid_soa_t; + +s_grid_soa_t generate_s_grid_soa(uint16_t height, uint16_t width, double y_range[2], double x_range[2]); +void free_s_grid_soa(s_grid_soa_t *grid); + +// for poles and zeros +typedef struct + +{ + double complex val; + uint8_t e; + double complex m; + double complex c; +} singularity_t; + +typedef struct +{ + singularity_t *data; + size_t size; + size_t capacity; +} singularity_array_t; + +// SoA structure for singularities (double precision) +typedef struct { + double *real; + double *imag; + double *m_re; + double *m_im; + double *c_re; + double *c_im; + uint8_t *e; + size_t size; +} singularity_soa_t; + +// SoA structure for SIMD (float precision) +typedef struct { + float *r; + float *i; + uint8_t *e; + float *m_r; + float *m_i; + float *c_r; + float *c_i; + size_t count; +} new_singularity_array_t; + +/* generate_s_grid declared above (line 16) */ + +void create_singularities(singularity_array_t *arr, size_t initial_capacity); +void add_singularity(singularity_array_t *arr, singularity_t s); +void remove_singularity(singularity_array_t *arr, size_t index); +void free_singularities(singularity_array_t *arr); + +// Normalize using SoA (real/imag arrays) +void normalize_H_soa(double *in_re, double *in_im, size_t size, double *out_re, double *out_im); +void normalize_H_log_soa(double *in_re, double *in_im, size_t size, double *out_re, double *out_im); +void normalize_H_log_steps_soa(double *in_re, double *in_im, size_t size, int steps, double *out_re, double *out_im); + +// Visualization using SoA +void H_g_img(double *re, double *im, size_t size, img_t H_img); +void H_c1_img(double *re, double *im, size_t size, img_t H_img); +void H_c2_img(double *re, double *im, size_t size, img_t H_img); +void H_c3_img(double *re, double *im, size_t size, img_t H_img); +void H_c4_img(double *re, double *im, size_t size, img_t H_img); +void H_c5_img(double *re, double *im, size_t size, img_t H_img); + +// Compute using SoA singularities +float **compute_H_simd_direct(s_grid_soa_t grid, singularity_soa_t zeros, singularity_soa_t poles, uint16_t height, uint16_t width, float **H); + +// SIMD data structure conversion +new_singularity_array_t convert_to_simd_array(const singularity_array_t *arr); +void free_simd_array(new_singularity_array_t *arr); + +double complex *compute_H(double complex *s_grid, singularity_array_t zeros_arr, singularity_array_t poles_arr, uint16_t height, uint16_t width, double complex *H); + +// Auto-dispatch function that selects best SIMD implementation +double complex *compute_H_auto(double complex *s_grid, singularity_array_t zeros_arr, singularity_array_t poles_arr, uint16_t height, uint16_t width, double complex *H); + + + #endif \ No newline at end of file diff --git a/include/ui.h b/include/ui.h new file mode 100644 index 0000000..90f3fe3 --- /dev/null +++ b/include/ui.h @@ -0,0 +1,109 @@ +#ifndef UI_H +#define UI_H + +#include "img_utils.h" +#include "transfer_function.h" +#include "simulation.h" + +/* ── Circle display mode (C key cycles) ─────────────────────────────── */ +typedef enum { + CIRCLE_SOLID = 0, /* filled dot + solid AA circle */ + CIRCLE_TRANSPARENT = 1, /* filled dot + low-alpha circle */ + CIRCLE_DASHED = 2, /* filled dot + dashed circle */ + CIRCLE_DOT_ONLY = 3, /* filled dot, no circle */ + CIRCLE_NONE = 4, /* nothing drawn */ + CIRCLE_MODE_COUNT = 5 +} CircleMode; + +/* ── Runtime display flags ───────────────────────────────────────────── */ +typedef struct { + int hud_visible; /* H key: grid + FPS + entities visible */ + CircleMode circle_mode; /* C key: entity rendering style */ + int screenshot_queued; /* 1 if screenshot requested */ + int orbit_mode; /* O key: auto rotate entities */ + int pulse_mode; /* U key: pulse entity magnitude */ + int fieldlines_visible; /* F key: draw phase trails */ + int domain_overlay; /* D key: domain coloring grid */ + int audio_enabled; /* A key: sonification */ +} UIFlags; + +/* ── Property popup panel ────────────────────────────────────────────── */ +typedef struct { + int is_open; + int x, y; + singularity_t *entity; + int entity_type; /* 1=zero, 2=pole */ + int edit_mode; /* 0=cartesian, 1=polar */ + int dragging_slider; /* -1 = none */ +} Popup; + +/* ── Settings panel ──────────────────────────────────────────────────── */ +typedef struct { + int is_open; + int drag_steps; /* -1 = not dragging */ +} SettingsPanel; + +/* ── Render config ───────────────────────────────────────────────────── */ +typedef struct { + int c_map; /* colour map 0-5 */ + int n_map; /* normalisation mode 0-2 */ + int steps; /* step count for log-steps */ +} RenderCfg; + +/* ── Drawing functions ───────────────────────────────────────────────── */ + +void ui_draw_grid(img_t *img, + double x0, double x1, double y0, double y1); + +/* Draw all entities; selected/move_target get a highlight ring */ +void ui_draw_entities(img_t *img, + Simulation *sim, + double x_range[2], double y_range[2], + singularity_t *selected, + singularity_t *move_target, + UIFlags flags); + +void ui_draw_tooltip(img_t *img, int mx, int my, + singularity_t *s, int type); + +void ui_draw_popup(img_t *img, Popup *p); + +void ui_draw_hud(img_t *img, float fps, const RenderCfg *cfg, + const UIFlags *flags, int cursor_mode, int is_fullscreen); + +/* settings_toggle_fs: if non-NULL and user presses fs button, set *=1 */ +void ui_draw_settings(img_t *img, SettingsPanel *sp, + RenderCfg *cfg, UIFlags *flags, + int screen_w, int screen_h, int is_fullscreen, + int *request_fs_toggle); + +/* ── Popup interaction ───────────────────────────────────────────────── */ + +/* Returns 1 and sets *out_idx if (mx,my) hits a slider area */ +int popup_hit_slider(const Popup *p, int mx, int my, int *out_idx); +/* Returns 1 if (mx,my) hits the mode toggle button */ +int popup_hit_mode_btn(const Popup *p, int mx, int my); +/* Apply slider drag at screen-x mx to entity */ +void popup_apply_slider(Popup *p, int mx); + +/* ── Settings interaction ────────────────────────────────────────────── */ +/* Returns action flags (bitmask) for what settings widget was clicked */ +#define SETTINGS_HIT_NONE 0 +#define SETTINGS_HIT_FS (1<<0) +#define SETTINGS_HIT_CMAP (1<<1) +#define SETTINGS_HIT_NMAP (1<<2) +#define SETTINGS_HIT_STEPS (1<<3) +#define SETTINGS_HIT_CLOSE (1<<4) +int settings_hit(SettingsPanel *sp, RenderCfg *cfg, UIFlags *flags, + int mx, int my, int screen_w, int screen_h, + int is_fullscreen, int *cmap_delta, int *nmap_delta); + +/* ── Vector Field / Phase Trails ─────────────────────────────────────── */ +void ui_draw_field_lines(img_t *img, Simulation *sim, + double x_range[2], double y_range[2]); + +/* ── Domain Coloring ─────────────────────────────────────────────────── */ +void ui_draw_domain_overlay(img_t *img, Simulation *sim, + double x_range[2], double y_range[2]); + +#endif /* UI_H */ diff --git a/include/undo.h b/include/undo.h new file mode 100644 index 0000000..37360a1 --- /dev/null +++ b/include/undo.h @@ -0,0 +1,31 @@ +#ifndef UNDO_H +#define UNDO_H + +#include "transfer_function.h" + +#define UNDO_MAX 128 + +typedef struct { + singularity_array_t zeros; + singularity_array_t poles; +} UndoFrame; + +typedef struct { + UndoFrame frames[UNDO_MAX]; + int pos; /* index of the current snapshot (-1 = none) */ + int total; /* number of valid snapshots */ +} UndoStack; + +void undo_init(UndoStack *u); +void undo_push(UndoStack *u, + const singularity_array_t *zeros, + const singularity_array_t *poles); +int undo_undo(UndoStack *u, + singularity_array_t *zeros, + singularity_array_t *poles); +int undo_redo(UndoStack *u, + singularity_array_t *zeros, + singularity_array_t *poles); +void undo_destroy(UndoStack *u); + +#endif /* UNDO_H */ diff --git a/include/utils.h b/include/utils.h new file mode 100644 index 0000000..8efa7e3 --- /dev/null +++ b/include/utils.h @@ -0,0 +1,22 @@ +#ifndef _UTILS_HH +#define _UTILS_HH + +#include +#include +#include +#include +#include + +int gcd(int a, int b); +double random_uniform(double min, double max); + +/* screen_w / screen_h as int to support large resolutions on resize */ +void screen_choords_to_complex(double x, double y, double complex *c, + double x_range[2], double y_range[2], + int width, int height); + +void complex_to_screen_choords(double complex c, uint32_t *x, uint32_t *y, + double x_range[2], double y_range[2], + int width, int height); + +#endif /* _UTILS_HH */ \ No newline at end of file diff --git a/makefile b/makefile new file mode 100644 index 0000000..6bebe82 --- /dev/null +++ b/makefile @@ -0,0 +1,60 @@ +# Configuration +CC = gcc +CFLAGS = -Wall -O3 -march=native -std=c11 -Iinclude -D_REENTRANT +LDFLAGS = -lm + +# OS-specific settings +ifeq ($(OS),Windows_NT) + CFLAGS += -I/mingw64/include/SDL2 + LDFLAGS += -L/mingw64/lib -lmingw32 -lSDL2main -lSDL2 + EXE_EXT = .exe +else + # Linux (Ubuntu) / macOS + CFLAGS += $(shell sdl2-config --cflags) + LDFLAGS += $(shell sdl2-config --libs) + EXE_EXT = +endif + +SRC_DIR = src +BUILD_DIR = build +BIN_DIR = bin +EXE_NAME = singularity +EXE = $(BIN_DIR)/$(EXE_NAME)$(EXE_EXT) + +# Get source files and corresponding object files +SRCS := $(wildcard $(SRC_DIR)/*.c) +OBJS := $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRCS)) + +# Default target +all: $(EXE) + +# Create bin and build directories if missing +$(EXE): $(OBJS) + @mkdir -p $(BIN_DIR) + $(CC) -o $@ $^ $(LDFLAGS) + @echo "Build complete: $(EXE)" + +# Compile source to object +$(BUILD_DIR)/%.o: $(SRC_DIR)/%.c + @mkdir -p $(BUILD_DIR) + $(CC) $(CFLAGS) -c $< -o $@ + +# Bundle SDL2.dll with executable (Windows only) +bundle-sdl2: all +ifeq ($(OS),Windows_NT) + @echo "Looking for SDL2.dll..." + @cp /mingw64/bin/SDL2.dll $(BIN_DIR)/ 2>/dev/null && echo "SDL2.dll copied to $(BIN_DIR)/" || echo "Warning: SDL2.dll not found in /mingw64/bin/" +else + @echo "Library bundling is not required on this OS (use dynamic linking)." +endif + +# Create distribution package +dist: bundle-sdl2 + @echo "Distribution package ready in $(BIN_DIR)/" + @ls -lh $(BIN_DIR)/ + +# Clean target +clean: + rm -rf $(BUILD_DIR) $(BIN_DIR) + +.PHONY: all clean bundle-sdl2 dist diff --git a/src/app.c b/src/app.c new file mode 100644 index 0000000..f706c30 --- /dev/null +++ b/src/app.c @@ -0,0 +1,143 @@ +#include "app.h" +#include "utils.h" +#include +#include +#include + +#define GRID_DIVISOR 4 +#define DEFAULT_WIN_W 1280 +#define DEFAULT_WIN_H 720 + +static void compute_ranges(App *a) { + /* Define a base scale: how many units the vertical axis covers in landscape */ + const double BASE_Y = 10.0; + double aspect = (double)a->screen_w / (double)a->screen_h; + + if (aspect >= 1.0) { + /* Landscape or square: fix vertical, expand horizontal */ + a->y_range[0] = a->cam_cy - (BASE_Y / a->cam_zoom); + a->y_range[1] = a->cam_cy + (BASE_Y / a->cam_zoom); + a->x_range[0] = a->cam_cx - (BASE_Y * aspect / a->cam_zoom); + a->x_range[1] = a->cam_cx + (BASE_Y * aspect / a->cam_zoom); + } else { + /* Portrait: fix horizontal, expand vertical */ + a->x_range[0] = a->cam_cx - (BASE_Y / a->cam_zoom); + a->x_range[1] = a->cam_cx + (BASE_Y / a->cam_zoom); + a->y_range[0] = a->cam_cy - (BASE_Y / (aspect * a->cam_zoom)); + a->y_range[1] = a->cam_cy + (BASE_Y / (aspect * a->cam_zoom)); + } +} + +static void alloc_framebuffers(App *a) { + free(a->H_img.data); + free(a->UI_img.data); + a->H_img.data = NULL; + a->UI_img.data = NULL; + + a->grid_w = a->screen_w / GRID_DIVISOR; + a->grid_h = a->screen_h / GRID_DIVISOR; + create_img((size_t)a->grid_w, (size_t)a->grid_h, &a->H_img); + create_img((size_t)a->screen_w, (size_t)a->screen_h, &a->UI_img); + a->UI_img.bg_ref = &a->H_img; +} + +static void rebuild_textures(App *a) { + if (a->H_tex) SDL_DestroyTexture(a->H_tex); + if (a->UI_tex) SDL_DestroyTexture(a->UI_tex); + + a->H_tex = SDL_CreateTexture(a->renderer, SDL_PIXELFORMAT_ARGB8888, + SDL_TEXTUREACCESS_STREAMING, + (int)a->H_img.width, (int)a->H_img.height); + a->UI_tex = SDL_CreateTexture(a->renderer, SDL_PIXELFORMAT_ARGB8888, + SDL_TEXTUREACCESS_STREAMING, + (int)a->UI_img.width, (int)a->UI_img.height); + SDL_SetTextureBlendMode(a->UI_tex, SDL_BLENDMODE_BLEND); +} + +int app_init(App *a) { + memset(a, 0, sizeof(*a)); + + if (SDL_Init(SDL_INIT_VIDEO | SDL_INIT_AUDIO) != 0) { + fprintf(stderr, "SDL_Init Error: %s\n", SDL_GetError()); + return 1; + } + + a->cursor_arrow = SDL_CreateSystemCursor(SDL_SYSTEM_CURSOR_ARROW); + a->cursor_hand = SDL_CreateSystemCursor(SDL_SYSTEM_CURSOR_HAND); + a->cursor_move = SDL_CreateSystemCursor(SDL_SYSTEM_CURSOR_SIZEALL); + + SDL_DisplayMode dm; + if (SDL_GetCurrentDisplayMode(0, &dm) != 0) { + fprintf(stderr, "SDL_GetCurrentDisplayMode Error: %s\n", SDL_GetError()); + SDL_Quit(); + return 1; + } + + a->screen_w = dm.w; + a->screen_h = dm.h; + a->windowed_w = DEFAULT_WIN_W; + a->windowed_h = DEFAULT_WIN_H; + a->is_fullscreen = 1; + + a->cam_zoom = 1.0; + a->cam_cx = 0.0; + a->cam_cy = 0.0; + + a->window = SDL_CreateWindow("SINGULARITY", + SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED, + a->screen_w, a->screen_h, + SDL_WINDOW_FULLSCREEN_DESKTOP | SDL_WINDOW_RESIZABLE); + if (!a->window) { + fprintf(stderr, "SDL_CreateWindow Error: %s\n", SDL_GetError()); + SDL_Quit(); + return 1; + } + + a->renderer = SDL_CreateRenderer(a->window, -1, SDL_RENDERER_ACCELERATED); + if (!a->renderer) { + fprintf(stderr, "SDL_CreateRenderer Error: %s\n", SDL_GetError()); + SDL_DestroyWindow(a->window); + SDL_Quit(); + return 1; + } + + compute_ranges(a); + alloc_framebuffers(a); + rebuild_textures(a); + return 0; +} + +void app_toggle_fullscreen(App *a) { + if (a->is_fullscreen) { + SDL_SetWindowFullscreen(a->window, 0); + SDL_SetWindowSize(a->window, a->windowed_w, a->windowed_h); + SDL_SetWindowPosition(a->window, SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED); + a->is_fullscreen = 0; + } else { + SDL_GetWindowSize(a->window, &a->windowed_w, &a->windowed_h); + SDL_SetWindowFullscreen(a->window, SDL_WINDOW_FULLSCREEN_DESKTOP); + a->is_fullscreen = 1; + } +} + +void app_handle_resize(App *a, int w, int h) { + if (w <= 0 || h <= 0) return; + a->screen_w = w; + a->screen_h = h; + compute_ranges(a); + alloc_framebuffers(a); + rebuild_textures(a); +} + +void app_destroy(App *a) { + free(a->H_img.data); + free(a->UI_img.data); + if (a->H_tex) SDL_DestroyTexture(a->H_tex); + if (a->UI_tex) SDL_DestroyTexture(a->UI_tex); + if (a->renderer) SDL_DestroyRenderer(a->renderer); + if (a->window) SDL_DestroyWindow(a->window); + SDL_FreeCursor(a->cursor_arrow); + SDL_FreeCursor(a->cursor_hand); + SDL_FreeCursor(a->cursor_move); + SDL_Quit(); +} diff --git a/src/audio.c b/src/audio.c new file mode 100644 index 0000000..686dbbb --- /dev/null +++ b/src/audio.c @@ -0,0 +1,77 @@ +#include "audio.h" +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +static SDL_AudioDeviceID audio_device = 0; +static double phase = 0.0; +static double current_freq = 440.0; +static double target_freq = 440.0; +static double current_amp = 0.0; +static double target_amp = 0.0; +static int is_enabled = 0; + +static void SDLCALL audio_callback(void *userdata, Uint8 *stream, int len) { + (void)userdata; + float *fstream = (float *)stream; + int samples = len / sizeof(float); + + for (int i = 0; i < samples; i++) { + /* Smooth parameter transitions to avoid clicking */ + current_freq += (target_freq - current_freq) * 0.005; + current_amp += (target_amp - current_amp) * 0.005; + + if (is_enabled && current_amp > 0.001) { + fstream[i] = (float)(sin(phase) * current_amp * 0.2); // 0.2 limits max volume safely + phase += 2.0 * M_PI * current_freq / 44100.0; + if (phase > 2.0 * M_PI) phase -= 2.0 * M_PI; + } else { + fstream[i] = 0.0f; + } + } +} + +void audio_init(void) { + SDL_AudioSpec want, have; + SDL_zero(want); + want.freq = 44100; + want.format = AUDIO_F32; + want.channels = 1; + want.samples = 1024; + want.callback = audio_callback; + + audio_device = SDL_OpenAudioDevice(NULL, 0, &want, &have, 0); + if (audio_device != 0) { + SDL_PauseAudioDevice(audio_device, 0); + } else { + SDL_Log("Failed to open audio: %s", SDL_GetError()); + } +} + +void audio_shutdown(void) { + if (audio_device != 0) { + SDL_CloseAudioDevice(audio_device); + audio_device = 0; + } +} + +void audio_set_params(double freq, double amp) { + target_freq = freq; + if (target_freq < 20.0) target_freq = 20.0; + if (target_freq > 2000.0) target_freq = 2000.0; + + target_amp = amp; + if (target_amp < 0.0) target_amp = 0.0; + if (target_amp > 1.0) target_amp = 1.0; +} + +void audio_set_enabled(int enabled) { + is_enabled = enabled; +} + +void audio_get_params(double *freq, double *amp) { + *freq = current_freq; + *amp = current_amp; +} diff --git a/src/cpu_detect.c b/src/cpu_detect.c new file mode 100644 index 0000000..7cd39a1 --- /dev/null +++ b/src/cpu_detect.c @@ -0,0 +1,64 @@ +#include "cpu_detect.h" +#include + +#ifdef _MSC_VER +#include +#else +#include +#endif + +static void run_cpuid(uint32_t eax, uint32_t ecx, uint32_t* regs) { +#ifdef _MSC_VER + __cpuidex((int*)regs, (int)eax, (int)ecx); +#else + __cpuid_count(eax, ecx, regs[0], regs[1], regs[2], regs[3]); +#endif +} + +cpu_features_t detect_cpu_features(void) { + cpu_features_t features = CPU_FEATURE_NONE; + uint32_t regs[4]; + + run_cpuid(0, 0, regs); + uint32_t max_level = regs[0]; + + if (max_level >= 1) { + run_cpuid(1, 0, regs); + + if (regs[3] & (1 << 26)) { + features |= CPU_FEATURE_SSE2; + } + + if (regs[2] & (1 << 28)) { + features |= CPU_FEATURE_AVX; + } + } + + if (max_level >= 7) { + run_cpuid(7, 0, regs); + + if (regs[1] & (1 << 5)) { + features |= CPU_FEATURE_AVX2; + } + + if (regs[1] & (1 << 16)) { + features |= CPU_FEATURE_AVX512F; + } + } + + return features; +} + +const char* get_simd_level_name(cpu_features_t features) { + if (features & CPU_FEATURE_AVX512F) { + return "AVX-512"; + } else if (features & CPU_FEATURE_AVX2) { + return "AVX2"; + } else if (features & CPU_FEATURE_AVX) { + return "AVX"; + } else if (features & CPU_FEATURE_SSE2) { + return "SSE2"; + } else { + return "Scalar (no SIMD)"; + } +} diff --git a/src/font_utils.c b/src/font_utils.c new file mode 100644 index 0000000..23c377d --- /dev/null +++ b/src/font_utils.c @@ -0,0 +1,168 @@ +#include "font_utils.h" +#include +#include + +// Standard 5x7 ASCII font data (Space to ~) +// Each byte represents a column of 7 pixels (LSB at top) +static const uint8_t font_data[][5] = { + {0x00, 0x00, 0x00, 0x00, 0x00}, // Space + {0x00, 0x00, 0x5F, 0x00, 0x00}, // ! + {0x00, 0x07, 0x00, 0x07, 0x00}, // " + {0x14, 0x7F, 0x14, 0x7F, 0x14}, // # + {0x24, 0x2A, 0x7F, 0x2A, 0x12}, // $ + {0x23, 0x13, 0x08, 0x64, 0x62}, // % + {0x36, 0x49, 0x55, 0x22, 0x50}, // & + {0x00, 0x05, 0x03, 0x00, 0x00}, // ' + {0x00, 0x1C, 0x22, 0x41, 0x00}, // ( + {0x00, 0x41, 0x22, 0x1C, 0x00}, // ) + {0x14, 0x08, 0x3E, 0x08, 0x14}, // * + {0x08, 0x08, 0x3E, 0x08, 0x08}, // + + {0x00, 0x50, 0x30, 0x00, 0x00}, // , + {0x08, 0x08, 0x08, 0x08, 0x08}, // - + {0x00, 0x60, 0x60, 0x00, 0x00}, // . + {0x20, 0x10, 0x08, 0x04, 0x02}, // / + {0x3E, 0x51, 0x49, 0x45, 0x3E}, // 0 + {0x00, 0x42, 0x7F, 0x40, 0x00}, // 1 + {0x42, 0x61, 0x51, 0x49, 0x46}, // 2 + {0x21, 0x41, 0x45, 0x4B, 0x31}, // 3 + {0x18, 0x14, 0x12, 0x7F, 0x10}, // 4 + {0x27, 0x45, 0x45, 0x45, 0x39}, // 5 + {0x3C, 0x4A, 0x49, 0x49, 0x30}, // 6 + {0x01, 0x71, 0x09, 0x05, 0x03}, // 7 + {0x36, 0x49, 0x49, 0x49, 0x36}, // 8 + {0x06, 0x49, 0x49, 0x29, 0x1E}, // 9 + {0x00, 0x36, 0x36, 0x00, 0x00}, // : + {0x00, 0x56, 0x36, 0x00, 0x00}, // ; + {0x08, 0x14, 0x22, 0x41, 0x00}, // < + {0x14, 0x14, 0x14, 0x14, 0x14}, // = + {0x00, 0x41, 0x22, 0x14, 0x08}, // > + {0x02, 0x01, 0x51, 0x09, 0x06}, // ? + {0x32, 0x49, 0x79, 0x41, 0x3E}, // @ + {0x7E, 0x11, 0x11, 0x11, 0x7E}, // A + {0x7F, 0x49, 0x49, 0x49, 0x36}, // B + {0x3E, 0x41, 0x41, 0x41, 0x22}, // C + {0x7F, 0x41, 0x41, 0x22, 0x1C}, // D + {0x7F, 0x49, 0x49, 0x49, 0x41}, // E + {0x7F, 0x09, 0x09, 0x09, 0x01}, // F + {0x3E, 0x41, 0x49, 0x49, 0x7A}, // G + {0x7F, 0x08, 0x08, 0x08, 0x7F}, // H + {0x00, 0x41, 0x7F, 0x41, 0x00}, // I + {0x20, 0x40, 0x41, 0x3F, 0x01}, // J + {0x7F, 0x08, 0x14, 0x22, 0x41}, // K + {0x7F, 0x40, 0x40, 0x40, 0x40}, // L + {0x7F, 0x02, 0x0C, 0x02, 0x7F}, // M + {0x7F, 0x04, 0x08, 0x10, 0x7F}, // N + {0x3E, 0x41, 0x41, 0x41, 0x3E}, // O + {0x7F, 0x09, 0x09, 0x09, 0x06}, // P + {0x3E, 0x41, 0x51, 0x21, 0x5E}, // Q + {0x7F, 0x09, 0x19, 0x29, 0x46}, // R + {0x46, 0x49, 0x49, 0x49, 0x31}, // S + {0x01, 0x01, 0x7F, 0x01, 0x01}, // T + {0x3F, 0x40, 0x40, 0x40, 0x3F}, // U + {0x1F, 0x20, 0x40, 0x20, 0x1F}, // V + {0x3F, 0x40, 0x38, 0x40, 0x3F}, // W + {0x63, 0x14, 0x08, 0x14, 0x63}, // X + {0x07, 0x08, 0x70, 0x08, 0x07}, // Y + {0x61, 0x51, 0x49, 0x45, 0x43}, // Z + {0x00, 0x7F, 0x41, 0x41, 0x00}, // [ + {0x02, 0x04, 0x08, 0x10, 0x20}, // \ (Backslash) + {0x00, 0x41, 0x41, 0x7F, 0x00}, // ] + {0x04, 0x02, 0x01, 0x02, 0x04}, // ^ + {0x40, 0x40, 0x40, 0x40, 0x40}, // _ + {0x00, 0x01, 0x02, 0x04, 0x00}, // ` + {0x20, 0x54, 0x54, 0x54, 0x78}, // a + {0x7F, 0x48, 0x44, 0x44, 0x38}, // b + {0x38, 0x44, 0x44, 0x44, 0x20}, // c + {0x38, 0x44, 0x44, 0x48, 0x7F}, // d + {0x38, 0x54, 0x54, 0x54, 0x18}, // e + {0x08, 0x7E, 0x09, 0x01, 0x02}, // f + {0x0C, 0x52, 0x52, 0x52, 0x3E}, // g + {0x7F, 0x08, 0x04, 0x04, 0x78}, // h + {0x00, 0x44, 0x7D, 0x40, 0x00}, // i + {0x20, 0x40, 0x44, 0x3D, 0x00}, // j + {0x7F, 0x10, 0x28, 0x44, 0x00}, // k + {0x00, 0x41, 0x7F, 0x40, 0x00}, // l + {0x7C, 0x04, 0x18, 0x04, 0x78}, // m + {0x7C, 0x08, 0x04, 0x04, 0x78}, // n + {0x38, 0x44, 0x44, 0x44, 0x38}, // o + {0x7C, 0x14, 0x14, 0x14, 0x08}, // p + {0x08, 0x14, 0x14, 0x18, 0x7C}, // q + {0x7C, 0x08, 0x04, 0x04, 0x08}, // r + {0x48, 0x54, 0x54, 0x54, 0x20}, // s + {0x04, 0x3F, 0x44, 0x40, 0x20}, // t + {0x3C, 0x40, 0x40, 0x20, 0x7C}, // u + {0x1C, 0x20, 0x40, 0x20, 0x1C}, // v + {0x3C, 0x40, 0x30, 0x40, 0x3C}, // w + {0x44, 0x28, 0x10, 0x28, 0x44}, // x + {0x0C, 0x50, 0x50, 0x50, 0x3C}, // y + {0x44, 0x64, 0x54, 0x4C, 0x44}, // z + {0x00, 0x08, 0x36, 0x41, 0x00}, // { + {0x00, 0x00, 0x7F, 0x00, 0x00}, // | + {0x00, 0x41, 0x36, 0x08, 0x00}, // } + {0x10, 0x08, 0x08, 0x10, 0x08} // ~ +}; + +// Removed unused draw_pixel_inverted + +void draw_char_colored(img_t *img, char c, uint32_t x, uint32_t y, uint32_t color) { + const uint8_t *glyph; + if (c >= ' ' && c <= '~') { + glyph = font_data[c - ' ']; + } else { + glyph = font_data[0]; // space + } + + for (int col = 0; col < 5; col++) { + uint8_t column = glyph[col]; + for (int row = 0; row < 7; row++) { + if (column & (1 << row)) { + if (x + col < img->width && y + row < img->height) { + img->data[(y + row) * img->width + (x + col)] = color; + } + } + } + } +} + +void draw_char(img_t *img, char c, uint32_t x, uint32_t y) { + draw_char_colored(img, c, x, y, 0xFFFFFFFF); // Default white +} + +void draw_char_dynamic(img_t *img, char c, uint32_t x, uint32_t y, int primary_channel) { + const uint8_t *glyph; + if (c >= ' ' && c <= '~') glyph = font_data[c - ' ']; + else glyph = font_data[0]; + + for (int col = 0; col < 5; col++) { + uint8_t column = glyph[col]; + for (int row = 0; row < 7; row++) { + if (column & (1 << row)) { + int px = x + col; + int py = y + row; + if ((uint32_t)px < img->width && (uint32_t)py < img->height) { + uint8_t tr, tg, tb; + channel_color(img, px, py, primary_channel, &tr, &tg, &tb); + img->data[py * img->width + px] = (0xFF << 24) | (tr << 16) | (tg << 8) | tb; + } + } + } + } +} + +void draw_string_colored(img_t *img, const char *str, uint32_t x, uint32_t y, uint32_t color) { + size_t len = strlen(str); + for (size_t i = 0; i < len; i++) { + draw_char_colored(img, str[i], x + (i * 6), y, color); + } +} + +void draw_string_dynamic(img_t *img, const char *str, uint32_t x, uint32_t y, int primary_channel) { + size_t len = strlen(str); + for (size_t i = 0; i < len; i++) { + draw_char_dynamic(img, str[i], x + (i * 6), y, primary_channel); + } +} + +void draw_string(img_t *img, const char *str, uint32_t x, uint32_t y) { + draw_string_colored(img, str, x, y, 0xFFFFFFFF); +} diff --git a/src/img_utils.c b/src/img_utils.c new file mode 100644 index 0000000..8747685 --- /dev/null +++ b/src/img_utils.c @@ -0,0 +1,335 @@ +#include "img_utils.h" +#include +#include +#include + +/* ── Buffer creation ─────────────────────────────────────────────────── */ + +img_t *create_img(size_t width, size_t height, img_t *img) { + if (img == NULL) { + img = malloc(sizeof(img_t)); + if (!img) return NULL; + } + img->width = width; + img->height = height; + img->data = malloc(width * height * sizeof(uint32_t)); + if (!img->data) return NULL; + return img; +} + +/* ── Pixel blending helpers ───────────────────────────────────────────── */ + +static inline void blend_pixel_fast(img_t *img, int32_t x, int32_t y, + uint32_t color, uint32_t alpha) { + if ((uint32_t)x >= img->width || (uint32_t)y >= img->height || alpha == 0) return; + uint32_t idx = (uint32_t)y * img->width + (uint32_t)x; + + if (alpha >= 255) { + img->data[idx] = (0xFF << 24) | (color & 0xFFFFFF); + return; + } + + uint32_t bg = img->data[idx]; + uint32_t ba = (bg >> 24) & 0xFF; + + if (ba == 0) { + /* Blending on empty pixel: just set color and alpha as-is */ + img->data[idx] = (alpha << 24) | (color & 0xFFFFFF); + return; + } + + /* Blending on existing UI element: simple linear blend */ + uint32_t br = (bg >> 16) & 0xFF, bg_g = (bg >> 8) & 0xFF, bb = bg & 0xFF; + uint32_t fr = (color >> 16) & 0xFF, fg = (color >> 8) & 0xFF, fb = color & 0xFF; + uint32_t ia = 256 - alpha; + + uint32_t out_r = (fr * alpha + br * ia) >> 8; + uint32_t out_g = (fg * alpha + bg_g * ia) >> 8; + uint32_t out_b = (fb * alpha + bb * ia) >> 8; + uint32_t out_a = alpha + (ba * ia >> 8); + if (out_a > 255) out_a = 255; + + img->data[idx] = (out_a << 24) | (out_r << 16) | (out_g << 8) | out_b; +} + +static inline void blend_pixel(img_t *img, int32_t x, int32_t y, + uint32_t color, float alpha) { + blend_pixel_fast(img, x, y, color, (uint32_t)(alpha * 255.0f)); +} + +/* Derive a channel colour from the primary_channel index used by dynamic functions */ +/* (Implementation moved to img_utils.h for inlining) */ + +/* ── Solid AA circle ─────────────────────────────────────────────────── */ + +void draw_circle_aa(img_t *img, int32_t cx, int32_t cy, int32_t radius, uint32_t color) { + int32_t x0 = cx-radius-2, y0 = cy-radius-2; + int32_t x1 = cx+radius+2, y1 = cy+radius+2; + if (x0 < 0) x0 = 0; + if (y0 < 0) y0 = 0; + if (x1 >= (int32_t)img->width) x1 = img->width - 1; + if (y1 >= (int32_t)img->height) y1 = img->height - 1; + for (int32_t y = y0; y <= y1; y++) { + for (int32_t x = x0; x <= x1; x++) { + float dist = sqrtf((float)((x-cx)*(x-cx)+(y-cy)*(y-cy))); + float diff = fabsf(dist - radius); + if (diff < 1.0f) blend_pixel(img, x, y, color, 1.0f - diff); + } + } +} + +void draw_circle_aa_dynamic(img_t *img, int32_t cx, int32_t cy, + int32_t radius, int primary_channel) { + int32_t x0 = cx-radius-2, y0 = cy-radius-2; + int32_t x1 = cx+radius+2, y1 = cy+radius+2; + if (x0 < 0) x0 = 0; + if (y0 < 0) y0 = 0; + if (x1 >= (int32_t)img->width) x1 = img->width - 1; + if (y1 >= (int32_t)img->height) y1 = img->height - 1; + for (int32_t y = y0; y <= y1; y++) { + for (int32_t x = x0; x <= x1; x++) { + float dist = sqrtf((float)((x-cx)*(x-cx)+(y-cy)*(y-cy))); + float diff = fabsf(dist - radius); + if (diff >= 1.0f) continue; + float alpha = 1.0f - diff; + uint8_t tr, tg, tb; + channel_color(img, x, y, primary_channel, &tr, &tg, &tb); + blend_pixel(img, x, y, (0xFF<<24)|(tr<<16)|(tg<<8)|tb, alpha); + } + } +} + +/* ── Transparent (low-alpha) AA circle ───────────────────────────────── */ + +void draw_circle_aa_alpha(img_t *img, int32_t cx, int32_t cy, + int32_t radius, int primary_channel, uint8_t alpha_scale) { + int32_t x0 = cx-radius-2, y0 = cy-radius-2; + int32_t x1 = cx+radius+2, y1 = cy+radius+2; + if (x0 < 0) x0 = 0; + if (y0 < 0) y0 = 0; + if (x1 >= (int32_t)img->width) x1 = img->width - 1; + if (y1 >= (int32_t)img->height) y1 = img->height - 1; + float ascale = alpha_scale / 255.0f; + for (int32_t y = y0; y <= y1; y++) { + for (int32_t x = x0; x <= x1; x++) { + float dist = sqrtf((float)((x-cx)*(x-cx)+(y-cy)*(y-cy))); + float diff = fabsf(dist - radius); + if (diff >= 1.0f) continue; + float alpha = (1.0f - diff) * ascale; + uint8_t tr, tg, tb; + channel_color(img, x, y, primary_channel, &tr, &tg, &tb); + blend_pixel(img, x, y, (0xFF<<24)|(tr<<16)|(tg<<8)|tb, alpha); + } + } +} + +/* ── Dashed AA circle ────────────────────────────────────────────────── */ + +void draw_circle_dashed_aa_dynamic(img_t *img, int32_t cx, int32_t cy, + int32_t radius, int primary_channel) { + int32_t x0 = cx-radius-2, y0 = cy-radius-2; + int32_t x1 = cx+radius+2, y1 = cy+radius+2; + if (x0 < 0) x0 = 0; + if (y0 < 0) y0 = 0; + if (x1 >= (int32_t)img->width) x1 = img->width - 1; + if (y1 >= (int32_t)img->height) y1 = img->height - 1; + /* 8 dashes: every two octants alternates on/off */ + for (int32_t y = y0; y <= y1; y++) { + for (int32_t x = x0; x <= x1; x++) { + float dx = (float)(x - cx), dy = (float)(y - cy); + float dist = sqrtf(dx*dx + dy*dy); + float diff = fabsf(dist - radius); + if (diff >= 1.0f) continue; + float angle = atan2f(dy, dx); /* -PI..PI */ + if (angle < 0) angle += 6.28318530f; + /* 48 segments total = 24 dashes */ + int seg = (int)(angle / (6.28318530f / 48.0f)); + if (seg % 2 == 0) continue; + float alpha = 1.0f - diff; + uint8_t tr, tg, tb; + channel_color(img, x, y, primary_channel, &tr, &tg, &tb); + blend_pixel(img, x, y, (0xFF<<24)|(tr<<16)|(tg<<8)|tb, alpha); + } + } +} + +/* ── Glow ring (highlight for selected/moving entity) ────────────────── */ + +void draw_circle_glow(img_t *img, int32_t cx, int32_t cy, + int32_t radius, int primary_channel) { + int32_t x0 = cx-radius-8, y0 = cy-radius-8; + int32_t x1 = cx+radius+8, y1 = cy+radius+8; + if (x0 < 0) x0 = 0; + if (y0 < 0) y0 = 0; + if (x1 >= (int32_t)img->width) x1 = img->width - 1; + if (y1 >= (int32_t)img->height) y1 = img->height - 1; + for (int32_t y = y0; y <= y1; y++) { + for (int32_t x = x0; x <= x1; x++) { + float dist = sqrtf((float)((x-cx)*(x-cx)+(y-cy)*(y-cy))); + float diff = dist - (float)radius; + if (diff < 0 || diff > 8.0f) continue; + float alpha = (1.0f - diff / 8.0f) * 0.5f; + uint8_t tr, tg, tb; + channel_color(img, x, y, primary_channel, &tr, &tg, &tb); + blend_pixel(img, x, y, (0xFF<<24)|(tr<<16)|(tg<<8)|tb, alpha); + } + } +} + +void draw_circle_glow_color(img_t *img, int32_t cx, int32_t cy, + int32_t radius, uint32_t color) { + int32_t x0 = cx-radius-8, y0 = cy-radius-8; + int32_t x1 = cx+radius+8, y1 = cy+radius+8; + if (x0 < 0) x0 = 0; + if (y0 < 0) y0 = 0; + if (x1 >= (int32_t)img->width) x1 = img->width - 1; + if (y1 >= (int32_t)img->height) y1 = img->height - 1; + + float base_alpha = ((color >> 24) & 0xFF) / 255.0f; + uint32_t rgb = color & 0xFFFFFF; // Strip alpha from color payload since blend_pixel takes 0xFF alpha in color arg + + for (int32_t y = y0; y <= y1; y++) { + for (int32_t x = x0; x <= x1; x++) { + float dist = sqrtf((float)((x-cx)*(x-cx)+(y-cy)*(y-cy))); + float diff = fabsf(dist - (float)radius); + if (diff > 8.0f) continue; + float alpha = (1.0f - diff / 8.0f) * base_alpha; + blend_pixel(img, x, y, (0xFF<<24)|rgb, alpha); + } + } +} + +/* ── Filled circle ───────────────────────────────────────────────────── */ + +void fill_circle(img_t *img, int32_t cx, int32_t cy, + int32_t radius, uint32_t color) { + int32_t r2 = radius * radius; + for (int32_t y = cy-radius; y <= cy+radius; y++) { + for (int32_t x = cx-radius; x <= cx+radius; x++) { + int32_t dx = x-cx, dy = y-cy; + if (dx*dx + dy*dy <= r2) + if (x >= 0 && x < (int32_t)img->width && + y >= 0 && y < (int32_t)img->height) + img->data[y*img->width+x] = color; + } + } +} + +void fill_circle_dynamic(img_t *img, int32_t cx, int32_t cy, + int32_t radius, int primary_channel) { + int32_t r2 = radius * radius; + for (int32_t y = cy-radius; y <= cy+radius; y++) { + for (int32_t x = cx-radius; x <= cx+radius; x++) { + int32_t dx = x-cx, dy = y-cy; + if (dx*dx + dy*dy > r2) continue; + if (x < 0 || x >= (int32_t)img->width || + y < 0 || y >= (int32_t)img->height) continue; + uint32_t idx = y*img->width + x; + uint8_t tr, tg, tb; + channel_color(img, x, y, primary_channel, &tr, &tg, &tb); + img->data[idx] = (0xFF<<24)|(tr<<16)|(tg<<8)|tb; + } + } +} + +void fill_circle_glow(img_t *img, int32_t cx, int32_t cy, + int32_t radius, int primary_channel) { + /* Soft filled dot with radial alpha falloff for highlight */ + int32_t ext = radius + 4; + for (int32_t y = cy-ext; y <= cy+ext; y++) { + for (int32_t x = cx-ext; x <= cx+ext; x++) { + float dist = sqrtf((float)((x-cx)*(x-cx)+(y-cy)*(y-cy))); + if (dist > (float)ext) continue; + if (x < 0 || x >= (int32_t)img->width || + y < 0 || y >= (int32_t)img->height) continue; + float alpha = (1.0f - dist / (float)ext) * 0.6f; + uint8_t tr, tg, tb; + channel_color(img, x, y, primary_channel, &tr, &tg, &tb); + blend_pixel(img, x, y, (0xFF<<24)|(tr<<16)|(tg<<8)|tb, alpha); + } + } +} + +/* ── Rectangle helpers ───────────────────────────────────────────────── */ + +void draw_rect_fill(img_t *img, int32_t x, int32_t y, + int32_t w, int32_t h, uint32_t color) { + int32_t x0 = (x < 0) ? 0 : x; + int32_t y0 = (y < 0) ? 0 : y; + int32_t x1 = (x+w > (int32_t)img->width) ? (int32_t)img->width : x+w; + int32_t y1 = (y+h > (int32_t)img->height) ? (int32_t)img->height : y+h; + for (int32_t cy = y0; cy < y1; cy++) { + uint32_t *p = &img->data[cy * img->width + x0]; + for (int32_t cx = x0; cx < x1; cx++) *p++ = color; + } +} + +void draw_rect_blend(img_t *img, int32_t x, int32_t y, + int32_t w, int32_t h, uint32_t color) { + if (x >= (int32_t)img->width || y >= (int32_t)img->height || + x+w <= 0 || y+h <= 0) return; + int32_t x0 = (x < 0) ? 0 : x; + int32_t y0 = (y < 0) ? 0 : y; + int32_t x1 = (x+w > (int32_t)img->width) ? (int32_t)img->width : x+w; + int32_t y1 = (y+h > (int32_t)img->height) ? (int32_t)img->height : y+h; + uint32_t sa = (color >> 24) & 0xFF, ia = 256 - sa; + uint32_t sr = (color >> 16) & 0xFF, sg = (color >> 8) & 0xFF, sb = color & 0xFF; + for (int32_t row = y0; row < y1; row++) { + uint32_t *ptr = &img->data[row * img->width + x0]; + for (int32_t col = x0; col < x1; col++) { + uint32_t dst = *ptr; + uint32_t r = (sr*sa + ((dst>>16)&0xFF)*ia) >> 8; + uint32_t g = (sg*sa + ((dst>> 8)&0xFF)*ia) >> 8; + uint32_t b = (sb*sa + (dst&0xFF)*ia) >> 8; + *ptr++ = (0xFF<<24)|(r<<16)|(g<<8)|b; + } + } +} + +void draw_rect_stroke(img_t *img, int32_t x, int32_t y, + int32_t w, int32_t h, uint32_t color) { + draw_rect_fill(img, x, y, w, 1, color); + draw_rect_fill(img, x, y+h-1, w, 1, color); + draw_rect_fill(img, x, y, 1, h, color); + draw_rect_fill(img, x+w-1, y, 1, h, color); +} + +/* ── Rounded rectangle ───────────────────────────────────────────────── */ + +void draw_rounded_rect_blend(img_t *img, int32_t x, int32_t y, + int32_t w, int32_t h, int32_t r, uint32_t color) { + if (r <= 0) { draw_rect_blend(img, x, y, w, h, color); return; } + int32_t x0 = (x < 0) ? 0 : x; + int32_t y0 = (y < 0) ? 0 : y; + int32_t x1 = (x+w > (int32_t)img->width) ? (int32_t)img->width : x+w; + int32_t y1 = (y+h > (int32_t)img->height) ? (int32_t)img->height : y+h; + uint32_t base_alpha = (color >> 24) & 0xFF; + for (int32_t row = y0; row < y1; row++) { + for (int32_t col = x0; col < x1; col++) { + int in_corner = 0, dx = 0, dy = 0; + if (col < x+r && row < y+r) { dx=x+r-col; dy=y+r-row; in_corner=1; } + else if (col >= x+w-r && row < y+r) { dx=col-(x+w-r-1); dy=y+r-row; in_corner=1; } + else if (col < x+r && row >= y+h-r) { dx=x+r-col; dy=row-(y+h-r-1); in_corner=1; } + else if (col >= x+w-r && row >= y+h-r) { dx=col-(x+w-r-1); dy=row-(y+h-r-1); in_corner=1; } + uint32_t final_alpha = base_alpha; + if (in_corner) { + float dist = sqrtf((float)(dx*dx+dy*dy)); + float diff = dist - r; + if (diff > 1.0f) continue; + if (diff > -1.0f) { + float aa = 0.5f - diff * 0.5f; + if (aa < 0) aa = 0; + if (aa > 1) aa = 1; + final_alpha = (uint32_t)(base_alpha * aa); + } + } + blend_pixel_fast(img, col, row, color, final_alpha); + } + } +} + +void draw_rounded_rect_fill(img_t *img, int32_t x, int32_t y, + int32_t w, int32_t h, int32_t r, uint32_t color) { + if (r <= 0) { draw_rect_fill(img, x, y, w, h, color); return; } + draw_rounded_rect_blend(img, x, y, w, h, r, color); +} \ No newline at end of file diff --git a/src/main.c b/src/main.c index 18f2e2e..f450b82 100644 --- a/src/main.c +++ b/src/main.c @@ -1,16 +1,687 @@ -#include -#include -#include -#include -#include "transfer_function.h" - -int main(int argc, char **argv) -{ - double complex *s = NULL; - singularity_array_t zeros, poles; - - create_singularities(&zeros, 5); - create_singularities(&poles, 5); - - return 0; +/* + * SINGULARITY — main.c (orchestrator) + * + * Key shortcuts: + * Alt+Enter Fullscreen / windowed toggle + * Escape Quit + * G Toggle Select ↔ Move cursor mode + * Z (Move mode) Add zero at mouse position + * X (Move mode) Add pole at mouse position + * Delete Delete selected entity + * Ctrl+Z Undo + * Ctrl+Y Redo (also Ctrl+Shift+Z) + * Ctrl+R Randomise all + * Tab Cycle colour map (0-5) + * Up/Down Increase / decrease step count + * H Toggle HUD visibility + * C Cycle circle display mode + * S / F1 Open / close settings panel + */ + +#include +#include +#include +#include +#include +#include +#include + +#include "app.h" +#include "simulation.h" +#include "transfer_function.h" +#include "ui.h" +#include "utils.h" +#include "presets.h" +#include "audio.h" +#include "undo.h" +#include "ui.h" +#include "transfer_function.h" +#include "utils.h" +#include "presets.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#define TARGET_FPS 60 +#define FRAME_DURATION_MS (1000 / TARGET_FPS) + +/* ── Helpers ─────────────────────────────────────────────────────────── */ + +static void randomise_sim(Simulation *sim, App *app) { + sim_randomise(sim, app->x_range, app->y_range); +} + +/* Invalidate any pointer-to-entity after array changes */ +static void clear_selection(Popup *popup, singularity_t **move_target, int *selected_count) { + popup->is_open = 0; + popup->entity = NULL; + popup->dragging_slider = -1; + *move_target = NULL; + *selected_count = 0; +} + +static void do_resize(App *app, Simulation *sim, int w, int h) { + app_handle_resize(app, w, h); + sim_rebuild_grid(sim, app->x_range, app->y_range, app->grid_w, app->grid_h); +} + +static void save_screenshot(App *app) { + static int shot_counter = 1; + char filename[64]; + snprintf(filename, sizeof(filename), "screenshot_%04d.bmp", shot_counter++); + + SDL_Surface *sshot = SDL_CreateRGBSurfaceWithFormat(0, app->screen_w, app->screen_h, 32, SDL_PIXELFORMAT_ARGB8888); + if (sshot) { + SDL_RenderReadPixels(app->renderer, NULL, SDL_PIXELFORMAT_ARGB8888, sshot->pixels, sshot->pitch); + SDL_SaveBMP(sshot, filename); + SDL_FreeSurface(sshot); + printf("Saved %s\n", filename); + } +} + +/* ═══════════════════════════════════════════════════════════════════════ + main + ═══════════════════════════════════════════════════════════════════════ */ +int main(int argc, char **argv) { + (void)argc; (void)argv; + + App app; + Simulation sim; + UndoStack undo; + + /* Runtime state */ + UIFlags flags = { .hud_visible = 1, .circle_mode = CIRCLE_SOLID }; + RenderCfg rcfg = { .c_map = 0, .n_map = 1, .steps = 5 }; + Popup popup = { .is_open = 0, .dragging_slider = -1 }; + SettingsPanel settings = { .is_open = 0, .drag_steps = -1 }; + + singularity_t *move_target = NULL; + singularity_t *hover_entity = NULL; + int hover_type = 0; + int cursor_mode = 0; /* 0=select 1=move */ + + /* FPS */ + uint32_t frame_count = 0; + uint32_t last_fps_time = 0; + float current_fps = 0.0f; + + /* Selection Box & Group Move State */ + int is_selecting = 0; + int sel_sx = 0, sel_sy = 0; + int sel_ex = 0, sel_ey = 0; + + /* We track which entities are 'selected'. Note: size_t since indices can change on delete */ + #define MAX_SELECTED 128 + singularity_t* selected_group[MAX_SELECTED]; + int selected_count = 0; + + /* Animation state */ + double orbit_angle = 0.0; + double pulse_t = 0.0; + + /* Init SDL, window, renderer */ + srand((unsigned int)time(NULL)); + if (app_init(&app) != 0) return 1; + audio_init(); + + sim_init(&sim, app.x_range, app.y_range, app.grid_w, app.grid_h); + undo_init(&undo); + undo_push(&undo, &sim.zeros, &sim.poles); /* initial snapshot */ + + int running = 1; + SDL_Event e; + + while (running) { + Uint32 frame_start = SDL_GetTicks(); + int mx, my; + SDL_GetMouseState(&mx, &my); + + /* ── Event loop ──────────────────────────────────────────────── */ + while (SDL_PollEvent(&e)) { + + if (e.type == SDL_MOUSEMOTION) { + /* No need to update local mx,my here as we pull at start of frame */ + } else if (e.type == SDL_MOUSEBUTTONDOWN || + e.type == SDL_MOUSEBUTTONUP) { + mx = e.button.x; my = e.button.y; + } + + /* ── Window events ───────────────────────────────────────── */ + if (e.type == SDL_QUIT) { + running = 0; + + } else if (e.type == SDL_WINDOWEVENT && + e.window.event == SDL_WINDOWEVENT_RESIZED) { + do_resize(&app, &sim, e.window.data1, e.window.data2); + clear_selection(&popup, &move_target, &selected_count); + + /* ── Keyboard ────────────────────────────────────────────── */ + } else if (e.type == SDL_KEYDOWN) { + SDL_Keycode key = e.key.keysym.sym; + SDL_Keymod mods = SDL_GetModState(); + int ctrl = (mods & KMOD_CTRL) != 0; + int shift = (mods & KMOD_SHIFT) != 0; + int alt = (mods & KMOD_ALT) != 0; + + if (key == SDLK_ESCAPE) { + running = 0; + + /* Alt+Enter — fullscreen toggle */ + } else if (key == SDLK_RETURN && alt) { + app_toggle_fullscreen(&app); + int w, h; + SDL_GetWindowSize(app.window, &w, &h); + do_resize(&app, &sim, w, h); + clear_selection(&popup, &move_target, &selected_count); + + /* Ctrl+Z — undo */ + } else if (key == SDLK_z && ctrl && !shift) { + if (undo_undo(&undo, &sim.zeros, &sim.poles)) + clear_selection(&popup, &move_target, &selected_count); + + /* Ctrl+Y or Ctrl+Shift+Z — redo */ + } else if ((key == SDLK_y && ctrl) || + (key == SDLK_z && ctrl && shift)) { + if (undo_redo(&undo, &sim.zeros, &sim.poles)) + clear_selection(&popup, &move_target, &selected_count); + + /* Ctrl+R — randomise */ + } else if (key == SDLK_r && ctrl) { + undo_push(&undo, &sim.zeros, &sim.poles); + randomise_sim(&sim, &app); + clear_selection(&popup, &move_target, &selected_count); + + /* G — toggle cursor mode */ + } else if (key == SDLK_g) { + cursor_mode = !cursor_mode; + popup.is_open = 0; + + /* 1-9 — Presets */ + } else if (key >= SDLK_1 && key <= SDLK_9) { + undo_push(&undo, &sim.zeros, &sim.poles); + load_preset(&sim, key - SDLK_0); + app.cam_zoom = 1.0; + app.cam_cx = 0.0; + app.cam_cy = 0.0; + do_resize(&app, &sim, app.screen_w, app.screen_h); + clear_selection(&popup, &move_target, &selected_count); + + /* Z — add zero */ + } else if (key == SDLK_z && !ctrl) { + double complex cpos; + screen_choords_to_complex(mx, my, &cpos, + app.x_range, app.y_range, app.screen_w, app.screen_h); + undo_push(&undo, &sim.zeros, &sim.poles); + sim_add_zero(&sim, cpos); + + /* X — add pole */ + } else if (key == SDLK_x) { + double complex cpos; + screen_choords_to_complex(mx, my, &cpos, + app.x_range, app.y_range, app.screen_w, app.screen_h); + undo_push(&undo, &sim.zeros, &sim.poles); + sim_add_pole(&sim, cpos); + + /* Delete — remove selected entity */ + } else if (key == SDLK_DELETE && popup.entity) { + undo_push(&undo, &sim.zeros, &sim.poles); + sim_delete_entity(&sim, popup.entity); + clear_selection(&popup, &move_target, &selected_count); + + /* Tab — cycle colour map */ + } else if (key == SDLK_TAB) { + rcfg.c_map = (rcfg.c_map + 1) % 6; + + /* Up/Down — step count */ + } else if (key == SDLK_UP) { + if (++rcfg.steps > 64) rcfg.steps = 64; + } else if (key == SDLK_DOWN) { + if (--rcfg.steps < 1) rcfg.steps = 1; + + /* Backspace — cycle normalisation */ + } else if (key == SDLK_BACKSPACE) { + rcfg.n_map = (rcfg.n_map + 1) % 3; + + /* H — toggle HUD */ + } else if (key == SDLK_h) { + flags.hud_visible = !flags.hud_visible; + + /* C — cycle circle mode */ + } else if (key == SDLK_c && !ctrl) { + flags.circle_mode = (CircleMode)((flags.circle_mode + 1) % CIRCLE_MODE_COUNT); + + /* P — save screenshot */ + } else if (key == SDLK_p && !ctrl) { + /* We actually save after the render frame finishes below, + so we just set a flag here */ + flags.screenshot_queued = 1; + + /* O — toggle orbit */ + } else if (key == SDLK_o && !ctrl) { + flags.orbit_mode = !flags.orbit_mode; + + /* U — toggle pulse */ + } else if (key == SDLK_u && !ctrl) { + flags.pulse_mode = !flags.pulse_mode; + + /* F — toggle field lines */ + } else if (key == SDLK_f && !ctrl) { + flags.fieldlines_visible = !flags.fieldlines_visible; + + /* D — toggle domain overlay */ + } else if (key == SDLK_d && !ctrl) { + flags.domain_overlay = !flags.domain_overlay; + + /* A — toggle audio */ + } else if (key == SDLK_a && !ctrl) { + flags.audio_enabled = !flags.audio_enabled; + audio_set_enabled(flags.audio_enabled); + + /* S / F1 — settings panel */ + } else if (key == SDLK_s || key == SDLK_F1) { + settings.is_open = !settings.is_open; + } + + /* ── Mouse Wheel (Zoom) ──────────────────────────────────── */ + } else if (e.type == SDL_MOUSEWHEEL) { + if (e.wheel.y != 0) { + double zoom_factor = (e.wheel.y > 0) ? 1.1 : (1.0 / 1.1); + + /* Find complex coord under mouse before zoom */ + double complex cpos_before; + screen_choords_to_complex(mx, my, &cpos_before, + app.x_range, app.y_range, app.screen_w, app.screen_h); + + /* Apply zoom */ + app.cam_zoom *= zoom_factor; + if (app.cam_zoom < 0.01) app.cam_zoom = 0.01; + if (app.cam_zoom > 1000.0) app.cam_zoom = 1000.0; + + /* Recompute ranges so we can map screen back to complex */ + do_resize(&app, &sim, app.screen_w, app.screen_h); + + /* Find complex coord under mouse after zoom (if center didn't move) */ + double complex cpos_after; + screen_choords_to_complex(mx, my, &cpos_after, + app.x_range, app.y_range, app.screen_w, app.screen_h); + + /* Offset camera so cpos_before == cpos_after */ + app.cam_cx += creal(cpos_before) - creal(cpos_after); + app.cam_cy += cimag(cpos_before) - cimag(cpos_after); + + /* Final recompute of ranges with shifted camera */ + do_resize(&app, &sim, app.screen_w, app.screen_h); + } + + /* ── Mouse button down ──────────────────────────────────── */ + } else if (e.type == SDL_MOUSEBUTTONDOWN && + e.button.button == SDL_BUTTON_LEFT) { + + /* --- Settings panel clicks first -------------------- */ + if (settings.is_open) { + int cd = 0, nd = 0, req_fs = 0; + int hit = settings_hit(&settings, &rcfg, &flags, + mx, my, app.screen_w, app.screen_h, + app.is_fullscreen, &cd, &nd); + if (hit == SETTINGS_HIT_FS) { + app_toggle_fullscreen(&app); + int w, h; SDL_GetWindowSize(app.window, &w, &h); + do_resize(&app, &sim, w, h); + clear_selection(&popup, &move_target, &selected_count); + } + (void)req_fs; (void)cd; (void)nd; + goto next_event; /* consumed */ + } + + /* --- Select mode ------------------------------------ */ + if (cursor_mode == 0) { + /* Clicking inside open popup */ + if (popup.is_open) { + if (mx >= popup.x && mx < popup.x + 340 && + my >= popup.y && my < popup.y + 320) { + if (popup_hit_mode_btn(&popup, mx, my)) { + popup.edit_mode = !popup.edit_mode; + } else { + int idx; + if (popup_hit_slider(&popup, mx, my, &idx)) { + popup.dragging_slider = idx; + popup_apply_slider(&popup, mx); + } + } + goto next_event; + } else { + popup.is_open = 0; popup.entity = NULL; + } + } + + /* Click on entity → open popup */ + double complex cpos; + screen_choords_to_complex(mx, my, &cpos, + app.x_range, app.y_range, app.screen_w, app.screen_h); + int etype = 0; + singularity_t *hit_e = sim_find_entity(&sim, cpos, 0.5, &etype); + if (hit_e) { + popup.entity = hit_e; + popup.entity_type = etype; + popup.is_open = 1; + popup.edit_mode = 0; + popup.dragging_slider = -1; + int px_try = mx + 40; + if (px_try + 340 > app.screen_w) px_try = mx - 340 - 40; + if (px_try < 0) px_try = 5; + popup.x = px_try; + popup.y = my - 20; + if (popup.y + 320 > app.screen_h) popup.y = app.screen_h - 320; + if (popup.y < 0) popup.y = 0; + } else { + popup.is_open = 0; popup.entity = NULL; + } + + /* --- Move mode -------------------------------------- */ + } else { + double complex cpos; + screen_choords_to_complex(mx, my, &cpos, + app.x_range, app.y_range, app.screen_w, app.screen_h); + move_target = sim_find_entity(&sim, cpos, 0.5, NULL); + + /* If we didn't click on a selected member but we hit someone else, make them the only selection */ + if (move_target) { + int in_group = 0; + for(int i=0; i= 0) { + undo_push(&undo, &sim.zeros, &sim.poles); + popup.dragging_slider = -1; + } + settings.drag_steps = -1; + + } else if (e.type == SDL_MOUSEBUTTONUP && + e.button.button == SDL_BUTTON_RIGHT && cursor_mode == 1) { + if (is_selecting) { + is_selecting = 0; + selected_count = 0; + + /* Build AABB in complex space */ + double complex c1, c2; + screen_choords_to_complex(sel_sx, sel_sy, &c1, app.x_range, app.y_range, app.screen_w, app.screen_h); + screen_choords_to_complex(sel_ex, sel_ey, &c2, app.x_range, app.y_range, app.screen_w, app.screen_h); + + double min_r = fmin(creal(c1), creal(c2)); + double max_r = fmax(creal(c1), creal(c2)); + double min_i = fmin(cimag(c1), cimag(c2)); + double max_i = fmax(cimag(c1), cimag(c2)); + + /* Find intersecting entities */ + for (int k = 0; k < 2; k++) { + singularity_array_t *arr = (k == 0) ? &sim.zeros : &sim.poles; + for (size_t i = 0; i < arr->size; i++) { + double r = creal(arr->data[i].val); + double im = cimag(arr->data[i].val); + if (r >= min_r && r <= max_r && im >= min_i && im <= max_i) { + if (selected_count < MAX_SELECTED) { + selected_group[selected_count++] = &arr->data[i]; + } + } + } + } + } + + /* ── Mouse motion ───────────────────────────────────────── */ + } else if (e.type == SDL_MOUSEMOTION) { + /* Drag slider in popup */ + if (popup.dragging_slider >= 0 && popup.entity) + popup_apply_slider(&popup, mx); + + /* Drag entity/group in move mode */ + if (move_target && cursor_mode == 1) { + double complex cpos; + screen_choords_to_complex(e.motion.x, e.motion.y, &cpos, + app.x_range, app.y_range, app.screen_w, app.screen_h); + + double complex diff = cpos - move_target->val; + for (int i = 0; i < selected_count; i++) { + selected_group[i]->val += diff; + } + } + + /* Selection box update */ + if (is_selecting && cursor_mode == 1) { + sel_ex = mx; + sel_ey = my; + } + + /* Middle-click pan */ + if (SDL_GetMouseState(NULL, NULL) & SDL_BUTTON(SDL_BUTTON_MIDDLE)) { + /* Convert pixel delta to complex delta using current range sizes */ + double dx = app.x_range[1] - app.x_range[0]; + double dy = app.y_range[1] - app.y_range[0]; + app.cam_cx -= (double)e.motion.xrel * (dx / app.screen_w); + app.cam_cy -= (double)e.motion.yrel * (dy / app.screen_h); + /* Important: Recompute grids immediately */ + do_resize(&app, &sim, app.screen_w, app.screen_h); + } + + /* Drag steps slider inside settings */ + if (settings.drag_steps >= 0) { + int srx = (app.screen_w - 340)/2 + 20; + double t = (double)(mx - srx) / 300.0; + if (t < 0) t = 0; + if (t > 1) t = 1; + rcfg.steps = 1 + (int)(t * 63 + 0.5); + } + } + + next_event:; + } /* end event loop */ + + /* ── Cursor update ───────────────────────────────────────────── */ + { + int cmx, cmy; SDL_GetMouseState(&cmx, &cmy); + double complex cpos; + screen_choords_to_complex(cmx, cmy, &cpos, + app.x_range, app.y_range, app.screen_w, app.screen_h); + hover_entity = sim_find_entity(&sim, cpos, 0.5, &hover_type); + + if (hover_entity) SDL_SetCursor(app.cursor_hand); + else if (cursor_mode == 1) SDL_SetCursor(app.cursor_move); + else SDL_SetCursor(app.cursor_arrow); + } + + /* ── Animations ──────────────────────────────────────────────── */ + double dt = 1.0 / 60.0; /* Approximate dt for smooth animation */ + if (current_fps > 10.0) dt = 1.0 / current_fps; + + if (flags.orbit_mode) { + double rotation_speed = 0.5; /* radians per second */ + double rot = rotation_speed * dt; + orbit_angle += rot; + double complex rot_factor = cexp(I * rot); + for (size_t i = 0; i < sim.zeros.size; i++) + sim.zeros.data[i].val *= rot_factor; + for (size_t i = 0; i < sim.poles.size; i++) + sim.poles.data[i].val *= rot_factor; + } + + if (flags.pulse_mode) { + double pulse_speed = 3.0; + pulse_t += pulse_speed * dt; + double pulse_factor = 1.0 + 0.3 * sin(pulse_t); + for (size_t i = 0; i < sim.zeros.size; i++) + sim.zeros.data[i].m = 1.0 * pulse_factor; /* Assuming base m is 1.0 for now, would be better to store base_m */ + for (size_t i = 0; i < sim.poles.size; i++) + sim.poles.data[i].m = 1.0 * pulse_factor; + } else if (pulse_t > 0.0) { + /* Reset to 1.0 when turned off */ + pulse_t = 0.0; + for (size_t i = 0; i < sim.zeros.size; i++) sim.zeros.data[i].m = 1.0; + for (size_t i = 0; i < sim.poles.size; i++) sim.poles.data[i].m = 1.0; + } + + if (flags.audio_enabled) { + int mx, my; + SDL_GetMouseState(&mx, &my); + double complex cpos; + screen_choords_to_complex(mx, my, &cpos, + app.x_range, app.y_range, app.screen_w, app.screen_h); + + double complex H = sim_eval_H(&sim, cpos); + double mag = cabs(H); + + double freq = 440.0; + if (mag > 0.001) { + // Pitch based on log-magnitude (octaves) + freq = 440.0 * pow(2.0, log10(mag) / 2.0); + } + if (freq < 20.0) freq = 20.0; + if (freq > 2000.0) freq = 2000.0; + + double min_dist = 1e9; + for (size_t i = 0; i < sim.zeros.size; i++) { + double d = cabs(cpos - sim.zeros.data[i].val); + if (d < min_dist) min_dist = d; + } + for (size_t i = 0; i < sim.poles.size; i++) { + double d = cabs(cpos - sim.poles.data[i].val); + if (d < min_dist) min_dist = d; + } + + double amp = 1.0 / (1.0 + min_dist * 5.0); // Fade out quickly + audio_set_params(freq, amp); + } + + /* ── Math render (low-res) ───────────────────────────────────── */ + sim_compute(&sim, app.grid_w, app.grid_h); + + switch (rcfg.n_map) { + case 0: normalize_H_soa(sim.nH_re, sim.nH_im, + (size_t)app.grid_w * app.grid_h, sim.nH_re, sim.nH_im); break; + case 1: normalize_H_log_soa(sim.nH_re, sim.nH_im, + (size_t)app.grid_w * app.grid_h, sim.nH_re, sim.nH_im); break; + case 2: normalize_H_log_steps_soa(sim.nH_re, sim.nH_im, + (size_t)app.grid_w * app.grid_h, rcfg.steps, + sim.nH_re, sim.nH_im); break; + } + + switch (rcfg.c_map) { + case 0: H_g_img (sim.nH_re, sim.nH_im, (size_t)app.grid_w*app.grid_h, app.H_img); break; + case 1: H_c1_img(sim.nH_re, sim.nH_im, (size_t)app.grid_w*app.grid_h, app.H_img); break; + case 2: H_c2_img(sim.nH_re, sim.nH_im, (size_t)app.grid_w*app.grid_h, app.H_img); break; + case 3: H_c3_img(sim.nH_re, sim.nH_im, (size_t)app.grid_w*app.grid_h, app.H_img); break; + case 4: H_c4_img(sim.nH_re, sim.nH_im, (size_t)app.grid_w*app.grid_h, app.H_img); break; + case 5: H_c5_img(sim.nH_re, sim.nH_im, (size_t)app.grid_w*app.grid_h, app.H_img); break; + } + + /* ── UI render (full-res) ────────────────────────────────────── */ + memset(app.UI_img.data, 0, app.screen_w * app.screen_h * sizeof(uint32_t)); + + if (flags.hud_visible) + ui_draw_grid(&app.UI_img, app.x_range[0], app.x_range[1], + app.y_range[0], app.y_range[1]); + + if (flags.domain_overlay) { + ui_draw_domain_overlay(&app.UI_img, &sim, app.x_range, app.y_range); + } + + if (flags.fieldlines_visible) { + ui_draw_field_lines(&app.UI_img, &sim, app.x_range, app.y_range); + } + + ui_draw_entities(&app.UI_img, &sim, + app.x_range, app.y_range, + popup.entity, move_target, flags); + + /* Highlight selected group members */ + for (int i = 0; i < selected_count; i++) { + uint32_t cx, cy; + complex_to_screen_choords(selected_group[i]->val, &cx, &cy, app.x_range, app.y_range, app.screen_w, app.screen_h); + draw_circle_glow_color(&app.UI_img, cx, cy, 20, 0xAAFFFF00); + } + + /* Draw active selection box */ + if (is_selecting) { + int x_min = (sel_sx < sel_ex) ? sel_sx : sel_ex; + int x_max = (sel_sx > sel_ex) ? sel_sx : sel_ex; + int y_min = (sel_sy < sel_ey) ? sel_sy : sel_ey; + int y_max = (sel_sy > sel_ey) ? sel_sy : sel_ey; + for (int x = x_min; x <= x_max; x+=4) { + if (x >= 0 && x < app.screen_w && y_min >= 0 && y_min < app.screen_h) app.UI_img.data[y_min * app.screen_w + x] = 0xFFFFFFFF; + if (x >= 0 && x < app.screen_w && y_max >= 0 && y_max < app.screen_h) app.UI_img.data[y_max * app.screen_w + x] = 0xFFFFFFFF; + } + for (int y = y_min; y <= y_max; y+=4) { + if (x_min >= 0 && x_min < app.screen_w && y >= 0 && y < app.screen_h) app.UI_img.data[y * app.screen_w + x_min] = 0xFFFFFFFF; + if (x_max >= 0 && x_max < app.screen_w && y >= 0 && y < app.screen_h) app.UI_img.data[y * app.screen_w + x_max] = 0xFFFFFFFF; + } + } + + if (hover_entity && !popup.is_open) + ui_draw_tooltip(&app.UI_img, mx, my, hover_entity, hover_type); + + ui_draw_hud(&app.UI_img, current_fps, &rcfg, &flags, + cursor_mode, app.is_fullscreen); + ui_draw_popup(&app.UI_img, &popup); + + int req_fs = 0; + ui_draw_settings(&app.UI_img, &settings, &rcfg, &flags, + app.screen_w, app.screen_h, app.is_fullscreen, &req_fs); + + /* ── Composite & present ─────────────────────────────────────── */ + SDL_UpdateTexture(app.H_tex, NULL, app.H_img.data, app.H_img.width * sizeof(uint32_t)); + SDL_UpdateTexture(app.UI_tex, NULL, app.UI_img.data, app.UI_img.width * sizeof(uint32_t)); + SDL_RenderClear(app.renderer); + SDL_RenderCopy(app.renderer, app.H_tex, NULL, NULL); + SDL_RenderCopy(app.renderer, app.UI_tex, NULL, NULL); + + if (flags.screenshot_queued) { + save_screenshot(&app); + flags.screenshot_queued = 0; + } + + SDL_RenderPresent(app.renderer); + + /* ── Frame timing & FPS ──────────────────────────────────────── */ + Uint32 elapsed = SDL_GetTicks() - frame_start; + frame_count++; + uint32_t now = SDL_GetTicks(); + if (now - last_fps_time >= 500) { + current_fps = frame_count * 1000.0f / (float)(now - last_fps_time); + frame_count = 0; + last_fps_time = now; + } + if (elapsed < FRAME_DURATION_MS) + SDL_Delay(FRAME_DURATION_MS - elapsed); + } + + /* ── Cleanup ─────────────────────────────────────────────────────── */ + undo_destroy(&undo); + sim_destroy(&sim); + audio_shutdown(); + app_destroy(&app); + printf("Bye.\n"); + return 0; } \ No newline at end of file diff --git a/src/presets.c b/src/presets.c new file mode 100644 index 0000000..3d0e74a --- /dev/null +++ b/src/presets.c @@ -0,0 +1,114 @@ +#include "presets.h" +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +void load_preset(Simulation *sim, int n) { + if (!sim) return; + + /* Clear current entities */ + sim->zeros.size = 0; + sim->poles.size = 0; + + switch (n) { + case 1: { + /* Preset 1: Butterworth Pattern (Poles on left unit semicircle) */ + int order = 5; + for (int k = 1; k <= order; k++) { + double angle = M_PI * (2 * k + order - 1) / (2 * order); + sim_add_pole(sim, cos(angle) + sin(angle) * I); + } + break; + } + case 2: { + /* Preset 2: Simple Resonator (Conjugate pole pair + zero at origin) */ + sim_add_zero(sim, 0.0 + 0.0 * I); + sim_add_pole(sim, -0.2 + 1.0 * I); + sim_add_pole(sim, -0.2 - 1.0 * I); + break; + } + case 3: { + /* Preset 3: Symmetric Cross (4 Zeros, 4 Poles) */ + sim_add_zero(sim, 1.0 + 0.0 * I); + sim_add_zero(sim, -1.0 + 0.0 * I); + sim_add_zero(sim, 0.0 + 1.0 * I); + sim_add_zero(sim, 0.0 - 1.0 * I); + + sim_add_pole(sim, 1.5 + 1.5 * I); + sim_add_pole(sim, -1.5 + 1.5 * I); + sim_add_pole(sim, 1.5 - 1.5 * I); + sim_add_pole(sim, -1.5 - 1.5 * I); + break; + } + case 4: { + /* Preset 4: Ring of Poles */ + int count = 8; + for (int k = 0; k < count; k++) { + double angle = 2.0 * M_PI * k / count; + sim_add_pole(sim, cos(angle) * 1.5 + sin(angle) * 1.5 * I); + } + break; + } + case 5: { + /* Preset 5: Alternating Ring (Zero/Pole) */ + int count = 6; + for (int k = 0; k < count * 2; k++) { + double angle = M_PI * k / count; + double complex pos = cos(angle) * 2.0 + sin(angle) * 2.0 * I; + if (k % 2 == 0) sim_add_zero(sim, pos); + else sim_add_pole(sim, pos); + } + break; + } + case 6: { + /* Preset 6: Elliptic Filter Pattern (Poles on ellipse, zeros on imaginary axis) */ + sim_add_pole(sim, -0.1 + 0.9 * I); + sim_add_pole(sim, -0.1 - 0.9 * I); + sim_add_pole(sim, -0.4 + 0.4 * I); + sim_add_pole(sim, -0.4 - 0.4 * I); + sim_add_pole(sim, -0.6 + 0.0 * I); + + sim_add_zero(sim, 0.0 + 1.5 * I); + sim_add_zero(sim, 0.0 - 1.5 * I); + sim_add_zero(sim, 0.0 + 2.5 * I); + sim_add_zero(sim, 0.0 - 2.5 * I); + break; + } + case 7: { + /* Preset 7: Multi-order pole at origin, surrounded by zeros */ + singularity_t *p = sim_add_pole(sim, 0.0 + 0.0 * I); + p->e = 4; /* Order 4 */ + + sim_add_zero(sim, 2.0 + 0.0 * I); + sim_add_zero(sim, -2.0 + 0.0 * I); + sim_add_zero(sim, 0.0 + 2.0 * I); + sim_add_zero(sim, 0.0 - 2.0 * I); + break; + } + case 8: { + /* Preset 8: Checkerboard (Grid of alternating entities) */ + for (int y = -1; y <= 1; y++) { + for (int x = -1; x <= 1; x++) { + if (x == 0 && y == 0) continue; + double complex pos = (double)x * 1.5 + (double)y * 1.5 * I; + if ((x + y) % 2 == 0) sim_add_pole(sim, pos); + else sim_add_zero(sim, pos); + } + } + break; + } + case 9: { + /* Preset 9: Unit Circle Poles + Center Zero (All-Pass-like structure) */ + sim_add_zero(sim, 0.0 + 0.0 * I); + for (int k = 0; k < 12; k++) { + double angle = 2.0 * M_PI * k / 12.0; + sim_add_pole(sim, cos(angle) + sin(angle) * I); + } + break; + } + default: + break; + } +} diff --git a/src/simd_avx.c b/src/simd_avx.c new file mode 100644 index 0000000..99d1072 --- /dev/null +++ b/src/simd_avx.c @@ -0,0 +1,241 @@ +#include "simd_avx.h" +#include + +// --- AVX family (256-bit) --- +static void compute_s8_from_index( + uint16_t height, + uint16_t width, + float y_range[2], + float x_range[2], + uint32_t start_index, // starting linear index + float out_re[8], // output real parts + float out_im[8] // output imaginary parts +) +{ + float dy = (y_range[1] - y_range[0]) / (height - 1); + float dx = (x_range[1] - x_range[0]) / (width - 1); + + for (int i = 0; i < 8; i++) + { + uint32_t idx = start_index + i; + if (idx >= height * width) + { + out_re[i] = 0.0f; + out_im[i] = 0.0f; + continue; + } + + uint16_t y = idx / width; + uint16_t x = idx % width; + + out_re[i] = (float)(x_range[0] + x * dx); + out_im[i] = (float)(y_range[0] + y * dy); + } +} + +// Multiply two complex SIMD numbers: (ar + i ai) * (br + i bi) +__attribute__((target("avx"))) static void avx_complex_mul( + __m256 ar, __m256 ai, + __m256 br, __m256 bi, + __m256 *out_r, __m256 *out_i) +{ + __m256 tmp_r = _mm256_sub_ps(_mm256_mul_ps(ar, br), _mm256_mul_ps(ai, bi)); + __m256 tmp_i = _mm256_add_ps(_mm256_mul_ps(ar, bi), _mm256_mul_ps(ai, br)); + *out_r = tmp_r; + *out_i = tmp_i; +} + +// Divide two complex SIMD numbers: (ar + i ai) / (br + i bi) +__attribute__((target("avx"))) static void avx_complex_div( + __m256 ar, __m256 ai, + __m256 br, __m256 bi, + __m256 *out_r, __m256 *out_i) +{ + // denominator = br^2 + bi^2 + __m256 br2 = _mm256_mul_ps(br, br); + __m256 bi2 = _mm256_mul_ps(bi, bi); + __m256 denom = _mm256_add_ps(br2, bi2); + + // real numerator: (ar*br + ai*bi) + __m256 num_r = _mm256_add_ps(_mm256_mul_ps(ar, br), _mm256_mul_ps(ai, bi)); + + // imag numerator: (ai*br - ar*bi) + __m256 num_i = _mm256_sub_ps(_mm256_mul_ps(ai, br), _mm256_mul_ps(ar, bi)); + + // divide both by denominator + *out_r = _mm256_div_ps(num_r, denom); + *out_i = _mm256_div_ps(num_i, denom); +} + +// Store 1-4 lanes from v into dst +__attribute__((target("avx"))) static void avx_store_n_floats(float *dst, __m256 v, int n) +{ + if (n <= 0 || n > 7) + return; // invalid + + float tmp[8]; + _mm256_storeu_ps(tmp, v); // store all 8 floats + for (int i = 0; i < n; i++) + { + dst[i] = tmp[i]; + } +} + +__attribute__((target("avx"))) float **avx_compute_H( + float x_range[2], float y_range[2], + new_singularity_array_t zeros_arr, new_singularity_array_t poles_arr, + uint16_t height, uint16_t width, + float **H) +{ + uint32_t length = height * width; + + // setting up H + if (H == NULL) + { + H = (float **)malloc(sizeof(float *) * 2); + if (!H) + { + return NULL; + } + H[0] = (float *)malloc(sizeof(float) * length); + H[1] = (float *)malloc(sizeof(float) * length); + if (!H[0] || !H[1]) + { + free(H[0]); + free(H[1]); + free(H); + return NULL; + } + } + + float sr8[8], si8[8]; + + for (uint32_t h_i = 0, h_di = 8; h_i < length;) + { + // load s (4 regs, total 4) + __m256 num_r = _mm256_set1_ps(1.0f); + __m256 num_i = _mm256_set1_ps(1.0f); + __m256 den_r = _mm256_set1_ps(1.0f); + __m256 den_i = _mm256_set1_ps(1.0f); + compute_s8_from_index(height, width, y_range, x_range, h_i, sr8, si8); + + // load s (2 regs, total 6) + __m256 sr8_v = _mm256_loadu_ps(sr8); // load 4 real parts + __m256 si8_v = _mm256_loadu_ps(si8); // load 4 imaginary parts + + for (size_t z_i = 0; z_i < zeros_arr.count; z_i += 1) + { + // terms (2 regs, total 8) + __m256 term_r, term_i; + + // load 1 zeros (2 regs, total 10) + __m256 zr_v = _mm256_set1_ps(zeros_arr.r[z_i]); + __m256 zi_v = _mm256_set1_ps(zeros_arr.i[z_i]); + + term_r = _mm256_sub_ps(sr8_v, zr_v); // term_r now holds s - z_r + term_i = _mm256_sub_ps(si8_v, zi_v); // term_i now holds s - z_i + + // --- do exponentiation per lane (binary loop) --- + __m256 res_r = _mm256_set1_ps(1.0f); + __m256 res_i = _mm256_set1_ps(0.0f); + __m256 base_r = term_r; + __m256 base_i = term_i; + int exp = zeros_arr.e[z_i]; + while (exp > 0) // O(logn) + { + if (exp & 1) + { + avx_complex_mul(res_r, res_i, base_r, base_i, &res_r, &res_i); + } + avx_complex_mul(base_r, base_i, base_r, base_i, &base_r, &base_i); + exp >>= 1; + } + term_r = res_r; + term_i = res_i; + + // (z.m * prod + z.c), with m,c complex + __m256 m_r = _mm256_set1_ps(zeros_arr.m_r[z_i]); + __m256 m_i = _mm256_set1_ps(zeros_arr.m_i[z_i]); + __m256 c_r = _mm256_set1_ps(zeros_arr.c_r[z_i]); + __m256 c_i = _mm256_set1_ps(zeros_arr.c_i[z_i]); + + __m256 tm_r, tm_i; + avx_complex_mul(m_r, m_i, term_r, term_i, &tm_r, &tm_i); + + term_r = _mm256_add_ps(tm_r, c_r); + term_i = _mm256_add_ps(tm_i, c_i); + + // num *= (that), applied per lane (2 regs, total 16) + avx_complex_mul(num_r, num_i, term_r, term_i, &num_r, &num_i); + } + // local regs released (-10 regs, total 6) + + for (size_t p_i = 0; p_i < poles_arr.count; p_i += 1) + { + // terms (2 regs, total 8) + __m256 term_r, term_i; + + // load 1 poles (2 regs, total 10) + __m256 pr_v = _mm256_set1_ps(poles_arr.r[p_i]); + __m256 pi_v = _mm256_set1_ps(poles_arr.i[p_i]); + + term_r = _mm256_sub_ps(sr8_v, pr_v); // term_r now holds s - p_r + term_i = _mm256_sub_ps(si8_v, pi_v); // term_i now holds s - p_i + + // --- do exponentiation per lane (binary loop) --- + __m256 res_r = _mm256_set1_ps(1.0f); + __m256 res_i = _mm256_set1_ps(0.0f); + __m256 base_r = term_r; + __m256 base_i = term_i; + int exp = poles_arr.e[p_i]; + while (exp > 0) // O(logn) + { + if (exp & 1) + { + avx_complex_mul(res_r, res_i, base_r, base_i, &res_r, &res_i); + } + avx_complex_mul(base_r, base_i, base_r, base_i, &base_r, &base_i); + exp >>= 1; + } + term_r = res_r; + term_i = res_i; + + // (p.m * prod + p.c), with m,c complex + __m256 m_r = _mm256_set1_ps(poles_arr.m_r[p_i]); + __m256 m_i = _mm256_set1_ps(poles_arr.m_i[p_i]); + __m256 c_r = _mm256_set1_ps(poles_arr.c_r[p_i]); + __m256 c_i = _mm256_set1_ps(poles_arr.c_i[p_i]); + + __m256 tm_r, tm_i; + avx_complex_mul(m_r, m_i, term_r, term_i, &tm_r, &tm_i); + + term_r = _mm256_add_ps(tm_r, c_r); + term_i = _mm256_add_ps(tm_i, c_i); + + // den *= (that), applied per lane (2 regs, total 16) + avx_complex_mul(den_r, den_i, term_r, term_i, &den_r, &den_i); + } + // local regs released (-10 regs, total 6) + + __m256 Hr, Hi; + avx_complex_div(num_r, num_i, den_r, den_i, &Hr, &Hi); + + if (h_i + 8 > length) + { + avx_store_n_floats(H[0] + h_i, Hr, length - h_i); + avx_store_n_floats(H[1] + h_i, Hi, length - h_i); + } + else + { + _mm256_storeu_ps(H[0] + h_i, Hr); + _mm256_storeu_ps(H[1] + h_i, Hi); + } + + if (h_di + h_i > length) + h_di = length - h_i; // adjust for last few elements + + h_i += h_di; + } + + return H; +} \ No newline at end of file diff --git a/src/simd_avx512.c b/src/simd_avx512.c new file mode 100644 index 0000000..881ebdb --- /dev/null +++ b/src/simd_avx512.c @@ -0,0 +1,221 @@ +#include "simd_avx512.h" +#include + +static void compute_s16_from_index( + uint16_t height, + uint16_t width, + float y_range[2], + float x_range[2], + uint32_t start_index, + float out_re[16], + float out_im[16] +) +{ + float dy = (y_range[1] - y_range[0]) / (height - 1); + float dx = (x_range[1] - x_range[0]) / (width - 1); + + for (int i = 0; i < 16; i++) + { + uint32_t idx = start_index + i; + if (idx >= height * width) + { + out_re[i] = 0.0f; + out_im[i] = 0.0f; + continue; + } + + uint16_t y = idx / width; + uint16_t x = idx % width; + + out_re[i] = (float)(x_range[0] + x * dx); + out_im[i] = (float)(y_range[0] + y * dy); + } +} + +__attribute__((target("avx512f"))) static inline void avx512_complex_mul( + __m512 ar, __m512 ai, + __m512 br, __m512 bi, + __m512 *out_r, __m512 *out_i) +{ + *out_r = _mm512_sub_ps(_mm512_mul_ps(ar, br), _mm512_mul_ps(ai, bi)); + *out_i = _mm512_add_ps(_mm512_mul_ps(ar, bi), _mm512_mul_ps(ai, br)); +} + +__attribute__((target("avx512f"))) static inline void avx512_complex_div( + __m512 ar, __m512 ai, + __m512 br, __m512 bi, + __m512 *out_r, __m512 *out_i) +{ + __m512 br2 = _mm512_mul_ps(br, br); + __m512 bi2 = _mm512_mul_ps(bi, bi); + __m512 denom = _mm512_add_ps(br2, bi2); + + __m512 num_r = _mm512_add_ps(_mm512_mul_ps(ar, br), _mm512_mul_ps(ai, bi)); + __m512 num_i = _mm512_sub_ps(_mm512_mul_ps(ai, br), _mm512_mul_ps(ar, bi)); + + *out_r = _mm512_div_ps(num_r, denom); + *out_i = _mm512_div_ps(num_i, denom); +} + +__attribute__((target("avx512f"))) static void avx512_store_n_floats(float *dst, __m512 v, int n) +{ + if (n <= 0 || n > 15) + return; + + float tmp[16]; + _mm512_storeu_ps(tmp, v); + for (int i = 0; i < n; i++) + { + dst[i] = tmp[i]; + } +} + +__attribute__((target("avx512f"))) float **avx512_compute_H( + float x_range[2], float y_range[2], + new_singularity_array_t zeros_arr, new_singularity_array_t poles_arr, + uint16_t height, uint16_t width, + float **H) +{ + uint32_t length = height * width; + + if (H == NULL) + { + H = (float **)malloc(sizeof(float *) * 2); + if (!H) + return NULL; + + H[0] = (float *)malloc(sizeof(float) * length); + H[1] = (float *)malloc(sizeof(float) * length); + if (!H[0] || !H[1]) + { + free(H[0]); + free(H[1]); + free(H); + return NULL; + } + } + + float sr16[16], si16[16]; + + for (uint32_t h_i = 0; h_i < length;) + { + uint32_t h_di = 16; + + __m512 num_r = _mm512_set1_ps(1.0f); + __m512 num_i = _mm512_set1_ps(0.0f); + __m512 den_r = _mm512_set1_ps(1.0f); + __m512 den_i = _mm512_set1_ps(0.0f); + + compute_s16_from_index(height, width, y_range, x_range, h_i, sr16, si16); + + __m512 sr16_v = _mm512_loadu_ps(sr16); + __m512 si16_v = _mm512_loadu_ps(si16); + + for (size_t z_i = 0; z_i < zeros_arr.count; z_i++) + { + __m512 zr_v = _mm512_set1_ps(zeros_arr.r[z_i]); + __m512 zi_v = _mm512_set1_ps(zeros_arr.i[z_i]); + + __m512 term_r = _mm512_sub_ps(sr16_v, zr_v); + __m512 term_i = _mm512_sub_ps(si16_v, zi_v); + + __m512 res_r = _mm512_set1_ps(1.0f); + __m512 res_i = _mm512_set1_ps(0.0f); + __m512 base_r = term_r; + __m512 base_i = term_i; + + int exp = zeros_arr.e[z_i]; + while (exp > 0) + { + if (exp & 1) + { + avx512_complex_mul(res_r, res_i, base_r, base_i, &res_r, &res_i); + } + avx512_complex_mul(base_r, base_i, base_r, base_i, &base_r, &base_i); + exp >>= 1; + } + + term_r = res_r; + term_i = res_i; + + __m512 m_r_v = _mm512_set1_ps(zeros_arr.m_r[z_i]); + __m512 m_i_v = _mm512_set1_ps(zeros_arr.m_i[z_i]); + __m512 c_r_v = _mm512_set1_ps(zeros_arr.c_r[z_i]); + __m512 c_i_v = _mm512_set1_ps(zeros_arr.c_i[z_i]); + + // Complex multiply term * m + __m512 tm_r, tm_i; + avx512_complex_mul(m_r_v, m_i_v, term_r, term_i, &tm_r, &tm_i); + + // Add c + term_r = _mm512_add_ps(tm_r, c_r_v); + term_i = _mm512_add_ps(tm_i, c_i_v); + + avx512_complex_mul(num_r, num_i, term_r, term_i, &num_r, &num_i); + } + + for (size_t p_i = 0; p_i < poles_arr.count; p_i++) + { + __m512 pr_v = _mm512_set1_ps(poles_arr.r[p_i]); + __m512 pi_v = _mm512_set1_ps(poles_arr.i[p_i]); + + __m512 term_r = _mm512_sub_ps(sr16_v, pr_v); + __m512 term_i = _mm512_sub_ps(si16_v, pi_v); + + __m512 res_r = _mm512_set1_ps(1.0f); + __m512 res_i = _mm512_set1_ps(0.0f); + __m512 base_r = term_r; + __m512 base_i = term_i; + + int exp = poles_arr.e[p_i]; + while (exp > 0) + { + if (exp & 1) + { + avx512_complex_mul(res_r, res_i, base_r, base_i, &res_r, &res_i); + } + avx512_complex_mul(base_r, base_i, base_r, base_i, &base_r, &base_i); + exp >>= 1; + } + + term_r = res_r; + term_i = res_i; + + __m512 m_r_v = _mm512_set1_ps(poles_arr.m_r[p_i]); + __m512 m_i_v = _mm512_set1_ps(poles_arr.m_i[p_i]); + __m512 c_r_v = _mm512_set1_ps(poles_arr.c_r[p_i]); + __m512 c_i_v = _mm512_set1_ps(poles_arr.c_i[p_i]); + + // Complex multiply term * m + __m512 tm_r, tm_i; + avx512_complex_mul(m_r_v, m_i_v, term_r, term_i, &tm_r, &tm_i); + + // Add c + term_r = _mm512_add_ps(tm_r, c_r_v); + term_i = _mm512_add_ps(tm_i, c_i_v); + + avx512_complex_mul(den_r, den_i, term_r, term_i, &den_r, &den_i); + } + + __m512 Hr, Hi; + avx512_complex_div(num_r, num_i, den_r, den_i, &Hr, &Hi); + + if (h_i + 16 > length) + { + avx512_store_n_floats(H[0] + h_i, Hr, length - h_i); + avx512_store_n_floats(H[1] + h_i, Hi, length - h_i); + } + else + { + _mm512_storeu_ps(H[0] + h_i, Hr); + _mm512_storeu_ps(H[1] + h_i, Hi); + } + + if (h_di + h_i > length) + h_di = length - h_i; + + h_i += h_di; + } + + return H; +} diff --git a/src/simd_sse.c b/src/simd_sse.c new file mode 100644 index 0000000..cd0aa6a --- /dev/null +++ b/src/simd_sse.c @@ -0,0 +1,344 @@ +#include "simd_sse.h" +#include + +// --- SSE family (128-bit) --- +static void compute_s4_from_index( + uint16_t height, + uint16_t width, + float y_range[2], + float x_range[2], + uint32_t start_index, // starting linear index + float out_re[4], // output real parts + float out_im[4] // output imaginary parts +) +{ + float dy = (y_range[1] - y_range[0]) / (height - 1); + float dx = (x_range[1] - x_range[0]) / (width - 1); + + for (int i = 0; i < 4; i++) + { + uint32_t idx = start_index + i; + if (idx >= height * width) + { + out_re[i] = 0.0f; + out_im[i] = 0.0f; + continue; + } + + uint16_t y = idx / width; + uint16_t x = idx % width; + + out_re[i] = (float)(x_range[0] + x * dx); + out_im[i] = (float)(y_range[0] + y * dy); + } +} + +// Multiply two complex SIMD numbers: (ar + i ai) * (br + i bi) +__attribute__((target("sse"))) static void sse_complex_mul( + __m128 ar, __m128 ai, + __m128 br, __m128 bi, + __m128 *out_r, __m128 *out_i) +{ + __m128 tmp_r = _mm_sub_ps(_mm_mul_ps(ar, br), _mm_mul_ps(ai, bi)); + __m128 tmp_i = _mm_add_ps(_mm_mul_ps(ar, bi), _mm_mul_ps(ai, br)); + *out_r = tmp_r; + *out_i = tmp_i; +} + +// Divide two complex SIMD numbers: (ar + i ai) / (br + i bi) +__attribute__((target("sse"))) static void sse_complex_div( + __m128 ar, __m128 ai, + __m128 br, __m128 bi, + __m128 *out_r, __m128 *out_i) +{ + // denominator = br^2 + bi^2 + __m128 br2 = _mm_mul_ps(br, br); + __m128 bi2 = _mm_mul_ps(bi, bi); + __m128 denom = _mm_add_ps(br2, bi2); + + // real numerator: (ar*br + ai*bi) + __m128 num_r = _mm_add_ps(_mm_mul_ps(ar, br), _mm_mul_ps(ai, bi)); + + // imag numerator: (ai*br - ar*bi) + __m128 num_i = _mm_sub_ps(_mm_mul_ps(ai, br), _mm_mul_ps(ar, bi)); + + // divide both by denominator + *out_r = _mm_div_ps(num_r, denom); + *out_i = _mm_div_ps(num_i, denom); +} + +// Store 1-4 lanes from v into dst +__attribute__((target("sse"))) static void sse_store_n_floats(float *dst, __m128 v, int n) +{ + if (n <= 0 || n > 3) + return; // invalid + + // lane 0 always exists + dst[0] = _mm_cvtss_f32(v); + + if (n >= 2) + dst[1] = _mm_cvtss_f32(_mm_shuffle_ps(v, v, _MM_SHUFFLE(1, 1, 1, 1))); + + if (n >= 3) + dst[2] = _mm_cvtss_f32(_mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 2, 2, 2))); +} + +__attribute__((target("sse2"))) float **sse2_compute_H(float x_range[2], float y_range[2], new_singularity_array_t zeros_arr, new_singularity_array_t poles_arr, uint16_t height, uint16_t width, float **H) +{ + uint32_t length = height * width; + + // setting up H + if (H == NULL) + { + H = (float **)malloc(sizeof(float *) * 2); + if (!H) + { + return NULL; + } + H[0] = (float *)malloc(sizeof(float) * length); + H[1] = (float *)malloc(sizeof(float) * length); + if (!H[0] || !H[1]) + { + free(H[0]); + free(H[1]); + free(H); + return NULL; + } + } + + float sr4[4], si4[4]; + + for (uint32_t h_i = 0, h_di = 4; h_i < length;) + { + // load s (4 regs, total 4) + __m128 num_r = _mm_set1_ps(1.0f); + __m128 num_i = _mm_set1_ps(1.0f); + __m128 den_r = _mm_set1_ps(1.0f); + __m128 den_i = _mm_set1_ps(1.0f); + compute_s4_from_index(height, width, y_range, x_range, h_i, sr4, si4); + + // load s (2 regs, total 6) + __m128 sr4_v = _mm_loadu_ps(sr4); // load 4 real parts + __m128 si4_v = _mm_loadu_ps(si4); // load 4 imaginary parts + + for (size_t z_i = 0; z_i < zeros_arr.count; z_i += 1) + { + // terms (2 regs, total 8) + __m128 term_r, term_i; + + // load 1 zeros (2 regs, total 10) + __m128 zr_v = _mm_set1_ps(zeros_arr.r[z_i]); + __m128 zi_v = _mm_set1_ps(zeros_arr.i[z_i]); + + term_r = _mm_sub_ps(sr4_v, zr_v); // term_r now holds s - z_r + term_i = _mm_sub_ps(si4_v, zi_v); // term_i now holds s - z_i + + // --- do exponentiation per lane (binary loop) --- + __m128 res_r = _mm_set1_ps(1.0f); + __m128 res_i = _mm_set1_ps(0.0f); + __m128 base_r = term_r; + __m128 base_i = term_i; + int exp = zeros_arr.e[z_i]; + while (exp > 0) // O(logn) + { + if (exp & 1) + { + sse_complex_mul(res_r, res_i, base_r, base_i, &res_r, &res_i); + } + sse_complex_mul(base_r, base_i, base_r, base_i, &base_r, &base_i); + exp >>= 1; + } + term_r = res_r; + term_i = res_i; + + // (z.m * prod + z.c), with m,c complex + __m128 m_r = _mm_set1_ps(zeros_arr.m_r[z_i]); + __m128 m_i = _mm_set1_ps(zeros_arr.m_i[z_i]); + __m128 c_r = _mm_set1_ps(zeros_arr.c_r[z_i]); + __m128 c_i = _mm_set1_ps(zeros_arr.c_i[z_i]); + + __m128 tm_r, tm_i; + sse_complex_mul(m_r, m_i, term_r, term_i, &tm_r, &tm_i); + + term_r = _mm_add_ps(tm_r, c_r); + term_i = _mm_add_ps(tm_i, c_i); + + // num *= (that), applied per lane (2 regs, total 16) + sse_complex_mul(num_r, num_i, term_r, term_i, &num_r, &num_i); + } + // local regs released (-10 regs, total 6) + + for (size_t p_i = 0; p_i < poles_arr.count; p_i += 1) + { + // terms (2 regs, total 8) + __m128 term_r, term_i; + + // load 1 poles (2 regs, total 10) + __m128 pr_v = _mm_set1_ps(poles_arr.r[p_i]); + __m128 pi_v = _mm_set1_ps(poles_arr.i[p_i]); + + term_r = _mm_sub_ps(sr4_v, pr_v); // term_r now holds s - p_r + term_i = _mm_sub_ps(si4_v, pi_v); // term_i now holds s - p_i + + // --- do exponentiation per lane (binary loop) --- + __m128 res_r = _mm_set1_ps(1.0f); + __m128 res_i = _mm_set1_ps(0.0f); + __m128 base_r = term_r; + __m128 base_i = term_i; + int exp = poles_arr.e[p_i]; + while (exp > 0) // O(logn) + { + if (exp & 1) + { + sse_complex_mul(res_r, res_i, base_r, base_i, &res_r, &res_i); + } + sse_complex_mul(base_r, base_i, base_r, base_i, &base_r, &base_i); + exp >>= 1; + } + term_r = res_r; + term_i = res_i; + + // (p.m * prod + p.c), with m,c complex + __m128 m_r = _mm_set1_ps(poles_arr.m_r[p_i]); + __m128 m_i = _mm_set1_ps(poles_arr.m_i[p_i]); + __m128 c_r = _mm_set1_ps(poles_arr.c_r[p_i]); + __m128 c_i = _mm_set1_ps(poles_arr.c_i[p_i]); + + __m128 tm_r, tm_i; + sse_complex_mul(m_r, m_i, term_r, term_i, &tm_r, &tm_i); + + term_r = _mm_add_ps(tm_r, c_r); + term_i = _mm_add_ps(tm_i, c_i); + + // den *= (that), applied per lane (2 regs, total 16) + sse_complex_mul(den_r, den_i, term_r, term_i, &den_r, &den_i); + } + // local regs released (-10 regs, total 6) + + __m128 Hr, Hi; + sse_complex_div(num_r, num_i, den_r, den_i, &Hr, &Hi); + + if (h_i + 4 > length) + { + sse_store_n_floats(H[0] + h_i, Hr, length - h_i); + sse_store_n_floats(H[1] + h_i, Hi, length - h_i); + } + else + { + _mm_storeu_ps(H[0] + h_i, Hr); + _mm_storeu_ps(H[1] + h_i, Hi); + } + + if (h_di + h_i > length) + h_di = length - h_i; // adjust for last few elements + + h_i += h_di; + } + + return H; +} + +__attribute__((target("sse2"))) float **sse_normalize_H(const float **H, size_t size, float **normalized) +{ + if (normalized == NULL) + { + normalized = (float **)malloc(sizeof(float *) * 2); + if (!normalized) + { + return NULL; + } + normalized[0] = malloc(sizeof(float) * size); + normalized[1] = malloc(sizeof(float) * size); + if (!normalized[0] || !normalized[1]) + { + free(normalized[0]); + free(normalized[1]); + free(normalized); + return NULL; + } + } + + // Step 1: find max magnitude + float max_mag = DBL_MIN; + size_t i = 0, idx = 4; + for (; i + 1 < size;) + { + // Load 2 complex values (4 floats: re0, im0, re1, im1) + __m128 term_r = _mm_loadu_ps(normalized[0] + i); + __m128 term_i = _mm_loadu_ps(normalized[1] + i); + + // Square: re^2 , im^2 + __m128 r_sq = _mm_mul_ps(term_r, term_r); + __m128 i_sq = _mm_mul_ps(term_i, term_i); + + // re^2 + im^2 + __m128 mag_sq = _mm_add_ps(r_sq, i_sq); + + // sqrt + __m128 mag = _mm_sqrt_ps(mag_sq); + + float buf[4]; + _mm_storeu_ps(buf, mag); + + if (buf[0] > max_mag) + max_mag = buf[0]; + if (buf[1] > max_mag) + max_mag = buf[1]; + if (buf[2] > max_mag) + max_mag = buf[2]; + if (buf[3] > max_mag) + max_mag = buf[3]; + + if (idx + i > size) + idx = size - i; // adjust for last few elements + + i += idx; + } + // Tail + for (size_t i = 0; i < size; i++) + { + float mag = sqrtf(H[0][i] * H[0][i] + H[1][i] * H[1][i]); + if (mag > max_mag) + max_mag = mag; + } + + float denom = max_mag + 1e-6f; + __m128 denom_v = _mm_set1_ps(denom); + + // Step 2: normalize magnitudes (preserve complex angle) + idx = 4; + for (i = 0; i + 1 < size;) + { + __m128 term_r = _mm_loadu_ps(normalized[0] + i); + __m128 term_i = _mm_loadu_ps(normalized[1] + i); + + // Square: re^2 , im^2 + __m128 r_sq = _mm_mul_ps(term_r, term_r); + __m128 i_sq = _mm_mul_ps(term_i, term_i); + + // re^2 + im^2 + __m128 hyp_sq = _mm_add_ps(r_sq, i_sq); + + // sqrt + __m128 hyp = _mm_sqrt_ps(hyp_sq); + + // scale = 1/denom + __m128 norm_r = _mm_mul_ps(_mm_div_ps(hyp, denom_v), term_r); + __m128 norm_i = _mm_mul_ps(_mm_div_ps(hyp, denom_v), term_i); + + _mm_storeu_ps(normalized[0] + i, norm_r); + _mm_storeu_ps(normalized[1] + i, norm_i); + + i += idx; + } + + // Tail + for (; i < size; i++) + { + float hyp = sqrtf(H[0][i] * H[0][i] + H[1][i] * H[1][i]); + normalized[0][i] = (hyp / denom) * H[0][i]; + normalized[1][i] = (hyp / denom) * H[1][i]; + } + + return normalized; +} diff --git a/src/simulation.c b/src/simulation.c new file mode 100644 index 0000000..83a761d --- /dev/null +++ b/src/simulation.c @@ -0,0 +1,225 @@ +#include "simulation.h" +#include "utils.h" +#include +#include +#include + +#define N_ZEROS_DEFAULT 3 +#define N_POLES_DEFAULT 3 + +/* ── SOA helpers ─────────────────────────────────────────────────────── */ + +static void ensure_soa_cap(singularity_soa_t *soa, size_t *cap, size_t need) { + if (need <= *cap) return; + size_t nc = (need < 8) ? 8 : need * 2; + soa->real = realloc(soa->real, nc * sizeof(double)); + soa->imag = realloc(soa->imag, nc * sizeof(double)); + soa->m_re = realloc(soa->m_re, nc * sizeof(double)); + soa->m_im = realloc(soa->m_im, nc * sizeof(double)); + soa->c_re = realloc(soa->c_re, nc * sizeof(double)); + soa->c_im = realloc(soa->c_im, nc * sizeof(double)); + soa->e = realloc(soa->e, nc * sizeof(uint8_t)); + *cap = nc; +} + +static void sync_soa(singularity_soa_t *soa, size_t *cap, + const singularity_array_t *arr) { + ensure_soa_cap(soa, cap, arr->size); + soa->size = arr->size; + for (size_t i = 0; i < arr->size; i++) { + soa->real[i] = creal(arr->data[i].val); + soa->imag[i] = cimag(arr->data[i].val); + soa->m_re[i] = creal(arr->data[i].m); + soa->m_im[i] = cimag(arr->data[i].m); + soa->c_re[i] = creal(arr->data[i].c); + soa->c_im[i] = cimag(arr->data[i].c); + soa->e[i] = arr->data[i].e; + } +} + +static void ensure_nH_cap(Simulation *s, size_t need) { + if (need <= s->nH_cap) return; + s->nH_re = realloc(s->nH_re, need * sizeof(double)); + s->nH_im = realloc(s->nH_im, need * sizeof(double)); + s->nH_cap = need; +} + +static void free_soa(singularity_soa_t *soa) { + free(soa->real); free(soa->imag); + free(soa->m_re); free(soa->m_im); + free(soa->c_re); free(soa->c_im); + free(soa->e); + memset(soa, 0, sizeof(*soa)); +} + +/* ── Random entity factory ───────────────────────────────────────────── */ + +static singularity_t make_random(double x_range[2], double y_range[2]) { + singularity_t s = {0}; + s.val = random_uniform(x_range[0], x_range[1]) + + random_uniform(y_range[0], y_range[1]) * I; + s.e = (uint8_t)random_uniform(1, 4); + s.m = random_uniform(0.5, 2.0) + random_uniform(-1, 1) * I; + s.c = random_uniform(0.5, 2.0) + random_uniform(-1, 1) * I; + return s; +} + +/* ── Public API ─────────────────────────────────────────────────────── */ + +void sim_init(Simulation *s, double x_range[2], double y_range[2], + int grid_w, int grid_h) { + memset(s, 0, sizeof(*s)); + create_singularities(&s->zeros, N_ZEROS_DEFAULT); + create_singularities(&s->poles, N_POLES_DEFAULT); + sim_randomise(s, x_range, y_range); + sim_rebuild_grid(s, x_range, y_range, grid_w, grid_h); +} + +void sim_rebuild_grid(Simulation *s, double x_range[2], double y_range[2], + int grid_w, int grid_h) { + free_s_grid_soa(&s->s_grid); + s->s_grid = generate_s_grid_soa((uint16_t)grid_h, (uint16_t)grid_w, + y_range, x_range); + ensure_nH_cap(s, (size_t)grid_w * grid_h); +} + +void sim_randomise(Simulation *s, double x_range[2], double y_range[2]) { + s->zeros.size = 0; + s->poles.size = 0; + for (int i = 0; i < N_ZEROS_DEFAULT; i++) + add_singularity(&s->zeros, make_random(x_range, y_range)); + for (int i = 0; i < N_POLES_DEFAULT; i++) + add_singularity(&s->poles, make_random(x_range, y_range)); +} + +singularity_t *sim_add_zero(Simulation *s, double complex pos) { + singularity_t z = {0}; + z.val = pos; z.e = 1; + z.m = 1.0 + 0.0 * I; + z.c = 0.5 + 0.0 * I; + add_singularity(&s->zeros, z); + return &s->zeros.data[s->zeros.size - 1]; +} + +singularity_t *sim_add_pole(Simulation *s, double complex pos) { + singularity_t p = {0}; + p.val = pos; p.e = 1; + p.m = 1.0 + 0.0 * I; + p.c = 0.5 + 0.0 * I; + add_singularity(&s->poles, p); + return &s->poles.data[s->poles.size - 1]; +} + +int sim_delete_entity(Simulation *s, singularity_t *entity) { + singularity_array_t *arrs[2] = {&s->zeros, &s->poles}; + for (int k = 0; k < 2; k++) { + for (size_t i = 0; i < arrs[k]->size; i++) { + if (&arrs[k]->data[i] == entity) { + remove_singularity(arrs[k], i); + return 1; + } + } + } + return 0; +} + +singularity_t *sim_find_entity(Simulation *s, double complex pos, + double threshold, int *out_type) { + singularity_array_t *arrs[2] = {&s->poles, &s->zeros}; + int types[2] = {2, 1}; + for (int k = 0; k < 2; k++) { + for (size_t i = 0; i < arrs[k]->size; i++) { + if (cabs(arrs[k]->data[i].val - pos) < threshold) { + if (out_type) *out_type = types[k]; + return &arrs[k]->data[i]; + } + } + } + return NULL; +} + +void sim_compute(Simulation *s, int grid_w, int grid_h) { + sync_soa(&s->zeros_soa, &s->soa_cap_z, &s->zeros); + sync_soa(&s->poles_soa, &s->soa_cap_p, &s->poles); + size_t n = (size_t)grid_w * grid_h; + ensure_nH_cap(s, n); + s->H_simd = compute_H_simd_direct(s->s_grid, s->zeros_soa, s->poles_soa, + (uint16_t)grid_h, (uint16_t)grid_w, s->H_simd); + for (size_t i = 0; i < n; i++) { + s->nH_re[i] = (double)s->H_simd[0][i]; + s->nH_im[i] = (double)s->H_simd[1][i]; + } +} + +void sim_destroy(Simulation *s) { + free_singularities(&s->zeros); + free_singularities(&s->poles); + free_soa(&s->zeros_soa); + free_soa(&s->poles_soa); + if (s->H_simd) { free(s->H_simd[0]); free(s->H_simd[1]); free(s->H_simd); } + free(s->nH_re); + free(s->nH_im); + free_s_grid_soa(&s->s_grid); +} + +double complex sim_eval_H(Simulation *s, double complex pos) { + double pr = creal(pos), pi = cimag(pos); + double Hr = 1.0, Hi = 0.0; + + /* Zeros */ + for (size_t i = 0; i < s->zeros.size; i++) { + singularity_t *z = &s->zeros.data[i]; + double zr = creal(z->val), zi = cimag(z->val); + double tr = pr - zr, ti = pi - zi; + + double fr = 1.0, fi = 0.0; + int e = z->e; + for (int j = 0; j < e; j++) { + double nfr = fr * tr - fi * ti; + double nfi = fr * ti + fi * tr; + fr = nfr; fi = nfi; + } + + double mr = creal(z->m), mi = cimag(z->m); + double cr = creal(z->c), ci = cimag(z->c); + + double tfr = fr * mr - fi * mi + cr; + double tfi = fr * mi + fi * mr + ci; + + double nHr = Hr * tfr - Hi * tfi; + double nHi = Hr * tfi + Hi * tfr; + Hr = nHr; Hi = nHi; + } + + /* Poles */ + for (size_t i = 0; i < s->poles.size; i++) { + singularity_t *p = &s->poles.data[i]; + double pr_pole = creal(p->val), pi_pole = cimag(p->val); + double tr = pr - pr_pole, ti = pi - pi_pole; + + double fr = 1.0, fi = 0.0; + int e = p->e; + for (int j = 0; j < e; j++) { + double nfr = fr * tr - fi * ti; + double nfi = fr * ti + fi * tr; + fr = nfr; fi = nfi; + } + + double mr = creal(p->m), mi = cimag(p->m); + double cr = creal(p->c), ci = cimag(p->c); + + double den_r = fr * mr - fi * mi + cr; + double den_i = fr * mi + fi * mr + ci; + + double den_mag2 = den_r * den_r + den_i * den_i; + if (den_mag2 > 1e-24) { + double nHr = (Hr * den_r + Hi * den_i) / den_mag2; + double nHi = (Hi * den_r - Hr * den_i) / den_mag2; + Hr = nHr; Hi = nHi; + } else { + Hr = 1e12; Hi = 1e12; + } + } + + return Hr + Hi * I; +} diff --git a/src/transfer_function.c b/src/transfer_function.c index 60da082..7c5606b 100644 --- a/src/transfer_function.c +++ b/src/transfer_function.c @@ -1,119 +1,401 @@ -#include "transfer_function.h" - -double complex *generate_s_grid(uint16_t height, uint16_t width, double y_range[2], double x_range[2]) -{ - double complex *s_grid = (double complex *)malloc(sizeof(double complex) * height * width); - - if (!s_grid) - return NULL; - - double y_start = y_range[0]; - double y_end = y_range[1]; - double x_start = x_range[0]; - double x_end = x_range[1]; - - // Calculate the step sizes - double dy = (y_end - y_start) / (height - 1); - double dx = (x_end - x_start) / (width - 1); - - for (uint16_t y = 0; y < height; ++y) - { - for (uint16_t x = 0; x < width; ++x) - { - double re = x_start + x * dx; - double im = y_start + y * dy; - s_grid[y * width + x] = re + im * I; - } - } - - return s_grid; -} - -void create_singularities(singularity_array_t *arr, size_t initial_capacity) -{ - arr->data = malloc(initial_capacity * sizeof(singularity_t)); - arr->size = 0; - arr->capacity = initial_capacity; -} - -void add_singularity(singularity_array_t *arr, singularity_t s) -{ - if (arr->size == arr->capacity) - { - arr->capacity = arr->capacity ? arr->capacity * 2 : 4; - arr->data = realloc(arr->data, arr->capacity * sizeof(singularity_t)); - } - arr->data[arr->size++] = s; -} - -void remove_singularity(singularity_array_t *arr, size_t index) -{ - if (index >= arr->size) - return; - for (size_t i = index; i < arr->size - 1; ++i) - { - arr->data[i] = arr->data[i + 1]; - } - arr->size--; -} - -void free_singularities(singularity_array_t *arr) -{ - free(arr->data); - arr->data = NULL; - arr->size = 0; - arr->capacity = 0; -} - -double complex *compute_H(double complex *s_grid, singularity_array_t zeros_arr, singularity_array_t poles_arr, uint16_t height, uint16_t width, double complex *H) -{ - uint32_t length = height * width; - - if (H == NULL) - { - H = malloc(sizeof(double complex) * length); - if (!H) - return NULL; - } - - for (uint32_t i = 0; i < length; i++) - { - double complex s = s_grid[i]; - double complex num = 1.0 + 0.0 * I; - double complex den = 1.0 + 0.0 * I; - - // Zeros - for (size_t z_i = 0; z_i < zeros_arr.size; z_i++) - { - singularity_t z = zeros_arr.data[z_i]; - double complex term = s - z.val; - double complex prod = 1.0 + 0.0 * I; - - for (uint8_t e = 0; e < z.e; e++) - { - prod *= term; - } - - num *= (z.m * prod + z.c); - } - - // Poles - for (size_t p_i = 0; p_i < poles_arr.size; p_i++) - { - singularity_t p = poles_arr.data[p_i]; - double complex term = s - p.val; - double complex prod = 1.0 + 0.0 * I; - - for (uint8_t e = 0; e < p.e; e++) - { - prod *= term; - } - - den *= (p.m * prod + p.c); - } - - H[i] = num / (den + 1e-15); // Avoid division by zero - } - - return H; -} +#include "transfer_function.h" +#include "cpu_detect.h" +#include "simd_sse.h" +#include "simd_avx.h" +#include "simd_avx512.h" +#include +#include +#include +#include +#include + +double complex *generate_s_grid(uint16_t height, uint16_t width, double y_range[2], double x_range[2]) +{ + double complex *s_grid = (double complex *)malloc(sizeof(double complex) * height * width); + if (!s_grid) return NULL; + + double y_start = y_range[0]; + double y_end = y_range[1]; + double x_start = x_range[0]; + double x_end = x_range[1]; + + double dy = (y_end - y_start) / (height - 1); + double dx = (x_end - x_start) / (width - 1); + + for (uint16_t y = 0; y < height; ++y) + { + for (uint16_t x = 0; x < width; ++x) + { + double re = x_start + x * dx; + double im = y_start + y * dy; + s_grid[y * width + x] = re + im * I; + } + } + return s_grid; +} + +s_grid_soa_t generate_s_grid_soa(uint16_t height, uint16_t width, double y_range[2], double x_range[2]) +{ + s_grid_soa_t grid; + uint32_t length = (uint32_t)height * width; + grid.length = length; + + grid.real = (float *)malloc(sizeof(float) * length); + grid.imag = (float *)malloc(sizeof(float) * length); + + if (!grid.real || !grid.imag) { + free(grid.real); free(grid.imag); + grid.real = NULL; grid.imag = NULL; grid.length = 0; + return grid; + } + + float y_start = (float)y_range[0]; + float y_end = (float)y_range[1]; + float x_start = (float)x_range[0]; + float x_end = (float)x_range[1]; + + float dy = (y_end - y_start) / (height - 1); + float dx = (x_end - x_start) / (width - 1); + + for (uint16_t y = 0; y < height; ++y) + { + for (uint16_t x = 0; x < width; ++x) + { + uint32_t idx = y * width + x; + grid.real[idx] = x_start + x * dx; + grid.imag[idx] = y_start + y * dy; + } + } + + return grid; +} + +void free_s_grid_soa(s_grid_soa_t *grid) +{ + if (grid) { + free(grid->real); free(grid->imag); + grid->real = NULL; grid->imag = NULL; grid->length = 0; + } +} + +void create_singularities(singularity_array_t *arr, size_t initial_capacity) +{ + arr->data = malloc(initial_capacity * sizeof(singularity_t)); + arr->size = 0; + arr->capacity = initial_capacity; +} + +void add_singularity(singularity_array_t *arr, singularity_t s) +{ + if (arr->size == arr->capacity) { + arr->capacity += 2; + arr->data = realloc(arr->data, arr->capacity * sizeof(singularity_t)); + } + arr->data[arr->size++] = s; +} + +void remove_singularity(singularity_array_t *arr, size_t index) +{ + if (index >= arr->size) return; + for (size_t i = index; i < arr->size - 1; ++i) { + arr->data[i] = arr->data[i + 1]; + } + arr->size--; +} + +void free_singularities(singularity_array_t *arr) +{ + free(arr->data); + arr->data = NULL; + arr->size = 0; + arr->capacity = 0; +} + +new_singularity_array_t convert_soa_to_simd(singularity_soa_t soa) +{ + new_singularity_array_t simd_arr; + simd_arr.count = soa.size; + + if (soa.size == 0) { + simd_arr.r = NULL; simd_arr.i = NULL; simd_arr.e = NULL; + simd_arr.m_r = NULL; simd_arr.m_i = NULL; simd_arr.c_r = NULL; simd_arr.c_i = NULL; + return simd_arr; + } + + simd_arr.r = (float *)malloc(sizeof(float) * soa.size); + simd_arr.i = (float *)malloc(sizeof(float) * soa.size); + simd_arr.e = (uint8_t *)malloc(sizeof(uint8_t) * soa.size); + simd_arr.m_r = (float *)malloc(sizeof(float) * soa.size); + simd_arr.m_i = (float *)malloc(sizeof(float) * soa.size); + simd_arr.c_r = (float *)malloc(sizeof(float) * soa.size); + simd_arr.c_i = (float *)malloc(sizeof(float) * soa.size); + + for (size_t j = 0; j < soa.size; j++) { + simd_arr.r[j] = (float)soa.real[j]; + simd_arr.i[j] = (float)soa.imag[j]; + simd_arr.e[j] = soa.e[j]; + simd_arr.m_r[j] = (float)soa.m_re[j]; + simd_arr.m_i[j] = (float)soa.m_im[j]; + simd_arr.c_r[j] = (float)soa.c_re[j]; + simd_arr.c_i[j] = (float)soa.c_im[j]; + } + + return simd_arr; +} + +void free_simd_array(new_singularity_array_t *arr) +{ + free(arr->r); free(arr->i); free(arr->e); + free(arr->m_r); free(arr->m_i); free(arr->c_r); free(arr->c_i); + arr->r = NULL; arr->i = NULL; arr->e = NULL; + arr->m_r = NULL; arr->m_i = NULL; arr->c_r = NULL; arr->c_i = NULL; + arr->count = 0; +} + +float **compute_H_simd_direct(s_grid_soa_t s_grid, singularity_soa_t zeros_soa, singularity_soa_t poles_soa, uint16_t height, uint16_t width, float **H) +{ + static cpu_features_t features = CPU_FEATURE_NONE; + static int detected = 0; + if (!detected) { + features = detect_cpu_features(); + detected = 1; + } + + uint32_t length = height * width; + new_singularity_array_t simd_zeros = convert_soa_to_simd(zeros_soa); + new_singularity_array_t simd_poles = convert_soa_to_simd(poles_soa); + + float x_range[2] = {s_grid.real[0], s_grid.real[length - 1]}; + float y_range[2] = {s_grid.imag[0], s_grid.imag[length - 1]}; + + if (features & CPU_FEATURE_AVX512F) { + H = avx512_compute_H(x_range, y_range, simd_zeros, simd_poles, height, width, H); + } else if (features & CPU_FEATURE_AVX) { + H = avx_compute_H(x_range, y_range, simd_zeros, simd_poles, height, width, H); + } else if (features & CPU_FEATURE_SSE2) { + H = sse2_compute_H(x_range, y_range, simd_zeros, simd_poles, height, width, H); + } + + free_simd_array(&simd_zeros); + free_simd_array(&simd_poles); + + return H; +} + +void normalize_H_soa(double *in_re, double *in_im, size_t size, double *out_re, double *out_im) +{ + double max_mag = DBL_MIN; + + for (size_t i = 0; i < size; ++i) { + double mag = sqrt(in_re[i]*in_re[i] + in_im[i]*in_im[i]); + if (mag > max_mag) max_mag = mag; + } + + double denom = max_mag + 1e-6; + + for (size_t i = 0; i < size; ++i) { + double mag = sqrt(in_re[i]*in_re[i] + in_im[i]*in_im[i]); + double scale = mag > 1e-9 ? (mag / denom) / mag : 0.0; + if (mag / denom > 1.0) scale = 1.0 / mag; + out_re[i] = in_re[i] * scale; + out_im[i] = in_im[i] * scale; + } +} + +// Fast Histogram-based 99th Percentile +static double find_99th_percentile(double *log_mag, size_t size) { + if (size == 0) return 1.0; + + double min_v = 1e30, max_v = -1e30; + for (size_t i = 0; i < size; i++) { + if (log_mag[i] < min_v) min_v = log_mag[i]; + if (log_mag[i] > max_v) max_v = log_mag[i]; + } + + if (max_v <= min_v) return max_v; + + #define HIST_BINS 1024 + int histogram[HIST_BINS] = {0}; + double range = max_v - min_v; + double bin_width = range / HIST_BINS; + + for (size_t i = 0; i < size; i++) { + int bin = (int)((log_mag[i] - min_v) / bin_width); + if (bin < 0) bin = 0; + if (bin >= HIST_BINS) bin = HIST_BINS - 1; + histogram[bin]++; + } + + size_t count = 0; + size_t target = (size_t)(size * 0.99); + for (int k = 0; k < HIST_BINS; k++) { + count += histogram[k]; + if (count >= target) { + return min_v + (k + 1) * bin_width; + } + } + return max_v; +} + +void normalize_H_log_soa(double *in_re, double *in_im, size_t size, double *out_re, double *out_im) +{ + double *log_mag = (double *)malloc(size * sizeof(double)); + if (!log_mag) return; + + for (size_t i = 0; i < size; ++i) { + log_mag[i] = log1p(sqrt(in_re[i]*in_re[i] + in_im[i]*in_im[i])); + } + + // Replace qsort with histogram approach + double max_val = find_99th_percentile(log_mag, size); + + double denom = max_val + 1e-6; + + for (size_t i = 0; i < size; ++i) { + double orig_mag = sqrt(in_re[i]*in_re[i] + in_im[i]*in_im[i]); + double new_mag = log_mag[i] / denom; + if (new_mag > 1.0) new_mag = 1.0; + if (new_mag < 0.0) new_mag = 0.0; + + double scale = orig_mag > 1e-9 ? new_mag / orig_mag : 0.0; + out_re[i] = in_re[i] * scale; + out_im[i] = in_im[i] * scale; + } + free(log_mag); +} + +void normalize_H_log_steps_soa(double *in_re, double *in_im, size_t size, int steps, double *out_re, double *out_im) +{ + if (steps <= 0) return; + double *log_mag = (double *)malloc(size * sizeof(double)); + if (!log_mag) return; + + for (size_t i = 0; i < size; ++i) { + log_mag[i] = log1p(sqrt(in_re[i]*in_re[i] + in_im[i]*in_im[i])); + } + + double max_val = find_99th_percentile(log_mag, size); + double denom = max_val + 1e-6; + double step_size = 1.0 / steps; + + for (size_t i = 0; i < size; ++i) { + double orig_mag = sqrt(in_re[i]*in_re[i] + in_im[i]*in_im[i]); + double raw_mag = log_mag[i] / denom; + if (raw_mag > 1.0) raw_mag = 1.0; + if (raw_mag < 0.0) raw_mag = 0.0; + + int step_index = (int)(raw_mag * steps + 0.5); + double quantized_mag = step_index * step_size; + + double scale = orig_mag > 1e-9 ? quantized_mag / orig_mag : 0.0; + out_re[i] = in_re[i] * scale; + out_im[i] = in_im[i] * scale; + } + free(log_mag); +} + +static inline uint8_t clamp(float x) { + return (uint8_t)(x < 0 ? 0 : (x > 255 ? 255 : x)); +} + +void hsv_to_rgb_uint8(uint8_t H, uint8_t S, uint8_t V, uint8_t *R, uint8_t *G, uint8_t *B) +{ + float h = H * 2.0f; + float s = S / 255.0f; + float v = V / 255.0f; + float c = v * s; + float x = c * (1 - fabsf(fmodf(h / 60.0f, 2) - 1)); + float m = v - c; + float r = 0, g = 0, b = 0; + + if (h < 60) { r = c; g = x; b = 0; } + else if (h < 120) { r = x; g = c; b = 0; } + else if (h < 180) { r = 0; g = c; b = x; } + else if (h < 240) { r = 0; g = x; b = c; } + else if (h < 300) { r = x; g = 0; b = c; } + else { r = c; g = 0; b = x; } + + *R = clamp((r + m) * 255.0f); + *G = clamp((g + m) * 255.0f); + *B = clamp((b + m) * 255.0f); +} + +void H_g_img(double *re, double *im, size_t size, img_t H_img) +{ + for (uint32_t i = 0; i < size; i++) + { + uint32_t c = (uint32_t)(sqrt(re[i]*re[i] + im[i]*im[i]) * 255); + if (c > 255) c = 255; + (H_img.data)[i] = (255 << 24) | (c << 16) | (c << 8) | c; + } +} + +void H_c1_img(double *re, double *im, size_t size, img_t H_img) +{ + for (uint32_t i = 0; i < size; i++) + { + uint32_t r = (uint32_t)(((re[i] + 1) / 2) * 255u); + uint32_t b = (uint32_t)(((im[i] + 1) / 2) * 255u); + if (r > 255) r = 255; + if (b > 255) b = 255; + (H_img.data)[i] = (255 << 24) | (r << 16) | b; + } +} + +void H_c2_img(double *re, double *im, size_t size, img_t H_img) +{ + for (uint32_t i = 0; i < size; i++) + { + double mag = sqrt(re[i]*re[i] + im[i]*im[i]); + uint32_t r = (uint32_t)(((re[i] + 1) / 2) * mag * 255u); + uint32_t b = (uint32_t)(((im[i] + 1) / 2) * mag * 255u); + if (r > 255) r = 255; + if (b > 255) b = 255; + (H_img.data)[i] = (255 << 24) | (r << 16) | b; + } +} + +void H_c3_img(double *re, double *im, size_t size, img_t H_img) +{ + for (uint32_t i = 0; i < size; i++) + { + double mag = sqrt(re[i]*re[i] + im[i]*im[i]); + uint32_t r = (uint32_t)(((re[i] + 1) / 2) * 255u); + uint32_t g = (uint32_t)(mag * 255u); + uint32_t b = (uint32_t)(((im[i] + 1) / 2) * 255u); + if (r > 255) r = 255; + if (b > 255) b = 255; + if (g > 255) g = 255; + (H_img.data)[i] = (255 << 24) | (r << 16) | (g << 8) | b; + } +} + +void H_c4_img(double *re, double *im, size_t size, img_t H_img) +{ + for (uint32_t i = 0; i < size; i++) + { + double mag = sqrt(re[i]*re[i] + im[i]*im[i]); + uint32_t r = (uint32_t)(((re[i] + 1) / 2) * mag * 255u); + uint32_t g = (uint32_t)((1 - mag) * 255u); + uint32_t b = (uint32_t)(((im[i] + 1) / 2) * mag * 255u); + if (r > 255) r = 255; + if (b > 255) b = 255; + if (g > 255) g = 255; + (H_img.data)[i] = (255 << 24) | (r << 16) | (g << 8) | b; + } +} + +void H_c5_img(double *re, double *im, size_t size, img_t H_img) +{ + for (uint32_t i = 0; i < size; i++) + { + double mag = sqrt(re[i]*re[i] + im[i]*im[i]); + uint8_t h = (uint8_t)(((re[i] + 1) / 2) * (2.f / 3.f) * 255); + uint8_t s = (uint8_t)(((im[i] + 1) / 2) * 255); + uint8_t v = (uint8_t)(mag * 255); + + uint8_t r, g, b; + hsv_to_rgb_uint8(h, s, v, &r, &g, &b); + (H_img.data)[i] = (255 << 24) | (uint32_t)r << 16 | (uint32_t)g << 8 | b; + } +} \ No newline at end of file diff --git a/src/ui.c b/src/ui.c new file mode 100644 index 0000000..5541554 --- /dev/null +++ b/src/ui.c @@ -0,0 +1,594 @@ +#include "ui.h" +#include "font_utils.h" +#include "img_utils.h" +#include "utils.h" +#include +#include +#include +#include +#include "audio.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +/* ── Layout constants ────────────────────────────────────────────────── */ +#define POPUP_W 340 +#define POPUP_H 320 +#define SLIDER_W 300 +#define SLIDER_X_OFF 20 +#define SLIDER_Y_START 50 +#define SLIDER_Y_STEP 45 +#define RADIUS_POPUP 12 +#define RADIUS_BTN 6 +#define RADIUS_HANDLE 6 +#define SETTINGS_W 340 + +/* ═══════════════════════════════════════════════════════════════════════ + Internal helpers + ═══════════════════════════════════════════════════════════════════════ */ + +static void clamp_rect(int *px, int *py, int w, int h, + int img_w, int img_h) { + if (*px + w > img_w) *px = img_w - w; + if (*py + h > img_h) *py = img_h - h; + if (*px < 0) *px = 0; + if (*py < 0) *py = 0; +} + +/* Slider range for a given index and edit_mode */ +static void slider_range(int idx, int edit_mode, + double *mn, double *mx, int *is_int) { + *mn = -5; *mx = 5; *is_int = 0; + if (edit_mode == 0) { + if (idx == 4) { *mn = 1; *mx = 10; *is_int = 1; } + } else { + if (idx == 0) { *mn = 0; *mx = 5; } + if (idx == 1) { *mn = -180; *mx = 180; } + if (idx == 2) { *mn = 1; *mx = 10; *is_int = 1; } + } +} + +/* Apply slider value at screen-x mx to entity (called from both click & drag) */ +static void apply_slider(singularity_t *e, int idx, int mx, int popup_x, + int edit_mode) { + double mn, mx_v; int is_int; + slider_range(idx, edit_mode, &mn, &mx_v, &is_int); + + int rx = popup_x + SLIDER_X_OFF; + double t = (double)(mx - rx) / (double)SLIDER_W; + if (t < 0) t = 0; + if (t > 1) t = 1; + double val = mn + t * (mx_v - mn); + if (is_int) val = (double)(int)(val + 0.5); + + if (edit_mode == 0) { + if (idx == 0) e->m = val + I * cimag(e->m); + if (idx == 1) e->m = creal(e->m) + I * val; + if (idx == 2) e->c = val + I * cimag(e->c); + if (idx == 3) e->c = creal(e->c) + I * val; + if (idx == 4) e->e = (uint8_t)val; + } else { + double complex ratio = (cabs(e->m) > 1e-6) ? (-e->c / e->m) : 0; + double cur_e = (double)e->e; + double cur_R = pow(cabs(ratio), 1.0 / cur_e); + double cur_Ang = carg(ratio) * 180.0 / M_PI; + if (idx == 0) cur_R = val; + if (idx == 1) cur_Ang = val; + if (idx == 2) cur_e = val; + e->m = 1.0 + 0.0 * I; + e->e = (uint8_t)cur_e; + e->c = -pow(cur_R, cur_e) * cexp(I * cur_Ang * M_PI / 180.0); + } +} + +/* ── Single slider widget ─────────────────────────────────────────────── */ +static void draw_slider_widget(img_t *img, int rx, int ry, + const char *label, double val, + double mn, double mx, int is_int) { + char buf[64]; + if (is_int) snprintf(buf, sizeof(buf), "%s: %d", label, (int)val); + else snprintf(buf, sizeof(buf), "%s: %.2f", label, val); + draw_string_colored(img, buf, rx, ry - 14, 0xFFE0E0E0); + + int ty = ry + 6; + draw_rounded_rect_fill(img, rx, ty, SLIDER_W, 4, 2, 0xFF404040); + + double t = (mx > mn) ? (val - mn) / (mx - mn) : 0; + if (t < 0) t = 0; + if (t > 1) t = 1; + int hx = rx + (int)(t * (SLIDER_W - 12)); + draw_rounded_rect_fill(img, hx, ty - 4, 12, 14, RADIUS_HANDLE, 0xFFC0C0C0); +} + +/* ═══════════════════════════════════════════════════════════════════════ + Grid + ═══════════════════════════════════════════════════════════════════════ */ +void ui_draw_grid(img_t *img, double x0, double x1, double y0, double y1) { + double dx = x1 - x0, dy = y1 - y0; + double sx = img->width / dx; + double sy = img->height / dy; + + int32_t cx = (int32_t)((0 - x0) * sx); + int32_t cy = (int32_t)((0 - y0) * sy); + + for (int i = (int)ceil(x0); i <= (int)floor(x1); i++) { + if (i == 0) continue; + int32_t x = (int32_t)((i - x0) * sx); + draw_rect_blend(img, x, 0, 1, (int32_t)img->height, 0x18FFFFFF); + } + for (int j = (int)ceil(y0); j <= (int)floor(y1); j++) { + if (j == 0) continue; + int32_t y = (int32_t)((j - y0) * sy); + draw_rect_blend(img, 0, y, (int32_t)img->width, 1, 0x18FFFFFF); + } + if (cx >= 0 && cx < (int32_t)img->width) + draw_rect_fill(img, cx - 1, 0, 2, (int32_t)img->height, 0x60FFFFFF); + if (cy >= 0 && cy < (int32_t)img->height) + draw_rect_fill(img, 0, cy - 1, (int32_t)img->width, 2, 0x60FFFFFF); +} + +/* ═══════════════════════════════════════════════════════════════════════ + Entities + ═══════════════════════════════════════════════════════════════════════ */ +void ui_draw_entities(img_t *img, Simulation *sim, + double x_range[2], double y_range[2], + singularity_t *selected, singularity_t *move_target, + UIFlags flags) { + if (!flags.hud_visible || flags.circle_mode == CIRCLE_NONE) return; + + double dx = x_range[1] - x_range[0]; + double dy = y_range[1] - y_range[0]; + double sx = img->width / dx; + double sy = img->height / dy; + + singularity_array_t *arrs[2] = {&sim->zeros, &sim->poles}; + int types[2] = {3, 1}; /* channel: 3=blue(zero), 1=red(pole) */ + + for (int k = 0; k < 2; k++) { + for (size_t i = 0; i < arrs[k]->size; i++) { + singularity_t *e = &arrs[k]->data[i]; + int32_t cx = (int32_t)((creal(e->val) - x_range[0]) * sx); + int32_t cy = (int32_t)((cimag(e->val) - y_range[0]) * sy); + + double math_R = (cabs(e->m) > 1e-6) + ? pow(cabs(e->c / e->m), 1.0 / e->e) : 0; + int32_t radius = (int32_t)(math_R * sx); + if (radius < 2) radius = 2; + + int ch = types[k]; + int is_highlighted = (e == selected || e == move_target); + + /* Highlight ring (outer glow) */ + if (is_highlighted) { + draw_circle_glow(img, cx, cy, radius + 4, ch); + fill_circle_glow(img, cx, cy, 4, ch); + } + + /* Circle itself */ + switch (flags.circle_mode) { + case CIRCLE_SOLID: + draw_circle_aa_dynamic(img, cx, cy, radius, ch); + fill_circle_dynamic(img, cx, cy, 3, ch); + break; + case CIRCLE_TRANSPARENT: + draw_circle_aa_alpha(img, cx, cy, radius, ch, 80); + fill_circle_dynamic(img, cx, cy, 3, ch); + break; + case CIRCLE_DASHED: + draw_circle_dashed_aa_dynamic(img, cx, cy, radius, ch); + fill_circle_dynamic(img, cx, cy, 3, ch); + break; + case CIRCLE_DOT_ONLY: + fill_circle_dynamic(img, cx, cy, 4, ch); + break; + default: + break; + } + } + } +} + +/* ═══════════════════════════════════════════════════════════════════════ + Tooltip + ═══════════════════════════════════════════════════════════════════════ */ +void ui_draw_tooltip(img_t *img, int mx, int my, singularity_t *s, int type) { + int w = 240, h = 100; + int x = mx + 20, y = my - 10; + if (x + w > (int)img->width) x = mx - w - 20; + if (y + h > (int)img->height) y = (int)img->height - h; + if (y < 0) y = 0; + + draw_rounded_rect_blend(img, x, y, w, h, RADIUS_POPUP, 0xE0101010); + + char buf[64]; + const char *typestr = (type == 1) ? "Zero" : "Pole"; + uint32_t tc = (type == 1) ? 0xFF9090FF : 0xFFFF9090; + draw_string_colored(img, typestr, x + 12, y + 10, tc); + snprintf(buf, sizeof(buf), "Pos: %.2f%+.2fi", creal(s->val), cimag(s->val)); + draw_string_colored(img, buf, x + 12, y + 30, 0xFFCCCCCC); + snprintf(buf, sizeof(buf), "e: %d", s->e); + draw_string_colored(img, buf, x + 12, y + 48, 0xFFCCCCCC); + snprintf(buf, sizeof(buf), "m: %.1f%+.1fi", creal(s->m), cimag(s->m)); + draw_string_colored(img, buf, x + 12, y + 66, 0xFFFFA0A0); + snprintf(buf, sizeof(buf), "c: %.1f%+.1fi", creal(s->c), cimag(s->c)); + draw_string_colored(img, buf, x + 130, y + 66, 0xFFA0A0FF); +} + +/* ═══════════════════════════════════════════════════════════════════════ + Popup panel + ═══════════════════════════════════════════════════════════════════════ */ +void ui_draw_popup(img_t *img, Popup *p) { + if (!p->is_open || !p->entity) return; + + int px = p->x, py = p->y; + clamp_rect(&px, &py, POPUP_W, POPUP_H, (int)img->width, (int)img->height); + + draw_rounded_rect_blend(img, px, py, POPUP_W, POPUP_H, RADIUS_POPUP, 0xE0181818); + + uint32_t hcol = (p->entity_type == 1) ? 0xFF005090 : 0xFF900000; + draw_rounded_rect_fill(img, px+2, py+2, POPUP_W-4, 32, RADIUS_POPUP-2, hcol|0xFF000000); + const char *title = (p->entity_type == 1) ? "ZERO" : "POLE"; + draw_string_colored(img, title, px + 12, py + 10, 0xFFFFFFFF); + + int btn_x = px + POPUP_W - 100, btn_y = py + 6; + draw_rounded_rect_fill(img, btn_x, btn_y, 80, 24, RADIUS_BTN, 0xFF303030); + draw_string_colored(img, (p->edit_mode == 0) ? "CART" : "POLAR", + btn_x + 10, btn_y + 8, 0xFFCCCCCC); + + int y_off = SLIDER_Y_START; + singularity_t *e = p->entity; + + if (p->edit_mode == 0) { + draw_slider_widget(img, px+SLIDER_X_OFF, py+y_off, "m.Re", creal(e->m), -5, 5, 0); y_off+=SLIDER_Y_STEP; + draw_slider_widget(img, px+SLIDER_X_OFF, py+y_off, "m.Im", cimag(e->m), -5, 5, 0); y_off+=SLIDER_Y_STEP; + draw_slider_widget(img, px+SLIDER_X_OFF, py+y_off, "c.Re", creal(e->c), -5, 5, 0); y_off+=SLIDER_Y_STEP; + draw_slider_widget(img, px+SLIDER_X_OFF, py+y_off, "c.Im", cimag(e->c), -5, 5, 0); y_off+=SLIDER_Y_STEP; + draw_slider_widget(img, px+SLIDER_X_OFF, py+y_off, "Exp", (double)e->e, 1, 10, 1); + } else { + double complex ratio = (cabs(e->m) > 1e-6) ? (-e->c / e->m) : 0; + double cur_e = (double)e->e; + double cur_R = pow(cabs(ratio), 1.0 / cur_e); + double cur_Ang = carg(ratio) * 180.0 / M_PI; + draw_slider_widget(img, px+SLIDER_X_OFF, py+y_off, "Rad", cur_R, 0, 5, 0); y_off+=SLIDER_Y_STEP; + draw_slider_widget(img, px+SLIDER_X_OFF, py+y_off, "Rot", cur_Ang, -180, 180, 0); y_off+=SLIDER_Y_STEP; + draw_slider_widget(img, px+SLIDER_X_OFF, py+y_off, "Exp", cur_e, 1, 10, 1); + } +} + +/* ═══════════════════════════════════════════════════════════════════════ + HUD overlay + ═══════════════════════════════════════════════════════════════════════ */ +static const char *circle_mode_name(CircleMode m) { + switch (m) { + case CIRCLE_SOLID: return "Solid"; + case CIRCLE_TRANSPARENT: return "Transparent"; + case CIRCLE_DASHED: return "Dashed"; + case CIRCLE_DOT_ONLY: return "Dot"; + case CIRCLE_NONE: return "None"; + default: return "?"; + } +} + +void ui_draw_hud(img_t *img, float fps, const RenderCfg *cfg, + const UIFlags *flags, int cursor_mode, int is_fullscreen) { + if (!flags->hud_visible) return; + + char buf[80]; + snprintf(buf, sizeof(buf), "FPS: %.1f", fps); + draw_string_dynamic(img, buf, 10, 10, 2); + + static const char *nmap_names[] = {"Linear", "Log", "Steps"}; + snprintf(buf, sizeof(buf), "Norm: %s", nmap_names[cfg->n_map]); + draw_string_dynamic(img, buf, 10, 22, 2); + + if (cfg->n_map == 2) { + snprintf(buf, sizeof(buf), "Steps: %d", cfg->steps); + draw_string_dynamic(img, buf, 10, 34, 2); + } + snprintf(buf, sizeof(buf), "CMap: %d | %s | %s", + cfg->c_map, + cursor_mode == 0 ? "SELECT" : "MOVE", + is_fullscreen ? "FS" : "WIN"); + draw_string_dynamic(img, buf, 10, 46, 2); + + snprintf(buf, sizeof(buf), "Circle: %s [H=HUD C=Circle S=Set]", + circle_mode_name(flags->circle_mode)); + draw_string_dynamic(img, buf, 10, 58, 2); + + snprintf(buf, sizeof(buf), "[P=Pic O=Orb U=Pls F=Fld D=Dom A=Aud 1-9=Pre]"); + draw_string_dynamic(img, buf, 10, 70, 2); + + if (flags->audio_enabled) { + double freq, amp; + audio_get_params(&freq, &); + snprintf(buf, sizeof(buf), "Audio: %.1f Hz | Amp: %.2f", freq, amp); + draw_string_dynamic(img, buf, 10, 82, 2); + } +} + +/* ═══════════════════════════════════════════════════════════════════════ + Settings panel + ═══════════════════════════════════════════════════════════════════════ */ +static const int SETTINGS_H = 380; + +void ui_draw_settings(img_t *img, SettingsPanel *sp, + RenderCfg *cfg, UIFlags *flags, + int screen_w, int screen_h, int is_fullscreen, + int *request_fs_toggle) { + (void)request_fs_toggle; + if (!sp->is_open) return; + + int px = (screen_w - SETTINGS_W) / 2; + int py = (screen_h - SETTINGS_H) / 2; + clamp_rect(&px, &py, SETTINGS_W, SETTINGS_H, (int)img->width, (int)img->height); + + draw_rounded_rect_blend(img, px, py, SETTINGS_W, SETTINGS_H, RADIUS_POPUP, 0xF0141414); + draw_rounded_rect_fill (img, px+2, py+2, SETTINGS_W-4, 30, RADIUS_POPUP-2, 0xFF202060); + draw_string_colored(img, "SETTINGS", px+12, py+10, 0xFFFFFFFF); + + /* Close button */ + draw_rounded_rect_fill(img, px+SETTINGS_W-36, py+6, 28, 22, 4, 0xFF602020); + draw_string_colored(img, "X", px+SETTINGS_W-26, py+12, 0xFFFFFFFF); + + int y = py + 44; + char buf[80]; + + /* Window section */ + draw_string_colored(img, "-- Window --", px+12, y, 0xFF8888CC); y += 14; + snprintf(buf, sizeof(buf), "Resolution: %d x %d", screen_w, screen_h); + draw_string_colored(img, buf, px+12, y, 0xFFCCCCCC); y += 14; + snprintf(buf, sizeof(buf), "Mode: %s", is_fullscreen ? "Fullscreen" : "Windowed"); + draw_string_colored(img, buf, px+12, y, 0xFFCCCCCC); y += 16; + draw_rounded_rect_fill(img, px+12, y, 120, 22, 4, 0xFF304060); + draw_string_colored(img, "Toggle FS [Alt+Enter]", px+16, y+7, 0xFFAAAAFF); y += 32; + + /* Render section */ + draw_string_colored(img, "-- Rendering --", px+12, y, 0xFF88CC88); y += 14; + + /* Colour map row */ + draw_string_colored(img, "Colour Map:", px+12, y, 0xFFCCCCCC); + for (int i = 0; i < 6; i++) { + int bx = px + 12 + i * 42; + uint32_t bc = (i == cfg->c_map) ? 0xFF506080 : 0xFF303030; + draw_rounded_rect_fill(img, bx, y+12, 36, 18, 3, bc); + snprintf(buf, sizeof(buf), "%d", i); + draw_string_colored(img, buf, bx+14, y+17, 0xFFCCCCCC); + } + y += 36; + + /* Normalisation row */ + draw_string_colored(img, "Normalise:", px+12, y, 0xFFCCCCCC); + static const char *nm[] = {"Lin","Log","Stp"}; + for (int i = 0; i < 3; i++) { + int bx = px + 12 + i * 60; + uint32_t bc = (i == cfg->n_map) ? 0xFF506080 : 0xFF303030; + draw_rounded_rect_fill(img, bx, y+12, 52, 18, 3, bc); + draw_string_colored(img, nm[i], bx+14, y+17, 0xFFCCCCCC); + } + y += 36; + + /* Steps slider */ + draw_slider_widget(img, px+SLIDER_X_OFF, y+14, "Steps", (double)cfg->steps, 1, 64, 1); + y += 40; + + /* Circle mode row */ + draw_string_colored(img, "Entity Style:", px+12, y, 0xFFCCCCCC); + static const char *cm[] = {"Sol","Trn","Dsh","Dot","None"}; + for (int i = 0; i < CIRCLE_MODE_COUNT; i++) { + int bx = px + 12 + i * 56; + uint32_t bc = (i == (int)flags->circle_mode) ? 0xFF506080 : 0xFF303030; + draw_rounded_rect_fill(img, bx, y+12, 50, 18, 3, bc); + draw_string_colored(img, cm[i], bx+10, y+17, 0xFFCCCCCC); + } + y += 36; + + /* HUD toggle */ + uint32_t hud_bc = flags->hud_visible ? 0xFF304050 : 0xFF502020; + draw_rounded_rect_fill(img, px+12, y, 100, 22, 4, hud_bc); + draw_string_colored(img, flags->hud_visible ? "HUD: ON [H]" : "HUD: OFF [H]", + px+16, y+7, 0xFFCCCCCC); +} + +/* ═══════════════════════════════════════════════════════════════════════ + Popup interaction + ═══════════════════════════════════════════════════════════════════════ */ +int popup_hit_slider(const Popup *p, int mx, int my, int *out_idx) { + if (!p->is_open || !p->entity) return 0; + int px = p->x, py = p->y; + int max_s = (p->edit_mode == 0) ? 5 : 3; + for (int idx = 0; idx < max_s; idx++) { + int y_off = SLIDER_Y_START + idx * SLIDER_Y_STEP; + int ry = py + y_off - 14; + int rx = px + SLIDER_X_OFF; + if (mx >= rx && mx <= rx + SLIDER_W && my >= ry && my <= ry + 30) { + *out_idx = idx; + return 1; + } + } + return 0; +} + +int popup_hit_mode_btn(const Popup *p, int mx, int my) { + if (!p->is_open) return 0; + int btn_x = p->x + POPUP_W - 100, btn_y = p->y + 6; + return (mx >= btn_x && mx <= btn_x + 80 && my >= btn_y && my <= btn_y + 24); +} + +void popup_apply_slider(Popup *p, int mx) { + if (!p->entity || p->dragging_slider < 0) return; + apply_slider(p->entity, p->dragging_slider, mx, p->x, p->edit_mode); +} + +/* ═══════════════════════════════════════════════════════════════════════ + Settings interaction + ═══════════════════════════════════════════════════════════════════════ */ +int settings_hit(SettingsPanel *sp, RenderCfg *cfg, UIFlags *flags, + int mx, int my, int screen_w, int screen_h, + int is_fullscreen, int *cmap_delta, int *nmap_delta) { + (void)is_fullscreen; + if (!sp->is_open) return SETTINGS_HIT_NONE; + *cmap_delta = 0; *nmap_delta = 0; + + int px = (screen_w - SETTINGS_W) / 2; + int py = (screen_h - SETTINGS_H) / 2; + if (px < 0) px = 0; + if (py < 0) py = 0; + + /* Close button */ + if (mx >= px+SETTINGS_W-36 && mx <= px+SETTINGS_W-8 && + my >= py+6 && my <= py+28) { + sp->is_open = 0; + return SETTINGS_HIT_CLOSE; + } + + int y = py + 44; + + /* Window section */ + y += 14; /* label */ + y += 14; /* resolution text */ + y += 14; /* mode text */ + /* Fullscreen toggle button */ + if (mx >= px+12 && mx <= px+132 && my >= y && my <= y+22) + return SETTINGS_HIT_FS; + y += 32; + + /* Render section */ + y += 14; /* label */ + + /* Colour map buttons */ + for (int i = 0; i < 6; i++) { + int bx = px + 12 + i * 42; + if (mx >= bx && mx <= bx+36 && my >= y+12 && my <= y+30) { + *cmap_delta = i - cfg->c_map; + cfg->c_map = i; + return SETTINGS_HIT_CMAP; + } + } + y += 36; + + /* Normalisation buttons */ + for (int i = 0; i < 3; i++) { + int bx = px + 12 + i * 60; + if (mx >= bx && mx <= bx+52 && my >= y+12 && my <= y+30) { + *nmap_delta = i - cfg->n_map; + cfg->n_map = i; + return SETTINGS_HIT_NMAP; + } + } + y += 36; + + /* Steps slider */ + int srx = px + SLIDER_X_OFF; + if (mx >= srx && mx <= srx + SLIDER_W && my >= y && my <= y+30) { + double t = (double)(mx - srx) / SLIDER_W; + if (t < 0) t = 0; + if (t > 1) t = 1; + cfg->steps = 1 + (int)(t * 63 + 0.5); + sp->drag_steps = 1; + return SETTINGS_HIT_STEPS; + } + y += 40; + + /* Circle mode buttons */ + for (int i = 0; i < CIRCLE_MODE_COUNT; i++) { + int bx = px + 12 + i * 56; + if (mx >= bx && mx <= bx+50 && my >= y+12 && my <= y+30) { + flags->circle_mode = (CircleMode)i; + return SETTINGS_HIT_NONE; + } + } + y += 36; + + /* HUD toggle button */ + if (mx >= px+12 && mx <= px+112 && my >= y && my <= y+22) { + flags->hud_visible = !flags->hud_visible; + return SETTINGS_HIT_NONE; + } + + return SETTINGS_HIT_NONE; +} + +/* ── Vector Field / Phase Trails ─────────────────────────────────────── */ +void ui_draw_field_lines(img_t *img, Simulation *sim, double x_range[2], double y_range[2]) { + int lines_x = 40; + int lines_y = 30; + + double step_size = (x_range[1] - x_range[0]) / img->width; + int max_steps = 150; + + for (int iy = 0; iy < lines_y; iy++) { + for (int ix = 0; ix < lines_x; ix++) { + double rx = x_range[0] + (x_range[1] - x_range[0]) * ((ix + 0.5) / lines_x); + double ry = y_range[0] + (y_range[1] - y_range[0]) * ((iy + 0.5) / lines_y); + double complex pos = rx + ry * I; + + for (int step = 0; step < max_steps; step++) { + double complex H = sim_eval_H(sim, pos); + double mag = cabs(H); + if (mag < 1e-6 || mag > 1e6) break; + + double complex dir = H / mag; + + double complex next_pos = pos + dir * step_size; + + uint32_t px, py; + complex_to_screen_choords(pos, &px, &py, x_range, y_range, img->width, img->height); + + if (px < img->width && py < img->height) { + uint32_t idx = py * img->width + px; + uint32_t bg = img->data[idx]; + uint32_t br = (bg >> 16) & 0xFF, bg_g = (bg >> 8) & 0xFF, bb = bg & 0xFF; + uint32_t a = 64; + uint32_t ia = 256 - a; + img->data[idx] = (0xFF << 24) | + (((255 * a + br * ia) >> 8) << 16) | + (((255 * a + bg_g * ia) >> 8) << 8) | + ((255 * a + bb * ia) >> 8); + } + + pos = next_pos; + } + } + } +} + +/* ── Domain Coloring ─────────────────────────────────────────────────── */ +void ui_draw_domain_overlay(img_t *img, Simulation *sim, + double x_range[2], double y_range[2]) { + double dx = (x_range[1] - x_range[0]) / img->width; + double dy = (y_range[1] - y_range[0]) / img->height; + + int step = 2; // evaluate every 2x2 pixels for performance + + for (int y = 0; y < img->height; y += step) { + double im = y_range[0] + y * dy; + for (int x = 0; x < img->width; x += step) { + double re = x_range[0] + x * dx; + double complex H = sim_eval_H(sim, re + im * I); + + double u = creal(H); + double v = cimag(H); + + double warp_scale = 5.0; + int grid_u = (int)floor(u * warp_scale); + int grid_v = (int)floor(v * warp_scale); + + if ((grid_u + grid_v) % 2 == 0) continue; + + uint32_t a = 32; + uint32_t ia = 256 - a; + for (int dy_bl = 0; dy_bl < step && y+dy_bl < img->height; dy_bl++) { + for (int dx_bl = 0; dx_bl < step && x+dx_bl < img->width; dx_bl++) { + uint32_t idx = (y+dy_bl) * img->width + (x+dx_bl); + uint32_t bg = img->data[idx]; + uint32_t br = (bg >> 16) & 0xFF, bg_g = (bg >> 8) & 0xFF, bb = bg & 0xFF; + img->data[idx] = (0xFF << 24) | + (((0 * a + br * ia) >> 8) << 16) | + (((0 * a + bg_g * ia) >> 8) << 8) | + ((0 * a + bb * ia) >> 8); + } + } + } + } +} diff --git a/src/undo.c b/src/undo.c new file mode 100644 index 0000000..586e305 --- /dev/null +++ b/src/undo.c @@ -0,0 +1,79 @@ +#include "undo.h" +#include +#include + +static void clone_array(singularity_array_t *dst, const singularity_array_t *src) { + dst->size = src->size; + dst->capacity = src->capacity ? src->capacity : 1; + dst->data = malloc(dst->capacity * sizeof(singularity_t)); + if (dst->data && src->data) + memcpy(dst->data, src->data, src->size * sizeof(singularity_t)); +} + +static void release_frame(UndoFrame *f) { + free(f->zeros.data); f->zeros.data = NULL; + free(f->poles.data); f->poles.data = NULL; +} + +static void restore_frame(const UndoFrame *f, + singularity_array_t *zeros, + singularity_array_t *poles) { + free(zeros->data); + free(poles->data); + clone_array(zeros, &f->zeros); + clone_array(poles, &f->poles); +} + +void undo_init(UndoStack *u) { + memset(u, 0, sizeof(*u)); + u->pos = -1; + u->total = 0; +} + +void undo_push(UndoStack *u, + const singularity_array_t *zeros, + const singularity_array_t *poles) { + /* Discard redo history above current position */ + for (int i = u->pos + 1; i < u->total; i++) + release_frame(&u->frames[i]); + u->total = u->pos + 1; + + if (u->total >= UNDO_MAX) { + /* Circular shift: drop oldest, slide everything left */ + release_frame(&u->frames[0]); + memmove(&u->frames[0], &u->frames[1], + (UNDO_MAX - 1) * sizeof(UndoFrame)); + u->total--; + u->pos--; + } + + u->pos++; + clone_array(&u->frames[u->pos].zeros, zeros); + clone_array(&u->frames[u->pos].poles, poles); + u->total = u->pos + 1; +} + +int undo_undo(UndoStack *u, + singularity_array_t *zeros, + singularity_array_t *poles) { + if (u->pos <= 0) return 0; + u->pos--; + restore_frame(&u->frames[u->pos], zeros, poles); + return 1; +} + +int undo_redo(UndoStack *u, + singularity_array_t *zeros, + singularity_array_t *poles) { + if (u->pos >= u->total - 1) return 0; + u->pos++; + restore_frame(&u->frames[u->pos], zeros, poles); + return 1; +} + +void undo_destroy(UndoStack *u) { + for (int i = 0; i < u->total; i++) + release_frame(&u->frames[i]); + u->pos = -1; + u->total = 0; +} diff --git a/src/utils.c b/src/utils.c new file mode 100644 index 0000000..f4667fc --- /dev/null +++ b/src/utils.c @@ -0,0 +1,47 @@ +#include "utils.h" + +int gcd(int a, int b) +{ + while (b != 0) { + int t = b; + b = a % b; + a = t; + } + return a; +} + +double random_uniform(double min, double max) +{ + return min + (rand() / (RAND_MAX + 1.0)) * (max - min); +} + +void screen_choords_to_complex(double x, double y, double complex *c, + double x_range[2], double y_range[2], + int width, int height) +{ + double re = x_range[0] + (x / (double)(width - 1)) * (x_range[1] - x_range[0]); + double im = y_range[0] + (y / (double)(height - 1)) * (y_range[1] - y_range[0]); + *c = re + im * I; +} + +void complex_to_screen_choords(double complex c, uint32_t *x, uint32_t *y, + double x_range[2], double y_range[2], + int width, int height) +{ + double re = creal(c); + double im = cimag(c); + + double dx = (re - x_range[0]) / (x_range[1] - x_range[0]); + double dy = (im - y_range[0]) / (y_range[1] - y_range[0]); + + int ix = (int)(dx * (width - 1) + 0.5); + int iy = (int)(dy * (height - 1) + 0.5); + + if (ix < 0) ix = 0; + if (ix >= width) ix = width - 1; + if (iy < 0) iy = 0; + if (iy >= height) iy = height - 1; + + *x = ix; + *y = iy; +} \ No newline at end of file