// these are shaping functions taken from various places of https://flong.com

// Polynominal shaping

float blinn_wyvill_cos_approx(float x) {
  float x2 = x*x;
  float x4 = x2*x2;
  float x6 = x4*x2;

  float fa = ( 4.0/9.0);
  float fb = (17.0/9.0);
  float fc = (22.0/9.0);

  return fa*x6 - fb*x4 + fc*x2;
}

// TODO: factor out clamping
float double_cubic_seat(float x, vec2 ab) {
  float epsilon = 0.00001;
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  float a = clamp(ab.x, min_param_a, max_param_a);
  float b = clamp(ab.y, min_param_b, max_param_b);

  if (x <= a){
    return b - b*pow(1.-x/a, 3.0);
  }
  return b + (1.-b)*pow((x-a)/(1.-a), 3.0);
}

float double_cubic_seat_linear_blend(float x, vec2 ab) {
  float epsilon = 0.00001;
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  float a = clamp(ab.x, min_param_a, max_param_a);
  float b = clamp(ab.y, min_param_b, max_param_b);
  b = 1.0 - b;

  if (x<=a){
    return b*x + (1.-b)*a*(1.-pow(1.-x/a, 3.));
  }
  return b*x + (1.-b)*(a + (1.-a)*pow((x-a)/(1.-a), 3.));
}

float double_odd_polynomial_seat(float x, float a, float b, float n) {
  float epsilon = 0.00001;
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0; 
  a = clamp(a, min_param_a, max_param_a);
  b = clamp(b, min_param_b, max_param_b);

  float p = 2.*n + 1.;
  if (x <= a){
    return b - b*pow(1.-x/a, p);
  }
  return b + (1.-b)*pow((x-a)/(1.-a), p);
}

float double_poly_sigmoid(float x, float a, float b, float n) {
  if(mod(n, 2.) == 0.) {
    if(x <= 0.5) {
      return pow(2.0*x, n)/2.0;
    }
    return 1.0 - pow(2.*(x-1.), n)/2.0;
  } else {
    if (x<=0.5){
      return pow(2.0*x, n)/2.0;
    }
    return 1.0 + pow(2.0*(x-1.), n)/2.0;
  }
}

float quad_through_point(float x, vec2 ab) {
  float epsilon = 0.00001;
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  float a = clamp(ab.x, min_param_a, max_param_a);
  float b = clamp(ab.y, min_param_b, max_param_b);

  float A = (1.-b)/(1.-a) - (b/a);
  float B = (A*(a*a)-b)/a;
  float y = A*(x*x) - B*(x);
  y = clamp(0., 1., y);

  return y;
}

// exponential shaping

float exp_easing(float x, float a) {
  float e = 0.00001;
  a = clamp(a, 0.+e, 1.-e);

  if (a < 0.5) {
    a = 2.0*(a);
    return pow(x, a);
  }
  a = 2.*(a-.5);
  return pow(x, 1./(1.-a));
}

float double_exp_seat(float x, float a) {
  float e = 0.00001;
  a = clamp(a, .0+e, 1.-e);

  if (x<=0.5) {
    return pow(2.*x, 1.-a)/2.0;
  }
  return 1. - pow(2.*(1.-x), 1.-a)/2.;
}

float double_exp_sigmoid(float x, float a) {
  float e = 0.00001;
  a = 1.-clamp(a, .0+e, 1.-e);

  if (x<=0.5){
    return pow(2.0*x, 1.0/a)/2.0;
  }
  return 1.0 - (pow(2.0*(1.0-x), 1.0/a))/2.0;
}

float normalized_log_sigmoid(float x, float a) {
  float e = 0.0001;
  a = 1./(1.-clamp(a, 0.+e, 1.-e)) - 1.;

  float b = 1.0 / (1.0 + exp(a));
  float c = 1.0 / (1.0 + exp(0.-a));
  a = 1. / (1. + exp(0. -((x-.5)*a*2.)));
  return (a-b)/(c-b);
}