File:TetSheldonTestImaT.png

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DB double
using namespace std;
#include <complex>
typedef complex<double> z_type;
//#include <complex.h>
//#define z_type complex<double>
#define Re(x) (x).real()
#define Im(x) (x).imag()
#define I z_type(0.,1.)
#define DO(x,y) for(x=0;x<y;x++)

#include "ado.cin"
#include "filog.cin"

#define M(x,y) fprintf(o,"%5.3f %5.3f M\n",1.*(x),1.*(y));
#define L(x,y) fprintf(o,"%5.3f %5.3f L\n",1.*(x),1.*(y));
#define o(x,y) fprintf(o,"%5.3f %5.3f o\n",1.*(x),1.*(y));

int main(){ int J,j,k,m,n; DB x,y, u, t; z_type z,c,d, cu,cd;
#include "GLxw2048.inc"

z_type b=z_type( 1.5259833851700000, 0.0178411853321000);
z_type a=log(b);
z_type Zo=Filog(a);
z_type Zc=conj(Filog(conj(a)));

int K=NPO; DB A=32.; printf("K=%3d A=%3.1f\n",K,A);
z_type E[2048],F[2048],G[2048],H[2048];

FILE *o;o=fopen("TetSheldonIma.eps","w");ado(o,2200,450);
fprintf(o,"1100 201 translate\n 100 100 scale\n");
for(j=-10;j<11;j+=1){M(j,-2)L(j,2);}
M(-10 , 2)L(10 , 2);
M(-10 , 1)L(10 , 1);
M(-10 ,-1)L(10 ,-1);
M(-10 ,-2)L(10 ,-2);
  fprintf(o,".006 W S\n");
M(-10.1,0)L(10.1,0); fprintf(o,".02 W S\n");
fprintf(o,".01 W S\n 1 setlinejoin 1\n");

DO(n,K){y=GLx[n]*A;
        if(y<-2.4) E[n]=F[n]=G[n]=Zc;
        else {if(y>2.4) E[n]=F[n]=G[n]=Zo;
                else { F[n]=c=1. ;//f3(z_type(0.,y));
                        E[n]=log(c)/a;
                        G[n]=exp(a*c);
                     }
             }
        }

DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Re(F[k]);if(j==0)M(y,u)else L(y,u)}
fprintf(o,".08 W 1 .2 1 RGB S\n");
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Im(F[k]);if(j==0)M(y,u)else L(y,u)}
//DO(k,K-512){y=GLx[k]*A; u=Im(F[k]);if(k==0)M(y,u)else L(y,u)}
fprintf(o,".08 W 1 .6 1 RGB S\n");

for(n=0;n<K;n+=2)
        { y=GLx[n]*A; z=z_type(0.,y); c=0.; //printf(" %3d",n);
         DO(k,K){t=A*GLx[k]; c+= GLw[k]*( G[k]/(z_type( 1.,t)-z) - E[k]/(z_type(-1.,t)-z) );}
         cu=.5-I/(2.*M_PI)*log( (z_type(1.,-A)+z)/(z_type(1., A)-z) );
         cd=.5-I/(2.*M_PI)*log( (z_type(1.,-A)-z)/(z_type(1., A)+z) );
         c=c*(A/(2.*M_PI)) +Zo*cu+Zc*cd;
         if(abs(n-K/2)<4)printf("%3d %6.3f %9.6f %9.6f %9.6f %9.6f\n",n,y,Re(F[n]),Im(F[n]),Re(c),Im(c));
         E[n]=log(c)/a; F[n]=c; G[n]=exp(c*a);}
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Re(F[k]);if(j==0)M(y,u)else L(y,u)}
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Im(F[k]);if(j==0)M(y,u)else L(y,u)}
//DO(k,K){y=GLx[k]*A; u=Re(F[k]);if(k==0)M(y,u)else L(y,u)} fprintf(o,".02 W .7 0 0 RGB S\n");
//DO(k,K){y=GLx[k]*A; u=Im(F[k]);if(k==0)M(y,u)else L(y,u)}
fprintf(o,".01 W 1 0 0 RGB S\n");

for(n=K-1;n>0;n-=2)
        { y=GLx[n]*A; z=z_type(0.,y); c=0.; //printf(" %3d",n);
         DO(k,K){t=A*GLx[k]; c+= GLw[k]*( G[k]/(z_type( 1.,t)-z) - E[k]/(z_type(-1.,t)-z) );}
         cu=.5-I/(2.*M_PI)*log( (z_type(1.,-A)+z)/(z_type(1., A)-z) );
         cd=.5-I/(2.*M_PI)*log( (z_type(1.,-A)-z)/(z_type(1., A)+z) );
         c=c*(A/(2.*M_PI)) +Zo*cu+Zc*cd;
         if(abs(n-K/2)<4)printf("%3d %6.3f %9.6f %9.6f %9.6f %9.6f\n",n,y,Re(F[n]),Im(F[n]),Re(c),Im(c));
         E[n]=log(c)/a; F[n]=c; G[n]=exp(c*a); }
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Re(F[k]);if(j==0)M(y,u)else L(y,u)}
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Im(F[k]);if(j==0)M(y,u)else L(y,u)}
//DO(k,K){y=GLx[k]*A; u=Re(F[k]);if(k==0)M(y,u)else L(y,u)} fprintf(o,".02 W 0 .6 0 RGB S\n");
//DO(k,K){y=GLx[k]*A; u=Im(F[k]);if(k==0)M(y,u)else L(y,u)}
fprintf(o,".02 W 0 1 0 RGB S\n");

for(n=0;n<K;n+=2)
        { y=GLx[n]*A; z=z_type(0.,y); c=0.; //printf(" %3d",n);
         DO(k,K){t=A*GLx[k]; c+= GLw[k]*( G[k]/(z_type( 1.,t)-z) - E[k]/(z_type(-1.,t)-z) );}
         cu=.5-I/(2.*M_PI)*log( (z_type(1.,-A)+z)/(z_type(1., A)-z) );
         cd=.5-I/(2.*M_PI)*log( (z_type(1.,-A)-z)/(z_type(1., A)+z) );
         c=c*(A/(2.*M_PI)) +Zo*cu+Zc*cd;
         if(abs(n-K/2)<4)printf("%3d %6.3f %9.6f %9.6f %9.6f %9.6f\n",n,y,Re(F[n]),Im(F[n]),Re(c),Im(c));
         E[n]=log(c)/a; F[n]=c; G[n]=exp(c*a); }
//DO(k,K){y=GLx[k]*A; u=Re(F[k]);if(k==0)M(y,u)else L(y,u)} fprintf(o,".01 W 0 0 .7 RGB S\n");
//DO(k,K){y=GLx[k]*A; u=Im(F[k]);if(k==0)M(y,u)else L(y,u)}
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Re(F[k]);if(j==0)M(y,u)else L(y,u)}
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Im(F[k]);if(j==0)M(y,u)else L(y,u)}
fprintf(o,".01 W 0 0 1 RGB S\n");

for(n=K-1;n>0;n-=2)
        { y=GLx[n]*A; z=z_type(0.,y); c=0.; //printf(" %3d",n);
         DO(k,K){t=A*GLx[k]; c+= GLw[k]*( G[k]/(z_type( 1.,t)-z) - E[k]/(z_type(-1.,t)-z) );}
         cu=.5-I/(2.*M_PI)*log( (z_type(1.,-A)+z)/(z_type(1., A)-z) );
         cd=.5-I/(2.*M_PI)*log( (z_type(1.,-A)-z)/(z_type(1., A)+z) );
         c=c*(A/(2.*M_PI)) +Zo*cu+Zc*cd;
         if(abs(n-K/2)<4)printf("%3d %6.3f %9.6f %9.6f %9.6f %9.6f\n",n,y,Re(F[n]),Im(F[n]),Re(c),Im(c));
         E[n]=log(c)/a; F[n]=c; G[n]=exp(c*a);}
//DO(k,K){y=GLx[k]*A; u=Re(F[k]);if(k==0)M(y,u)else L(y,u)} fprintf(o,".01 W 0 0 0 RGB S\n");
//DO(k,K){y=GLx[k]*A; u=Im(F[k]);if(k==0)M(y,u)else L(y,u)} fprintf(o,".01 W 0 0 0 RGB S\n");

DO(j,4)
{
for(n=0;n<K;n+=2)
        { y=GLx[n]*A; z=z_type(0.,y); c=0.; //printf(" %3d",n);
         DO(k,K){t=A*GLx[k]; c+= GLw[k]*( G[k]/(z_type( 1.,t)-z) - E[k]/(z_type(-1.,t)-z) );}
         cu=.5-I/(2.*M_PI)*log( (z_type(1.,-A)+z)/(z_type(1., A)-z) );
         cd=.5-I/(2.*M_PI)*log( (z_type(1.,-A)-z)/(z_type(1., A)+z) );
         c=c*(A/(2.*M_PI)) +Zo*cu+Zc*cd;
         if(abs(n-K/2)<2)printf("%3d %6.3f %9.6f %9.6f %9.6f %9.6f\n",n,y,Re(F[n]),Im(F[n]),Re(c),Im(c));
         E[n]=log(c)/a; F[n]=c; G[n]=exp(c*a); }
for(n=K-1;n>0;n-=2)
        { y=GLx[n]*A; z=z_type(0.,y); c=0.; //printf(" %3d",n);
         DO(k,K){t=A*GLx[k]; c+= GLw[k]*( G[k]/(z_type( 1.,t)-z) - E[k]/(z_type(-1.,t)-z) );}
         cu=.5-I/(2.*M_PI)*log( (z_type(1.,-A)+z)/(z_type(1., A)-z) );
         cd=.5-I/(2.*M_PI)*log( (z_type(1.,-A)-z)/(z_type(1., A)+z) );
         c=c*(A/(2.*M_PI)) +Zo*cu+Zc*cd;
         if(abs(n-K/2)<2)printf("%3d %6.3f %9.6f %9.6f %9.6f %9.6f\n",n,y,Re(F[n]),Im(F[n]),Re(c),Im(c));
         E[n]=log(c)/a; F[n]=c; G[n]=exp(c*a);}
}
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Re(F[k]);if(j==0)M(y,u)else L(y,u)}
DO(j,440){k=K/2+j-220; y=GLx[k]*A; u=Im(F[k]);if(j==0)M(y,u)else L(y,u)}
//DO(k,K){y=GLx[k]*A; u=Re(F[k]);if(k==0)M(y,u)else L(y,u)} fprintf(o,".01 W 0 0 0 RGB S\n");
//DO(k,K){y=GLx[k]*A; u=Im(F[k]);if(k==0)M(y,u)else L(y,u)}
fprintf(o,".01 W 0 0 0 RGB S\n");
fprintf(o,"showpage\n\%c%cTrailer",'%','%'); fclose(o);
system("epstopdf TetSheldonIma.eps");
system( "open TetSheldonIma.pdf");

DO(j,40)
{
for(n=0;n<K;n+=2)
        { y=GLx[n]*A; z=z_type(0.,y); c=0.; //printf(" %3d",n);
         DO(k,K){t=A*GLx[k]; c+= GLw[k]*( G[k]/(z_type( 1.,t)-z) - E[k]/(z_type(-1.,t)-z) );}
         cu=.5-I/(2.*M_PI)*log( (z_type(1.,-A)+z)/(z_type(1., A)-z) );
         cd=.5-I/(2.*M_PI)*log( (z_type(1.,-A)-z)/(z_type(1., A)+z) );
         c=c*(A/(2.*M_PI)) +Zo*cu+Zc*cd;
         if(abs(n-K/2)<2)printf("%3d %6.3f %9.6f %9.6f %9.6f %9.6f\n",n,y,Re(F[n]),Im(F[n]),Re(c),Im(c));
         E[n]=log(c)/a; F[n]=c; G[n]=exp(c*a); }
for(n=K-1;n>0;n-=2)
        { y=GLx[n]*A; z=z_type(0.,y); c=0.; //printf(" %3d",n);
         DO(k,K){t=A*GLx[k]; c+= GLw[k]*( G[k]/(z_type( 1.,t)-z) - E[k]/(z_type(-1.,t)-z) );}
         cu=.5-I/(2.*M_PI)*log( (z_type(1.,-A)+z)/(z_type(1., A)-z) );
         cd=.5-I/(2.*M_PI)*log( (z_type(1.,-A)-z)/(z_type(1., A)+z) );
         c=c*(A/(2.*M_PI)) +Zo*cu+Zc*cd;
         if(abs(n-K/2)<2)printf("%3d %6.3f %9.6f %9.6f %9.6f %9.6f\n",n,y,Re(F[n]),Im(F[n]),Re(c),Im(c));
         E[n]=log(c)/a; F[n]=c; G[n]=exp(c*a);}
printf("%3d %19.16lf %19.16lf %19.16lf %19.16lf\n",j,Re(F[1023]),Im(F[1023]),Re(F[1024]),Im(F[1024]) );
}
o=fopen("TetSheldonIma.inc","w");
fprintf(o,"z_type F[%4d]={\n",K);
DO(k,K-1) fprintf(o,"z_type(%19.16lf,%19.16lf),\n",Re(F[k]),Im(F[k]));
          fprintf(o,"z_type(%19.16lf,%19.16lf)};\n",Re(F[K-1]),Im(F[K-1]));
fclose(o);
}