Tetration/Code/Polynomial25power

// Generator of the upper right part of the figure 
// In order to generate the figure above, you need 
// the plotter routines
// ContourPlot/code/ado.cin , than makes the header of the eps file, and
// ContourPlot/code/conto.cin, that draws there the contour lines.
// With these routines, you may compile and execute the source below
//
// Copyleft 2008 by Dmitrii Kouznetsov
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
#include <complex.h>
#define z_type complex<double>
#define Re(x) x.real()
#define Im(x) x.imag()
#define I z_type(0.,1.)
z_type F4natu(z_type z) // polybonial approcimation
{ int K=51,k;
DB coef[51]={ 1.
, 1.091767351258322138
, 0.271483212901696469
, 0.212453248176258214
, 0.069540376139988952
, 0.044291952090474256
, 0.014736742096390039
, 0.008668781817225539
, 0.002796479398385586
, 0.001610631290584341
, 0.000489927231484419
, 0.000288181071154065
, 0.000080094612538551
, 0.000050291141793809
, 0.000012183790344901
, 0.000008665533667382
, 0.000001687782319318
, 0.000001493253248573
, 0.000000198760764204
, 0.000000260867356004
, 0.000000014709954143
, 0.000000046834497327
,-0.000000001549241666
, 0.000000008741510781
,-0.000000001125787310
, 0.000000001707959267
,-0.000000000377858315
, 0.000000000349577877
,-0.000000000105377012
, 0.000000000074590971
,-0.000000000027175982
, 0.000000000016460766
,-0.000000000006741873
, 0.000000000003725329
,-0.000000000001639087
, 0.000000000000858364
,-0.000000000000394374
, 0.000000000000200252
,-0.000000000000094420
, 0.000000000000047121
,-0.000000000000022563
, 0.000000000000011155
,-0.000000000000005391
, 0.000000000000002652
,-0.000000000000001289
, 0.000000000000000633
,-0.000000000000000309
, 0.000000000000000151
,-0.000000000000000074
, 0.000000000000000036
,-0.000000000000000018
};
z_type s=1.,t=z; 
for(k=1;k<=25;k++) { s+=coef[k]*t; t*=z; }	//only 25 terms are used
return s;
}
//
#include "conto.cin" 
main(){ int j,k,m,n; DB x,y, p,q, t; z_type z,c,d, cu,cd;
z_type Zo=z_type(.31813150520476413, 1.3372357014306895);
z_type Zc=z_type(.31813150520476413,-1.3372357014306895);
int M=150,M1=M+1;
int N=151,N1=N+1;
DB X[M1],Y[N1], g[M1*N1],f[M1*N1], w[M1*N1]; // w is working array.
char v[M1*N1]; // v is working array
FILE *o;o=fopen("TetrationPolynom25terms.eps","w");ado(o,0,0,124,124);
fprintf(o,"62 62 translate\n 20 20 scale\n");
//
DO(m,M1) X[m]=-3+.04*m;
DO(n,N1) Y[n]=-3+.04*(n-.5);
//
for(m=-3;m<4;m++) {	if(m==0){M(m,-3.1)L(m,3.1)}
			else	{M(m,-3)L(m,3)}			}
for(n=-3;n<4;n++) {M(  -3,n)L(3,n)} fprintf(o,".006 W 0 0 0 RGB S\n");
//
z_type tm,tp,F[M1*N1];
//
DO(m,M1)DO(n,N1){	g[m*N1+n]=9999;
			f[m*N1+n]=9999;
		}
for(m=0;m<M1;m++){x=X[m];
			DO(n,N1){y=Y[n]; z=z_type(x,y);
				c=F4natu(z); F[m*N1+n]=c;
				}
		}
//
DO(m,M1)
DO(n,N1){
		c=(F[m*N1+n]); p=Re(c); q=Im(c);
		if(p>-999 && p<999) g[m*N1+n]=p;
		if(q>-999 && q<999) f[m*N1+n]=q;
	}
//
p=1;q=99;
                 conto(o,f,w,v,X,Y,M,N, (-3*M_PI ),-q,q); fprintf(o,".04 W 1 0 1 RGB S\n");
                 conto(o,f,w,v,X,Y,M,N, ( -M_PI  ),-q,q); fprintf(o,".04 W 1 0 1 RGB S\n");
for(n=2;n<10;n+=2)conto(o,f,w,v,X,Y,M,N, (-3.+.1*n),-p, p); fprintf(o,".01 W 0 1 0 RGB S\n");
                 conto(o,f,w,v,X,Y,M,N, (-2.     ),-q ,q); fprintf(o,".03 W 1 0 0 RGB S\n");
for(n=2;n<10;n+=2)conto(o,f,w,v,X,Y,M,N, (-2.+.1*n),-p ,p); fprintf(o,".01 W 0 1 0 RGB S\n");
                 conto(o,f,w,v,X,Y,M,N, (-1.     ),-q,q); fprintf(o,".03 W 1 0 0 RGB S\n");
for(n=2;n<10;n+=2)conto(o,f,w,v,X,Y,M,N, (-1.+.1*n),-p ,p); fprintf(o,".01 W 0 1 0 RGB S\n");
                 conto(o,f,w,v,X,Y,M,N, (0.      ),-q,q); fprintf(o,".02 W 0 0 0 RGB S\n");
for(n=2;n<10;n+=2)conto(o,f,w,v,X,Y,M,N, (    .1*n),-p ,p); fprintf(o,".01 W 0 1 0 RGB S\n");
                 conto(o,f,w,v,X,Y,M,N, ( 1.     ),-q ,q); fprintf(o,".03 W 0 0 1 RGB S\n");
for(n=2;n<10;n+=2)conto(o,f,w,v,X,Y,M,N, ( 1.+.1*n),-p ,p); fprintf(o,".01 W 0 1 0 RGB S\n"); 
                conto(o,f,w,v,X,Y,M,N, ( 2.     ),-q ,q);  fprintf(o,".03 W 0 0 1 RGB S\n");
for(n=2;n<10;n+=2)conto(o,f,w,v,X,Y,M,N, ( 2.+.1*n),-p, p); fprintf(o,".01 W 0 1 0 RGB S\n");
                 conto(o,f,w,v,X,Y,M,N, ( M_PI   ),-q,q); fprintf(o,".04 W 1 0 1 RGB S\n");
                 conto(o,f,w,v,X,Y,M,N, (3*M_PI  ),-q,q); fprintf(o,".04 W 1 0 1 RGB S\n");
//
                 conto(o,g,w,v,X,Y,M,N, (-2.     ),-p,p); fprintf(o,".03 W 0 0 0 RGB S\n");
for(n=2;n<10;n+=2)conto(o,g,w,v,X,Y,M,N, (-2.+.1*n),-p,p); fprintf(o,".01 W 1 0 0 RGB S\n");
                 conto(o,g,w,v,X,Y,M,N, (-1.     ),-q,q); fprintf(o,".03 W 0 0 0 RGB S\n");
for(n=2;n<10;n+=2)conto(o,g,w,v,X,Y,M,N, (-1.+.1*n),-p,p); fprintf(o,".01 W 1 0 0 RGB S\n");
                 conto(o,g,w,v,X,Y,M,N, (0.      ),-q,q); fprintf(o,".03 W 0 0 0 RGB S\n");
for(n=2;n<10;n+=2)conto(o,g,w,v,X,Y,M,N, (    .1*n),-p,p); fprintf(o,".01 W 0 0 1 RGB S\n");
                 conto(o,g,w,v,X,Y,M,N, ( 1.     ),-q,q); fprintf(o,".03 W 0 0 0 RGB S\n");
for(n=2;n<10;n+=2)conto(o,g,w,v,X,Y,M,N, ( 1.+.1*n),-p,p); fprintf(o,".01 W 0 0 1 RGB S\n");
                 conto(o,g,w,v,X,Y,M,N, ( 2.     ),-p,p); fprintf(o,".03 W 0 0 0 RGB S\n");
                 conto(o,g,w,v,X,Y,M,N, ( 3.     ),-p,p); fprintf(o,".03 W 0 0 0 RGB S\n");
                 conto(o,g,w,v,X,Y,M,N, ( 4.     ),-p,p); fprintf(o,".03 W 0 0 0 RGB S\n");
//
//M(-10,0)L(-2,0)fprintf(o,".04 W 1 0 1 RGB S\n");
fprintf(o,"showpage\n\%\%\%Trailer"); fclose(o);
//system( "ggv TetrationPolynom25Terms.eps");    // for linux
//system(  "open TetrationPolynom25Terms.eps"); // for macintosh
system("ps2pdf TetrationPolynom25terms.eps");  
//getchar(); system("killall Preview"); //for macintosh
}
//
// end of plotting program
//