program isophotes

c     contour plotting program
 
      PARAMETER (MAXSIZE=1024)

      character image*80, errmsg*80, title*40, yn, ok, chl
      character hcopy*80
      integer  	im, dtype, ier, hard
      integer*2 ipix(MAXSIZE)
      dimension pix(MAXSIZE)
      logical   there
      
      dimension GR(MAXSIZE,MAXSIZE)
      DIMENSION clev(30), TR(6)
      DATA TR/0.,1.,0.,0.,0.,1./

      lin=11

      call opninp (1,im,1,image,title,dtype,iext,jext,ier,*999)
      call clargc (2, hcopy, ier)
	 if ( ier.ne.0 ) hcopy='no'

      CALL MATRD(im,GR,dtype,iext,jext,pix,ipix,ier)


C     *** SELECT PROCESS ***
								       
      hard=0

C     Terminal
  777 continue
      call PGBEG(0,'/xserve',1,1)

C---- Set the user coordinates and make a coordinate box
 3481 close (lin)
      open (lin,file='cntpar',access='SEQUENTIAL')
      read (lin,*,end=888) clev(1)
      write(*,*)
      do i=2,30
      y   read (lin,*,end=529) clev(i)
      end do
  888 continue
      close(lin)
      write(*,*) ' Edit the file for contour levels '
      write(*,*) '   Give 1 level in 1 line '
      write(*,'(''   Press ENTER-key when ready'',$)')
      read (*,'(a)') yn
      call system('vi cntpar')
      chl='n'
      goto 3481
  529 close(lin)
      nlev=i-1
      if ( chl.eq.'y' ) then
	 call system('vi cntpar')
	 chl='n'
	 goto 3481
      end if

      call PGWNAD(0.,float(iext),0.,float(jext))
      call PGBOX('BCST',0.0,-1,'BCSTV',0.0,-1)
 1000 continue
      call PGLAB(title,' ',' ')
      call PGCONT(GR,iext,jext,1,iext,1,jext,clev,nlev,TR)

      if ( hard.eq.0 )then
	 write(*,'(a,$)') 'Change levels?  (y/) ===> '
	 read (*,'(a)') chl
	 if ( chl.eq.'y' )then
	    call PGEND
	    goto 777
	 end if
	 call PGEND

 2000    if ( hcopy.ne.'no' ) hard=1
	 if ( hard.eq.1 ) then
	    x2=float(iext)
	    y2=float(jext)
	    l=len(hcopy)
	    do i=l,1,-1
	       if ( hcopy(i:i).ne.' ' ) goto 1
	    end do
	    i=0
    1       LNBLNK=i
	    call PGBEG(0,hcopy(1:LNBLNK)//'/vps',1,1)
	    call PGSCF(2)
	    call PGWNAD(0.,x2,0.,y2,-1,-1)
	    call PGSLW(1)
	    call PGBOX('BCST',0.0,-1,'BCSTV',0.0,-1)
	    go to 1000
	 end if
      end if

      if ( hard.eq.1 ) then
C----    Finish the plot
	 call PGEND				! Finish the plot
	 write(*,*) '### Output(PostScript) was written in '//hcopy
      end if

 9000 call PGEND
      STOP
  999 call imemsg (ier, errmsg)
      write(*,'('' Error: '',a80)') errmsg
      stop
      END