#include <iostream>
#include <fstream>
#include <cmath>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define DEGTORAD(x) ((x) * 0.0174532925)
#define PI 3.14159

using namespace std;

using std::cout;
using std::cin;
using std::endl;
using std::ofstream;



///////////////////////// INITIALISATION \\\\\\\\\\\\\\\\\\\\\\\\




int a = 0;
long i;
long j;
int l = 0;

// Create Target Data File

ofstream myfile ("example7.csv", ios::out | ios::app | ios::binary);


// Initialise Fractal Properties

float reduction_factor;
float compound_reduction_factor;
int divergence_angle;
int generations;
int k = 1000;


// Initialise Line Properties

struct line{
double length;
double radius;
int colour;
int generation;
long identifier;
double start_coord[3];
double end_coord[3];
double vector[4];
double reference_geometry[4];
double x;
bool Child1;
bool Child2;
double angle1;
double angle2;
}DataList[20002];

//Function Prototypes

line fractal_branch1 (line DataList1);
line fractal_branch2 (line DataList2);
line fractal_branch3 (line DataList3);




///////////////////////// FUNCTION MAIN \\\\\\\\\\\\\\\\\\\\\\\\




int main()
{

// Generate Initial Line


DataList[0].length = 10.0;
DataList[0].radius = 10;
DataList[0].colour = 255;
DataList[0].generation = 0;
DataList[0].identifier = 0;
DataList[0].start_coord[0] = 0;
DataList[0].start_coord[1] = 0;
DataList[0].start_coord[2] = 0;
DataList[0].end_coord[0] = 0;
DataList[0].end_coord[1] = 5.0;
DataList[0].end_coord[2] = 0;
DataList[0].vector[0] = 0;
DataList[0].vector[1] = 5.0;
DataList[0].vector[2] = 0;
DataList[0].vector[3] = 1;
DataList[0].x = 5;
DataList[0].reference_geometry[0] = 0;
DataList[0].reference_geometry[1] = 0;
DataList[0].reference_geometry[2] = 1;
DataList[0].reference_geometry[2] = 1;
DataList[0].Child1 = false;
DataList[0].Child2 = false;
DataList[0].angle1 = 25;
DataList[0].angle2 = 70;



// Initialise pointers


i = 0;
j = 1;
srand ( time(NULL) );

//cout<< "Enter number of generations"<<endl<<endl;
//cin>> generations;

//k = (4*(2^generations));


// USE IF STATEMENT TO CALL FUNCTION CHILD1 OR CHILD2

// has parent already been used for child1?
//if no then call Child1 function, specify that parent has been used for Child1 and increment exchange2 by one
//if yes then ask has parent already been used for child 2?
//if no then call Child2 function, specify that parent has been used for Child2 and increment exchange2 by one
//if yes then increment exchange1 by one (so that it is assigned to Datalist[1])


cout<< "check 1"<<endl;

do {



if ((DataList[i].Child1 == false) && (DataList[i].Child2 == false))
{

DataList[i].angle1 = ((rand() % 100)+30)/3;
DataList[i].angle2 = (rand() % 10)+30;

DataList[j] = fractal_branch1 (DataList[i]);
DataList[j].identifier = j;
DataList[i].Child1 = true;

//cout<<j<<" child 1"<<endl;

j++;

DataList[j] = fractal_branch3 (DataList[j-1]);
DataList[j].identifier = j;
j++;

DataList[j] = fractal_branch3 (DataList[j-1]);
DataList[j].identifier = j;
j++;

DataList[j] = fractal_branch3 (DataList[j-1]);
DataList[j].identifier = j;
j++;


}

else if ((DataList[i].Child1 == true) && (DataList[i].Child2 == false))
{

DataList[i].angle1 = -((rand() % 100)+30)/3;
DataList[i].angle2 = (rand() % 50)+30;

DataList[j] = fractal_branch2 (DataList[i]);
DataList[j].identifier = j;
DataList[i].Child2 = true;


j++;

DataList[j] = fractal_branch3 (DataList[j-1]);
DataList[j].identifier = j;
j++;

DataList[j] = fractal_branch3 (DataList[j-1]);
DataList[j].identifier = j;
j++;

DataList[j] = fractal_branch3 (DataList[j-1]);
DataList[j].identifier = j;
j++;

}

else if ((DataList[i].Child1 == true) && (DataList[i].Child2 == true))
{
i+=4;
//i++;


}

}
while (i <= k);

// // // END OF TEST FUNCTION \\ \\ \\

cout<< "check 2"<<endl;


do
{
cout<<DataList[l].identifier<<"\t"<<DataList[l].x<<"\t"<<DataList[l].radius<<endl<<endl;

l++;
}

while (l <= 50);



cout<< "check 3"<<endl;


// Write data to .csv file


do{

myfile <<DataList[a].identifier<<", "<<DataList[a].x<<", "<<DataList[a].radius<<endl;

a++;

}

while ( a < 2002);


myfile.close();


// End main

return 0;

}





///////////////////////// FUNCTION DEFINITIONS \\\\\\\\\\\\\\\\\\\\\\\\







line fractal_branch1 (line DataList1)
{

// Define Newline

typedef struct line newline;

newline exchange1;
newline exchange2;
exchange1 = DataList1;
double vector1[4][1] = {exchange1.vector[0],
exchange1.vector[1],
exchange1.vector[2],
exchange1.vector[3]};

double vector2[4][1];
double vector3[4][1];
double vector4[4][1] = {exchange1.reference_geometry[0],
exchange1.reference_geometry[1],
exchange1.reference_geometry[2],
exchange1.reference_geometry[3]};

double theta1 = exchange1.angle1;
double theta2 = exchange1.angle2;
double c, s, t, x, y, z, magnitude;


// Initialise Start Coordinates

exchange2.start_coord[0] = exchange1.end_coord[0];
exchange2.start_coord[1] = exchange1.end_coord[1];
exchange2.start_coord[2] = exchange1.end_coord[2];


// Step 1: Rotate vector1 (line vector) around reference geometry to generate vector2

c = cos(DEGTORAD(theta1));
s = sin(DEGTORAD(theta1));
t = (1-cos(DEGTORAD(theta1)));
x = exchange1.reference_geometry[0]/(sqrt((exchange1.reference_geometry[0]*exchange1.reference_geometry[0])+(exchange1.reference_geometry[1]*exchange1.reference_geometry[1])+(exchange1.reference_geometry[2]*exchange1.reference_geometry[2])));
y = exchange1.reference_geometry[1]/(sqrt((exchange1.reference_geometry[0]*exchange1.reference_geometry[0])+(exchange1.reference_geometry[1]*exchange1.reference_geometry[1])+(exchange1.reference_geometry[2]*exchange1.reference_geometry[2])));
z = exchange1.reference_geometry[2]/(sqrt((exchange1.reference_geometry[0]*exchange1.reference_geometry[0])+(exchange1.reference_geometry[1]*exchange1.reference_geometry[1])+(exchange1.reference_geometry[2]*exchange1.reference_geometry[2])));

double rotation_matrix1[4][4] = {((t*x*x)+c), ((t*x*y)+(s*z)), ((t*x*z)-(s*y)), 0,
((t*x*y)-(s*z)), ((t*y*y)+c), ((t*y*z)+(s*x)), 0,
((t*x*z)+(s*y)), ((t*y*z)-(s*x)), ((t*z*z)+c), 0,
0, 0, 0, 1};


vector2[0][0] = (rotation_matrix1[0][0]*vector1[0][0])+(rotation_matrix1[0][1]*vector1[1][0])+(rotation_matrix1[0][2]*vector1[2][0])+(rotation_matrix1[0][3]*vector1[3][0]);
vector2[1][0] = (rotation_matrix1[1][0]*vector1[0][0])+(rotation_matrix1[1][1]*vector1[1][0])+(rotation_matrix1[1][2]*vector1[2][0])+(rotation_matrix1[1][3]*vector1[3][0]);
vector2[2][0] = (rotation_matrix1[2][0]*vector1[0][0])+(rotation_matrix1[2][1]*vector1[1][0])+(rotation_matrix1[2][2]*vector1[2][0])+(rotation_matrix1[2][3]*vector1[3][0]);
vector2[3][0] = (rotation_matrix1[3][0]*vector1[0][0])+(rotation_matrix1[3][1]*vector1[1][0])+(rotation_matrix1[3][2]*vector1[2][0])+(rotation_matrix1[3][3]*vector1[3][0]);



// Step 2: Rotate vector2 around line to generate vector3

c = cos(DEGTORAD(theta2));
s = sin(DEGTORAD(theta2));
t = (1-cos(DEGTORAD(theta2)));
x = exchange1.vector[0]/(sqrt((exchange1.vector[0]*exchange1.vector[0])+(exchange1.vector[1]*exchange1.vector[1])+(exchange1.vector[2]*exchange1.vector[2])));
y = exchange1.vector[1]/(sqrt((exchange1.vector[0]*exchange1.vector[0])+(exchange1.vector[1]*exchange1.vector[1])+(exchange1.vector[2]*exchange1.vector[2])));
z = exchange1.vector[2]/(sqrt((exchange1.vector[0]*exchange1.vector[0])+(exchange1.vector[1]*exchange1.vector[1])+(exchange1.vector[2]*exchange1.vector[2])));

double rotation_matrix2[4][4] = {((t*x*x)+c), ((t*x*y)+(s*z)), ((t*x*z)-(s*y)), 0,
((t*x*y)-(s*z)), ((t*y*y)+c), ((t*y*z)+(s*x)), 0,
((t*x*z)+(s*y)), ((t*y*z)-(s*x)), ((t*z*z)+c), 0,
0, 0, 0, 1};


vector3[0][0] = (rotation_matrix2[0][0]*vector2[0][0])+(rotation_matrix2[0][1]*vector2[1][0])+(rotation_matrix2[0][2]*vector2[2][0])+(rotation_matrix2[0][3]*vector2[3][0]);
vector3[1][0] = (rotation_matrix2[1][0]*vector2[0][0])+(rotation_matrix2[1][1]*vector2[1][0])+(rotation_matrix2[1][2]*vector2[2][0])+(rotation_matrix2[1][3]*vector2[3][0]);
vector3[2][0] = (rotation_matrix2[2][0]*vector2[0][0])+(rotation_matrix2[2][1]*vector2[1][0])+(rotation_matrix2[2][2]*vector2[2][0])+(rotation_matrix2[2][3]*vector2[3][0]);
vector3[3][0] = (rotation_matrix2[3][0]*vector2[0][0])+(rotation_matrix2[3][1]*vector2[1][0])+(rotation_matrix2[3][2]*vector2[2][0])+(rotation_matrix2[3][3]*vector2[3][0]);


// Step 3: Rotate reference geometry around line

exchange2.reference_geometry[0] = (rotation_matrix2[0][0]*vector4[0][0])+(rotation_matrix2[0][1]*vector4[1][0])+(rotation_matrix2[0][2]*vector4[2][0])+(rotation_matrix2[0][3]*vector4[3][0]);
exchange2.reference_geometry[1] = (rotation_matrix2[1][0]*vector4[0][0])+(rotation_matrix2[1][1]*vector4[1][0])+(rotation_matrix2[1][2]*vector4[2][0])+(rotation_matrix2[1][3]*vector4[3][0]);
exchange2.reference_geometry[2] = (rotation_matrix2[2][0]*vector4[0][0])+(rotation_matrix2[2][1]*vector4[1][0])+(rotation_matrix2[2][2]*vector4[2][0])+(rotation_matrix2[2][3]*vector4[3][0]);
exchange2.reference_geometry[3] = (rotation_matrix2[3][0]*vector4[0][0])+(rotation_matrix2[3][1]*vector4[1][0])+(rotation_matrix2[3][2]*vector4[2][0])+(rotation_matrix2[3][3]*vector4[3][0]);

// Step 4: Write resultant vector to exchange2

exchange2.vector[0] = 0.7*vector3[0][0];
exchange2.vector[1] = 0.7*vector3[1][0];
exchange2.vector[2] = 0.7*vector3[2][0];
exchange2.vector[3] = vector3[3][0];

// Step 5: Apply Vector Transformation to Start Coordinates

exchange2.end_coord[0] = (exchange1.end_coord[0] + 0.7*vector3[0][0]);
exchange2.end_coord[1] = (exchange1.end_coord[1] + 0.7*vector3[1][0]);
exchange2.end_coord[2] = (exchange1.end_coord[2] + 0.7*vector3[2][0]);

// Step 6: Calculate Radius

exchange2.radius = 0.7*exchange1.radius;
//exchange2.angle1 = exchange1.angle1;
//exchange2.angle2 = exchange1.angle2;

// Step 7: Calculate x
magnitude = sqrt((exchange2.vector[0]*exchange2.vector[0])+(exchange2.vector[1]*exchange2.vector[1])+(exchange2.vector[2]*exchange2.vector[2]));
exchange2.x = (exchange1.x + magnitude);

return (exchange2);

}


line fractal_branch2 (line DataList2)
{


typedef struct line newline;

newline exchange1;
newline exchange2;
exchange1 = DataList2;
double vector1[4][1] = {exchange1.vector[0],
exchange1.vector[1],
exchange1.vector[2],
exchange1.vector[3]};

double vector2[4][1];
double vector3[4][1];
double vector4[4][1] = {exchange1.reference_geometry[0],
exchange1.reference_geometry[1],
exchange1.reference_geometry[2],
exchange1.reference_geometry[3]};

double theta1 = exchange1.angle1;
double theta2 = exchange1.angle2;
double c, s, t, x, y, z, magnitude;



// Initialise Start Coordinates

exchange2.start_coord[0] = exchange1.end_coord[0];
exchange2.start_coord[1] = exchange1.end_coord[1];
exchange2.start_coord[2] = exchange1.end_coord[2];


// Step 1: Rotate vector1 (line vector) around reference geometry to generate vector2

c = cos(DEGTORAD(theta1));
s = sin(DEGTORAD(theta1));
t = (1-cos(DEGTORAD(theta1)));
x = exchange1.reference_geometry[0]/(sqrt((exchange1.reference_geometry[0]*exchange1.reference_geometry[0])+(exchange1.reference_geometry[1]*exchange1.reference_geometry[1])+(exchange1.reference_geometry[2]*exchange1.reference_geometry[2])));
y = exchange1.reference_geometry[1]/(sqrt((exchange1.reference_geometry[0]*exchange1.reference_geometry[0])+(exchange1.reference_geometry[1]*exchange1.reference_geometry[1])+(exchange1.reference_geometry[2]*exchange1.reference_geometry[2])));
z = exchange1.reference_geometry[2]/(sqrt((exchange1.reference_geometry[0]*exchange1.reference_geometry[0])+(exchange1.reference_geometry[1]*exchange1.reference_geometry[1])+(exchange1.reference_geometry[2]*exchange1.reference_geometry[2])));

double rotation_matrix1[4][4] = {((t*x*x)+c), ((t*x*y)+(s*z)), ((t*x*z)-(s*y)), 0,
((t*x*y)-(s*z)), ((t*y*y)+c), ((t*y*z)+(s*x)), 0,
((t*x*z)+(s*y)), ((t*y*z)-(s*x)), ((t*z*z)+c), 0,
0, 0, 0, 1};


vector2[0][0] = (rotation_matrix1[0][0]*vector1[0][0])+(rotation_matrix1[0][1]*vector1[1][0])+(rotation_matrix1[0][2]*vector1[2][0])+(rotation_matrix1[0][3]*vector1[3][0]);
vector2[1][0] = (rotation_matrix1[1][0]*vector1[0][0])+(rotation_matrix1[1][1]*vector1[1][0])+(rotation_matrix1[1][2]*vector1[2][0])+(rotation_matrix1[1][3]*vector1[3][0]);
vector2[2][0] = (rotation_matrix1[2][0]*vector1[0][0])+(rotation_matrix1[2][1]*vector1[1][0])+(rotation_matrix1[2][2]*vector1[2][0])+(rotation_matrix1[2][3]*vector1[3][0]);
vector2[3][0] = (rotation_matrix1[3][0]*vector1[0][0])+(rotation_matrix1[3][1]*vector1[1][0])+(rotation_matrix1[3][2]*vector1[2][0])+(rotation_matrix1[3][3]*vector1[3][0]);



// Step 2: Rotate vector2 around line to generate vector3

c = cos(DEGTORAD(theta2));
s = sin(DEGTORAD(theta2));
t = (1-cos(DEGTORAD(theta2)));
x = exchange1.vector[0]/(sqrt((exchange1.vector[0]*exchange1.vector[0])+(exchange1.vector[1]*exchange1.vector[1])+(exchange1.vector[2]*exchange1.vector[2])));
y = exchange1.vector[1]/(sqrt((exchange1.vector[0]*exchange1.vector[0])+(exchange1.vector[1]*exchange1.vector[1])+(exchange1.vector[2]*exchange1.vector[2])));
z = exchange1.vector[2]/(sqrt((exchange1.vector[0]*exchange1.vector[0])+(exchange1.vector[1]*exchange1.vector[1])+(exchange1.vector[2]*exchange1.vector[2])));

double rotation_matrix2[4][4] = {((t*x*x)+c), ((t*x*y)+(s*z)), ((t*x*z)-(s*y)), 0,
((t*x*y)-(s*z)), ((t*y*y)+c), ((t*y*z)+(s*x)), 0,
((t*x*z)+(s*y)), ((t*y*z)-(s*x)), ((t*z*z)+c), 0,
0, 0, 0, 1};


vector3[0][0] = (rotation_matrix2[0][0]*vector2[0][0])+(rotation_matrix2[0][1]*vector2[1][0])+(rotation_matrix2[0][2]*vector2[2][0])+(rotation_matrix2[0][3]*vector2[3][0]);
vector3[1][0] = (rotation_matrix2[1][0]*vector2[0][0])+(rotation_matrix2[1][1]*vector2[1][0])+(rotation_matrix2[1][2]*vector2[2][0])+(rotation_matrix2[1][3]*vector2[3][0]);
vector3[2][0] = (rotation_matrix2[2][0]*vector2[0][0])+(rotation_matrix2[2][1]*vector2[1][0])+(rotation_matrix2[2][2]*vector2[2][0])+(rotation_matrix2[2][3]*vector2[3][0]);
vector3[3][0] = (rotation_matrix2[3][0]*vector2[0][0])+(rotation_matrix2[3][1]*vector2[1][0])+(rotation_matrix2[3][2]*vector2[2][0])+(rotation_matrix2[3][3]*vector2[3][0]);


// Step 3: Rotate reference geometry around line

exchange2.reference_geometry[0] = (rotation_matrix2[0][0]*vector4[0][0])+(rotation_matrix2[0][1]*vector4[1][0])+(rotation_matrix2[0][2]*vector4[2][0])+(rotation_matrix2[0][3]*vector4[3][0]);
exchange2.reference_geometry[1] = (rotation_matrix2[1][0]*vector4[0][0])+(rotation_matrix2[1][1]*vector4[1][0])+(rotation_matrix2[1][2]*vector4[2][0])+(rotation_matrix2[1][3]*vector4[3][0]);
exchange2.reference_geometry[2] = (rotation_matrix2[2][0]*vector4[0][0])+(rotation_matrix2[2][1]*vector4[1][0])+(rotation_matrix2[2][2]*vector4[2][0])+(rotation_matrix2[2][3]*vector4[3][0]);
exchange2.reference_geometry[3] = (rotation_matrix2[3][0]*vector4[0][0])+(rotation_matrix2[3][1]*vector4[1][0])+(rotation_matrix2[3][2]*vector4[2][0])+(rotation_matrix2[3][3]*vector4[3][0]);

// Step 4: Write resultant vector to exchange2

exchange2.vector[0] = 0.92*vector3[0][0];
exchange2.vector[1] = 0.92*vector3[1][0];
exchange2.vector[2] = 0.92*vector3[2][0];
exchange2.vector[3] = vector3[3][0];

// Step 5: Apply Vector Transformation to Start Coordinates

exchange2.end_coord[0] = (exchange1.end_coord[0] + 0.92*vector3[0][0]);
exchange2.end_coord[1] = (exchange1.end_coord[1] + 0.92*vector3[1][0]);
exchange2.end_coord[2] = (exchange1.end_coord[2] + 0.92*vector3[2][0]);

// Step 6: Calculate Radius

exchange2.radius = 0.92*exchange1.radius;
//exchange2.angle1 = exchange1.angle1;
//exchange2.angle2 = exchange1.angle2;

// Step 7: Calculate x
magnitude = sqrt((exchange2.vector[0]*exchange2.vector[0])+(exchange2.vector[1]*exchange2.vector[1])+(exchange2.vector[2]*exchange2.vector[2]));
exchange2.x = (exchange1.x + magnitude);

return (exchange2);

}



line fractal_branch3 (line DataList3)
{

// Define Newline

typedef struct line newline;

newline exchange1;
newline exchange2;
double magnitude;
exchange1 = DataList3;

// Initialise Start Coordinates

exchange2.start_coord[0] = exchange1.end_coord[0];
exchange2.start_coord[1] = exchange1.end_coord[1];
exchange2.start_coord[2] = exchange1.end_coord[2];

// Calculate Transformation Vector

// x component
exchange2.vector[0] = exchange1.vector[0];
// y component
exchange2.vector[1] = exchange1.vector[1];
// z component
exchange2.vector[2] = exchange1.vector[2];


// Apply Vector Transformation to Start Coordinates

exchange2.end_coord[0] = (exchange1.end_coord[0]+ exchange2.vector[0]);
exchange2.end_coord[1] = (exchange1.end_coord[1]+ exchange2.vector[1]);
exchange2.end_coord[2] = (exchange1.end_coord[2]+ exchange2.vector[2]);

// Step 7: Calculate x
magnitude = sqrt((exchange2.vector[0]*exchange2.vector[0])+(exchange2.vector[1]*exchange2.vector[1])+(exchange2.vector[2]*exchange2.vector[2]));
exchange2.x = (exchange1.x + magnitude);

return (exchange2);

}