/* bibilioteca complexos.h  
	Assunto:  	Álgebra de números complexos e funções complexas
	Programa 	sem erros
	

   Descrição: operações com números complexos e funções complexas
         

   por Tarcisio Praciano Pereira - 10 lições para aprender C
   Sobral, Novembro de 2002 - na era Lula	- UVA

	Palavras-chave: imagem de curvas por polinômios complexos
	palavras chave: struct, estrutura
	palavras chave: fronteira de uma superficie de Riemann 	
   por Tarcisio Praciano Pereira - 10 lições para aprender C
   Sobral,  de 2008	- UeVA    	      		
*/

# include <stdio.h>  	// (1) leitura da biblioteca padrão
# include <stdlib.h>
# include <math.h>
# include "/home/tarcisio/tex/c/ufg/ambiente.h"

			//	(2)	O tipo, número complexo
typedef	struct { float a; float b; } Complexo;

// funções complexas
Complexo	P(Complexo z,  Complexo a0, Complexo a1,Complexo a2);
Complexo	Q(Complexo z, Complexo a0, Complexo a1,Complexo a2, Complexo a3);
Complexo Exp(Complexo z);
Complexo	mais(Complexo z, Complexo w);
Complexo	vezes(Complexo z, Complexo w);
float		Re(Complexo z);
float		Im(Complexo z);
// *******************************
// utilitários  ******************
void		escreve(Complexo z); // escreve  a + b i
Complexo	recebe(Complexo z, Complexo w); // atribuição  w --> z
Complexo	faz_complexo(float a, float b); // cria  a,b  --> z = (a,b)
// *********  gnuplot  ***********************************
// produzem  arquivos:  dados, transfere    para gnuplot
int	curva_retangulo(int selecao, Complexo z, float raio,
         Complexo a0, Complexo a1,Complexo a2, Complexo a3); // imagem de quadrado de centro em z
int	curva_circulo(int selecao, Complexo z,float raio,
         Complexo a0, Complexo a1,Complexo a2, Complexo a3); // imagem de quadrado de centro em z
void		prepara_gnuplot(char msg[80]); // transfere




// ******************  funções complexas  *************************
//      a0 + z*(a1 + z*( a2 + z*a3))) 
//     mais(a0, vezes(z, mais(a1, vezes(z, mais(a2, vezes(a3,z))))))
Complexo	P(Complexo z,  Complexo a0, Complexo a1,Complexo a2)   // segundo grau
{
	Complexo temp = vezes(a2,z);
	temp = mais(temp, a1);
	temp = vezes(temp, z);
	temp = mais(temp, a0);
	return(temp);
}
// a0 + z(a1 + z(a2 + z(a3 ))))
Complexo	Q(Complexo z,  Complexo a0, Complexo a1,Complexo a2, Complexo a3)   // terceiro  grau
{
	Complexo temp = vezes(a3,z);
	temp = mais(temp, a2);
	temp = vezes(temp, z);
	temp = mais(temp, a1);
	temp = vezes(temp, z);
	temp = mais(temp, a0);	
	return(temp);
}



//  exp(z) = exp(x+iy) = exp(x)*exp(iy) = exp(x)( cos(y) + i sin(y) )
Complexo	Exp(Complexo z)
{
	Complexo temp;
	temp.a = exp(z.a)*cos(z.b); // exp(x)cos(y)
	temp.b = exp(z.a)*sin(z.b); // exp(x)sin(y)
	return(temp);
}

Complexo	mais(Complexo z, Complexo w)
{
	Complexo temp;
	temp.a = z.a + w.a;
	temp.b = z.b + w.b;
	return(temp);
}

Complexo vezes(Complexo z, Complexo w)
{
	Complexo temp;
	temp.a = z.a*w.a - z.b*w.b;
	temp.b = z.b*w.a + z.a*w.b;
	return(temp);
}


float		Re(Complexo z)
	{
		return z.a;
	}

float		Im(Complexo z)
	{
		return z.b;
	}




// ***************  utilitários  ******************

void	escreve(Complexo z)  // saida de dados
{
	printf(" %2.2f + %2.2f i", Re(z),Im(z));
}

Complexo	recebe(Complexo z, Complexo w) // atribuição  w -> z
{
	z.a = w.a;
	z.b = w.b;
	return(z);
}

Complexo	faz_complexo(float a, float b) // cria a,b -> (a,b)
{
	Complexo temp;
	temp.a= a; temp.b = b;
	return(temp);
}



// ********************  gnuplot    ******************************
// produzem  arquivos:  dados, transfere    para gnuplot
// percorre um quadrado unitário centrado em z 
int	curva_retangulo(int selecao, Complexo z, float raio,
           Complexo a0, Complexo a1,Complexo a2, Complexo a3) 
{
  float a=z.a ,b=z.b;
  float x,y;
  float delta = 0.001;
  //  Complexo Z; // local  - inutil ??
  FILE	*dados;
  dados	= fopen("dados", "w");
  // primeiro lado do quadrado - fronteira de [0,1]x[-3.1415,3.1415]
  switch(selecao)
  {
  case 4:
  x = -0.5;
  y = -0.5;//-3.1415;
  while(y <= 0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = Exp(z); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 y+=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }
  // terceiro  lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = -0.5;
  y = 0.5;
  while(x <= 0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = Exp(z); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 x+=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }  
  // segundo  lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = 0.5;
  y = 0.5;
  while(y >= -0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b);  // translada para o centro z
			 z = Exp(z); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 y-=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }    
  // quarto  lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = 0.5;
  y = -0.5;// 3.1415;
  while(x >=-0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = Exp(z); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 x-=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }    
  break;
  
  case 5:
  x = -0.5;
  y = -0.5;//-3.1415; // primeiro lado
  while(y <= 0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = P(z, a0,  a1,  a2); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 y+=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }
  // segundo  lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = -0.5;
  y = 0.5;
  while(x <= 0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = P(z, a0,  a1,  a2); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 x+=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }  
  // terceiro   lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = 0.5;
  y = 0.5;
  while(y >= -0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b);  // translada para o centro z
			 z = P(z, a0,  a1,  a2); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 y-=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }    
  // quarto  lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = 0.5;
  y = -0.5;// 3.1415;
  while(x >=-0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = P(z, a0,  a1,  a2); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 x-=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }    
  break;  
  
  case 6:
  x = -0.5;
  y = -0.5;//-3.1415; // primeiro lado
  while(y <= 0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = Q(z, a0,  a1,  a2, a3); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 y+=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }
  // segundo  lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = -0.5;
  y = 0.5;
  while(x <= 0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = Q(z, a0,  a1,  a2, a3); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 x+=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }  
  // terceiro   lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = 0.5;
  y = 0.5;
  while(y >= -0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b);  // translada para o centro z
			 z = Q(z, a0,  a1,  a2, a3); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 y-=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }    
  // quarto  lado do quadrado - fronteira de [-3,3]x[-3,3]  
  x = 0.5;
  y = -0.5;// 3.1415;
  while(x >=-0.5)
  {
			 z = faz_complexo(raio*x + a, raio*y + b); // translada para o centro z
			 z = Q(z, a0,  a1,  a2, a3); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z));
			 x-=delta; //  y+= delta; ou delta2  // (10)  troque aqui para alterar a ênfase...
  }    
  break;  
  } // end of switch
  fclose(dados);
  return(0);
}

int	curva_circulo(int	selecao, Complexo z, float raio,
     Complexo a0, Complexo a1,Complexo a2, Complexo a3) // imagem de quadrado de centro em z
{
  float a=z.a ,b=z.b;
  float rho=raio,theta;
  float delta = 0.001;
  float doispi	= 6.28318530717958647696;
  // Complexo Z; // local  inutil 
  FILE	*dados;
  dados	= fopen("dados", "w");
  // círculo de raio rho  e centro em z = (a,b)
  switch(selecao)
  {
  case 1:
  theta = 0;
  while( theta <= doispi)
  {
          float x = a + rho*cos(theta),y = b + rho*sin(theta);
			 z = faz_complexo(x, y); // translada para o centro z
			 z = Exp(z); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z)); // imprime em dados
			 theta+=delta; 
  }
  break;
  case 2:
  theta = 0;
  while( theta <= doispi)
  {
          float x = a + rho*cos(theta),y = b + rho*sin(theta);
			 z = faz_complexo(x, y); // translada para o centro z
			 z = P(z, a0,  a1,  a2); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z)); // imprime em dados
			 theta+=delta; 
  }
  break;
  case 3:
  theta = 0;
  while( theta <= doispi)
  {
          float x = a + rho*cos(theta),y = b + rho*sin(theta);
			 z = faz_complexo(x, y); // translada para o centro z
			 z = Q(z, a0,  a1,  a2, a3); // exp(x+iy)
			 fprintf(dados, "%f  %f\n", Re(z), Im(z)); // imprime em dados
			 theta+=delta; 
  }
  break;  
  } // fim do switch 
  fclose(dados);
  return(0);
}


void	prepara_gnuplot(char msg[80])
{
  FILE  *transfere;
  transfere =  fopen("transfere", "w");
  fprintf(transfere, "set title \" %s \" \n", msg); // qual é a figura 
  fprintf(transfere, "set pointsize 0.1 \n");
  fprintf(transfere, "set grid  \n");  
  fprintf(transfere, "plot \"dados\" with points \n");
  fprintf(transfere, "pause -2\n");
  fclose(transfere);
}




/* Comentários:  A numeração dos comentários não
		  é continuada, pode dar saltos...para
		  facilitar a reutilização de programas.

	(2)	O tipo, número complexo. Com "typedef" podemos criar novos tipos 
	de dados. Aqui foi criado o tipo "Complexo" com auxílio de "struct".
	Um número complexo tem dois campos indicados por "a" e por "b".

(10) na função curva_imagem()  a escolha de uma malha mais fina em
    um dos eixos determina a ênfase - quer dizer, se quero a imagem
    de segmentos de reta paralelo ao eixo OX  ou  ao eixo OY 
*/