c      --------------------------------------------------
c      PostScript file name  /home/g3/ohta/sample/sub0.f
c      Created by T. Ogino and K. Ohta
c      STELab Nagoya University at April 1992
c      f77 -c filename.f   f77 fname1.o fname2.o   a.out
c
       subroutine xyopen
       write(10,10)
   10  format( "%!" )
       write(10,20)
   20  format( "%%Page: 1 1" )
       write(10,*) "/m { moveto } def"
       write(10,*) "/l { lineto } def"
       write(10,*) "/s { show } def"
       write(10,*) "/aws { awidthshow } def"
       write(10,*) "/sl { setlinewidth } def"
       write(10,*) "/sn { stroke newpath } def"
       write(10,*) "/tr { /Times-Roman findfont } def"
       write(10,*) "/sb { /Symbol findfont } def"
       write(10,*) "/sf { scalefont } def"
       write(10,*) "/se { setfont } def"
       write(10,*) "/ro { rotate } def"
       write(10,*) "/sd { setdash } def"
       write(10,*) "/st { stroke } def"
       write(10,*) "/sp { stroke showpage } def"
       write(10,*) "/tl { translate } def"
       write(10,*) "/sc { scale } def"
       write(10,*) "/np { newpath } def"
       write(10,*) "/cp { closepath } def"
       write(10,*) "/gs { gsave } def"
       write(10,*) "/gr { grestore } def"
       write(10,*) "/grc { grestore 0.0 0.0 0.0 setrgbcolor } def"
       write(10,*) "/sg { setgray } def"
       write(10,*) "/chs { sethsbcolor } def"
       write(10,*) "/crg { setrgbcolor } def"
       write(10,*) "/f { fill } def"
       write(10,*) "/cr { 0 360 arc } def"
       write(10,*) "/rt { arcto } def"
       write(10,*) "/cup { currentpoint } def"
       write(10,*) "/ctr { /Times-Roman findfont } def"
       write(10,*) "/cti { /Times-Italic findfont } def"
       write(10,*) "/ctb { /Times-Bold findfont } def"
       write(10,*) "/ctbi { /Times-BoldItalic findfont } def"
       write(10,*) "/chr { /Helvetica findfont } def"
       write(10,*) "/cho { /Helvetica-Oblique findfont } def"
       write(10,*) "/chb { /Helvetica-Bold findfont } def"
       write(10,*) "/chbo { /Helvetica-BoldOblique findfont } def"
       write(10,*) "/ccr { /Courier findfont } def"
       write(10,*) "/cco { /Courier-Oblique findfont } def"
       write(10,*) "/ccb { /Courier-Bold findfont } def"
       write(10,*) "/ccbo { /Courier-BoldOblique findfont } def"
       write(10,*) "/c {setrgbcolor 1 0 translate newpath 0 0 moveto
     &         0 1 lineto 1 1 lineto 1 0 lineto closepath fill} def"
       write(10,*) "/ch {sethsbcolor 1 0 translate newpath 0 0 moveto
     &         0 1 lineto 1 1 lineto 1 0 lineto closepath fill} def"
       write(10,*) "/krh { /Ryumin-Light-EUC-H findfont } def"
       write(10,*) "/krv { /Ryumin-Light-EUC-V findfont } def"
       write(10,*) "/kgh { /GothicBBB-Medium-EUC-H findfont } def"
       write(10,*) "/kgv { /GothicBBB-Medium-EUC-V findfont } def"
       return
       end
       subroutine plots(ipci,ity,isz,ino)
c      ipci=1 point, ipci=2 centimeter, ipci=3 inch
c      ity=1 vertical, ity=2 horizontal
c      isz=1 a4 size, isz=2 b4 size,  ino=page number
       common /cpci/lpci
       lpci=ipci
       if(ino.ge.2) write(10,10) ino,ino
   10  format("%%Page:",1x,i2,1x,i2)
c
       if(ipci.eq.1) then
       if(ity.eq.1.and.(isz.eq.1.or.isz.eq.2)) then
       write(10,*) " 1.00 1.00 sc"
       write(10,*) " 39.0 44.0 tl"
       else if(ity.eq.2.and.isz.eq.1) then
       write(10,*) " 1.00 1.00 sc"
       write(10,*) " 90.0 ro"
       write(10,*) " 45.0 -557.0 tl"
       else if(ity.eq.2.and.isz.eq.2) then
       write(10,*) " 1.00 1.00 sc"
       write(10,*) " 90.0 ro"
       write(10,*) " 42.0 -689.0 tl"
       end if
       else if(ipci.eq.2) then
       if(ity.eq.1.and.(isz.eq.1.or.isz.eq.2)) then
       write(10,*) " 28.35  28.35 sc"
       write(10,*) "  1.30  1.65 tl"
       else if(ity.eq.2.and.isz.eq.1) then
       write(10,*) " 90.0 ro"
       write(10,*) " 28.35  28.35 sc"
       write(10,*) " 1.70  -19.75 tl"
       else if(ity.eq.2.and.isz.eq.2) then
       write(10,*) " 90.0 ro"
       write(10,*) " 28.35  28.35 sc"
       write(10,*) " 1.70  -24.30 tl"
       end if
       else if(ipci.eq.3) then
       if(ity.eq.1.and.(isz.eq.1.or.isz.eq.2)) then
       write(10,*) " 72.00 72.00 sc"
       write(10,*) " 0.52 0.60 tl"
       else if(ity.eq.2.and.isz.eq.1) then
       write(10,*) " 90.0 ro"
       write(10,*) " 72.00 72.00 sc"
       write(10,*) " 0.62 -7.75 tl"
       else if(ity.eq.2.and.isz.eq.2) then
       write(10,*) " 90.0 ro"
       write(10,*) " 72.00 72.00 sc"
       write(10,*) " 0.63 -9.58 tl"
       end if
       end if
       call newpen(1)
       return
       end
       subroutine plote
       write(10,*) "sp"
       return
       end
       subroutine linee
       write(10,*) "st"
       return
       end
       subroutine dashe
       write(10,*) "sn [ ] 0 sd st"
       return
       end
       subroutine plot(x1,y1,im)
       common /cpci/lpci
       if(lpci.eq.1) then
       if(im.eq.3) write(10,11) x1,y1
       if(im.eq.2) write(10,12) x1,y1
       if(im.eq.-3) write(10,13) x1,y1
       if(im.eq.-2) write(10,14) x1,y1,x1,y1
   11  format(f6.1,1x,f6.1," m")
   12  format(f6.1,1x,f6.1," l")
   13  format(f6.1,1x,f6.1," tl")
   14  format(f6.1,1x,f6.1," l sn",2(1x,f6.1)," tl")
       else if(lpci.ge.2) then
       if(im.eq.3) write(10,21) x1,y1
       if(im.eq.2) write(10,22) x1,y1
       if(im.eq.-3) write(10,23) x1,y1
       if(im.eq.-2) write(10,24) x1,y1,x1,y1
   21  format(f6.2,1x,f6.2," m")
   22  format(f6.2,1x,f6.2," l")
   23  format(f6.2,1x,f6.2," tl")
   24  format(f6.2,1x,f6.2," l sn",2(1x,f6.2)," tl")
       end if
       return
       end
       subroutine factor(fct)
       write(10,10) fct,fct
   10  format(f6.3,1x,f6.3," sc")
       return
       end
       subroutine newpen(ipn)
       common /cpci/lpci
       if(lpci.eq.1) pn=0.350
       if(lpci.eq.2) pn=0.012
       if(lpci.eq.3) pn=0.005
       ip=ipn
       ip=max0(ip,1)
       ip=min0(ip,10)
       fp=float(ip)*pn
       write(10,30) fp
   30  format("sn",1x,f6.3," sl")
       return
       end
       subroutine newdsh(ipn)
       common /cpci/lpci
       dimension xa(7),ya(7)
       data xa/2.0,4.0,7.0,10.0,6.0,9.0,12.0/
       data ya/2.0,4.0,7.0,10.0,3.0,3.0, 3.0/
       if(lpci.eq.1) fl=1.00
       if(lpci.eq.2) fl=1.0/28.35
       if(lpci.eq.3) fl=1.0/72.0
       ip=max0(ipn,1)
       ip=min0(ipn,10)
       is=ip
       if(ip.gt.7) is=ip-3
       x1=xa(is)*fl
       y1=ya(is)*fl
       if(lpci.eq.1) then
       if(ip.le.7) write(10,40) x1,y1
   40  format("[",2(1x,f4.1),"] 0 sd")
       if(ip.gt.7) write(10,42) x1,y1,y1,y1
   42  format("[",4(1x,f4.1),"] 0 sd")
       else if(lpci.eq.2) then
       if(ip.le.7) write(10,50) x1,y1
   50  format("[",2(1x,f4.2),"] 0 sd")
       if(ip.gt.7) write(10,52) x1,y1,y1,y1
   52  format("[",4(1x,f4.2),"] 0 sd")
       else if(lpci.eq.3) then
       if(ip.le.7) write(10,50) x1,y1
       if(ip.gt.7) write(10,52) x1,y1,y1,y1
       end if
       return
       end
       subroutine symbol(x,y,h,isymb,ang,n)
       character isymb*80,ica*80,ich(80)*1
       equivalence (ica,ich)
       ica=isymb
       write(10,*) "gs"
       h0=h*1.50
       write(10,10) h0
   10  format("tr ",f6.2," sf se")
       if(x.gt.998.0.or.y.gt.998.0) then
       write(10,*) "cup m"
       else
       write(10,20) x,y
   20  format(f6.2,1x,f6.2," m")
       end if
       write(10,30) ang
   30  format(f6.1," ro")
       h1=h*0.37
       h2=h*0.06
       write(10,40) h1,h2
   40  format(f6.3," 0.0 8#040 ",f6.3," 0.0")
c       write(10,*) "(",(ich(i),i=1,n),") s"
       write(10,*) "(",(ich(i),i=1,n),") aws"
       write(10,*) "gr"
       return
       end
       subroutine stypet(ichr)
       common /cchr/lchr
       lchr=min0(ichr,4)
       lchr=max0(ichr,1)
       return
       end
       subroutine stypeh(ichr)
       common /cchr/lchr
       lchr=min0(ichr,4)+4
       lchr=max0(ichr,1)+4
       return
       end
       subroutine stypec(ichr)
       common /cchr/lchr
       lchr=min0(ichr,4)+8
       lchr=max0(ichr,1)+8
       return
       end
       subroutine stypek(ichr)
       common /kchr/lchr
       lchr=min0(ichr,4)
       lchr=max0(ichr,1)
       return
       end
       subroutine symblb(x,y,h,isymb,ang,n)
       character isymb*80,ica*80,ich(80)*1,icb(12)*4
       equivalence (ica,ich)
       common /cchr/lchr
       data icb/"ctr ","cti ","ctb ","ctbi","chr ","cho ",
     &          "chb ","chbo","ccr ","cco ","ccb ","ccbo"/
       ica=isymb
       write(10,*) "gs"
       if(lchr.eq.0) write(10,*) "tr"
       if(lchr.ge.1) write(10,*) icb(lchr)
       h0=h*1.50
       write(10,10) h0
   10  format(f5.2," sf se")
       if(x.gt.998.0.or.y.gt.998.0) then
       write(10,*) "cup m"
       else
       write(10,20) x,y
   20  format(f6.2,1x,f6.2," m")
       end if
       write(10,30) ang
   30  format(f6.1," ro")
       if(lchr.le.8) then
       h1=h*0.37
       h2=h*0.06
       write(10,40) h1,h2
   40  format(f6.3," 0.0 8#040 ",f6.3," 0.0")
       write(10,*) "(",(ich(i),i=1,n),") aws"
       else if(lchr.ge.9) then 
       write(10,*) "(",(ich(i),i=1,n),") s"
       end if
       write(10,*) "gr"
       return
       end
       subroutine symblc(x,y,h,isymb,ang,n)
       character isymb*80,ica*80,ich(80)*1
       equivalence (ica,ich)
       ica=isymb
       write(10,*) "gs"
       h0=h*1.50
       write(10,10) h0
   10  format("sb ",f6.2," sf se")
       if(x.gt.998.0.or.y.gt.998.0) then
       write(10,*) "cup m"
       else
       write(10,20) x,y
   20  format(f6.2,1x,f6.2," m")
       end if
       write(10,30) ang
   30  format(f6.1," ro")
       h1=h*0.37
       h2=h*0.06
       write(10,40) h1,h2
   40  format(f6.3," 0.0 8#040 ",f6.3," 0.0")
c       write(10,*) "(",(ich(i),i=1,n),") s"
       write(10,*) "(",(ich(i),i=1,n),") aws"
       write(10,*) "gr"
       return
       end
       subroutine symblk(x,y,h,isymb,ang,n)
       common /kchr/lchr
       character isymb*80,ica*80,ich(80)*1,icb(4)*3
       equivalence (ica,ich)
       data icb/"krh","kgh","krv","kgv"/
       ica=isymb
       write(10,*) "gs"
       if(lchr.eq.0) write(10,*) "krh"
       if(lchr.ge.1) write(10,*) icb(lchr)
       h0=h*1.50
       write(10,10) h0
   10  format(f6.2," sf se")
       write(10,20) x,y
   20  format(f6.2,1x,f6.2," m")
       write(10,30) ang
   30  format(f6.1," ro")
       h1=h*0.37
       h2=h*0.06
c       write(10,40) h1,h2
c   40  format(f6.3," 0.0 8#040 ",f6.3," 0.0")
       write(10,*) "(",(ich(i),i=1,n),") s"
c       write(10,*) "(",(ich(i),i=1,n),") aws"
       write(10,*) "gr"
       return
       end
       subroutine circ1(x,y,r)
       x1=x+r
       write(10,10) x1,y
   10  format(f6.2,1x,f6.2," m")
       write(10,20) x,y,r
   20  format(3(1x,f6.2)," 0 360 arc")
       return
       end
       subroutine circ2(x,y,r)
       x1=x+r
       write(10,10) x1,y
   10  format(f6.2,1x,f6.2," m")
       write(10,20) x,y,r
   20  format(2(f6.2),f6.3," cr")
       return
       end
       subroutine circl(x,y,angs,ange,r,r,dir)
       common /cpci/lpci
       write(10,10) x,y
   10  format(f6.2,1x,f6.2," m")
       if(dir.lt.0.05) then
       write(10,20) x,y,r,angs,ange
   20  format(3(1x,f6.2),2(1x,f6.1)," arc")
       write(10,30) x,y
   30  format(f6.2,1x,f6.2," l")
       else if(dir.gt.0.45) then
       if(lpci.eq.1) write(10,*) "sn [1.50 1.50] 0 sd"
       if(lpci.eq.2) write(10,*) "sn [0.05 0.05] 0 sd"
       if(lpci.eq.3) write(10,*) "sn [0.02 0.02] 0 sd"
       write(10,20) x,y,r,angs,ange
       write(10,30) x,y
       rad=3.1415926/180.0
       x1=x+r*cos(angs*rad)
       y1=y+r*sin(angs*rad)
       write(10,30) x1,y1
       write(10,*) "sn [ ] 0 sd"
       end if
       return
       end
       subroutine frame(xl,yl,xu,yu)
       call plot(xl,yl,3)
       call plot(xl,yu,2)
       call plot(xu,yu,2)
       call plot(xu,yl,2)
       call plot(xl,yl,2)
       return
       end
       subroutine framer(xl,yl,xu,yu)
       common /cpci/lpci
       if(lpci.eq.1) pn=0.350
       if(lpci.eq.2) pn=0.012
       if(lpci.eq.3) pn=0.005
       xw=0.5*(xu-xl)
       yw=0.5*(yu-yl)
       r1=pn*(xw+yw)*5.0
       x1=xl+xw
       y1=yl+yw
       call plot(x1,yu,3)
       write(10,10) xl,yu,xl,y1,r1
       write(10,10) xl,yl,x1,yl,r1
       write(10,10) xu,yl,xu,y1,r1
       write(10,10) xu,yu,x1,yu,r1
   10  format(4f6.1,f6.2' rt')
       write(10,*) "cp"
       return
       end
       subroutine rect(x,y,h,w,ang,ipn)
       ar=ang*3.1415926/180.0
       as=sin(ar)
       ac=cos(ar)
       xw=w*ac
       yw=w*as
       xh=h*as
       yh=h*ac
       x1=x
       y1=y
       x2=x-xh
       y2=y+yh
       x3=x+xw-xh
       y3=y+yw+yh
       x4=x+xw
       y4=y+yw
       call plot(x1,y1,ipn)
       call plot(x2,y2,2)
       call plot(x3,y3,2)
       call plot(x4,y4,2)
       call plot(x1,y1,2)
       return
       end
       subroutine rectr(x,y,h,w,ang,ipn)
       common /cpci/lpci
       if(lpci.eq.1) pn=0.350
       if(lpci.eq.2) pn=0.012
       if(lpci.eq.3) pn=0.005
       r1=pn*(w+h)*5.0
       ar=ang*3.1415926/180.0
       as=sin(ar)
       ac=cos(ar)
       xw=w*ac
       yw=w*as
       xh=h*as
       yh=h*ac
       x1=x
       y1=y
       x2=x-xh
       y2=y+yh
       x3=x+xw-xh
       y3=y+yw+yh
       x4=x+xw
       y4=y+yw
       x23=0.5*(x2+x3)
       y23=0.5*(y2+y3)
       x12=0.5*(x1+x2)
       y12=0.5*(y1+y2)
       x14=0.5*(x1+x4)
       y14=0.5*(y1+y4)
       x34=0.5*(x3+x4) 
       y34=0.5*(y3+y4) 
       call plot(x23,y23,ipn)
       write(10,10) x2,y2,x12,y12,r1
       write(10,10) x1,y1,x14,y14,r1
       write(10,10) x4,y4,x34,y34,r1
       write(10,10) x3,y3,x23,y23,r1
   10  format(4f7.2,f7.3' rt')
       write(10,*) "cp"
       return
       end
       subroutine grid(x0,y0,dx,dy,nx,ny)
       dimension dlx(10),dly(10)
       nxa=abs(nx)
       nya=abs(ny)
       do 10 i=1,10
       dlx(i)=dx*log10(float(i))
       dly(i)=dy*log10(float(i))
   10  continue
       x3=x0+dx*float(nxa)
       y3=y0+dy*float(nya)
       nx1=nxa+1
       ny1=nya+1
       if(nx.gt.0) then
       do 20 i=1,nx1
       x1=x0+dx*float(i-1)
       call plot(x1,y0,3)
       call plot(x1,y3,2)
   20  continue
       else if(nx.lt.0) then
       do 22 i=1,nxa
       do 22 j=1,10
       x1=x0+dx*float(i-1)+dlx(j)
       call plot(x1,y0,3)
       call plot(x1,y3,2)
   22  continue
       endif
       if(ny.gt.0) then
       do 30 i=1,ny1
       y1=y0+dy*float(i-1)
       call plot(x0,y1,3)
       call plot(x3,y1,2)
   30  continue
       else if(ny.lt.0) then
       do 32 i=1,nya
       do 32 j=1,10
       y1=y0+dy*float(i-1)+dly(j)
       call plot(x0,y1,3)
       call plot(x3,y1,2)
   32  continue
       endif
       return
       end
       subroutine colb
       write(10,*) "gs"
       return
       end
       subroutine cole
       write(10,*) "gr"
       return
       end
       subroutine newhsb(h,s,b)
       ah=amax1(h,0.0)
       as=amax1(s,0.0)
       ab=amax1(b,0.0)
       ah=amin1(h,1.0)
       as=amin1(s,1.0)
       ab=amin1(b,1.0)
       write(10,10) ah,as,ab
   10  format(3(1x,f4.2)," chs")
       return
       end
       subroutine newrgb(r,g,b)
       ar=amax1(r,0.0)
       ag=amax1(g,0.0)
       ab=amax1(b,0.0)
       ar=amin1(r,1.0)
       ag=amin1(g,1.0)
       ab=amin1(b,1.0)
       write(10,10) ar,ag,ab
   10  format(3(1x,f4.2)," crg")
       return
       end
       subroutine newgry(g)
c      0.0 black, 1.0 white
       fg=amax1(g,0.0)
       fg=amin1(g,1.0)
       write(10,10) fg
   10  format(f4.2," sg")
       return
       end
       subroutine paint
       write(10,*) "f"
       return
       end
       subroutine number(x,y,h,fnp,ang,n)
       character fmt1*6,fmt2*7,fmt3*4,fmt4*5,fna*20
       fnp1=fnp
       fnpa=abs(fnp1)+1.0e-30
       keta=1
       if(fnpa.gt.1.0) keta=ifix(log10(fnpa)*1.00001)+1
       if(fnpa.lt.0.01) keta=-(ifix(log10(fnpa)*1.00001))-1
       if(fnpa.lt.1.0e-9) keta=1
       if(abs(fnp).lt.1.0e-9) fnp1=1.0e-9
       minus=0
       if(fnp1.lt.0.0) minus=1
       if(n.lt.0) go to 1
       n1=keta+n+1+minus
       if(n1.le.0) n1=1
       if(n1.ge.10) go to 2
       write(fmt1,100) n1,n
  100  format('(f',i1,'.',i1,')')
       write(fna,fmt1) fnp1
       call symbol(x,y,h,fna,ang,n1)
    2  if(n1.lt.10) go to 1
       write(fmt2,110) n1,n
  110  format('(f',i2,'.',i1,')')
       write(fna,fmt2) fnp1
       call symbol(x,y,h,fna,ang,n1)
c
    1  if(n.ge.0) go to 3
       if(fnp1.gt.0.0) ifnp=ifix(fnp1*(10.0**(n+1))+0.5)
       if(fnp1.lt.0.0) ifnp=ifix(fnp1*(10.0**(n+1))-0.5)
       n2=keta+n+1+minus
       if(n2.le.0) n2=1
       if(n2.ge.10) go to 4
       write(fmt3,200) n2
  200  format('(i',i1,')')
       write(fna,fmt3) ifnp
       call symbol(x,y,h,fna,ang,n2)
    4  if(n2.lt.10) go to 3
       write(fmt4,210) n2
  210  format('(i',i2,')')
       write(fna,fmt4) ifnp
       call symbol(x,y,h,fna,ang,n2)
    3  return
       end
       subroutine color1(nx,ny,xb,yb,xl,yl,u,ico)
       dimension u(1)
       dx=xl/float(nx)
       dy=yl/float(ny)
c
       write(10,*) "gs"
       write(10,20) dx,dy
   20  format(f7.4,1x,f7.4," sc")
c
       ya=yb/dy-1.0
       xa=xb/dx-1.0+float(nx)
       write(10,22) xa,ya
c
       do 100 j=1,ny
       y=1.0
       x=-float(nx)
       write(10,22) x,y
   22  format(f7.2,1x,f7.2," tl")
       do 100 i=1,nx
       i1=i+nx*(j-1)
       a1=u(i1)
c
       call colora(a1,r,g,b,ico)
c
       write(10,24) r,g,b
   24  format(f4.2,1x,f4.2,1x,f4.2," c")
c*   24  format(f4.2,1x,f4.2,1x,f4.2," ch")
  100  continue
c
       x1=-float(nx)
       y1=-float(ny)
       write(10,22) x1,y1
c
       x1=1.0
       y1=1.0
       x2=x1+float(nx)
       y2=y1+float(ny)
       call newrgb(1.0,0.0,0.0)
       call plot(x1,y1,3)
       call plot(x1,y2,2)
       call plot(x2,y2,2)
       call plot(x2,y1,2)
       call plot(x1,y1,2)
       write(10,*) "st"
       write(10,*) "grc"
       return
       end
       subroutine imagc1(nx,ny,xb,yb,xl,yl,u,ico)
c      image color 4bit  limits 360*360
       parameter (nc=60)
       dimension u(1),ird(nc),igd(nc),ibd(nc)
c
       write(10,*) "gs"
       write(10,22) xb,yb
   22  format(f7.2,1x,f7.2," tl")
       write(10,20) xl,yl
   20  format(f7.4,1x,f7.4," sc")
c
       do 100 k=1,3
       do 100 j=1,ny
       do 100 i=1,nx
       i1=i+nx*(j-1)
       a1=u(i1)
c
       call colora(a1,r,g,b,ico)
c
       i2=mod(i1-1,nc)+1
       ird(i2)=r*15.0
       igd(i2)=g*15.0
       ibd(i2)=b*15.0
c
       if(k.eq.1) then
       if(j.eq.1.and.i.eq.1) write(10,51)
   51  format("/Red <")
       if(mod(i1,nc).eq.0) write(10,52) ird
   52  format(60z1)
       if(j.eq.ny.and.i.eq.nx) write(10,53) ird 
   53  format(60z1,/,"> def")
       else if(k.eq.2) then
       if(j.eq.1.and.i.eq.1) write(10,55)
   55  format("/Green <")
       if(mod(i1,nc).eq.0) write(10,52) igd
       if(j.eq.ny.and.i.eq.nx) write(10,53) igd
       else if(k.eq.3) then
       if(j.eq.1.and.i.eq.1) write(10,57)
   57  format("/Blue <")
       if(mod(i1,nc).eq.0) write(10,52) ibd
       if(j.eq.ny.and.i.eq.nx) write(10,53) ibd
       endif
  100  continue
c
       write(10,61) nx,ny
   61  format(i3,1x,i3," 4")
c      write(10,62) nx,-ny,ny
c  62  format("[",i3," 0 0 ",i4," 0 ",i3,"]")
       write(10,62) nx,ny
   62  format("[",i3," 0 0 ",i3," 0 0]")
       write(10,63)
   63  format("{Red} {Green} {Blue}",/,"true 3 colorimage")      
c
       x1=0.0
       y1=0.0
       x2=1.0
       y2=1.0
c      call newpen(1)
c      call newrgb(1.0,0.0,0.0)
c      call plot(x1,y1,3)
c      call plot(x1,y2,2)
c      call plot(x2,y2,2)
c      call plot(x2,y1,2)
c      call plot(x1,y1,2)
c      write(10,*) "st"
       write(10,*) "grc"
       return
       end
       subroutine imagc2(nx,ny,xb,yb,xl,yl,u,ico)
c      image color 8bit  limits 360*360
       parameter (nc=30)
       dimension u(1),ird(nc),igd(nc),ibd(nc)
c
       write(10,*) "gs"
       write(10,22) xb,yb
   22  format(f7.2,1x,f7.2," tl")
       write(10,20) xl,yl
   20  format(f7.4,1x,f7.4," sc")
c
       do 100 k=1,3
       do 100 j=1,ny
       do 100 i=1,nx
       i1=i+nx*(j-1)
       a1=u(i1)
c
       call colora(a1,r,g,b,ico)
c
       i2=mod(i1-1,nc)+1
       ird(i2)=r*255.0
       igd(i2)=g*255.0
       ibd(i2)=b*255.0
c
       if(k.eq.1) then
       if(j.eq.1.and.i.eq.1) write(10,51)
   51  format("/Red <")
       if(mod(i1,nc).eq.0) write(10,52) ird
   52  format(30z2.2)
       if(j.eq.ny.and.i.eq.nx) write(10,53) ird 
   53  format(30z2.2,/,"> def")
       else if(k.eq.2) then
       if(j.eq.1.and.i.eq.1) write(10,55)
   55  format("/Green <")
       if(mod(i1,nc).eq.0) write(10,52) igd
       if(j.eq.ny.and.i.eq.nx) write(10,53) igd
       else if(k.eq.3) then
       if(j.eq.1.and.i.eq.1) write(10,57)
   57  format("/Blue <")
       if(mod(i1,nc).eq.0) write(10,52) ibd
       if(j.eq.ny.and.i.eq.nx) write(10,53) ibd
       endif
  100  continue
c
       write(10,61) nx,ny
   61  format(i3,1x,i3," 8")
c      write(10,62) nx,-ny,ny
c  62  format("[",i3," 0 0 ",i4," 0 ",i3,"]")
       write(10,62) nx,ny
   62  format("[",i3," 0 0 ",i3," 0 0]")
       write(10,63)
   63  format("{Red} {Green} {Blue}",/,"true 3 colorimage")      
c
       x1=0.0
       y1=0.0
       x2=1.0
       y2=1.0
c      call newpen(1)
c      call newrgb(1.0,0.0,0.0)
c      call plot(x1,y1,3)
c      call plot(x1,y2,2)
c      call plot(x2,y2,2)
c      call plot(x2,y1,2)
c      call plot(x1,y1,2)
c      write(10,*) "st"
       write(10,*) "grc"
       return
       end
       subroutine imagc3(nx,ny,xb,yb,xl,yl,u,ico)
c      image color 4bit  limits 720*720
       parameter (nc=60,n1=360)
       dimension u(1),ird(nc),igd(nc),ibd(nc)
c
       nxn=nx/n1+1
       nyn=ny/n1+1
c
       do 200 jj=1,nyn
       if(jj.eq.1) then
       ny2=min0(ny,n1)
       nya=1
       nyb=ny2
       yba=yb
       yla=yl*float(ny2)/float(ny)
       else if(jj.eq.2) then
       ny2=ny-n1
       nya=n1+1
       nyb=ny
       yba=yb+yl*float(n1)/float(ny)
       yla=yl*float(ny-nya)/float(ny)
       endif
c
       do 200 ii=1,nxn
       if(ii.eq.1) then
       nx2=min0(nx,n1)
       nxa=1
       nxb=nx2
       xba=xb
       xla=xl*float(nx2)/float(nx)
       else if(ii.eq.2) then
       nx2=nx-n1
       nxa=n1+1
       nxb=nx
       xba=xb+xl*float(n1)/float(nx)
       xla=xl*float(nx-nxa)/float(nx)
       endif
c
       write(10,*) "gs"
       write(10,22) xba,yba
   22  format(f7.2,1x,f7.2," tl")
       write(10,20) xla,yla
   20  format(f7.4,1x,f7.4," sc")
c
       do 100 k=1,3
       num=0
       do 100 j=nya,nyb
       do 100 i=nxa,nxb
       i1=i+nx*(j-1)
       a1=u(i1)
       num=num+1
c
       call colora(a1,r,g,b,ico)
c
       i2=mod(num-1,nc)+1
       ird(i2)=r*15.0
       igd(i2)=g*15.0
       ibd(i2)=b*15.0
c
       if(k.eq.1) then
       if(j.eq.nya.and.i.eq.nxa) write(10,51)
   51  format("/Red <")
       if(mod(num,nc).eq.0) write(10,52) ird
   52  format(60z1)
       if(j.eq.nyb.and.i.eq.nxb) write(10,53) ird 
   53  format(60z1,/,"> def")
       else if(k.eq.2) then
       if(j.eq.nya.and.i.eq.nxa) write(10,55)
   55  format("/Green <")
       if(mod(num,nc).eq.0) write(10,52) igd
       if(j.eq.nyb.and.i.eq.nxb) write(10,53) igd
       else if(k.eq.3) then
       if(j.eq.nya.and.i.eq.nxa) write(10,57)
   57  format("/Blue <")
       if(mod(num,nc).eq.0) write(10,52) ibd
       if(j.eq.nyb.and.i.eq.nxb) write(10,53) ibd
       endif
  100  continue
c
       write(10,61) nx2,ny2
   61  format(i3,1x,i3," 4")
c      write(10,62) nx2,-ny2,ny2
c  62  format("[",i3," 0 0 ",i4," 0 ",i3,"]")
       write(10,62) nx2,ny2
   62  format("[",i3," 0 0 ",i3," 0 0]")
       write(10,63)
   63  format("{Red} {Green} {Blue}",/,"true 3 colorimage")      
       write(10,*) "grc"
  200  continue
c
       x1=0.0
       y1=0.0
       x2=1.0
       y2=1.0
c      call newpen(1)
c      call newrgb(1.0,0.0,0.0)
c      call plot(x1,y1,3)
c      call plot(x1,y2,2)
c      call plot(x2,y2,2)
c      call plot(x2,y1,2)
c      call plot(x1,y1,2)
c      write(10,*) "st"
c      write(10,*) "grc"
       return
       end
       subroutine image1(nx,ny,xb,yb,xl,yl,u,ico)
c      image gray 4bit  limits 360*360
       parameter (nc=60)
       dimension u(1),iwd(nc)
c
       write(10,*) "gs"
       write(10,22) xb,yb
   22  format(f7.2,1x,f7.2," tl")
       write(10,20) xl,yl
   20  format(f7.4,1x,f7.4," sc")
c
       a0=1.05
       b0=0.05
c
       write(10,*) "<"
       do 100 j=1,ny
       do 100 i=1,nx
       i1=i+nx*(j-1)
       x1=u(i1)
c
       i2=mod(i1-1,nc)+1
       w=a0-(a0-b0)*x1
       w=amax1(w,0.0)
       w=amin1(w,1.0)
       iwd(i2)=w*15.0 
       if(mod(i1,nc).eq.0) write(10,52) iwd
       if(j.eq.ny.and.i.eq.nx) write(10,52) iwd
   52  format(60z1)
c
  100  continue
c
       write(10,*) ">"
       write(10,61) nx,ny
   61  format(i3,1x,i3," 4")
c      write(10,62) nx,-ny,ny
c  62  format("[",i3," 0 0 ",i4," 0 ",i3,"]")
       write(10,62) nx,ny
   62  format("[",i3," 0 0 ",i3," 0 0]")
       write(10,63)
   63  format("{ } image")
       write(10,*) "grc"
c
       return
       end
       subroutine image2(nx,ny,xb,yb,xl,yl,u,ico)
c      image gray 8bit  limits 360*360
       parameter (nc=30)
       dimension u(1),iwd(nc)
c
       write(10,*) "gs"
       write(10,22) xb,yb
   22  format(f7.2,1x,f7.2," tl")
       write(10,20) xl,yl
   20  format(f7.4,1x,f7.4," sc")
c
       a0=1.05
       b0=0.05
c
       write(10,*) "<"
       do 100 j=1,ny
       do 100 i=1,nx
       i1=i+nx*(j-1)
       x1=u(i1)
c
       i2=mod(i1-1,nc)+1
       w=a0-(a0-b0)*x1
       w=amax1(w,0.0)
       w=amin1(w,1.0)
       iwd(i2)=w*255.0 
       if(mod(i1,nc).eq.0) write(10,52) iwd
       if(j.eq.ny.and.i.eq.nx) write(10,52) iwd
   52  format(30z2.2)
c
  100  continue
c
       write(10,*) ">"
       write(10,61) nx,ny
   61  format(i3,1x,i3," 8")
c      write(10,62) nx,-ny,ny
c  62  format("[",i3," 0 0 ",i4," 0 ",i3,"]")
       write(10,62) nx,ny
   62  format("[",i3," 0 0 ",i3," 0 0]")
       write(10,63)
   63  format("{ } image")
       write(10,*) "grc"
c
       return
       end
       subroutine image3(nx,ny,xb,yb,xl,yl,u,ico)
c      image gray 4bit  limits 720*720
       parameter (nc=60,n1=360)
       dimension u(1),iwd(nc)
c
       nxn=nx/n1+1
       nyn=ny/n1+1
c
       do 200 jj=1,nyn
       if(jj.eq.1) then
       ny2=min0(ny,n1)
       nya=1
       nyb=ny2
       yba=yb
       yla=yl*float(ny2)/float(ny)
       else if(jj.eq.2) then
       ny2=ny-n1
       nya=n1+1
       nyb=ny
       yba=yb+yl*float(n1)/float(ny)
       yla=yl*float(ny-nya)/float(ny)
       endif
c
       do 200 ii=1,nxn
       if(ii.eq.1) then
       nx2=min0(nx,n1)
       nxa=1
       nxb=nx2
       xba=xb
       xla=xl*float(nx2)/float(nx)
       else if(ii.eq.2) then
       nx2=nx-n1
       nxa=n1+1
       nxb=nx
       xba=xb+xl*float(n1)/float(nx)
       xla=xl*float(nx-nxa)/float(nx)
       endif
c
       write(10,*) "gs"
       write(10,22) xba,yba
   22  format(f7.2,1x,f7.2," tl")
       write(10,20) xla,yla
   20  format(f7.4,1x,f7.4," sc")
c
       a0=1.05
       b0=0.05
c
       nom=0
       write(10,*) "<"
       do 100 j=nya,nyb
       do 100 i=nxa,nxb
       i1=i+nx*(j-1)
       x1=u(i1)
       num=num+1
c
       i2=mod(num-1,nc)+1
       w=a0-(a0-b0)*x1
       w=amax1(w,0.0)
       w=amin1(w,1.0)
       iwd(i2)=w*15.0 
       if(mod(num,nc).eq.0) write(10,52) iwd
       if(j.eq.nyb.and.i.eq.nxb) write(10,52) iwd
   52  format(60z1)
c
  100  continue
c
       write(10,*) ">"
       write(10,61) nx2,ny2
   61  format(i3,1x,i3," 4")
c      write(10,62) nx2,-ny2,ny2
c  62  format("[",i3," 0 0 ",i4," 0 ",i3,"]")
       write(10,62) nx2,ny2
   62  format("[",i3," 0 0 ",i3," 0 0]")
       write(10,63)
   63  format("{ } image")
       write(10,*) "grc"
  200  continue
c
       return
       end
       subroutine colora(a1,r,g,b,ico)
c
       a0=1.20
       if(ico.eq.1) go to 31
       if(ico.eq.2) go to 32
       if(ico.eq.3) go to 33
       if(ico.eq.4) go to 34
       if(ico.eq.5) go to 35
       if(ico.eq.6) go to 36
       if(ico.eq.7) go to 37
       if(ico.eq.8) go to 38
       if(ico.ge.10) go to 50
       r=a0-(a0-0.5)*a1
       g=a0-(a0-0.5)*a1
       b=a0-(a0-0.5)*a1
       go to 40
   31  continue
       r=1.0
       g=a0-a0*a1
       b=a0-a0*a1
       go to 40
   32  continue
       r=a0-a0*a1
       g=1.0
       b=a0-a0*a1
       go to 40
   33  continue
       r=a0-a0*a1
       g=a0-a0*a1
       b=1.0
       go to 40
   34  continue
       a0=2.0
       a1=a1*(a0+4.0)
       r=a0-a1
       if(a1.ge.a0) r=a1-(a0+1.0)
       g=a0+3.0-a1
       b=a0+1.0-a1
       if(a1.ge.a0+2.0) b=a1-(a0+3.0)
       go to 40
   35  continue
       a0=2.0
       a1=a1*(a0+4.0)
       r=a0+3.0-a1
       g=a0+1.0-a1
       if(a1.ge.a0+2.0) g=a1-(a0+3.0)
       b=a0-a1
       if(a1.ge.a0) b=a1-(a0+1.0)
       go to 40
   36  continue
       a0=2.0
       a1=a1*(a0+4.0)
       r=a0-a1
       if(a1.ge.a0) r=a1-(a0+1.0)
       g=a0+1.0-a1
       if(a1.ge.a0+2.0) g=a1-(a0+3.0)
       b=a0+3.0-a1
       go to 40
   37  continue
       a1=10.0*a1
       r=a1
       if(a1.ge.1.0) r=2.0-a1
       if(a1.ge.4.0) r=a1-4.0
       g=a1-2.0
       if(a1.ge.5.0) g=6.0-a1
       if(a1.ge.6.0) g=0.25*(a1-6.0)
       b=a1
       if(a1.ge.3.0) b=4.0-a1
       if(a1.ge.6.0) b=0.25*(a1-6.0)
       go to 40
   38  continue
       a1=7.0*a1
       r=a1-3.0
       g=a1-1.0
       if(a1.ge.4.0) g=5.0-a1
       if(a1.ge.5.0) g=0.5*(a1-5.0)
       b=a1
       if(a1.ge.2.0) b=3.0-a1
       if(a1.ge.5.0) b=0.5*(a1-5.0)
       go to 40
c
c
   50  continue
       ic1=ico/10
       ic2=mod(ico,10)
       data a02,a03,a04,a05,a06,a07/2.0,3.0,4.0,5.0,6.0,7.0/
c
       if(ic1.eq.1) go to 51
       if(ic1.eq.2) go to 52
       if(ic1.eq.3) go to 53
       if(ic1.eq.4) go to 54
       if(ic1.eq.5) go to 55
       if(ic1.eq.6) go to 56
       if(ic1.eq.7) go to 57
       if(ic1.eq.8) go to 58
       go to 40
c
   51  continue
c      color1
c
       a0=1.0
       a2=a1*a0*3.0
       r=0.0
       if(a2.gt.a0)  r=(a2-a0)
       if(a2.gt.a02) r=1.0
       g=0.0
       b=3.0*a1
       if(a2.gt.a0)  b=1.0
       if(a2.gt.a02) b=1.0-(a2-a02)
       if(ic2.eq.3) then
         g=3.0*a1
         if(a2.gt.a0)  g=1.0
         if(a2.gt.a02) g=1.0-(a2-a02)
         b=0.0
       endif
       go to 40
c
   52  continue
c      color2
c
       r=a1
       g=a1
       b=a1
       if(ic2.eq.1) then
       r=a1
       g=0.0
       b=0.0
       else if(ic2.eq.2) then
       r=0.0
       g=0.0
       b=a1
       else if(ic2.eq.3) then
       r=0.0
       g=a1
       b=0.0
       else if(ic2.eq.4) then
       r=a1
       g=0.0
       b=a1
       else if(ic2.eq.5) then
       r=a1
       g=a1
       b=0.0
       else if(ic2.eq.6) then
       r=0.0
       g=a1
       b=a1
       end if
       go to 40
c
   53  continue
c      color3
c
       a0=1.0
       a2=a1*a0*2.0
       r=0.0
       if(a2.gt.a0) r=(a2-a0)
       g=0.0
       b=0.0
       if(a2.le.a0) b=1.0-a2
       if(ic2.eq.3) then
         g=0.0
         if(a2.le.a0) g=1.0-a2
         b=0.0
       end if
       go to 40
c
   54  continue
c      color4
c
       r=a1
       g=0.0
       b=1.0-a1
       if(ic2.eq.1) then
         g=1.0-a1
         b=0.0
       else if(ic2.eq.2) then
         g=0.0
         b=1.0
       else if(ic2.eq.3) then
         g=1.0
         b=0.0
       else if(ic2.eq.4) then
         r=1.0
         g=0.0
         b=a1
       else if(ic2.eq.5) then
         r=1.0
         g=a1
         b=0.0
       end if
       go to 40
c
   55  continue
c      color5
c
       a0=1.0
       a2=a1*a0*3.0
       if(a2.le.a0) r=3.0*a1
       if(a2.gt.a0) r=1.0
       g=0.0
       if(a2.gt.a0)  g=3.0*a1-a0 
       if(a2.gt.a02) g=1.0
       b=0.0
       if(a2.gt.a02) b=3.0*a1-a02
       if(ic2.gt.0) then
         r=2.0*a1
         if(a2.gt.a0*1.5) r=1.0
         if(a2.gt.a0*1.5) g=2.0*a1-1.0
         b=0.0
       end if
       go to 40
c
   56  continue
c      color6
c
       a0=1.0
       a2=a1*a0*5.0
       r=0.0
       if(a2.gt.a02) r=a2-a02
       g=0.0
       if(a2.gt.a0 ) g=a2-a0
       if(a2.gt.a04) g=a05-a2
       b=a2
       if(a2.gt.a03) b=a04-a2
c
       if(ic2.eq.1) then
         r=0.0
         if(a2.gt.a03) r=a2-a03
         b=a2
         if(a2.gt.a02) b=a03-a2
       endif
       go to 40
c
   57  continue
c      color7
c
       a0=1.0
       a2=a1*a0*4.0
       r=0.0
       if(a2.gt.a02) r=a2-a02
       if(a2.gt.a03) r=1.0
       g=1.0-a2
       if(a2.gt.a0 ) g=0.0
       if(a2.gt.a03) g=a2-a03
       b=1.0
       if(a2.gt.a0 ) b=a02-a2
       if(a2.gt.a02) b=0.0
       if(ic2.eq.3) then
         g=1.0
         if(a2.gt.a0 ) g=a02-a2
         if(a2.gt.a02) g=0.0
         b=1.0-a2
         if(a2.gt.a0 ) b=0.0
         if(a2.gt.a03) b=a2-a03
       end if
       go to 40
c
   58  continue
c      color8
c
       a0=1.0
       a2=a1*a0*7.0
       r=a2-a03
       if(a2.gt.a04) r=a2-a04
       if(a2.gt.a05) r=a2-a05
       if(a2.gt.a06) r=a2-a06
       g=a2-a0
       if(a2.gt.a02) g=a2-a02
       if(a2.gt.a03) g=a2-a03
       if(a2.gt.a04) g=a2-a06
       b=a2
       if(a2.gt.a0 ) b=a2-a0
       if(a2.gt.a02) b=a2-a04
       if(a2.gt.a05) b=a2-a06
c
       am=1.20
       r=am*r
       g=am*g
       b=am*b
   40  continue
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       return
       end