/* See intg_problem.h */   
/* Last edited on 2004-08-17 23:11:47 by stolfi */


#include <intg_gen.h>
#include <intg_problem.h>

#include <js.h>
#include <vec.h>
#include <math.h>

/* INTERNAL PROTOTYPES */

void Quadr_Sol(Time t, State *s);
bool Quadr_RHS(Time t, State *s, Velocity *v);

void Quart_Sol(Time t, State *s);
bool Quart_RHS(Time t, State *s, Velocity *v);

void Quint_Sol(Time t, State *s);
bool Quint_RHS(Time t, State *s, Velocity *v);

void Sinus_Sol(Time t, State *s);
bool Sinus_RHS(Time t, State *s, Velocity *v);

void Twang_Sol(Time t, State *s);
bool Twang_RHS(Time t, State *s, Velocity *v);

void Gauss_Sol(Time t, State *s);
bool Gauss_RHS(Time t, State *s, Velocity *v);

void Hyper_Sol(Time t, State *s);
bool Hyper_RHS(Time t, State *s, Velocity *v);

void Notch_Sol(Time t, State *s);
bool Notch_RHS(Time t, State *s, Velocity *v);

/* PACK THEM ALL: */

Intg_Problem_vec_t Intg_Problem_Sample(void)
  {
    Intg_Problem_vec_t p = Intg_Problem_vec_new(0);
    int np = 0;
    
    auto void mkproblem(Intg_Solution *sol, Intg_RHS *rhs, char *tag, int n);
    void mkproblem(Intg_Solution *sol, Intg_RHS *rhs, char *tag, int n)
      {
        Intg_Problem_vec_expand(&p, np);
        Intg_Problem *q = &(p.el[np]);
        q->tag = tag;
        q->n = n;
        q->rhs = rhs;
        q->sol = sol;
        q->t0 = -1.0;
        q->t1 = +1.0;
        q->dt = 0.1;
        q->dtMin = 1.0e-6;
        q->dtMax = 1.0;
        q->sol = sol;
        np++;
      }
    
    mkproblem(&Quadr_Sol, &Quadr_RHS, "quadr", 2);
    mkproblem(&Quart_Sol, &Quart_RHS, "quart", 2);
    mkproblem(&Quint_Sol, &Quint_RHS, "quint", 2);
    mkproblem(&Sinus_Sol, &Sinus_RHS, "sinus", 2);
    mkproblem(&Gauss_Sol, &Gauss_RHS, "gauss", 2);
    mkproblem(&Twang_Sol, &Twang_RHS, "twang", 2);
    mkproblem(&Hyper_Sol, &Hyper_RHS, "hyper", 2);
    mkproblem(&Notch_Sol, &Notch_RHS, "notch", 1);
    
    Intg_Problem_vec_trim(&p, np);
    return p;
  }

Intg_Problem_vec_t Intg_Problem_vec_new(nat nel)
  { /* This is not a macro only because gcc does not allow cast of struct: */
    vec_t v = vec_new(nel, sizeof(Intg_Problem));
    Intg_Problem_vec_t r;
    r.nel = v.nel; r.el = (Intg_Problem *)v.el;
    return r;
  }

/* THE PROBLEMS */

/* Quadratic */

void Quadr_Sol(Time t, State *s)
  { s->el[0] = t*t;
    s->el[1] = 2.0 * t;
  }

bool Quadr_RHS(Time t, State *s, Velocity *v)
  { v->el[0] = s->el[1];
    v->el[1] = 2.0;
    return FALSE;
  }

/* Quartic */

void Quart_Sol(Time t, State *s)
  { s->el[0] = t*(t*t*t - 1.0);
    s->el[1] = 4.0 * t*t*t - 1.0;
  }

bool Quart_RHS(Time t, State *s, Velocity *v)
  { v->el[0] = s->el[1];
    v->el[1] = 12.0*t*t;
    return FALSE;
  }

/* Quintic */

void Quint_Sol(Time t, State *s)
  { s->el[0] = t*(t*t*t*t - 1.0);
    s->el[1] = 5.0 * t*t*t*t - 1.0;
  }

bool Quint_RHS(Time t, State *s, Velocity *v)
  { v->el[0] = s->el[1];
    v->el[1] = 20.0*t*t*t;
    return FALSE;
  }

/* Sinusoid */

#define SinusW (20.0)

void Sinus_Sol(Time t, State *s)
  { double W = SinusW;
    double ct = cos(W*t);
    double st = sin(W*t);
    s->el[0] = ct;
    s->el[1] = -W*st;
  }

bool Sinus_RHS(Time t, State *s, Velocity *v)
  { double W = SinusW;
    v->el[0] = s->el[1];
    v->el[1] = - W*W*s->el[0];
    return FALSE;
  }

/* Decaying sinusoid */

#define TwangW (20.0)
#define TwangK (2.0)

void Twang_Sol(Time t, State *s)
  { double W = TwangW;
    double K = TwangK;
    double et = exp(-K*t);
    double ct = et * cos(W*t);
    double st = et * sin(W*t);
    s->el[0] = ct;
    s->el[1] = st;
  }

bool Twang_RHS(Time t, State *s, Velocity *v)
  { double W = TwangW;
    double K = TwangK;
    v->el[0] = - K * s->el[0] - W * s->el[1];
    v->el[1] = - K * s->el[1] + W * s->el[0];
    return FALSE;
  }

/* Gaussian bell */

#define GaussK (3.0)

void Gauss_Sol(Time t, State *s)
  { double K = GaussK;
    double kt = K*t;
    double et = exp(-kt*kt);
    s->el[0] = et;
    s->el[1] = -2.0*K*K*t*et;
  }

bool Gauss_RHS(Time t, State *s, Velocity *v)
  { double K = GaussK;
    double r = s->el[1]/s->el[0];
    v->el[0] = s->el[1];
    v->el[1] = s->el[0]*(r*r - 2.0*K*K);
    return FALSE;
  }

/* Hyperbola */

#define HyperE (0.05)

void Hyper_Sol(Time t, State *s)
  { double E = HyperE;
    double et = sqrt(t*t + E*E);
    s->el[0] = et;
    s->el[1] = t/et;
  }

bool Hyper_RHS(Time t, State *s, Velocity *v)
  { double E = HyperE;
    double r = 1.0/s->el[0];
    v->el[0] = s->el[1];
    v->el[1] = E*E * r*r*r;
    return FALSE;
  }

/* Notches: log(2 + sin(t)) */

#define NotchK (1.41421356237309504880)
#define NotchW (4.0)

void Notch_Sol(Time t, State *s)
  { double K = NotchK;
    double W = NotchW;
    double et = log(K + sin(W*t));
    s->el[0] = et;
  }

bool Notch_RHS(Time t, State *s, Velocity *v)
  { double W = NotchW;
    v->el[0] = W*cos(W*t)/exp(s->el[0]);
    return FALSE;
  }