!plotting routines for general purpose using postscript
   subroutine annulus(x1, y1, r1, r2)
      real(4) x1, y1, r1, r2
      call plot(x1+r2, y1, 3)
      call arc(x1, y1, r2, 1.0, 360.0)
      call plot(x1+r1, y1, 2)
      call arcn(x1, y1, r1, 360.0, 0.0)
      call plot(x1+r2, y1, 2)
      call closepath
      return
      end

   subroutine arc(x1, y1, r1, t1, t2)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, r1, t1, t2
      real(4) 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
      write(1,100) xx1,yy1,rr1,t1,t2
      100 format(1x,3(1pe12.5,' cm '),1pe12.5,1pe12.5,' arc ')
      return
      end

   subroutine arcn(x1, y1, r1, t1, t2)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, r1, t1, t2
      real(4) 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
      write(1,100) xx1,yy1,rr1,t1,t2
      100 format(1x,3(1pe12.5,' cm '),1pe12.5,1pe12.5,' arcn ')
      return
      end

   subroutine arcto(x1, y1, x2, y2, r1)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, x2, y2, r1
      real(4) 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
      write(1,100) xx1,yy1,xx2,yy2,rr1
      100 format(1x,5(1pe12.5,' cm '),' arcto ')
      return
      end

   subroutine arctorot(x1, y1, x2, y2, r1, t1)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, x2, y2, r1, t1
      real(4) 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
      write(1,100) xx1,yy1,xx2,yy2,rr1
      100 format(1x,5(1pe12.5,' cm '),' arcto ')
      return
      end

   subroutine arrow(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4),parameter::alpha=2.2
      real(4) 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
      call plot(x1, y1, 3)
      call plot(x2, y2, 2)
      call plot(x3, y3, 3)
      call plot(x2, y2, 2)
      call plot(x4, y4, 2)
      return
      end

   subroutine arrowa(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4),parameter::alpha=2.2
      real(4) 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
      call plot(x1, y1, 3)
      call plot((x3+x4)/2.0, (y3+y4)/2.0, 2)
      call stroke
      call plot(x3, y3, 3)
      call plot(x2, y2, 2)
      call plot(x4, y4, 2)
      call closepath
      call fill
      return
      end

   subroutine arrowb(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4),parameter::alpha=2.2
      real(4) 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
      call plot(x1, y1, 3)
      call plot(x43, y43, 2)
      call stroke
      call plot(x3, y3, 3)
      call plot(x32, y32, 2)
      call plot(x2, y2, 2)
      call plot(x24, y24, 2)
      call plot(x4, y4, 2)
      call plot(x43,y43,2)
      call closepath
      call fill
      return
      end

   subroutine arrowc(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4),parameter::alpha=2.2
      real(4) 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
      call plot(x1, y1, 3)
      call plot(x5, y5, 2)
      call stroke
      call plot(x3, y3, 3)
      call plot(x2, y2, 2)
      call plot(x4, y4, 2)
      call plot(x5, y5, 2)
      call closepath
      call fill
      return
      end

   subroutine arrowfill(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4),parameter::alpha=2.2
      real(4) 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
      call plot(x1, y1, 3)
      call plot(x2, y2, 2)
      call stroke
      call plot(x3, y3, 3)
      call plot(x2, y2, 2)
      call plot(x4, y4, 2)
      call closepath
      call fill
      return
      end

   subroutine arrowrot(x1, y1, x2, y2, t1, d)
      real(4) x1, y1, x2, y2, t1, d
      real(4),parameter::alpha=2.2
      real(4) 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
      call plotrot(x1, y1, t1, 3)
      call plotrot(x2, y2, t1, 2)
      call plotrot(x3, y3, t1, 3)
      call plotrot(x2, y2, t1, 2)
      call plotrot(x4, y4, t1, 2)
      return
      end

   subroutine arrowwide(x1, y1, x2, y2, d, al)
      real(4) x1, y1, x2, y2, d, al
      real(4) 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
      call plot(x6, y6, 3)
      call plot(x8, y8, 2)
      call plot(x4, y4, 2)
      call plot(x2, y2, 2)
      call plot(x3, y3, 2)
      call plot(x7, y7, 2)
      call plot(x5, y5, 2)
      call closepath()
      call stroke()
      return
      end

   subroutine battery(x1, y1, x2, y2, d1, d2)
      real(4) x1, y1, x2, y2, d1, d2
      real(4) dx, dy, ex, ey, xx, yy
      dx=x2-x1
      dy=y2-y1
      ex=-dy
      ey=dx
      xx=x1+d1*ex
      yy=y1+d1*ey
      call plot(xx, yy, 3)
      xx=x1-d1*ex
      yy=y1-d1*ey
      call plot(xx, yy, 2)
      xx=x2+d2*ex
      yy=y2+d2*ey
      call plot(xx, yy, 3)
      xx=x2-d2*ex
      yy=y2-d2*ey
      call plot(xx, yy, 2)
      return
      end

   subroutine brokenlines(x1, y1, n)
      real(4),dimension(0:n)::x1, y1
      integer n
      integer i
      call newpath
      call plot(x1(0),y1(0),3)
      do i=1,n
         call plot(x1(i),y1(i),2)
      end do
      return
      end

   subroutine circ(x1, y1, r1)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, r1
      real(4) 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
      call newpath
      write(1,100) xx1,yy1,rr1
      100 format(1x,3(1pe12.5,' cm '), ' 0 360 arc' )
      return
      end

   subroutine circn(x1, y1, r1)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, r1
      real(4) 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
      call newpath
      write(1,100) xx1,yy1,rr1
      100 format(1x,3(1pe12.5,' cm '), ' 0 360 arcn' )
      return
      end

   subroutine clipoff
      write(1,*) ' initclip'
      return
      end

   subroutine clipon
      write(1,*) ' clip'
      return
      end

   subroutine cliponrec(x1, y1, x2, y2)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, x2, y2
      call newpath
      call plot(x1, y1, 3)
      call plot(x2, y1, 2)
      call plot(x2, y2, 2)
      call plot(x1, y2, 2)
      call closepath
      write(1,*) ' clip'
      return
      end

   subroutine closepath
      write(1,*) ' closepath'
      return
      end

   subroutine coil(x1, y1, x2, y2, d, n)
      real(4) x1, y1, x2, y2, d
      integer n
      real(4) s,ex,ey,ds,a,xx,yy,t
      integer 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)/(real(n)+0.5)
      call plot(x1, y1, 3)
      do i=1,40*n+20
         t=2.0*3.14159*real(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
         call plot(xx, yy, 2)
      end do
      return
      end

   subroutine coilnew(x1, y1, x2, y2, d, n)
      real(4) x1, y1, x2, y2, d
      integer n
      real(4) s,ex,ey,xx,yy,t
      integer i,j
      s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
      ex=(x2-x1)/s
      ey=(y2-y1)/s
      call plot(x1, y1, 3)
      do i=1,n
      do j=1,20
         t=3.14159*real(j)/20.0
         xx=x1+(real(i)*2.0-1.0-cos(t))*d*ex-d*sin(t)*ey
         yy=y1+(real(i)*2.0-1.0-cos(t))*d*ey+d*sin(t)*ex
         call plot(xx, yy, 2)
      end do
      end do
      return
      end

   subroutine curvto(x1, y1, x2, y2, x3, y3)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, x2, y2, x3, y3
      real(4) 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
      write(1,100) xx1,yy1,xx2,yy2,xx3,yy3
      100 format(1x,6(1pe12.5,' cm '),' curveto ')
      return
      end

   subroutine curvton(x1, y1, n)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4),dimension(0:n)::x1, y1
      integer n
      real(4) xx1, yy1
      integer i
      do i=0,n
         xx1=((x1(i)-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize
         yy1=((y1(i)-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize
         write(1,100) xx1,yy1
         100 format(1x,2(1pe12.5,' cm '))
      end do
      write(1,*) ' curveto'
      return
      end

   subroutine curvtona(x1, y1, t, n)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4),dimension(0:n)::x1, y1
      real(4) t
      integer n
      real(4),dimension(0:n)::x12, y12, x21, y21
      integer i
      do i=0,n-1
         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)
      end do
      call plot(x21(0),y21(0), 2)
      do i=0,n-2
         call curvto(x21(i),y21(i),x1(i+1),y1(i+1),x12(i+1),y12(i+1))
      end do
      call plot(x1(n),y1(n),2)
      return
      end

   subroutine dims(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4) x3, y3, x4, y4, s, t
      integer 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
      call newpath
      call plot(x4,y4,3)
      do i=1,n
         t=0.7*real(i)/real(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
         call plot(x4,y4,2)
      end do
      call stroke
      x4=x2
      y4=y2
      call newpath
      call plot(x4,y4,3)
      do i=1,n
         t=0.7*real(i)/real(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
         call plot(x4,y4,2)
      end do
      call stroke
      return
      end

   subroutine dimsrot(x1, y1, x2, y2, d, t1)
      real(4) x1, y1, x2, y2, d, t1
      real(4) x3, y3, x4, y4, s, t
      integer 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
      call newpath
      call plotrot(x4,y4,t1,3)
      do i=1,n
         t=0.7*real(i)/real(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
         call plotrot(x4,y4,t1,2)
      end do
      call stroke
      x4=x2
      y4=y2
      call newpath
      call plotrot(x4,y4,t1,3)
      do i=1,n
         t=0.7*real(i)/real(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
         call plotrot(x4,y4,t1,2)
      end do
      call stroke
      return
      end

   subroutine ellipse(x1, y1, rx, ry)
      real(4) x1, y1, rx, ry
      real(4) x, y, dt
      integer i
      call newpath
      x=x1+rx
      y=y1
      dt=3.14156*2.0/40.0
      call plot(x,y,3)
      do i=1,40
         x=x1+rx*cos(dt*i)
         y=y1+ry*sin(dt*i)
         call plot(x,y,2)
      end do
      return
      end

   subroutine ellipserot(x1, y1, rx, ry, t)
      real(4) x1, y1, rx, ry, t
      real(4) x, y, dt
      integer i
      call newpath
      x=x1+rx
      y=y1
      dt=3.14156*2.0/40.0
      call plotrot(x,y,t,3)
      do i=1,40
         x=x1+rx*cos(dt*i)
         y=y1+ry*sin(dt*i)
         call plotrot(x,y,t,2)
      end do
      return
      end

   subroutine eoclipon
      write(1,*) ' eoclip'
      return
      end

   subroutine eocliponrec(x1, y1, x2, y2)
      real(4) x1, y1, x2, y2
      call newpath
      call plot(x1, y1, 3)
      call plot(x2, y1, 2)
      call plot(x2, y2, 2)
      call plot(x1, y2, 2)
      call closepath
      write(1,*) ' eoclip'
      return
      end

   subroutine fill
      write(1,*) ' fill'
      return
      end

   subroutine fin
      call stroke
      write(1,*) ' showpage'
      close(1)
      return
      end

   subroutine grestore
      write(1,*) ' grestore'
      return
      end

   subroutine gsave
      write(1,*) ' gsave'
      return
      end

   subroutine init
      common /paper/ xsize, ysize
      open(1,file='temp1.ps')
!aspect ratio of the paper is 1.4143
      write(1,*) '%! Created by Jiro Mizushima'
      write(1,*) '/Times-Roman findfont 18 scalefont setfont'
      write(1,*) '/cm {28.35 mul } def '
      xsize=21.0
!      ysize=29.7
      ysize=21.0
      call viewport(0.2,0.2,0.8,0.8)
      call xyworld(0.0,0.0,1.0,1.0)
      call linety(1)
      call linewidth(1.0)
      return
      end

   subroutine line(x, y, xx, yy)
      real(4) x, y, xx, yy
      call plot(x, y, 3)
      call plot(xx,yy,2)
      return
      end

   subroutine linerot(x, y, xx, yy, t)
      real(4) x, y, xx, yy, t
      call plotrot(x, y, t, 3)
      call plotrot(xx,yy,t, 2)
      return
      end

   subroutine linecap(ip)
      integer ip
      if (ip==0) then 
         write(1,*) '0 setlinecap'
      else if (ip==1) then
         write(1,*) '1 setlinecap'
      else if (ip==2) then
         write(1,*) '2 setlinecap'
      end if
      return
      end

   subroutine linety(ichar)
      integer ichar
      if (ichar==1) then 
         write(1,*) '[] 0 setdash'
         else if (ichar==2) then
         write(1,*) '[2 2] 0 setdash'
         else if (ichar==3) then
         write(1,*) '[4 2] 0 setdash'
         else if (ichar==4) then
         write(1,*) '[8 2] 0 setdash'
         else if (ichar==5) then
         write(1,*) '[16 4] 0 setdash'
         else if (ichar==6) then
         write(1,*) '[32 8] 0 setdash'
         else if (ichar==7) then
         write(1,*) '[4 1 1 1] 0 setdash'
         else if (ichar==8) then
         write(1,*) '[8 1 2 1] 0 setdash'
         else if (ichar==9) then
         write(1,*) '[16 1 4 1] 0 setdash'
      end if
      return
      end

   subroutine linewidth(w)
      write(1,*) w,' setlinewidth'
      return
      end

   subroutine lgcurves(x1, y1, n)
      real(4),dimension(0:n)::x1, y1
      integer n
      integer,parameter::ns=10
      integer i,j
      real(4) u,x,y
      do j=0,n-1
         u=real(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)
         call plot(x,y,3)
         do i=1,ns
            u=real(j-1)+real(i)/real(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)
            call plot(x,y,2)
         end do
      end do
      return
      end

   subroutine newpath
      write(1,*) ' newpath'
      return
      end

   subroutine parabola(x1, y1, x2, y2, x3, y3)
      real(4) x1, y1, x2, y2, x3, y3
      real(4) u,x,y
      integer i
      integer,parameter::ns=20
      u=0.0
      call 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
      call plot(x,y,3)
      do i=1,2*ns
         u=real(i)/real(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
         call plot(x,y,2)
      end do
      call stroke
      return
      end

   subroutine plot(x1,y1,ipen)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1
      integer ipen
      real(4) x,y
      x=((x1-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize
      y=((y1-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize
      if(ipen==3) write(1,100) x,y
      if(ipen==2) write(1,200) x,y
      100 format(1x,1pe11.4,' cm ',1pe11.4,' cm moveto')
      200 format(1x,1pe11.4,' cm ',1pe11.4,' cm lineto')
      return
      end

   subroutine plotrot(x1,y1,t1,ipen)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x1, y1, t1
      integer ipen
      real(4) x,y,t
      t=t1*3.14159/180.0;
      x=x1*cos(t)-y1*sin(t);
      y=x1*sin(t)+y1*cos(t);
      x=((x-xw1)/(xw2-xw1)*(xv2-xv1) + xv1)*xsize;
      y=((y-yw1)/(yw2-yw1)*(yv2-yv1) + yv1)*ysize;
      if(ipen==3) write(1,100) x,y
      if(ipen==2) write(1,200) x,y
      100 format(1x,1pe11.4,' cm ',1pe11.4,' cm moveto')
      200 format(1x,1pe11.4,' cm ',1pe11.4,' cm lineto')
      return
      end

   subroutine polygon(x1, y1, r, n)
      real(4) x1, y1, r
      integer n
      real(4) x,y,dt
      integer i
      call newpath
      x=x1
      y=y1+r
      dt=3.14156*2.0/real(n)
      call plot(x,y,3)
      do i=1,n
         x=x1+r*sin(dt*i)
         y=y1+r*cos(dt*i)
         call plot(x,y,2)
      end do
      return
      end

   subroutine polygonrot(x1, y1, r, t1, n)
      real(4) x1, y1, r, t1
      integer n
      real(4) x,y,t,dt
      integer i
      call newpath
      t=3.14156*t1/180.0
      x=x1+r*sin(t1)
      y=y1+r*cos(t1)
      dt=3.14156*2.0/real(n)
      call plot(x,y,3)
      do i=1,n
         x=x1+r*sin(dt*i+t)
         y=y1+r*cos(dt*i+t)
         call plot(x,y,2)
      end do
      return
      end

   subroutine rect(x1, y1, x2, y2)
      real(4) x1, y1, x2, y2
      call newpath
      call plot(x1,y1,3)
      call plot(x2,y1,2)
      call plot(x2,y2,2)
      call plot(x1,y2,2)
      call closepath
      return
      end

   subroutine rectrot(x1, y1, x2, y2, t)
      real(4) x1, y1, x2, y2, t
      call newpath
      call plotrot(x1,y1,t,3)
      call plotrot(x2,y1,t,2)
      call plotrot(x2,y2,t,2)
      call plotrot(x1,y2,t,2)
      call closepath
      return
      end

   subroutine rectround(x1, y1, x2, y2, r1)
      real(4) x1, y1, x2, y2, r1
      call newpath
      call plot((x1+x2)/2.0,y1,3)
      call arcto(x2,y1,x2,y2,r1)
      call arcto(x2,y2,x1,y2,r1)
      call arcto(x1,y2,x1,y1,r1)
      call arcto(x1,y1,x2,y1,r1)
      call closepath
      return
      end

   subroutine rectroundrot(x1, y1, x2, y2, r1, t1)
      real(4) x1, y1, x2, y2, r1, t1
      call newpath
      call plotrot((x1+x2)/2.0,y1,t1,3)
      call arctorot(x2,y1,x2,y2,r1,t1)
      call arctorot(x2,y2,x1,y2,r1,t1)
      call arctorot(x1,y2,x1,y1,r1,t1)
      call arctorot(x1,y1,x2,y1,r1,t1)
      call closepath
      return
      end

   subroutine resist(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4) dx,dy,ex,ey,xx,yy        !,s
      integer 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
      call plot(x1, y1, 3)
      xx=x1+dx+ex
      yy=y1+dy+ey
      call plot(xx, yy, 2)
      do i=1,5
         xx=x1+(i*2+1)*dx+(-1)**i*ex
         yy=y1+(i*2+1)*dy+(-1)**i*ey
         call plot(xx, yy, 2)
      end do
      call plot(x2, y2, 2)
      return
      end

   subroutine rotate(theta)
      real(4) theta
      call translateo(theta)
      write(1,100) theta
      100 format(1x,1pe12.5,' rotate ')
      return
      end

   subroutine rplot(dx,dy,ipen)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) dx,dy
      integer ipen
      real(4) dx1,dy1
      dx1=((dx-xw1)/(xw2-xw1)*(xv2-xv1))*xsize
      dy1=((dy-yw1)/(yw2-yw1)*(yv2-yv1))*ysize
      if(ipen==3) write(1,100) dx1,dy1
      if(ipen==2) write(1,200) dx1,dy1
      100 format(1x,1pe11.4,' cm ',1pe11.4,' cm rmoveto')
      200 format(1x,1pe11.4,' cm ',1pe11.4,' cm rlineto')
      return
      end

   subroutine rplotrot(dx, dy, t, ipen)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) dx, dy, t
      integer ipen
      real(4) 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
      if(ipen==3) write(1,100) dx1,dy1
      if(ipen==2) write(1,200) dx1,dy1
      100 format(1x,1pe11.4,' cm ',1pe11.4,' cm rmoveto')
      200 format(1x,1pe11.4,' cm ',1pe11.4,' cm rlineto')
      return
      end

   subroutine scale(tx, ty)
      real(4) tx, ty
      write(1,100) tx,ty
      100 format(1x,2(1pe12.5),' scale ')
      return
      end

   subroutine setchar(ichar,ip)
      integer ichar, ip
      if (ichar==1) then 
         write(1,*) '/Times-Roman findfont ',ip,' scalefont setfont'
      else if (ichar==2) then
         write(1,*) '/Times-Bold findfont ',ip,' scalefont setfont'
      else if (ichar==3) then
         write(1,*) '/Times-Italic findfont ',ip,' scalefont setfont'
      else if (ichar==4) then
         write(1,*) '/Times-BoldItalic findfont ',ip,' scalefont setfont'
      else if (ichar==5) then
         write(1,*) '/Helvetica findfont ',ip,' scalefont setfont'
      else if (ichar==6) then
         write(1,*) '/Helvetica-Bold findfont ',ip,' scalefont setfont'
      else if (ichar==7) then
         write(1,*) '/Helvetica-Oblique findfont ',ip,' scalefont setfont'
      else if (ichar==8) then
         write(1,*) '/Helvetica-BoldOblique findfont ',ip,' scalefont setfont'
      else if (ichar==9) then
         write(1,*) '/Courier findfont ',ip,' scalefont setfont'
      else if (ichar==10) then
         write(1,*) '/Courier-Bold findfont ',ip,' scalefont setfont'
      else if (ichar==11) then
         write(1,*) '/Courier-Oblique findfont ',ip,' scalefont setfont'
      else if (ichar==12) then
         write(1,*) '/Courier-BoldOblique findfont ',ip,' scalefont setfont'
      else if (ichar==13) then
         write(1,*) '/Symbol findfont ',ip,' scalefont setfont'
      end if
      return
      end

   subroutine setgray(g)
      real(4) g
      write(1,*) g,' setgray'
      return
      end

   subroutine setlinejoin(ljoin)
      integer ljoin
      write(1,*) ljoin,' setlinejoin'
      return
      end

   subroutine setrgb(r, g, b)
      real(4) r, g, b
      write(1,*) r,g,b,' setrgbcolor'
      return
      end

   subroutine setcmyk(c, m, y, k)
      real(4) c, m, y, k
      write(1,*) c,m,y,k,' setcmykcolor'
      return
      end

   subroutine spline(x1,y1,x2,y2,x3,y3,x4,y4,ipart)
      real(4) x1,y1,x2,y2,x3,y3,x4,y4
      integer ipart
      integer,parameter::ns=10
      integer i
      real(4) u,x,y
      call newpath
      u=real(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
      call plot(x,y,3)
      do i=1,ns
         u=real(ipart)+real(i)/real(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
         call plot(x,y,2)
      end do
      call stroke
      return
      end

   subroutine square(x, y, ell)
      real(4) x, y, ell
      call plot(x-ell/2.0,y-ell/2.0,3)
      call plot(x+ell/2.0,y-ell/2.0,2)
      call plot(x+ell/2.0,y+ell/2.0,2)
      call plot(x-ell/2.0,y+ell/2.0,2)
      call closepath
      return
      end

   subroutine squarerot(x, y, ell, t)
      real(4) x, y, ell, t
      call plotrot(x-ell/2.0,y-ell/2.0,t,3)
      call plotrot(x+ell/2.0,y-ell/2.0,t,2)
      call plotrot(x+ell/2.0,y+ell/2.0,t,2)
      call plotrot(x-ell/2.0,y+ell/2.0,t,2)
      call closepath
      return
      end

   subroutine star(x0, y0, r, t, n)
      real(4) x0, y0, r, t
      integer n
      real(4) pi,dth,dthi,alpha
      integer i
      if(n<3) then
         write(*,*) 'can''t make a polygon'
      else if(n==6)then
         alpha=0.575
      else if(n<=4)then
	 alpha=0.24
      else 
         alpha=0.36
      end if
      pi=3.14159
      dth=pi/real(n)
      call newpath
      call plot(x0+r*cos(t+pi/2.0),y0+r*sin(t+pi/2.0),3)
      do i=1,2*n-1
         if(mod(i,2)==0)then
            dthi=dth*real(i)
            call plot(x0+r*cos(t+pi/2.0+dthi),&
                     &y0+r*sin(t+pi/2.0+dthi),2)
         else 
            dthi=dth*real(i)
            call plot(x0+alpha*r*cos(t+pi/2.0+dthi),&
                     &y0+alpha*r*sin(t+pi/2.0+dthi),2)
         end if
      end do
      call closepath
      return
      end

   subroutine starx(x0, y0, r, t, n)
      real(4) x0, y0, r, t
      integer n
      real(4) pi,dth,dthi
      integer i
      if(n<3) write(*,*) 'can''t make a polygon'
      pi=3.14159
      dth=pi/real(n)*2.0
      call newpath
      call plot(x0+r*cos(t+pi/2.0),y0+r*sin(t+pi/2.0),3)
      if(mod(n,2)==1)then
	 do i=2,2*n,2
            dthi=dth*real(i)
            call plot(x0+r*cos(t+pi/2.0+dthi),&
                     &y0+r*sin(t+pi/2.0+dthi),2)
         end do
         call closepath
      else
	 do i=2,n,2
            dthi=dth*real(i)
            call plot(x0+r*cos(t+pi/2.0+dthi),&
                     &y0+r*sin(t+pi/2.0+dthi),2)
         end do
         call closepath
         i=1
         dthi=dth*real(i)
         call plot(x0+r*cos(t+pi/2.0+dthi),&
                  &y0+r*sin(t+pi/2.0+dthi),3)
         do i=3,n-1,2
            dthi=dth*real(i)
            call plot(x0+r*cos(t+pi/2.0+dthi),&
                     &y0+r*sin(t+pi/2.0+dthi),2)
         end do
         call closepath
      end if
      return
      end

   subroutine stroke
      write(1,*) ' stroke'
      return
      end

   subroutine text(x, y, str1)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x, y
      character str1*(*)
      real(4) x1,y1
      x1=((x-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize
      y1=((y-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize
      write(1,100) x1,y1
      100 format(1x,1pe12.5,' cm ',1pe12.5,' cm moveto')
      write(1,*) '(',str1,') show'
      return
      end

   subroutine text1(x, y, s1)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x, y
      character s1
      real(4) x1,y1
      x1=((x-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize
      y1=((y-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize
      write(1,100) x1,y1
      100 format(1x,1pe12.5,' cm ',1pe12.5,' cm moveto')
      write(1,*) '(',s1,') show'
      return
      end

   subroutine textx(x, y, str1)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x, y
      character str1*(*)
      integer n
      n=len_trim(str1)
      write(*,*) 'length',n
      call plot(x, y, 3)
      write(1,100) -n/2.0*9.0,-3.0/2.0*15.0
      100 format(1x,1pe12.5,1x,1pe12.5,' rmoveto')
      write(1,*) '(',str1,') show'
      return
      end

   subroutine textrot(x, y, t, str1)
      real(4) x, y, t
      character str1*(*)
      call rotate(t)
      call text(x,y,str1)
      call rotate(-t)
      return
      end

   subroutine texty(x, y, str1)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x, y
      character str1*(*)
      integer n
      n=len_trim(str1)
      write(*,*) 'length',n
      call plot(x, y, 3)
      write(1,100) -(n+1)*9.0,-1.0/2.0*12.0
      100 format(1x,1pe12.5,1x,1pe12.5,' rmoveto')
      write(1,*) '(',str1,') show'
      return
      end

   subroutine translateto(x, y)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x, y
      real(4) x1, y1
      x1=((x-xw1)/(xw2-xw1)*(xv2-xv1)+xv1)*xsize
      y1=((y-yw1)/(yw2-yw1)*(yv2-yv1)+yv1)*ysize
      write(1,100) x1,y1
      100 format(1x,1pe12.5,' cm ',1pe12.5,' cm translate')
      return
      end

   subroutine translaterel(x, y)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) x, y
      real(4) x1,y1
      x1=(x/(xw2-xw1)*(xv2-xv1))*xsize
      y1=(y/(yw2-yw1)*(yv2-yv1))*ysize
      write(1,100) x1,y1
      100 format(1x,1pe12.5,' cm ',1pe12.5,' cm translate')
      return
      end

   subroutine translateo(theta)
      common /paper/ xsize, ysize
      common /viewp/ xv1, yv1, xv2, yv2
      common /world/ xw1, yw1, xw2, yw2
      real(4) theta
      real(4) 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)
      write(1,100) xx-x1,yy-y1
      100 format(1x,1pe12.5,' cm ',1pe12.5,' cm translate')
      return
      end

   subroutine triangle(x1, y1, d1, a1)
      real(4) x1, y1, d1, a1
      real(4) 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
      call newpath
      call plot(x1+r1*cos(aa1),y1+r1*sin(aa1),3)
      call plot(x1+r1*cos(aa2),y1+r1*sin(aa2),2)
      call plot(x1+r1*cos(aa3),y1+r1*sin(aa3),2)
      call closepath
      return
      end

   subroutine damper(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4) s,ex1,ey1,ex2,ey2
      real(4) Ax1,Ay1,Ax2,Ay2,Ax3,Ay3,Ax4,Ay4,Ax5,Ay5
      real(4) Bx1,By1,Bx2,By2,Bx3,By3
      s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
      ex1=(x2-x1)/s
      ey1=(y2-y1)/s
      ex2=-(y2-y1)/s
      ey2=(x2-x1)/s
      Ax1=x1+0.45*s*ex1
      Ay1=y1+0.45*s*ey1
      Ax2=Ax1+0.5*d*ex2
      Ay2=Ay1+0.5*d*ey2
      Ax3=Ax1-0.5*d*ex2
      Ay3=Ay1-0.5*d*ey2
      Ax4=Ax2+0.15*s*ex1
      Ay4=Ay2+0.15*s*ey1
      Ax5=Ax3+0.15*s*ex1
      Ay5=Ay3+0.15*s*ey1
      Bx1=x2-0.45*s*ex1
      By1=y2-0.45*s*ey1
      Bx2=Bx1+0.4*d*ex2
      By2=By1+0.4*d*ey2
      Bx3=Bx1-0.4*d*ex2
      By3=By1-0.4*d*ey2
      call plot(x1,y1,3)
      call plot(Ax1,Ay1,2)
      call plot(Ax2,Ay2,2)
      call plot(Ax4,Ay4,2)
      call plot(Ax1,Ay1,3)
      call plot(Ax3,Ay3,2)
      call plot(Ax5,Ay5,2)
      call plot(x2,y2,3)
      call plot(Bx1,By1,2)
      call plot(Bx2,By2,3)
      call plot(Bx3,By3,2)
      return
      end

   subroutine walll(x1, y1, x2, y2, n)
      real(4) x1, y1, x2, y2
      integer n
      integer i
      real(4) pi
      real(4) s,ex1,ey1,ex2,ey2
      real(4) x11, y11, x22, y22
      pi=acos(-1.0)
      s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
      ex1=(x2-x1)/s
      ey1=(y2-y1)/s
      ex2=ex1*cos(135.0/180.0*pi)-ey1*sin(135.0/180.0*pi)
      ey2=ex1*sin(135.0/180.0*pi)+ey1*cos(135.0/180.0*pi)
      call plot(x1,y1,3)
      call plot(x2,y2,2)
      x11=x1
      y11=y1
      do i=1,n
         x11=x11+0.01*ex1
         y11=y11+0.01*ey1
         x22=x11+0.01*ex2
         y22=y11+0.01*ey2
         call plot(x11,y11,3)
         call plot(x22,y22,2)
      end do
      return
      end

   subroutine wallr(x1, y1, x2, y2, n)
      real(4) x1, y1, x2, y2
      integer n
      integer i
      real(4) pi
      real(4) s,ex1,ey1,ex2,ey2
      real(4) x11, y11, x22, y22
      pi=acos(-1.0)
      s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
      ex1=(x2-x1)/s
      ey1=(y2-y1)/s
      ex2=ex1*cos(-135.0/180.0*pi)-ey1*sin(-135.0/180.0*pi)
      ey2=ex1*sin(-135.0/180.0*pi)+ey1*cos(-135.0/180.0*pi)
      call plot(x1,y1,3)
      call plot(x2,y2,2)
      x11=x1
      y11=y1
      do i=1,n
         x11=x11+0.01*ex1
         y11=y11+0.01*ey1
         x22=x11+0.01*ex2
         y22=y11+0.01*ey2
         call plot(x11,y11,3)
         call plot(x22,y22,2)
      end do
      return
      end

   subroutine spring(x1, y1, x2, y2, d, n)
      real(4) x1, y1, x2, y2, d
      integer n
      integer i
      real(4) s,ex1,ey1,ex2,ey2,xx,yy
      s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
      ex1=(x2-x1)/s
      ey1=(y2-y1)/s
      ex2=-(y2-y1)/s
      ey2=(x2-x1)/s
      call plot(x1,y1,3)
      xx=x1+0.25*s*ex1
      yy=y1+0.25*s*ey1
      call plot(xx,yy,2)
      do i=1,n
         xx=xx+1.0/(4.0*real(n))*0.5*s*ex1+0.5*d*ex2
         yy=yy+1.0/(4.0*real(n))*0.5*s*ey1+0.5*d*ey2
         call plot(xx,yy,2)
         xx=xx+1.0/(4.0*real(n))*s*ex1-d*ex2
         yy=yy+1.0/(4.0*real(n))*s*ey1-d*ey2
         call plot(xx,yy,2)
         xx=xx+1.0/(4.0*real(n))*0.5*s*ex1+0.5*d*ex2
         yy=yy+1.0/(4.0*real(n))*0.5*s*ey1+0.5*d*ey2
         call plot(xx,yy,2)
      end do
      call plot(x2,y2,2)
      return
      end

   subroutine msmi(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4),parameter::arrowd=0.0075,bard=0.02
      real(4) s,ex1,ey1,ex2,ey2,xx,yy
      s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
      ex1=(x2-x1)/s
      ey1=(y2-y1)/s
      ex2=-(y2-y1)/s
      ey2=(x2-x1)/s
      xx=x1+d*ex1
      yy=y1+d*ey1
      call arrow(xx,yy,x1+0.025*s*ex1,y1+0.025*s*ey1,arrowd)
      xx=x1+bard*ex2
      yy=y1+bard*ey2
      call plot(xx,yy,3)
      call plot(x1,y1,2)
      xx=x1-bard*ex2
      yy=y1-bard*ey2
      call plot(xx,yy,3)
      call plot(x1,y1,2)
      xx=x2-d*ex1
      yy=y2-d*ey1
      call arrow(xx,yy,x2-0.025*s*ex1,y2-0.025*s*ey1,arrowd)
      xx=x2+bard*ex2
      yy=y2+bard*ey2
      call plot(xx,yy,3)
      call plot(x2,y2,2)
      xx=x2-bard*ex2
      yy=y2-bard*ey2
      call plot(xx,yy,3)
      call plot(x2,y2,2)
      return
      end

   subroutine msmo(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      real(4),parameter::arrowd=0.0075,bard=0.02
      real(4) s,ex1,ey1,ex2,ey2,xx,yy
      s=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
      ex1=(x2-x1)/s
      ey1=(y2-y1)/s
      ex2=-(y2-y1)/s
      ey2=(x2-x1)/s
      xx=x1-d*ex1
      yy=y1-d*ey1
      call arrow(xx,yy,x1-0.025*s*ex1,y1-0.025*s*ey1,arrowd)
      xx=x1+bard*ex2
      yy=y1+bard*ey2
      call plot(xx,yy,3)
      call plot(x1,y1,2)
      xx=x1-bard*ex2
      yy=y1-bard*ey2
      call plot(xx,yy,3)
      call plot(x1,y1,2)
      xx=x2+d*ex1
      yy=y2+d*ey1
      call arrow(xx,yy,x2+0.025*s*ex1,y2+0.025*s*ey1,arrowd)
      xx=x2+bard*ex2
      yy=y2+bard*ey2
      call plot(xx,yy,3)
      call plot(x2,y2,2)
      xx=x2-bard*ex2
      yy=y2-bard*ey2
      call plot(xx,yy,3)
      call plot(x2,y2,2)
      return
      end

   subroutine hatch(x1, y1, x2, y2, d)
      real(4) x1, y1, x2, y2, d
      integer i,n,m
      real(4) x,y
      call newpath
      call plot(x1,y1,3)
      call plot(x2,y1,2)
      call plot(x2,y2,2)
      call plot(x1,y2,2)
      call closepath
      call clipon
      call newpath
      n=(x2-x1)/d
      m=(y2-y1)/d
      do i=-m,n
         x=real(i)*d+x1
         y=y1
         call plot(x,y,3)
         call plot(x+(y2-y1),y+(y2-y1),2)
         call stroke
      end do
      call clipoff
      return
      end

   subroutine viewport(xv1d, yv1d, xv2d, yv2d)
      common /viewp/ xv1, yv1, xv2, yv2
      real(4) xv1d, yv1d, xv2d, yv2d
      xv1=xv1d
      yv1=yv1d
      xv2=xv2d
      yv2=yv2d
      return
      end

   subroutine xyworld(xw1d, yw1d, xw2d, yw2d)
      common /world/ xw1, yw1, xw2, yw2
      real(4) xw1d, yw1d, xw2d, yw2d
      xw1=xw1d
      yw1=yw1d
      xw2=xw2d
      yw2=yw2d
      return
      end