// feigenbaum.cpp -- attempts to calculate the feigenbaum number based on the behavior of the logistic population growth map, using the supercycle ratio method

#include <iostream.h>
#include "apvector.h"

const int MAXNUM = 5000;
const int INITNTH = 5000;

double f(double A, double x) {  // logistic map
  return A*x*(1-x);
}

double nthf(double A, double x, int n) { // gives the nth f(x)
  if (n == 0) return x;
  return nthf(A, A*x*(1-x), n-1);
}

int appeq(double a, double b) { // approximately equal for same-sign #s
  double diff = a - b;
  if (diff < 0) diff *= -1;
  if (diff < 0.000001) return 1;
  else return 0;
}

int pernum(double a) {
  double x = 0.1, xinit;
  apvector<double> perx(0);
  int size = 0, i;
  x = nthf(a, 0.1, INITNTH);
  xinit = x;
  for (i = 0; i < MAXNUM; i++) {
    perx.resize(size+1);
    size++;
    perx[size-1] = x;
    x = f(a, x);
    if (appeq(x, xinit)) goto don;
  }
 don:
  return size;
}

int show_perx() {
  double x = 0.1, xinit;
  apvector<double> perx(0);
  int size = 0, i;
  double A;

  cout << endl << "What constant A would you like to use? ";
  cin >> A;
  cout << endl;

  x = nthf(A, 0.1, INITNTH);
  xinit = x;
  cout << "Initial x = " << xinit << endl;
  for (i = 0; i < MAXNUM; i++) {
    perx.resize(size+1);
    size++;
    perx[size-1] = x;
    x = f(A, x);
    if (appeq(x, xinit)) goto don;
  }

 don:

  cout << endl << "Periodic values of x (" << size << "):" << endl;
  for (i = 0; i < size; i++) cout << perx[i] << " ";

  return 0;
}

apvector<double> calc_feigenbaum() { // returns apvector with 3 values, each an approximation of Feigenbaum's number
  const double inc = 0.001;
  double i = 2.5;
  int tmp = 0, tmpnext = pernum(i);
  apvector<double> fei(3, 0.0);
  double num12 = -1.0, num24 = -1.0, num48 = -1.0, num816 = -1.0, num1632 = -1.0;
  while (num12 == -1.0 || num24 == -1.0 || num48 == -1.0 || num816 == -1.0 || num1632 == -1.0) {
    tmp = tmpnext;
    tmpnext = pernum(i+inc);
    if (tmp == 1 && tmpnext == 2) num12 = i;
    if (tmp == 2 && tmpnext == 4) num24 = i;
    if (tmp == 4 && tmpnext == 8) num48 = i;
    if (tmp == 8 && tmpnext == 16) num816 = i;
    if (tmp == 16 && tmpnext == 32) num1632 = i;
    if (tmp == tmpnext == MAXNUM) {
      cout << "Error finding Feigenbaum number!";
      break;
    }
    i += inc;
    if (i > 4.5) {
      cout << "A > 4.5 when calculating Feigenbaum delta, Error." << endl;
      break;
    }
  }
  cout << "1 -> 2 Bifurcation: A = " << num12 << endl;
  cout << "2 -> 4 Bifurcation: A = " << num24 << endl;
  cout << "4 -> 8 Bifurcation: A = " << num48 << endl;
  cout << "8 -> 16 Bifurcation: A = " << num816 << endl;
  cout << "16 -> 32 Bifurcation: A = " << num1632 << endl;
  fei[0] = (num12-num24)/(num24-num48);
  fei[1] = (num24-num48)/(num48-num816);
  fei[2] = (num48-num816)/(num816-num1632);
  return fei;
}

int main() {
  int sel, n, end = 0; 
  double x0, A1, A2, inc, i, tmp;
  apvector<double> feig(3);
  while (end == 0) {
    sel = n = 0;
    x0 = A1 = A2 = inc = i = tmp = 0.0;
    cout << endl << endl;
    cout << "-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-" << endl;
    cout << "Logistic Map Simulator" << endl << "by Alex Ellis" << endl << endl;
    cout << "x[n+1] = A * x[n] * (1-x[n])" << endl;
    cout << "-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-" << endl << endl;
    cout << "Make your selection:" << endl;
    cout << "1. Find the nth x for given A and x[0]" << endl;;
    cout << "2. See period-n behavior for various values of A" << endl;
    cout << "3. Calculate the Feigenbaum delta" << endl;
    cout << "4. Exit" << endl;
    cout << "?";
    cin >> sel;

    switch (sel) {
    case 1:
      cout << endl << endl << "n = ?";
      cin >> n;
      cout << endl << "A = ?";
      cin >> A1;
      cout << endl << "x[0] = ?";
      cin >> x0;
      cout << endl << endl << "f(x[n]) = " << nthf(A1, x0, n) << endl << endl;
      break;
    case 2:
      cout << endl << endl << "Starting A = ?";
      cin >> A1;
      cout << endl << "Ending A = ?";
      cin >> A2;
      cout << endl << "A-Incrementation = ?";
      cin >> inc;
      cout << endl;
      for (i = A1; i <= A2; i += inc) {
	tmp = pernum(i);
	cout << "A = " << i << ", " << tmp << " periods";
	if (tmp == MAXNUM) cout << " (max detectable)";
	cout << endl;
      }
      cout << endl;
      break;
    case 3:
      feig = calc_feigenbaum();
      cout << "Try #1: " << feig[0] << endl;
      cout << "Try #2: " << feig[1] << endl;
      cout << "Try #3: " << feig[2] << endl;
      cout << endl;
      break;
    case 4:
      end = 1;
      break;
    default:
      cout << endl;
      break;
    }
  }
  cout << endl;
}