define copyleft()
{
printf("Este programa é distribuido sob GPL. \n");
printf("Se você não souber o que é a GPL, envie um e-mail ao autor \n");
printf("tarcisio@member.ams.org \n");
printf("mas em suma isto quer dizer que você pode usar livremente o sistema\n");
printf("de programas aqui mencionados, desde que esta copyleftção seja mantida\n");
printf("em todas as alterações que você fizer do programa e das cópias que \n");
printf("você dele fizer - é a proteção do direito autoral que não lhe impede\n");
printf("de usar livremente o conhecimento mas que mantém a memória de quem \n");
printf("participou de sua produção.\n");
}


define  apeteco2()
{
   local temporario=0;
   printf("%s","================================================= \n");
   printf(" Aperte <enter> para continuar \n");
   scanf("%c",temporario);
}

define  apetecof()
{
   local temporario=0;
   printf("%s","================================================= \n");
   printf(" Aperte <enter> para terminar  \n");
   scanf("%c",temporario);
}

define informa(){
print "===============  A lei do coseno para áreas de superfícies em 3D.  ============= ";
print "Este programa constrói um círculo unitário centrado no ponto (a,b) no plano XOY e";
print "sua imagem sobre o plano tangente ao gráfico de uma função z = F(x,y) no ponto";
print "(a,b,F(a,b)). O objetivo é mostrar que o resultado é uma elipse contida num ";
print "cilindro unitário mas no plano onde o círculo é projetado. Isto pode ser";
print "observado se você rodar a imagem usando a janela do gnuplot - escolha a posição";
print "em que a projeção fique perpendicular ao plano XOY e verá que a elipse ficará ";
print "projetada de volta dentro do círculo unitário no plano de saída, mas ela se expande";
print "e portanto sua área aumenta além da área do círculo.";
print "O real objetivo é mostrar-lhe que a área da elipse-imagem de um círculo projetado ";
print "em um plano é maior do que a área do círculo e consequentemente que o coeficiente";
print "de distorção provido pela regra da cadéia tem que ser maior do que 1. ";
print "Isto vai justificar a regra do coseno para áreas, ver Apostol, Calculus II página";
print "285, ou consulte o índice remissivo de qualquer livro de Cálculo a duas variáveis.";
print "Para este caso, área da elipse, o seu valor é (área do círculo)/cos(t) em que t é ";
print "o ângulo entre o plano tangente e o plano XOY - logo o ângulo do vetor perpendicular";
print "ao plano o que nos leva ao gradiente da função z = F(x,y). ";
print "Este programa cria dois arquivos, \"transfere\" - com comandos para gnuplot e o arquivo";
print "\"dados\" contendo as coordenadas do círculo e da elipse-imagem. Você poderá desejar";
print "apagar estes arquivos quando esta experiência estiver encerrada. Também você poderá ";
print "visualizar a imagem sem usar este programa usando, na linha de comando (em Linux funciona!)";
print "                        gnuplot transfere";
print "observe, sem aspas, gnuplot vai entender que você está lhe passando um arquivo de comandos.";
print "Se você desejar experimentar com outras funções diferentes da usada no programa, procure com";
print "um editor de textos F(x,y), DF_x(x,y), DF_y(x,y) e corrija as equações - observe que ";
print "          DF_x(x,y), DF_y(x,y)  ";
print "é a notação que estou usando para as derivadas parciais de F.  ";
} 



define  F(x, y)
{
  return power(x,2) + 3*x*y + power(y,2);
}
define DF_x(x, y)
{
  return 2*x + 3*y;
}

define   DF_y(x, y)
{
  return  3*x + 2*y;
}

define P(x,y,a,b){
 return F(a,b) + DF_x(a,b)*(x-a) + DF_y(a,b)*(y-b);
}

define cria_dados(a,b,rho)
{
  local dados = fopen("dados", "w");
  local theta = 0, passo = 0.01;
  local dois_pi = 8*atan(1);
  local x,y;
  while(theta < dois_pi){
  x = a + rho*cos(theta); y = b +rho*sin(theta);  ## círculo unitário com centro em  (a,b)
  fprintf(dados, "%f %f %f \n",x, y,P(x,y,a,b));
  theta +=passo;
  }
  theta = 0;
  fprintf(dados, "\n \n");  ## para separar os gráficos - gnuplot muda as cores
  while(theta < dois_pi){
  x = a + rho*cos(theta); y = b +rho*sin(theta);  ## círculo unitário com centro em  (a,b)
  fprintf(dados, "%f %f %f \n",x, y, 0.0); ## no plano XOY pontos para o splot 
  theta +=passo;
  }
  fclose(dados)
}	

define cria_transfere(a,b)
{
local transfere  = fopen("transfere", "w");
local c = F(a,b);
fprintf(transfere, "a = %f; b = %f \n",a,b );
fprintf(transfere, "pow(x,n) = x**n; \n" );
fprintf(transfere, "F(x,y) = pow(x,2) + 3.0*x*y + pow(y,2);\n" );  ##  z = F(x,y)
fprintf(transfere, "DF_x(x,y) = 2.0*x + 3.0*y ;\n" );  ##  Derivada parcial relativamente a x                  
fprintf(transfere, "DF_y(x,y) = 3.0*x + 2.0*y ;\n" );   ##  Derivada parcial relativamente a y
fprintf(transfere, "P(x,y) = F(a,b) + DF_x(a,b)*(x- %f)+ DF_y(a,b)*(y- %f);\n",a,b);
fprintf(transfere, "set pointsize 0.1 \n"); 
fprintf(transfere, "set xrange [%f-5:%f+5]\n",a);
fprintf(transfere, "set yrange [%f-5:%f+5]\n",b);
fprintf(transfere, "set zrange [%f-15:%f+15]\n",c);
fprintf(transfere, "splot F(x,y), P(x,y), 0, \"dados\" with points\n" );
fprintf(transfere, "print \"Aperte uma tecla para terminar\";\n");
fprintf(transfere, "pause -2 \n");
fclose(transfere);
}

define main(x,y){
	local a=x,b=y;
	system("clear"); ## sob outros sistemas substitua clear por cls ....pode ser que funcione!
	informa(); apeteco2();
	print "Forneça o ponto (a,b)  onde a imagem será centrada ";
	print " a = "; scanf("%f",a);
	print " b = "; scanf("%f",b);
	print "Você deseja a imagem de um círculo unitário centrado no ponto (",a,",",b,")";
	print "O programa rodou com valores default! ";
	print "Rode novamente o programa digitando: ";
	print "  main()                                                      \n";
	print "......   é só copiar até fim da linha e colar  ";
	print "e dando enter ...!   no terminal do calc. ";
	apeteco2();
	cria_dados(a,b,1);
	cria_transfere(a,b);
	system("gnuplot -persist transfere "); 
        apetecof();
}


main(-2,-1);