#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define pi 3.14159265358979323844

#define nt 60

FILE *psfile;
char fontname[30];
long IDUM=-1234;double ran2(void);

void psheader(double xmin,double ymin,double xmax,double ymax);
void psfooter(void);
void rectdraw(double xcenter,double ycenter,
		  double width,double height,double angle);
void rectdrawfill(double xcenter,double ycenter,
		  double width,double height,double angle);
void linedraw(double xa,double ya,double xb, double yb);
void textwrite(double xx, double yy,double fontsize,double angle,char *s);
void spheredraw(double xx,double yy,double radius,int it,int nsppheredraw);
void circledraw(double radius,double xx,double yy);
void coildraw(int ncoils,double radius,double xxi,double yyi,
	      double xxf,double yyf);
void xyfind(double xx,double yy,double zz,double *x0,double *y0);
void rgb(double red,double green,double blue);
void shell_sort(int n);
void euler(int ip);

int main(){
  double x,y,fontsize=24.0,angle;
  double xtext,ytext,delxtext,delytext,theta;
  double redchar[20],greenchar[20],bluechar[20];
  double amag[20],bmag[20],omega[20];
  double xchar[20],ychar[20];
  int ichar,it,nchars;
  char filename[30];
  char title[50];
  char letter[1];

  nchars=4;
  sprintf(title," LIZ");
  sprintf(fontname,"Courier-Bold");

  delxtext=fontsize/3.5;
  delytext=fontsize/3.5;
  x=280;
  y=550;
  for(ichar=1;ichar<nchars;ichar++){
    amag[ichar]=100.0*(1.0-2.0*ran2());
    bmag[ichar]=100.0*(1.0-2.0*ran2());
    xchar[ichar]=x;
    ychar[ichar]=y;
    x=x+fontsize;
PICK_OMEGA:
    omega[ichar]=(360.0/nt)*floor(3.0-5.0*ran2());
    if(fabs(omega[ichar])<0.001) goto PICK_OMEGA;
    redchar[ichar]=(floor(6.0*ran2()))/5.0;
    bluechar[ichar]=(floor(6.0*ran2()))/5.0;
    greenchar[ichar]=(floor(6.0*ran2()))/5.0;
  }
  for(it=0;it<nt;it++){
    sprintf(filename,"toy%02d.eps",it);
    psfile=fopen(filename,"w");
    psheader(100,400,550,700);
    rgb(0.0,0.0,0.0);
    for(ichar=1;ichar<nchars;ichar++){
      angle=omega[ichar]*it;
      x=xchar[ichar]+amag[ichar]*sin(2.0*pi*it/nt);
      y=ychar[ichar]+bmag[ichar]*sin(2.0*pi*it/nt);
      coildraw(6,0.3*fontsize,xchar[ichar],ychar[ichar],x,y);
      rgb(redchar[ichar],greenchar[ichar],bluechar[ichar]);
      rectdrawfill(x,y,fontsize,fontsize,angle);
      circledraw(0.25*fontsize,xchar[ichar],ychar[ichar]);
      rgb(0.0,0.0,0.0);
      sprintf(letter,"%c",title[ichar]);
      printf(letter,"%c",title[ichar]);
      theta=angle*pi/180.0;
      xtext=x-delxtext*cos(theta)+delytext*sin(theta);
      ytext=y-delytext*cos(theta)-delxtext*sin(theta);
      textwrite(xtext,ytext,fontsize,angle,letter);
    }
    printf("\n");
    psfooter();
    fclose(psfile);
  }
  return(0);
}


/* long IDUM=-123 */
/* double ran2(void); */

/* To use this routine, call the ran2 with the form */
/*   x=ran2();  */
/* and x will be a random number between 0 and 1, which won't repeat */
/* before every 10^18 numbers */
/* To initialize, call with IDUM being a NEGATIVE number */
/* This came out of the ran2 function from numerical recipes */

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double ran2() {
  long J;
  long K;
  static long IDUM2=123456789;
  static long IY=0;
  static long IV[NTAB];
  double TEMP;
  if(IDUM <=0){
    if(-(IDUM)<1) IDUM=1;
    else IDUM = -(IDUM);
    IDUM2=(IDUM);
    for (J=NTAB+7;J>=0;J--) {
      K=(IDUM)/IQ1;
      IDUM=IA1*(IDUM-K*IQ1)-K*IR1;
      if(IDUM<0) IDUM +=IM1;
      if(J<NTAB) IV[J]=IDUM;
    }
    IY=IV[0];
    //printf("Initialization Completed %d\n",IDUM);
  }
  K=(IDUM)/IQ1;
  IDUM=IA1*(IDUM-K*IQ1)-K*IR1;
  if(IDUM<0) IDUM+=IM1;
  K=IDUM2/IQ2;
  IDUM2=IA2*(IDUM2-K*IQ2)-K*IR2;
  if(IDUM2<0) IDUM2+=IM2;
  J=IY/NDIV;
  IY=IV[J]-IDUM2;
  IV[J]=IDUM;
  if(IY<1) IY+=IMM1;
  TEMP=AM*IY;
  if(TEMP>RNMX) return RNMX;
  else return TEMP;
}
double gauss(){
  double x,y,r2,r,c,randy;
TRY_AGAIN:
  x=1.0-2.0*ran2();
  y=1.0-2.0*ran2();
  r2=x*x+y*y;
  if(r2>1.0) goto TRY_AGAIN;
  r=sqrt(r2);
  c=x/r;
  randy=c*sqrt(-2.0*log(r2));
  return randy;
}

void psheader(double xmin,double ymin,double xmax,double ymax){
  fprintf(psfile,"%%!PS-Adobe-2.0 EPSF-2.0 \n");
  fprintf(psfile,
	  "%%%%BoundingBox: %5.1f %5.1f %5.1f %5.1f \n",xmin,ymin,xmax,ymax);
  fprintf(psfile,"%%%%Creator: S. Pratt, M.S.U. \n");
  fprintf(psfile,"%%%%EndComments \n");
}

void psfooter(void){
  fprintf(psfile,"showpage\n");
}

void xyfind(double xx,double yy,double zz,double *x0,double *y0){
  double factor;
  double xperspective,yperspective,zperspective;
  factor=-zperspective/(-zperspective+zz);
  *x0=factor*xx+(1.0-factor)*xperspective;
  *y0=factor*yy+(1.0-factor)*yperspective;
}

void rgb(double red,double green,double blue){
  fprintf(psfile,"%4.2f %4.2f %4.2f setrgbcolor \n",red,green,blue);
}

void linedraw(double xa,double ya,double xb,double yb){
  fprintf(psfile,"%6.2f %6.2f moveto %6.2f %6.2f lineto stroke \n",xa,ya,xb,yb);
}

void rectdrawfill(double xcenter,double ycenter,double width,double height,double angle){
  double xa,ya,xb,yb,xc,yc,xd,yd;
  double delx1,delx2,dely1,dely2;
  if(fabs(angle)<0.001){
    xa=xcenter-0.5*width;
    xb=xcenter+0.5*width;
    xc=xb;
    xd=xa;
    ya=ycenter-0.5*height;
    yc=ycenter+0.5*height;
    yb=ya;
    yd=yc;
  }
  else{
    angle=angle*pi/180.0;
    delx1=0.5*width*cos(angle)-0.5*height*sin(angle);
    dely1=0.5*width*sin(angle)+0.5*height*cos(angle);
    delx2=0.5*width*cos(angle)+0.5*height*sin(angle);
    dely2=-0.5*width*sin(angle)+0.5*height*cos(angle);
    xa=xcenter-delx1;
    xc=xcenter+delx1;
    xb=xcenter+delx2;
    xd=xcenter-delx2;
    ya=ycenter-dely1;
    yc=ycenter+dely1;
    yb=ycenter-dely2;
    yd=ycenter+dely2;
  }
  
  fprintf(psfile,"newpath \n%6.2f %6.2f moveto %6.2f %6.2f lineto \n",
	  xa,ya,xb,yb);
fprintf(psfile,"%6.2f %6.2f lineto %6.2f %6.2f lineto \nclosepath fill \n",
xc,yc,xd,yd);
}

void rectdraw(double xcenter,double ycenter,double width,double height,double angle){
  double xa,ya,xb,yb,xc,yc,xd,yd;
  double delx1,delx2,dely1,dely2;
  if(fabs(angle)<0.001){
    xa=xcenter-0.5*width;
    xb=xcenter+0.5*width;
    xc=xb;
    xd=xa;
    ya=ycenter-0.5*height;
    yc=ycenter+0.5*height;
    yb=ya;
    yd=yc;
  }
  else{
    angle=angle*pi/180.0;
    delx1=0.5*width*cos(angle)-0.5*height*sin(angle);
    dely1=0.5*width*sin(angle)+0.5*height*cos(angle);
    delx2=0.5*width*cos(angle)+0.5*height*sin(angle);
    dely2=-0.5*width*sin(angle)+0.5*height*cos(angle);
    xa=xcenter-delx1;
    xc=xcenter+delx1;
    xb=xcenter+delx2;
    xd=xcenter-delx2;
    ya=ycenter-dely1;
    yc=ycenter+dely1;
    yb=ycenter-dely2;
    yd=ycenter+dely2;
  }
  
  fprintf(psfile,"newpath \n%6.2f %6.2f moveto %6.2f %6.2f lineto \n",
	  xa,ya,xb,yb);
fprintf(psfile,"%6.2f %6.2f lineto %6.2f %6.2f lineto \nclosepath stroke\n",
xc,yc,xd,yd);
}

void spheredraw(double xx,double yy,double radius,int it,int nspheredraw){
  int i;
  double theta0,thetaf,graylevel,graylevel0;
  double r,r2,yx,xci,yci;
  theta0=0.0;
  thetaf=360.0;
  for(i=1;i<=nspheredraw;i++){
    if(it==1)graylevel0=.9;
    if(it==2)graylevel0=.7;
    if(it==3)graylevel0=.8;
    graylevel=graylevel0;
    if(nspheredraw>1)
      graylevel=graylevel-graylevel0*0.3*(double)(i-1)/(double)(nspheredraw-1);
    r2=pow(radius,2)*(1.0-(double)(i-1)/(double)(nspheredraw));
    r=sqrt(r2);
    fprintf(psfile,"newpath \n");
    xci=xx+((double)(i-1)/(double)(nspheredraw))*r*0.15;
    yci=yy-((double)(i-1)/(double)(nspheredraw))*r*0.1;
    fprintf(psfile," %6.2f %6.2f %6.2f %6.2f %6.2f arc closepath \n",
	    xci,yci,r,theta0,thetaf);
    if(it==1) rgb(graylevel,0.3*graylevel,0.3*graylevel);
    if(it==2) rgb(0.3*graylevel,graylevel,0.3*graylevel);
    if(it==3) rgb(0.3*graylevel,graylevel,0.3*graylevel);
    fprintf(psfile," fill stroke \n");
  }
}

void textwrite(double xtext,double ytext, double fontsize,double angle,char *textstring){
  int n;
  fprintf(psfile,"/%s findfont %5.1f scalefont setfont ",fontname,fontsize);
  fprintf(psfile,"%6.2f %6.2f moveto \n",xtext,ytext);
  if(fabs(angle)>0.0001){
    n=(int)(angle/360.0);
    angle=angle-360.0*n;
    fprintf(psfile,"%6.1f rotate (%c) show %6.1f rotate \n",
	    angle,textstring[0],-angle);
  }
  else{
    fprintf(psfile,"(%s)show \n",textstring);
  }
}

void coildraw(int ncoils,double radius,double xxi,double yyi,double xxf,double yyf){
  int icoil;
  double x1,x2,x3,x4,y1,y2,y3,y4,xperp,yperp,rforward,delx,dely,length;
  delx=xxf-xxi;
  dely=yyf-yyi;
  length=sqrt(0.0000001+pow(delx,2)+pow(dely,2));
  if(fabs(length)>0.001){
    delx=delx/length;
    dely=dely/length;
    xperp=-dely*radius;
    yperp=delx*radius;
    rforward=length/(0.25*ncoils+0.75);
    delx=delx*rforward;
    dely=dely*rforward;
    /* print first half semi-circle */
    x1=xxi;
    y1=yyi;
    x4=x1+delx;
    y4=y1+dely;
    x2=x1+xperp;
    y2=y1+yperp;
    x3=x4+xperp;
    y3=y4+yperp;
    fprintf(psfile,"%6.1f %6.1f moveto \n",x1,y1);
    fprintf(psfile,"%6.1f %6.1f %6.1f %6.1f %6.1f %6.1f curveto \n",
	    x2,y2,x3,y3,x4,y4);
    fprintf(psfile,"stroke \n");
    for(icoil=1;icoil<ncoils;icoil++){
      /* print backward half coil */
      x1=x4;
      y1=y4;
      x4=x1-0.75*delx;
      y4=y1-0.75*dely;
      x2=x1-xperp;
      y2=y1-yperp;
      x3=x4-xperp;
      y3=y4-yperp;
      fprintf(psfile,"%6.1f %6.1f moveto \n",x1,y1);
      fprintf(psfile,"%6.1f %6.1f %6.1f %6.1f %6.1f %6.1f curveto \n",
	      x2,y2,x3,y3,x4,y4);
      fprintf(psfile,"stroke \n");
      /* print forward half coil */
      x1=x4;
      y1=y4;
      x4=x1+delx;
      y4=y1+dely;
      x2=x1+xperp;
      y2=y1+yperp;
      x3=x4+xperp;
      y3=y4+yperp;
      fprintf(psfile,"%6.1f %6.1f moveto \n",x1,y1);
      fprintf(psfile,"%6.1f %6.1f %6.1f %6.1f %6.1f %6.1f curveto \n",
	      x2,y2,x3,y3,x4,y4);
      fprintf(psfile,"stroke \n");
    }
  }
}

void circledraw(double radius,double x,double y){
  double theta0=0.0;
  double thetaf=360.0;
  fprintf(psfile,"newpath \n");
  fprintf(psfile," %6.2f %6.2f %6.2f %6.2f %6.2f arc closepath \n",
	  x,y,radius,theta0,thetaf);
  fprintf(psfile,"fill\n");
}