#include <math.h>
#include <stdlib.h>   // For calloc

// add functions missing in math.h???
double fmax(double a, double b) { return a > b ? a : b; }
double fmin(double a, double b) { return a < b ? a : b; }

// This file contains the C++ code for a QSPICE C-block that implements
// Space Vector PWM (SVPWM). The block takes a modulation index and angle
// as inputs and outputs the six PWM signals (high and low side for each phase)
// with dead time.

// QSPICE C-block specific union for data exchange
// This definition is crucial for the 'data' argument in the svpwm_logic
// function.
union uData {
  double d;
  int    i;
  char  *s;
};

// Constants
#define M_PI  3.14159265358979323846
#define SQRT3 1.7320508075688772

// This structure holds per-instance data that persists between calls.
// We use this to store the previous simulation time, which is necessary
// for calculating the elapsed time (Ts) and checking if a new PWM cycle
// has started.
struct sSVPWM_X1 {
  double last_t;
  double T_s;
};

// Function prototypes
double max3(double a, double b, double c);
double min3(double a, double b, double c);

// Main C-block function
// extern "C" __declspec(dllexport) void svpwm_logic(void **opaque, double t,
// union uData *data)
//
// changed to match expected evaluation function name -- untested beyond that
//
extern "C" __declspec(dllexport) void spwn_x5(
    void **opaque, double t, union uData *data) {

  // Initialize per-instance data on the first call
  if (*opaque == 0) {
    *opaque                = calloc(1, sizeof(struct sSVPWM_X1));
    struct sSVPWM_X1 *inst = (struct sSVPWM_X1 *)*opaque;
    inst->last_t           = t;
  }
  struct sSVPWM_X1 *inst = (struct sSVPWM_X1 *)*opaque;

  // Inputs from the schematic
  double mod_index = data[0].d;   // Modulation index (0 to 1)
  double angle     = data[1].d;   // Angle of the reference vector (radians)
  double freq_sw   = data[2].d;   // Switching frequency (Hz)
  double dead_time = data[3].d;   // Dead time (seconds)

  // Outputs to the schematic
  double &pwm_a_high = data[4].d;
  double &pwm_a_low  = data[5].d;
  double &pwm_b_high = data[6].d;
  double &pwm_b_low  = data[7].d;
  double &pwm_c_high = data[8].d;
  double &pwm_c_low  = data[9].d;

  // The simulation time step can be variable. We need to define a fixed
  // switching period for the PWM calculation.
  inst->T_s = 1.0 / freq_sw;

  // Calculate the normalized angle for a 360-degree range
  double theta = fmod(angle, 2.0 * M_PI);
  if (theta < 0) { theta += 2.0 * M_PI; }

  // Scale modulation index
  double modulation_index = mod_index * 2.0 / SQRT3;

  // --- SVPWM Core Logic ---
  double T_a, T_b, T_c;   // Dwell times (normalized to T_s)

  // Get the voltages in the alpha-beta frame
  double v_alpha = modulation_index * cos(theta);
  double v_beta  = modulation_index * sin(theta);

  // Clarke transformation inverse to get phase voltages (Va, Vb, Vc)
  double v_a_ref = v_alpha;
  double v_b_ref = -0.5 * v_alpha + SQRT3 / 2.0 * v_beta;
  double v_c_ref = -0.5 * v_alpha - SQRT3 / 2.0 * v_beta;

  // Find the min and max reference voltages
  double v_max = max3(v_a_ref, v_b_ref, v_c_ref);
  double v_min = min3(v_a_ref, v_b_ref, v_c_ref);

  // Calculate the common-mode offset voltage
  double v_offset = -0.5 * (v_max + v_min);

  // Apply common-mode offset to center the modulation
  double v_a_svpwm = v_a_ref + v_offset;
  double v_b_svpwm = v_b_ref + v_offset;
  double v_c_svpwm = v_c_ref + v_offset;   // Corrected typo here previously

  // Calculate the duty cycles for each phase leg
  double Ta_on = (v_a_svpwm + 1.0) / 2.0;
  double Tb_on = (v_b_svpwm + 1.0) / 2.0;
  double Tc_on = (v_c_svpwm + 1.0) / 2.0;

  // --- PWM Signal Generation with Dead Time ---
  double t_local = fmod(t, inst->T_s);
  double t_carrier =
      2.0 * t_local / inst->T_s;   // Triangle carrier, 0 to 2, centered at 1
  if (t_carrier > 1.0) {
    t_carrier = 2.0 - t_carrier;   // Make it a triangle wave
  }

  double t_dead_norm = dead_time / inst->T_s;

  // Phase A
  if (t_carrier > Ta_on + t_dead_norm) {
    pwm_a_high = 1.0;
    pwm_a_low  = 0.0;
  } else if (t_carrier < Ta_on - t_dead_norm) {
    pwm_a_high = 0.0;
    pwm_a_low  = 1.0;
  } else {
    pwm_a_high = 0.0;
    pwm_a_low  = 0.0;   // Dead time
  }

  // Phase B
  if (t_carrier > Tb_on + t_dead_norm) {
    pwm_b_high = 1.0;
    pwm_b_low  = 0.0;
  } else if (t_carrier < Tb_on - t_dead_norm) {
    pwm_b_high = 0.0;
    pwm_b_low  = 1.0;
  } else {
    pwm_b_high = 0.0;
    pwm_b_low  = 0.0;   // Dead time
  }

  // Phase C
  if (t_carrier > Tc_on + t_dead_norm) {
    pwm_c_high = 1.0;
    pwm_c_low  = 0.0;
  } else if (t_carrier < Tc_on - t_dead_norm) {
    pwm_c_high = 0.0;
    pwm_c_low  = 1.0;
  } else {
    pwm_c_high = 0.0;
    pwm_c_low  = 0.0;   // Dead time
  }
}

// Helper function to find the maximum of three doubles
double max3(double a, double b, double c) { return fmax(fmax(a, b), c); }

// Helper function to find the minimum of three doubles
double min3(double a, double b, double c) { return fmin(fmin(a, b), c); }