#include "pssub.h"

//papersize
double xsize, ysize;

//viewpoint
double xv1, yv1, xv2, yv2;

//world coordinate
double xw1, yw1, xw2, yw2;

FILE* stream;

void annulus(double x1,double y1,double r1,double r2)
{
    plot(x1+r2, y1, 3);
    arc(x1, y1, r2, 1.0, 360.0);
    plot(x1+r1, y1, 2);
    arcn(x1, y1, r1, 360.0, 0.0);
    plot(x1+r2, y1, 2);
    closepath();
}

void arc(double x1, double y1, double r1, double t1, double t2)
{ 
   double xx1,yy1,rr1;
   
   xx1=((x1-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize; 
   yy1=((y1-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize; 
   rr1=(r1/(xw2-xw1)*(xv2-xv1))*xsize; 
   fprintf(stream,"%lf cm %lf cm %lf cm %lf %lf arc \n",xx1,yy1,rr1,t1,t2);
   return;
}

void arcn(double x1, double y1, double r1, double t1, double t2)
{ 
   double xx1,yy1,rr1;
   
   xx1=((x1-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize; 
   yy1=((y1-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize; 
   rr1=(r1/(xw2-xw1)*(xv2-xv1))*xsize; 
   fprintf(stream,"%lf cm %lf cm %lf cm %lf %lf arcn \n",xx1,yy1,rr1,t1,t2);
   return;
}

void arcto(double x1, double y1, double x2, double y2, double r1)
{ 
   double xx1,yy1,xx2,yy2,rr1;
   
   xx1=((x1-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize; 
   yy1=((y1-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize; 
   xx2=((x2-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize; 
   yy2=((y2-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize; 
   rr1=(r1/(xw2-xw1)*(xv2-xv1))*xsize; 
   fprintf(stream,"%lf cm %lf cm %lf cm %lf cm %lf cm arcto \n",xx1,yy1,xx2,yy2,rr1);
   return;
}

void arctorot(double x1, double y1, double x2, double y2, double r1, double t1)
{ 
   double xx1,yy1,xx2,yy2,rr1,tt1;
     tt1 = t1*3.14159/180.0;
     xx1 = x1*cos(tt1)-y1*sin(tt1);
     yy1 = x1*sin(tt1)+y1*cos(tt1);
     xx2 = x2*cos(tt1)-y2*sin(tt1);
     yy2 = x2*sin(tt1)+y2*cos(tt1);
     xx1=((xx1-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize; 
     yy1=((yy1-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize; 
     xx2=((xx2-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize; 
     yy2=((yy2-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize; 
     rr1=(r1/(xw2-xw1)*(xv2-xv1))*xsize; 
     fprintf(stream,"%lf cm %lf cm %lf cm %lf cm %lf cm arcto \n",xx1,yy1,xx2,yy2,rr1);
   return;
}

void arrow(double x1, double y1, double x2, double y2, double d)
{
   double alpha=2.2,s,x3,y3,x4,y4;
   s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); 
   x3=x2-(x2-x1)/s*d*alpha-(y2-y1)/s*d;
   y3=y2-(y2-y1)/s*d*alpha+(x2-x1)/s*d;
   x4=x2-(x2-x1)/s*d*alpha+(y2-y1)/s*d;
   y4=y2-(y2-y1)/s*d*alpha-(x2-x1)/s*d;
   plot(x1, y1, 3); plot(x2, y2, 2); 
   plot(x3, y3, 3); plot(x2, y2, 2); plot(x4, y4, 2); 
   return; 
}

void arrowa(double x1, double y1, double x2, double y2, double d)
{
   double alpha=2.2,s,x3,y3,x4,y4;
   s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); 
   x3=x2-(x2-x1)/s*d*alpha-(y2-y1)/s*d;
   y3=y2-(y2-y1)/s*d*alpha+(x2-x1)/s*d;
   x4=x2-(x2-x1)/s*d*alpha+(y2-y1)/s*d;
   y4=y2-(y2-y1)/s*d*alpha-(x2-x1)/s*d;
   plot(x1, y1, 3); plot((x3+x4)/2.0, (y3+y4)/2.0, 2); stroke();
   plot(x3, y3, 3); plot(x2, y2, 2); plot(x4, y4, 2); 
   closepath();fill();
   return; 
}

void arrowb(double x1, double y1, double x2, double y2, double d)
{
   double alpha=2.2,s,x3,y3,x4,y4,x32,y32,x24,y24,x43,y43;
   s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); 
   x3=x2-(x2-x1)/s*d*alpha-(y2-y1)/s*d;
   y3=y2-(y2-y1)/s*d*alpha+(x2-x1)/s*d;
   x4=x2-(x2-x1)/s*d*alpha+(y2-y1)/s*d;
   y4=y2-(y2-y1)/s*d*alpha-(x2-x1)/s*d;
   x32=(x3+x2)/2.0-(y3-y2)/12.0;
   y32=(y3+y2)/2.0+(x3-x2)/12.0;
   x24=(x2+x4)/2.0-(y2-y4)/12.0;
   y24=(y2+y4)/2.0+(x2-x4)/12.0;
   x43=(x4+x3)/2.0-(y4-y3)/12.0;
   y43=(y4+y3)/2.0+(x4-x3)/12.0;
   plot(x1, y1, 3); plot(x43, y43, 2); stroke();
   plot(x3, y3, 3); plot(x32, y32, 2); plot(x2, y2, 2);
   plot(x24, y24, 2); plot(x4, y4, 2); plot(x43,y43,2);
   closepath();fill();
   return; 
}

void arrowc(double x1, double y1, double x2, double y2, double d)
{
   double alpha=2.2,s,x3,y3,x4,y4,x5,y5;
   s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); 
   x3=x2-(x2-x1)/s*d*alpha-(y2-y1)/s*d;
   y3=y2-(y2-y1)/s*d*alpha+(x2-x1)/s*d;
   x4=x2-(x2-x1)/s*d*alpha+(y2-y1)/s*d;
   y4=y2-(y2-y1)/s*d*alpha-(x2-x1)/s*d;
   x5=x2-(x2-x1)/s*d*alpha*0.8;
   y5=y2-(y2-y1)/s*d*alpha*0.8;
   plot(x1, y1, 3); plot(x5, y5, 2); stroke();
   plot(x3, y3, 3); plot(x2, y2, 2); plot(x4, y4, 2); plot(x5, y5, 2);
   closepath();fill();
   return; 
}

void arrowfill(double x1, double y1, double x2, double y2, double d)
{
   double alpha=2.2,s,x3,y3,x4,y4,x5,y5;
   s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
   x3=x2-(x2-x1)/s*d*alpha-(y2-y1)/s*d;
   y3=y2-(y2-y1)/s*d*alpha+(x2-x1)/s*d;
   x4=x2-(x2-x1)/s*d*alpha+(y2-y1)/s*d;
   y4=y2-(y2-y1)/s*d*alpha-(x2-x1)/s*d;
   plot(x1, y1, 3); plot(x2, y2, 2); stroke();
   plot(x3, y3, 3); plot(x2, y2, 2); plot(x4, y4, 2); closepath(); 
   fill();
   return;
}

void arrowrot(double x1, double y1, double x2, double y2, double t1, double d)
{
   double alpha=2.2,s,x3,y3,x4,y4;
   s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); 
   x3=x2-(x2-x1)/s*d*alpha-(y2-y1)/s*d;
   y3=y2-(y2-y1)/s*d*alpha+(x2-x1)/s*d;
   x4=x2-(x2-x1)/s*d*alpha+(y2-y1)/s*d;
   y4=y2-(y2-y1)/s*d*alpha-(x2-x1)/s*d;
   plotrot(x1, y1, t1, 3); plotrot(x2, y2, t1, 2); 
   plotrot(x3, y3, t1, 3); plotrot(x2, y2, t1, 2); plotrot(x4, y4, t1, 2); 
   return; 
}

void arrowwide(double x1,double y1, double x2,double y2, double d, double al)
{
   double s,dx,dy,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7,x8,y8;
   s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); 
   dx=al; dy=al;
   x3=x2-(x2-x1)*d-(y2-y1)/s*dx;
   y3=y2-(y2-y1)*d+(x2-x1)/s*dy;
   x4=x2-(x2-x1)*d+(y2-y1)/s*dx;
   y4=y2-(y2-y1)*d-(x2-x1)/s*dy;
   x7=x2-(x2-x1)*d-(y2-y1)/s*dx*0.4;
   y7=y2-(y2-y1)*d+(x2-x1)/s*dy*0.4;
   x8=x2-(x2-x1)*d+(y2-y1)/s*dx*0.4;
   y8=y2-(y2-y1)*d-(x2-x1)/s*dy*0.4;
   x5=x1-(y2-y1)/s*dx*0.4;
   y5=y1+(x2-x1)/s*dy*0.4;
   x6=x1+(y2-y1)/s*dx*0.4;
   y6=y1-(x2-x1)/s*dy*0.4;
   plot(x6, y6, 3); plot(x8, y8, 2); plot(x4, y4, 2);
   plot(x2, y2, 2); plot(x3, y3, 2); plot(x7, y7, 2); plot(x5, y5, 2); 
   closepath(); stroke();
  return;
}

void battery(double x1, double y1, double x2, double y2, double d1, double d2)
{
    double s,dx,dy,ex,ey,xx,yy; 
  
    s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); 
     dx=x2-x1; dy=y2-y1; 
     ex=-dy; ey=dx; 
     xx=x1+d1*ex; yy=y1+d1*ey; plot(xx, yy, 3); 
     xx=x1-d1*ex; yy=y1-d1*ey; plot(xx, yy, 2); 
     xx=x2+d2*ex; yy=y2+d2*ey; plot(xx, yy, 3); 
     xx=x2-d2*ex; yy=y2-d2*ey; plot(xx, yy, 2); 
    return;
}

void brokenlines(double *x1, double *y1, int n)
{
	 int i;
	 newpath();
     plot(x1[0],y1[0],3); 
     for (i=1; i<n;i++) {
     plot(x1[i], y1[i], 2); 
      }
    return;
}

void circ(double x1, double y1, double r1)
{ 
  double xx1,yy1,rr1;
  xx1 = ((x1-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
  yy1 = ((y1-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
  rr1 = r1/(xw2-xw1)*(xv2-xv1)*xsize;
  newpath();
  fprintf(stream, "%lf cm %lf cm %lf cm 0 360 arc\n", xx1, yy1, rr1);
  return;
}

void circn(double x1, double y1, double r1)
{ 
  double xx1,yy1,rr1;
  xx1 = ((x1-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
  yy1 = ((y1-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
  rr1 = r1/(xw2-xw1)*(xv2-xv1)*xsize;
  newpath();
  fprintf(stream, "%lf cm %lf cm %lf cm 0 360 arcn\n", xx1, yy1, rr1);
  return;
}

void clipoff(void)
{
  fprintf(stream, "clip\n");
}

void clipon(void)
{ 
  fprintf(stream, "clip\n");
}

void cliponrec(double x1, double y1, double x2, double y2)
{ 
  newpath();
  plot(x1, y1, 3);plot(x2, y1, 2);plot(x2, y2,2); plot(x1, y2, 2);
  closepath();
  fprintf(stream, "clip\n");
}

void closepath(void)
{
  fprintf(stream, "closepath\n"); return;
}

void coil(double x1, double y1, double x2, double y2, double d, int n)
{
   double s,ex,ey,ds,a,xx,yy,t;
   int i;

   s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
   ex=(x2-x1)/s;
   ey=(y2-y1)/s;
   a=2.0; ds=(s-2.0*d)/((double) n +0.5);
     plot(x1, y1, 3);
     for(i=1;i<=(40*n+20);i++) {
      t=2.0*3.14159*(double) i/40.0; 
      xx=x1+((d-d*cos(t))+t/2.0/3.14159*ds)*ex-a*d*sin(t)*ey;
      yy=y1+((d-d*cos(t))+t/2.0/3.14159*ds)*ey+a*d*sin(t)*ex;
      plot(xx, yy, 2);
      } 
    return;
}

void curvto(double x1, double y1, double x2, double y2, double x3, double y3)
{
      double xx1,yy1,xx2,yy2,xx3,yy3;

     xx1=((x1-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize;
     yy1=((y1-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize;
     xx2=((x2-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize;
     yy2=((y2-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize;
     xx3=((x3-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize;
     yy3=((y3-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize;
     fprintf(stream, "%lf cm %lf cm %lf cm %lf cm %lf cm %lf cm curveto\n",xx1,yy1,xx2,yy2,xx3,yy3); 
     return;
}

void curvton(double * x1, double * y1, int n)
{
      double xx1[100],yy1[100];int i;

     for (i=0;i<n;i++) {
		xx1[i]=((x1[i]-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize;
        yy1[i]=((y1[i]-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize; 
        fprintf(stream, "%lf cm %lf cm ",xx1[i],yy1[i]); }
     fprintf(stream, " curveto\n");
     return;
}

void curvtona(double * x1, double * y1, double t, int n)
{
     double x12[100],y12[100],x21[100],y21[100];int i;

     for (i=0;i<n-1;i++) {
		x12[i]=t*x1[i]+(1.0-t)*x1[i+1];
		x21[i]=(1.0-t)*x1[i]+t*x1[i+1];
		y12[i]=t*y1[i]+(1.0-t)*y1[i+1];
		y21[i]=(1.0-t)*y1[i]+t*y1[i+1];}
     plot(x21[0],y21[0], 2); 
     for (i=0;i<n-2;i++) {
		curvto(x21[i],y21[i],x1[i+1],y1[i+1],x12[i+1],y12[i+1]); }
	 plot(x1[n-1],y1[n-1],2);
     return;
}

void dims(double x1, double y1, double x2, double y2,double d)
{     double x3,y3,x4,y4,s,t; 
      int i,n;
       n=10;
       s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
       x3=(x1+x2)*0.5-d*(y2-y1)/s;
       y3=(y1+y2)*0.5+d*(x2-x1)/s;
       x4=x1;
       y4=y1;
       newpath(); plot(x4,y4,3);
       for (i=1;i<=n;i++) {
         t=0.7*(double)i/(double)n; 
         x4=x1*(t-1.0)*(t-2.0)*0.5-x3*t*(t-2.0)+x2*t*(t-1.0)*0.5;
         y4=y1*(t-1.0)*(t-2.0)*0.5-y3*t*(t-2.0)+y2*t*(t-1.0)*0.5;
         plot(x4,y4,2);
       };     stroke(); 
       x4=x2;
       y4=y2;
       newpath(); plot(x4,y4,3);
       for (i=1;i<=n;i++) {
         t=0.7*(double)i/(double)n; 
         x4=x2*(t-1.0)*(t-2.0)*0.5-x3*t*(t-2.0)+x1*t*(t-1.0)*0.5;
         y4=y2*(t-1.0)*(t-2.0)*0.5-y3*t*(t-2.0)+y1*t*(t-1.0)*0.5;
         plot(x4,y4,2);
       };     stroke(); 
      return;
}

void dimsrot(double x1, double y1, double x2, double y2,double d,double t1)
{     double x3,y3,x4,y4,s,t; 
      int i,n;
       n=10;
       s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
       x3=(x1+x2)*0.5-d*(y2-y1)/s;
       y3=(y1+y2)*0.5+d*(x2-x1)/s;
       x4=x1;
       y4=y1;
       newpath(); plotrot(x4,y4,t1,3);
       for (i=1;i<=n;i++) {
         t=0.7*(double)i/(double)n; 
         x4=x1*(t-1.0)*(t-2.0)*0.5-x3*t*(t-2.0)+x2*t*(t-1.0)*0.5;
         y4=y1*(t-1.0)*(t-2.0)*0.5-y3*t*(t-2.0)+y2*t*(t-1.0)*0.5;
         plotrot(x4,y4,t1,2);
       };     stroke(); 
       x4=x2;
       y4=y2;
       newpath(); plotrot(x4,y4,t1,3);
       for (i=1;i<=n;i++) {
         t=0.7*(double)i/(double)n; 
         x4=x2*(t-1.0)*(t-2.0)*0.5-x3*t*(t-2.0)+x1*t*(t-1.0)*0.5;
         y4=y2*(t-1.0)*(t-2.0)*0.5-y3*t*(t-2.0)+y1*t*(t-1.0)*0.5;
         plotrot(x4,y4,t1,2);
       };     stroke(); 
      return;
}

void ellipse(double x1, double y1, double rx, double ry)
{     
      double x, y, dt; int i;
      newpath();
      x=x1+rx; y=y1; dt=3.14156*2.0/40.0;
      plot(x,y,3); 
      for (i=1;i<=40;i++) { 
             x=x1+rx*cos(dt*i);y=y1+ry*sin(dt*i); plot(x,y,2);
             }
      return; 
}

void ellipserot(double x1, double y1, double rx, double ry, double t)
{     
      double x, y, dt; int i;
      newpath();
      x=x1+rx; y=y1; dt=3.14156*2.0/40.0;
      plotrot(x,y,t,3); 
      for (i=1;i<=40;i++) { 
             x=x1+rx*cos(dt*i);y=y1+ry*sin(dt*i); plotrot(x,y,t,2);
             }
      return; 
}

void eoclipon()
{ 
  fprintf(stream, "eoclip\n");
}

void eocliponrec(double x1, double y1, double x2, double y2)
{ 
  newpath();
  plot(x1, y1, 3);plot(x2, y1, 2);plot(x2, y2,2); plot(x1, y2, 2);
  closepath();
  fprintf(stream, "eoclip\n");
}

void fill(void)
{
  fprintf(stream, "fill\n"); return;
}

void fin(void)
{
  stroke(); fprintf(stream, "showpage\n"); fclose(stream);
  return;
}

void grestore(void)
{
  stroke(); fprintf(stream, "grestore\n"); fclose(stream);
  return;
}

void gsave(void)
{
  stroke(); fprintf(stream, "gsave\n"); fclose(stream);
  return;
}

void init(void)
{ stream = fopen("temp1.ps", "w");

  // aspect ratio of the paper is 1.4143
    fprintf(stream, "%%! Created by Jiro Mizushima \n");
    fprintf(stream, "/Times-Roman findfont 18 scalefont setfont \n");
    fprintf(stream, "/cm {28.35 mul} def \n");

  xsize = 21.0;/* ysize = 29.7;*/ ysize = 21.0;
  viewport(0.2, 0.2, 0.8, 0.8);
  xyworld(0.0, 0.0, 1.0, 1.0);
  linety(1); linewidth(1);
  return;
}

void init1(int n)
{   
    char str[12];
    sprintf(str,"temp%3d.ps",n);
    stream = fopen(str, "w");

  // aspect ratio of the paper is 1.4143
    fprintf(stream, "%%! Created by Jiro Mizushima \n");
    fprintf(stream, "/Times-Roman findfont 18 scalefont setfont \n");
    fprintf(stream, "/cm {28.35 mul} def \n");

  xsize = 21.0;/* ysize = 29.7;*/ ysize = 21.0;
  viewport(0.2, 0.2, 0.8, 0.8);
  xyworld(0.0, 0.0, 1.0, 1.0);
  linety(1); linewidth(1);
  return;
}

void line(double x, double y, double xx, double yy)
{
  
  plot(x, y, 3); 
  plot(xx,yy,2);
  return;
}


void lgcurves(double * x1, double *y1, int n)
{ 
  int ns=10,i,j; double u,x,y; 

  for (j=0; j<n-1; j++) {
    u=(double) j-1; 
    x=-u*(u-1.)*(u-2.)/6.*x1[0]+(u+1.)*(u-1.)*(u-2.)/2.*x1[1]-(u+1.)*u*(u-2.)/2.*x1[2] 
       +(u+1.)*u*(u-1.)/6.*x1[3]; 
    y=-u*(u-1.)*(u-2.)/6.*y1[0]+(u+1.)*(u-1.)*(u-2.)/2.*y1[1]-(u+1.)*u*(u-2.)/2.*y1[2] 
       +(u+1.)*u*(u-1.)/6.*y1[3]; 
   plot(x,y,3); 
   for (i=1;i<=ns;i++){
    u=(double) j-1 +(double) i/(double) ns;
    x=-u*(u-1.)*(u-2.)/6.*x1[0]+(u+1.)*(u-1.)*(u-2.)/2.*x1[1]-(u+1.)*u*(u-2.)/2.*x1[2] 
       +(u+1.)*u*(u-1.)/6.*x1[3]; 
    y=-u*(u-1.)*(u-2.)/6.*y1[0]+(u+1.)*(u-1.)*(u-2.)/2.*y1[1]-(u+1.)*u*(u-2.)/2.*y1[2] 
      +(u+1.)*u*(u-1.)/6.*y1[3]; 
    plot(x,y,2); 
    }
   }
   return;
}

void linerot(double x, double y, double xx, double yy, double t)
{
  plotrot(x, y, t, 3); 
  plotrot(xx,yy,t, 2);
  return;
}

void linecap(int ip)
{
  fprintf(stream, "%d setlinecap\n", ip);
  return;
}

void linety(int ichar)
{
  switch (ichar) {
    case 1: fprintf(stream, "[] 0 setdash \n"); break;
    case 2: fprintf(stream, "[2 2] 0 setdash \n");break;
    case 3: fprintf(stream, "[4 2] 0 setdash \n");break;
    case 4: fprintf(stream, "[8 2] 0 setdash \n");break;
    case 5: fprintf(stream, "[16 4] 0 setdash \n");break;
    case 6: fprintf(stream, "[32 8] 0 setdash \n");break;
    case 7: fprintf(stream, "[4 1 1 1] 0 setdash \n");break;
    case 8: fprintf(stream, "[8 1 2 1] 0 setdash \n");break;
    case 9: fprintf(stream, "[16 1 4 1] 0 setdash \n");break;
  }
  return;
}

void linewidth(double w)
{
  fprintf(stream, "%lf setlinewidth\n", w); return;
}

void newpath(void)
{
  fprintf(stream, "newpath\n"); return;
}

void parabola(double x1, double y1, double x2, double y2, double x3, double y3)
{
      double u,x,y; int i,ns=20; u=0.0;
      newpath(); 
      x=(u-1.)*(u-2.)/2.*x1-u*(u-2.)*x2+u*(u-1.)/2.*x3; 
      y=(u-1.)*(u-2.)/2.*y1-u*(u-2.)*y2+u*(u-1.)/2.*y3; 
      plot(x,y,3); 
      for (i=1; i<=2*ns; i++)
        { u=(double) i/(double) ns; 
          x=(u-1.)*(u-2.)/2.*x1-u*(u-2.)*x2+u*(u-1.)/2.*x3; 
          y=(u-1.)*(u-2.)/2.*y1-u*(u-2.)*y2+u*(u-1.)/2.*y3; 
          plot(x,y,2); 
         } stroke(); 
      return; 
}

void plot(double x, double y, int ipen)
{
  double x1, y1;

  x1 = ((x-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
  y1 = ((y-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
  switch (ipen) {
    case 3: fprintf(stream, "%lf cm %lf cm moveto\n", x1,y1);
      break;
    case 2: fprintf(stream, "%lf cm %lf cm lineto\n", x1,y1);
      break;
     }
  return;
}

void plotrot(double x, double y, double t, int ipen)
{
  double x1, y1, t1;

  t1 = t*3.14159/180.0;
  x1 = x*cos(t1)-y*sin(t1);
  y1 = x*sin(t1)+y*cos(t1);
  x1 = ((x1-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
  y1 = ((y1-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
  switch (ipen) {
    case 3: fprintf(stream, "%lf cm %lf cm moveto\n", x1,y1);
      break;
    case 2: fprintf(stream, "%lf cm %lf cm lineto\n", x1,y1);
      break;
     }
  return;
}


void polygon(double x1, double y1, double r, int n)
{     
      double x, y, dt; int i;
      newpath();
      x=x1; y=y1+r; dt=3.14156*2.0/(double)n;
      plot(x,y,3); 
      for (i=1;i<=n;i++) { 
             x=x1+r*sin(dt*i);y=y1+r*cos(dt*i); plot(x,y,2);
             }
      return; 
}

void polygonrot(double x1, double y1, double r, double t1, int n)
{     
      double x, y, t, dt; int i;
      newpath();
      t=3.14156*t1/180.0;
      x=x1+r*sin(t1); y=y1+r*cos(t1); 
      dt=3.14156*2.0/(double)n; plot(x,y,3); 
      for (i=1;i<=n;i++) { 
             x=x1+r*sin(dt*i+t);y=y1+r*cos(dt*i+t); plot(x,y,2);
             }
      return; 
}

void rect(double x1, double y1, double x2, double y2)
{
  newpath();plot(x1,y1,3);plot(x2,y1,2);plot(x2,y2,2);plot(x1,y2,2);
  closepath();
  return;
}

void rectrot(double x1, double y1, double x2, double y2,double t)
{
  newpath();plotrot(x1,y1,t,3);plotrot(x2,y1,t,2);plotrot(x2,y2,t,2);
  plotrot(x1,y2,t,2);
  closepath();
  return;
}


void rectround(double x1, double y1, double x2, double y2,double r1)
{
  newpath();plot((x1+x2)/2.0,y1,3);arcto(x2,y1,x2,y2,r1);
  arcto(x2,y2,x1,y2,r1);arcto(x1,y2,x1,y1,r1);arcto(x1,y1,x2,y1,r1);
  closepath();
  return;
}

void rectroundrot(double x1, double y1, double x2, double y2, double r1, double t1)
{
  newpath();plotrot((x1+x2)/2.0,y1,t1,3);arctorot(x2,y1,x2,y2,r1,t1);
  arctorot(x2,y2,x1,y2,r1,t1);arctorot(x1,y2,x1,y1,r1,t1);arctorot(x1,y1,x2,y1,r1,t1);
  closepath();
  return;
}

void resist(double x1, double y1, double x2, double y2, double d)
{
   double s, dx, dy, ex, ey, xx, yy;int i; 
      s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); 
      dx=(x2-x1)/12.0; dy=(y2-y1)/12.0; 
      ex=-d*dy; ey=d*dx; 
       plot(x1, y1, 3); xx=x1+dx+ex; yy=y1+dy+ey; plot(xx, yy, 2); 
      for(i=1;i<=5; i++){
        xx=x1+(i*2+1)*dx+pow(-1,i)*ex; yy=y1+(i*2+1)*dy+pow(-1,i)*ey; 
        plot(xx, yy, 2);
        }
      plot(x2, y2, 2); 
     return;
}

void rotate(double theta)
{ 
  translateo(theta);
  fprintf(stream, "%lf rotate\n",theta); 
}

void rplot(double dx, double dy, int ipen)
{
  double dx1, dy1;

  dx1 = ((dx-xw1)/(xw2-xw1)*(xv2-xv1) )*xsize;
  dy1 = ((dy-yw1)/(yw2-yw1)*(yv2-yv1) )*ysize;
  switch (ipen) {
    case 3: fprintf(stream, "%lf cm %lf cm rmoveto\n", dx1,dy1);
      break;
    case 2: fprintf(stream, "%lf cm %lf cm rlineto\n", dx1,dy1);
      break;
     }
  return;
}

void rplotrot(double dx, double dy, double t, int ipen)
{
  double dx1, dy1, t1;

  t1 = t*3.14159/180.0;
  dx1 = dx*cos(t1)-dy*sin(t1);
  dy1 = dx*sin(t1)+dy*cos(t1);
  dx1 = ((dx1-xw1)/(xw2-xw1)*(xv2-xv1) )*xsize;
  dy1 = ((dy1-yw1)/(yw2-yw1)*(yv2-yv1) )*ysize;
  switch (ipen) {
    case 3: fprintf(stream, "%lf cm %lf cm rmoveto\n", dx1,dy1);
      break;
    case 2: fprintf(stream, "%lf cm %lf cm rlineto\n", dx1,dy1);
      break;
     }
  return;
}

void scale(double tx, double ty)
{ 
  fprintf(stream, "%lf %lf scale\n",tx, ty); 
}

void setchar(int ichar, int ip)
{
 switch (ichar) {
  case 1:fprintf(stream, "/Times-Roman findfont %d scalefont setfont \n",ip);
      break;
  case 2:fprintf(stream, "/Times-Bold findfont %d scalefont setfont \n",ip);
      break;
  case 3:fprintf(stream, "/Times-Italic findfont %d scalefont setfont \n",ip);
      break;
  case 4:fprintf(stream, "/Times-BoldItalic findfont %d scalefont setfont \n",ip);
      break;
  case 5:fprintf(stream, "/Helvetica findfont %d scalefont setfont \n",ip);
      break;
  case 6:fprintf(stream, "/Helvetica-Bold findfont %d scalefont setfont \n",ip);
      break;
  case 7:fprintf(stream, "/Helvetica-Oblique findfont %d scalefont setfont \n",ip);
      break;
  case 8:fprintf(stream, "/Helvetica-BoldOblique findfont %d scalefont setfont \n",ip);
      break;
  case 9:fprintf(stream, "/Courier findfont %d scalefont setfont \n",ip);
      break;
  case 10:fprintf(stream, "/Courier-Bold findfont %d scalefont setfont \n",ip);
      break;
  case 11:fprintf(stream, "/Courier-Oblique findfont %d scalefont setfont \n",ip);
      break;
  case 12:fprintf(stream, "/Courier-BoldOblique findfont %d scalefont setfont \n",ip);
      break;
  case 13:fprintf(stream, "/Symbol findfont %d scalefont setfont \n",ip);
      break;
   }
  return;
}

void setgray(double g)
{
  fprintf(stream, "%lf setgray\n", g); return;
}

void setlinejoin(int ljoin)
{
  fprintf(stream, "%d setlinejoin\n",ljoin); return;
}

void setrgb(double r,double g,double b)
{
  fprintf(stream, "%lf %lf %lf setrgbcolor\n",r,g,b); return;
}

void setcmyk(double c,double m,double y,double k)
{
  fprintf(stream, "%lf %lf %lf %lf setcmykcolor\n",c,m,y,k); return;
}


void spline(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, int ipart)
{ int ns=10,i; double u,x,y; 
   
  newpath(); 
  u=(double) ipart; 
  x=-u*(u-1.)*(u-2.)/6.*x1+(u+1.)*(u-1.)*(u-2.)/2.*x2-(u+1.)*u*(u-2.)/2.*x3+(u+1.)*u*(u-1.)/6.*x4; 
  y=-u*(u-1.)*(u-2.)/6.*y1+(u+1.)*(u-1.)*(u-2.)/2.*y2-(u+1.)*u*(u-2.)/2.*y3+(u+1.)*u*(u-1.)/6.*y4; 
  plot(x,y,3); 
   for (i=1;i<=ns;i++){
    u=(double) ipart +(double) i/(double) ns;
    x=-u*(u-1.)*(u-2.)/6.*x1+(u+1.)*(u-1.)*(u-2.)/2.*x2-(u+1.)*u*(u-2.)/2.*x3+(u+1.)*u*(u-1.)/6.*x4; 
    y=-u*(u-1.)*(u-2.)/6.*y1+(u+1.)*(u-1.)*(u-2.)/2.*y2-(u+1.)*u*(u-2.)/2.*y3+(u+1.)*u*(u-1.)/6.*y4; 
    plot(x,y,2); 
    }
   stroke();
   return;
}

void square(double x, double y, double ell)
{
  plot(x-ell/2.0,y-ell/2.0,3);
  plot(x+ell/2.0,y-ell/2.0,2);
  plot(x+ell/2.0,y+ell/2.0,2);
  plot(x-ell/2.0,y+ell/2.0,2);
  closepath();
  return;
}

void squarerot(double x, double y, double ell, double t)
{
  plotrot(x-ell/2.0,y-ell/2.0,t,3);
  plotrot(x+ell/2.0,y-ell/2.0,t,2);
  plotrot(x+ell/2.0,y+ell/2.0,t,2);
  plotrot(x-ell/2.0,y+ell/2.0,t,2);
  closepath();
  return;
}

void star(double x0, double y0, double r, double t, int n)
{
  double pi,dth,dthi,alpha;int i;
  if(n<3){
	printf("can't make a polygon");
	return;
  }
  if (n==6){
     alpha=0.575;
    } else if(n<=4) {
	 alpha=0.24;
    } else {
     alpha=0.36;
  }
  pi=3.14159; dth=pi/(double)n;
  newpath();plot(x0+r*cos(t+pi/2.0),y0+r*sin(t+pi/2.0),3);
  for(i=1; i<=2*n-1 ; i++){
	if(i%2==0){
     dthi=dth*(double)i;
     plot(x0+r*cos(t+pi/2.0+dthi),y0+r*sin(t+pi/2.0+dthi),2);
	} else {
	 dthi=dth*(double)i;
	 plot(x0+alpha*r*cos(t+pi/2.0+dthi),y0+alpha*r*sin(t+pi/2.0+dthi),2);
    }
  }
  closepath();
  return;
}

void starx(double x0, double y0, double r, double t, int n)
{
  double pi,dth,dthi,alpha;int i;
  if(n<3){
	printf("can't make a polygon");
	return;
  }
  pi=3.14159;
  dth=pi/(double)n*2.0;
  newpath();plot(x0+r*cos(t+pi/2.0),y0+r*sin(t+pi/2.0),3);
  if(n%2==1){
	for(i=2 ; i<=2*n ; i=i+2){
	  dthi=dth*(double)i; 
	  plot(x0+r*cos(t+pi/2.0+dthi),y0+r*sin(t+pi/2.0+dthi),2);
	}
	closepath();
  } else {
	for(i=2 ; i<=n ; i=i+2){
	  dthi=dth*(double)i;
	  plot(x0+r*cos(t+pi/2.0+dthi),y0+r*sin(t+pi/2.0+dthi),2);
	}
	closepath();
	i=1;dthi=dth*(double)i;
    plot(x0+r*cos(t+pi/2.0+dthi),y0+r*sin(t+pi/2.0+dthi),3);
	for(i=3 ; i<=n-1 ; i=i+2){
	  dthi=dth*(double)i;
	  plot(x0+r*cos(t+pi/2.0+dthi),y0+r*sin(t+pi/2.0+dthi),2);
	}
	closepath();
  }
  return;
}


void stroke(void)
{ 
  fprintf(stream, "stroke\n"); 
}

void text(double x, double y, char* str1)
{
  double x1, y1;

  x1 = ((x-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
  y1 = ((y-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
    fprintf(stream, "%lf cm %lf cm moveto\n", x1, y1);
    fprintf(stream, "(%s) show\n", str1);
  return;
}

void text1(double x, double y, char s1)
{
  double x1, y1;

  x1 = ((x-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
  y1 = ((y-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
    fprintf(stream, "%lf cm %lf cm moveto\n", x1, y1);
    fprintf(stream, "(%c) show\n", s1);
  return;
}

void textx(double x, double y, char* str1)
{
	  int n; n=strlen(str1);
	  printf("%d length\n",n);
      plot(x, y, 3); 
      fprintf(stream, "%lf %lf rmoveto\n", -n/2.0*9.0,-3.0/2.0*15.0);
      fprintf(stream, "(%s) show\n", str1);
      return;
}

void textrot(double x, double y, double t, char* str1)
{
  rotate(t);
  text(x, y,str1);
  rotate(-t);
  return;
}


void  texty(double x, double y, char* str1)
{ 
	  int n; n=strlen(str1);
      plot(x, y, 3); 
      fprintf(stream, "%lf %lf rmoveto\n", -(n+1)*9.0,-1.0/2.0*12.0);
      fprintf(stream, "(%s) show\n", str1);
      return;
}

void translate(double x, double y)
{ double x1,y1;
  x1 = ((x-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
  y1 = ((y-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
  fprintf(stream, "%lf cm %lf cm translate\n", x1, y1);
}

void translateo(double theta)
{ double xx,yy,x1,y1,t;
  xx = ((0.0-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
  yy = ((0.0-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
  t=theta*3.14159/180.0;
  x1 = xx*cos(t)-yy*sin(t);
  y1 = xx*sin(t)+yy*cos(t);
  fprintf(stream, "%lf cm %lf cm translate\n", xx-x1, yy-y1);
}

void triangle(double x1, double y1, double d1, double a1)
{
    double r1,aa1,aa2,aa3;

     r1=d1/sqrt(3.0);
     aa1=(a1-30.0)*3.14159/180.0;
     aa2=aa1+2.0*3.14159/3.0;
     aa3=aa2+2.0*3.14159/3.0;
     newpath(); plot(x1+r1*cos(aa1),y1+r1*sin(aa1),3); 
     plot(x1+r1*cos(aa2),y1+r1*sin(aa2),2); plot(x1+r1*cos(aa3),y1+r1*sin(aa3),2);
     closepath();
     return; 
}

void viewport(double xv1d, double yv1d, double xv2d, double yv2d)
{
  xv1 = xv1d; yv1 = yv1d; xv2 = xv2d; yv2 = yv2d;
  return;
}

void xyworld(double xw1d, double yw1d, double xw2d, double yw2d)
{
  xw1 = xw1d; yw1 = yw1d; xw2 = xw2d; yw2 = yw2d;
  return;
}