#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
//#include <mpi.h>


// fractal computing routine of Nexus project
// one computer


// globals


	#define ELEMENT_SIZE	4                       // bytes per sended depth array elements
							// 2 = depth_value < 64000
	int l_apu= 0;





int char2int(const char *cline)
{

// converts character table to the double


	int value;
        int i,j,number, negflag;

	char c;
	

	c = *cline;  
	i = 1;
	value = 0;            
	
	negflag = 0;
	if (c=='-') {
		negflag = 1;
		cline++;
		c = *cline;
		}

        while ((c != 0) && (c != '.')) { 

	   
	   number = c - '0';			// i:s number
           value += number;
	   
	   cline++;
	   c = *cline;
		
	   if ((c != 0) && (c != '.'))	{
		value *= 10;	// more values to come
		} 
 	}


	if (negflag != 0) value = 0-value;
	return value;
}


double char2double(const char *cline)
{

// converts character table to the double


	double value,decimal,divided;
        int i,j,number, negflag;

	char c;       
	

	c = *cline;  
	i = 1;
	value = 0;            
	
	negflag = 0;
	if (c=='-') {
		negflag = 1;
		cline++;
		c = *cline;
		}	

        while ((c != 0) && (c != '.')) { 

	   
	   number = c - '0';			// i:s number
           value += number; 
	   
	   cline++;
	   c = *cline;
		
	   if ((c != 0) && (c != '.'))	{
		value *= 10;	// more values to come
		}
 	}
	
// decimals:


	if (c != '.') { 			// no decimals
		if (negflag != 0) value = 0 - value;		//negativity		
		return value; 
		}

	cline++;		// jump the decimalpoint

	
	divided = 10;
	c = *cline;
	while (c != 0) {
	   number = c - '0';		
           
  	   if (number != 0) {	
	   decimal = number; 
	   decimal /= divided;	
	   value += decimal;
	}	
	 
	   divided *= 10;	  
	   cline++;
	   c = *cline;
	}


	if (negflag != 0) value = 0 - value;		// negativity	
        return value;
   
}


char make_char(int value)
{
        int i,help,hand;
        char ret;
        ret = 0x00;
        hand = 0x01;

        for (i=0; i < 8; i++) {
                help = (value & hand);
                if (help == hand) {
                        ret = (ret | hand);
                        }
                hand = (hand << 1);
        }
        ret = (ret | 0x01);


        return ret;
}



int code_integer(int value)
{
        int i,help,mult,ret,hor;
        hor = 0x01;

        ret = 0x00;
        //ret = 1;
        //ret = (ret | !ret);              // bits off
        //ret = (ret & 0x01);


        mult = 1;

        if (value == 0) {
                return ret;
                }

        if (value == 1) {
                ret = (ret | hor);
                return ret;
                }

        for (i=0; (i < 31 || value > 0); i++) {
                help = (value % 2);

                if (help != 0) {
                        ret = (ret | hor);
                        value -= 1;
                }
                value /= 2;
                hor = (hor << 1);

        }

        return ret;

}


// encoding/decoding-functions


void encode_line(const int *depth_line, char *depth_line_char, int steps_x)
{
// encodes line from integer array to char array, which can be sended
// upper bytes from begin
// in each byte there will be seven bits of integer, the first bit needs
// to be 1 to avoid the char to be zero, wich means end of a char-table


	int i,j,k,l,shiftbits;
	int element, help1;
	char help2;
	k = 0;
        l= 0;
	for (i=0; i < steps_x; i++) {
		element = depth_line[i];
		shiftbits = (ELEMENT_SIZE * 6) - 6;


                //element = l_apu;
                //l_apu++;
                //if (l_apu >= 100000) l_apu = 0;


		for (j=0; j < ELEMENT_SIZE; j++) {

			help2 = 0xff;

			if (element >= 0) {
                                element = code_integer(element);

                                help1 = element;
				help1 = (help1 >> shiftbits);

                                help1 = (help1 << 1);
                                help1 = (help1 & 0xff);
                                //help1 = (help1 & 0x7f);

                                help1 = (help1 | 0x01);
                                help1 = (help1 & 0x7f);
				help2 = (char)help1;
                                //help2 = make_char(help1);
                                //help2 = (help2 & 0x7f);
				shiftbits -= 6;
			}
			depth_line_char[k + j] = help2;
			}
		k += ELEMENT_SIZE;
		}
		depth_line_char[k] = 0;			// the end of a line

}

void decode_line(int *depth_line, const char *depth_line_char, int steps_x)
{

	int i,j,k,flag;
	int help1;
	int value;
	char element;
	k = 0;
	for (i=0; i < steps_x+1; i++) {
		value = 0;
                value = (value & ~value);
		flag = 1;
		for (j=0; j < ELEMENT_SIZE; j++)
		{
			element = depth_line_char[k+j];
			help1 = 0;
			help1 = (int) element;
			help1 = (help1 & 0xff);
			if (help1 != 0xff) flag = 0; 	

			//help1 = (help1 & 0x7f);   // first 1 bit off
                        //help1 = (help1 & 0x80);

			value = (value | help1);
			if (j < ELEMENT_SIZE-1) value = value << 7;
		}
		if (flag != 0) value = -1;			// negative value
		depth_line[i] = value;
		k += ELEMENT_SIZE;
	   }				
		
}



// fractal calculating class


	class CalculateFractal
{


private:
	double Start_x;
	double Start_y; 
	double Start_z; 
	double End_x;
	double End_y; 
	double End_z; 
	double C_x; 
	double C_y; 
	double C_z; 
	double C_w; 
	double W_value;
	double Sin_x;
	double Sin_y; 
	double Sin_z;
	double Cos_x; 
	double Cos_y; 
	double Cos_z;
	double Origo_x;
	double Origo_y; 
	double Origo_z;



	int Steps_x;
	int Steps_y;
	int Steps_z; 
	int Steplength_z; 
	int Iterations; 
	int IBound;
	int ID;
	int Tasks;




private:
	
	

 	int is_in_julia(double x,double y,double z)
{

// calculating routine

	 double xt,yt,zt,w;
	 double length, temp;
	 xt = x;
	 yt = y;
	 zt = z;	
	 w=W_value;
 	 int i=0;

       	do {
  	  temp=xt+xt;
  	  xt=xt*xt-yt*yt-zt*zt-w*w;
	    yt=temp*yt;
  	    zt=temp*zt;
  	    w=temp*w;

  	    length=xt*xt+yt*yt+zt*zt+w*w;
		

	    xt += C_x;
            yt += C_y;
	    zt += C_z;
	    w += C_w;
	    i++;

 	 } while ((i < Iterations) && (length < IBound));
	
	if (length >= IBound) return 0;		// out of set
	return 1;

}



 void rotate(double *x, double *y, double *z)
{
// rotates a point describted in x,y,z
 	
	
	 double xt,yt, zt;
	
	 double	 xx,xy,xz,yx,yy,yz,zx,zy,zz;
	 xt = *x;
	 yt = *y;
  	 zt = *z;
	 xt-=Origo_x;yt-=Origo_y;zt-=Origo_z;


	
 	 xy=yt*Cos_x-zt*Sin_x;
 	 xz=yt*Sin_x+zt*Cos_x;
 	 xx=xt;

 	 xt=xx;
 	 yt=xy;
 	 zt=xz;
 	 yx=xt*Cos_y+zt*Sin_y;
 	 yz=-xt*Sin_y+zt*Cos_y;
 	 xt=yx;
 	 zt=yz;
 	 zx=xt*Cos_z-yt*Sin_z;
 	 zy=xt*Sin_z+yt*Cos_z;
 	 xt=zx;
 	 yt=zy;
 	
	 xt+=Origo_x;yt+=Origo_y;zt+=Origo_z;
	 	
	 *x = xt;
	 *y = yt;
	 *z = zt;
}


 public:
	

	int get_parameters	(char *argv[], int argc)
{
// take parameters

	double angle_x, angle_y,angle_z;

        

	Start_x = -1.2;
	Start_y = -1.2;
	Start_z	= -1.2;
	End_x = 1.2;
	End_y = 1.2;
	End_z = 1.2;
	C_x = -0.1;
	C_y = 0.2;
	C_z = 0.3;
	C_w = -0.1;
	W_value = 0;
	angle_x = 0;
	angle_y = 0;
	angle_z = 0;
	Steps_x = 300;
	Steps_y = 300;
	Steps_z = 100;
	Steplength_z = 10;
	Iterations = 10;
	IBound = 16;
	

        if (argc >= 19) {
/*	   	float p[20];


		sscanf(argv[1], "%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f",
	  		&p[0], &p[1], &p[2]
	  		, &p[3], &p[4], &p[5]
	  		, &p[6], &p[7], &p[8]
	  		, &p[9], &p[10], &p[11]
	  		, &p[12], &p[13], &p[14]
	  		, &p[15], &p[16], &p[17]
	  		, &p[18], &p[19], &p[20]

	  		);	   
	   
*/
	        Start_x = char2double(argv[1]);
	        Start_y = char2double(argv[2]);
	        Start_z	= char2double(argv[3]);
	        End_x = char2double(argv[4]);
	        End_y = char2double(argv[5]);
	        End_z = char2double(argv[6]);
	        C_x =char2double(argv[7]);
	        C_y = char2double(argv[8]);
	        C_z = char2double(argv[9]);
	        C_w = char2double(argv[10]);
	        W_value = char2double(argv[11]);
	        angle_x = char2double(argv[12]);
	        angle_y = char2double(argv[13]);
	        angle_z = char2double(argv[14]);
	        Steps_x = char2int(argv[15]);
	        Steps_y = char2int(argv[16]);
	        Steps_z = char2int(argv[17]);
	        Steplength_z = char2int(argv[18]);
	        Iterations = char2int(argv[19]);
	        IBound = char2int(argv[20]);
	
        }

	/*cout << " Steps_x: ";
	cout << Steps_x;
	cout << " Steps_y: ";
	cout << Steps_y;*/
	
	Sin_x = sin(angle_x);
	Sin_y = sin(angle_y);
	Sin_z = sin(angle_z);

	Cos_x = cos(angle_x);
	Cos_y = cos(angle_y);
	Cos_z = cos(angle_z);

	Origo_x = Start_x + ((End_x - Start_x) / 2);
	Origo_y = Start_y + ((End_y - Start_y) / 2);
	Origo_z = Start_z + ((End_z - Start_z) / 2);

}



	int take_depth_value(int xi, int yi)
{
// gives a depth value to the point

	double w, value;
	double tx,ty,tz;
	double x,y,z;
	int depthvalue = -1;
	int k;
	double temp;
	int zi = 0;
	do  {
		x = xi;
		x /= Steps_x;
		x *= (End_x - Start_x);
		x += Start_x;

		y = yi;
		y /= Steps_y;
		y *= (End_y - Start_y);
		y += Start_y;

		z = zi;
		z /= (Steps_z * Steplength_z);
		z *= (End_z - Start_z);
		z += Start_z;
		
		tx = x; ty = y; tz = z;
		rotate(&tx,&ty,&tz);
		if (is_in_julia(tx,ty,tz) != 0) {

// search "Steplength" steps backwards to get a smoother value:

			depthvalue = zi;
			

			if (zi >= 1) {		
				zi--;

				for (k=1; k < Steplength_z; k++) {
					z = zi;
					z /= (Steps_z * Steplength_z);
					z *= (End_z - Start_z);
					z += Start_z;

					tx = x; ty = y; tz = z;
					rotate(&tx,&ty,&tz);

					if (is_in_julia(tx,ty,tz) != 0) {
						depthvalue = zi;
						}
					zi--;
					}

				}
				}
		
				zi += Steplength_z;

	} while ((zi <= (Steps_z * Steplength_z)) && (depthvalue < 0)); 
	
	return depthvalue;	

}


// get methods;

	int get_Steps_x() {
		return Steps_x;
	
	}

	int get_Steps_y() {
		return Steps_y;
	}

	int get_Id() {
		return ID;
	}
	int get_Tasks() {
		return Tasks;
	}

};




int main(int argc, char *argv[])

{

// parameter order:
//  1. Start_x 	    (double)
//  2. Start_y 	    (double)
//  3. Start_z	    (double)
//  4. End_x 	    (double)
//  5. End_y 	    (double)
//  6. End_z 	    (double)
//  7. C_x 	    (double)
//  8. C_y	    (double)
//  9. C_z	    (double)
// 10. C_w	    (double)
// 11. W_value		(double)
// 12. angle_x		(double)
// 13. angle_y 		(double)
// 14. angle_z		(double)
// 15. Steps_x	    (integer)
// 16. Steps_y      (integer)
// 17. Steps_z	    (integer)
// 18. Steplength_z (integer)
// 19. Iterations	(integer)
// 20. IBound		(integer)
// 21. ID		(integer)
// 22. Tasks		(integer)

	int steps_x,steps_y;




	CalculateFractal *cl= new CalculateFractal();

	int *depth_line,*depth_test_line;
	char *depth_line_char,apuchar;
	int i,j,k;



        cl->get_parameters(argv,argc);

	steps_x = cl->get_Steps_x();
	steps_y = cl->get_Steps_y();
	

	depth_line = NULL;
	depth_line_char = NULL;
	depth_test_line = (int*)calloc(steps_x+1, sizeof(int));
	depth_line = (int*)calloc(steps_x+1, sizeof(int));
	depth_line_char = (char*)calloc((steps_x+1) * ELEMENT_SIZE + 1, sizeof(char));
	if ((depth_line == NULL) || (depth_line_char == NULL)) { cout<<("\nCan't allocate memory \n");
					return 1; }


	for (i=0; i < steps_y; i++) {
		depth_line[0] = i;					// line number
		for (j=0; j < steps_x; j++) {
			depth_line[j+1] = cl->take_depth_value(j,i);
			}

            
		encode_line(depth_line, depth_line_char,steps_x+1);

                //decode_line(depth_line,depth_line_char,steps_x);
                //for (k=0; k < steps_x;k++) {
                //        cout <<depth_line[k];
                //        cout << " ";
                //        }


		//cout << depth_line_char;
                fwrite(depth_line_char, 1, (steps_x+1)*ELEMENT_SIZE, stdout);

		//fflush(stdout);

		}

	free(depth_line);
	free(depth_line_char);
	return 0;

}