#include "libwasm.h"
#include "drawings.h"
#include "complex.h"


static constexpr double sincos (const double x, const bool even) {
  double result (0.0), element(1.0), divider(0.0);
  if (even) { element *= x; divider += 1.0; }
  constexpr double eps = 1.0e-9;      // maximální chyba výpočtu
  const double aa = - (x * x);
  for (;;) {
    result  += element;
    if (fabs  (element) < eps) break;
    divider += 1.0;
    double fact = divider;
    divider += 1.0;
    fact    *= divider;
    element *= aa / fact;
  }
  return result;
}

/*********************************************************************************/
void BackGround::drawings() {
  const real maxx = getMaxX(),  maxy = getMaxY();
  fill (Color(0, 0, 0, 0xFF));
  
  const real hx = 0.5 * maxx, margin = 5.0;
  const Matrix m  (hx-2.0*margin, 0.0, 0.0, -(maxy-2*margin), hx, maxy - margin);
  setMatrix (m);
  setColor (Color (0x60, 0x60, 0x60));
  const real st = 0.1;
  for (real x=-1.0; x<+1.0; x+=st) {
    line (FPoint(+x,-1.0), FPoint(+x,+1.0), true);
  }
  for (real y=st; y<1.0; y+=st) {
    line (FPoint(-1.0,+y), FPoint(+1.0,+y), true);
  }
  setColor (Color (0x80, 0x80, 0));
  line (FPoint(-1.0, 0.0), FPoint(1.0, 0.0));
  line (FPoint( 0.0, 0.0), FPoint(0.0, 1.0));
  FPoint o (-1.0, 1.0);
  setColor (Color (0xFF, 0xFF, 0));
  /*
  for (real x=-1.0; x<+1.0; x+=0.05) {
    const real y = x * x;
    FPoint n (x, y);
    line (o, n);
    o = n;
  }
  */
  line (FPoint (-1.0, 0.0), FPoint (+1.0, 0.0));
  line (FPoint (-1.0, 0.0), FPoint (-1.0, 1.0));
  line (FPoint (+1.0, 0.0), FPoint (+1.0, 1.0));
}
///////////////////////////////////////////////////////////////////////////////////
void ForeGround::drawings() {
  const real maxx = getMaxX(),  maxy = getMaxY();
  fill (Color(0xFF, 0xFF, 0xFF, 0));
  const real hx = 0.5 * maxx, margin = 5.0, zx = hx - 2.0 * margin, ix = 1.0 / zx;
  const Matrix m (zx, 0.0, 0.0, -(maxy - 2 * margin), hx, maxy - margin);
  setMatrix (m);
  for (real a=-1.0; a<+1.0; a+=ix) {
    const double  y = probality (a);
    const double cz = y * 255.0;
    int c = 0;
    if      (cz > 255.0) c = 0xFF;
    else if (cz < 0.0)   c = 0;
    else                 c = cz;
    setColor (Color(c, 0xFF -c, 0xFF - c, 0xE0));
    line (FPoint (a, 0.0), FPoint (a, y));
  }
}
/**
 * Výpočty vlnových funkcí viz
 * https://en.wikipedia.org/wiki/Quantum_well,
 * časový vývoj https://www.aldebaran.cz/studium/tf.pdf str. 182, výraz (2.183)
 * */
static const unsigned EN [NumWeight] = {
//1,  4,  9, 16, 25, 36, 49, 64, 81, 100,   // energie roste kvadraticky, ale je to moc divoké
  1,  2,  3,  4,  5,  7,  9, 11, 13,  16,   // tohle sice není úplně správně, ale zase se na to dá koukat
};
double ForeGround::probality (const double x) const {
  const double z = 0.5 * (x + 1.0);                               // jáma je od -1 do 1
  double fz [NumWeight];
  for (unsigned n=0; n<NumWeight; n++) {
    const double a = (double) (n+1) * M_PI * z;
    fz [n] = sincos (a, true);
  }
  complex psi;
  for (unsigned n=0; n<NumWeight; n++) {
    complex y;
    y.exp (index * EN[n]);                                        // otočení komplexní jednotky - proměnná index představuje čas
    const double fac = weights [n] * fz [n];                      // |n>
    y   *= fac;                                                   // |Ψ>
    psi += y;                                                     // sčítají se komplexní vlnové funkce
  }
  return psi.abs();                                               // vrací se <Ψ*|Ψ>
}
void ForeGround::step() {
  drawings();
  phase();
  index += speed;
}
// jen čárka ukazující fázi základního stavu
void ForeGround::phase() {
  complex c, o (-0.8, 0.9);
  c.exp (index);
  c *= 0.1;
  c += o;
  const FPoint a (o.re, o.im), b (c.re, c.im);
  setColor (Color(0xFF, 0x80, 0, 0xFF));
  line (a, b);
}
void ForeGround::normalize() {
  index = 0; speed = 1u;
  for (unsigned n=0; n<NumWeight; n++) {
    weights [n] = 0.0;
  }
}