subroutine gra_init

	call cho_choose('G')
	end


	subroutine gra_open(hardcopy, ok)

	logical hardcopy, ok

	character htype*32, dtype*32, file*128, buf*256, popt*1, dest*128
	character num*10
	integer lf, l, i, j, iostat
	integer hrddev/0/
	integer gradev/0/
	integer pgcnt/0/
	integer pan/0/
	integer filcnt/1/
	logical init/.true./
	logical xwin/.true./
	save htype, dtype, file, lf, hrddev, gradev, pgcnt, pan,init,xwin

	integer pgopen
	external gra_close, gra_print_exit

	ok=.false.
	if (hardcopy) then
	  if (gradev .ne. 0) then
	    call pgslct(gradev)
	    call pgclos
	    gradev=0
	  endif
	  if (hrddev .eq. 0) then
	    call cho_choose('P')
	    call sys_getenv('CHOOSER_POPT', popt)
	    if (popt .eq. 'L' .or. popt .eq. 'I') then
	      call sys_temp_name('gra', file)
	      call str_trim(file, file, lf)
	    else
	      call sys_getenv('CHOOSER_FILE', file)
	      call str_trim(file, file, lf)
	      if (popt .eq. 'F') then
	        if (file .eq. '?') then
	          file=' '
	          print '(x,a,$)', 'filename for graphics: '
	          read(*,'(a)', iostat=iostat) file
	          if (iostat .ne. 0) file=' '
	          call str_trim(file, file, lf)
	        else
	          buf=' '
	          j=0
	          do i=lf,1,-1
	            if (file(i:i) .eq. '.') then
	              j=lf-i+1
	              lf=i-1
	              buf=file(i:)
	              exit
	            endif
	          enddo
	          write(num, '(i10)') -filcnt
	          i=index(num,'-')
	          num(i:i)='_'
	          write(file(lf+1:), '(a,a)') num(i:), buf(1:max(1,j))
	          lf=lf+11-i+j
	        endif
	      endif
	      filcnt=filcnt+1
	      if (file(1:lf) .eq. ' ') then
	        lf=10
	        file(1:lf)='pgplot.dat'
	      endif
	    endif
	    call sys_getenv('CHOOSER_PAN', num)
	    pan=1
	    read(num, *, err=10, end=10) pan
10	    continue

	    call sys_getenv('CHOOSER_PDEV', htype)
	    if (htype .eq. ' ') then
	      print *,'unknown printer device'
	      return
	    endif
	    hrddev=pgopen(file(1:lf)//'/'//htype)
	    if (hrddev .le. 0) return

	    print *,'Opened print file ',file(1:lf)

	    call cho_calc_pan(pan, i, j)
	    pan=i*j
	    call pgsubp(i,j)
	    call pgpanl(1,1)
	    pgcnt=1
	  else
	    pgcnt=pgcnt+1
	    call pgslct(hrddev)
	    call pgpage
	  endif
	elseif (gradev .eq. 0) then
	  dtype(1:1)='/'
	  call sys_getenv('CHOOSER_GDEV', dtype(2:))
	  if (dtype .eq. '/') then
	    print *,'unknown display device'
	    return
	  endif
	  if (xwin .and. dtype(1:2) .eq. '/X') then
	    xwin=.false.
	    call sys_getenv('PG_START', buf)
	    if (buf .ne. ' ') call sys_cmd(buf)
	  endif
	  gradev=pgopen(dtype)
	  if (gradev .le. 0) then
	    print *,'cannot open display'
	    gradev=0
	    return
	  endif
	  call sys_err_hdl(gra_close)
	else
	  call pgpage	!cgt
	endif
	ok=.true.
	return


	entry gra_end_hcopy

	if (init) then
	  call sys_exit_hdl(gra_print_exit)
	  init=.false.
	endif
	if (mod(pgcnt,pan) .ne. 0) return
	call sys_getenv('CHOOSER_POPT', popt)
	if (popt .ne. 'F' .and. popt .ne. 'I') return
	goto 90


	entry gra_print

	call sys_getenv('CHOOSER_POPT', popt)
90	if (popt .ne. 'I' .and. popt .ne. 'L') then
	  if (hrddev .ne. 0) then
	    call pgslct(hrddev)
	    call pgclos
	    print *,'saved ',pgcnt,' graphs on ',(pgcnt-1)/pan+1,' pages'
	    print *,'on file ',file(1:lf)
	    hrddev=0
	  endif
	  return
	endif	  

	if (hrddev .ne. 0) then
	  call sys_getenv('CHOOSER_PCMD', buf)
	  call str_trim(buf, buf, l)
	  i=index(buf,'*')
	  call sys_getenv('CHOOSER_DEST', dest)
	  if (dest .eq. ' ' .and. i .ne. 0 .or. buf .eq. ' ') then
	    print *,'missing print command or destination for file '
	1           ,file(1:lf)
	  else
c	    if (htype .ne. 'PS' .and.
c	1       htype .ne. 'VPS' .and.
c	1       htype .ne. 'CPS' .and.
c	1       htype .ne. 'VCPS') return
	    if (i .ne. 0) then
	      call str_trim(dest, dest, j)
	      buf(i:)=dest(1:j)//buf(i+1:)
	      call str_trim(buf, buf, l)
	    endif
	    l=l+1
	    buf(l+1:)=file
	    l=l+lf
	    call pgslct(hrddev)
	    call pgclos
	    hrddev=0
	    print *,'print ',pgcnt,' graphs on ',(pgcnt-1)/pan+1,' pages'
	    print *,buf(1:l)
	    print *
	    call sys_cmd(buf(1:l))
	  endif
	endif
	return


	entry gra_close_it

	if (gradev .ne. 0) then
	  call pgslct(gradev)
	  call pgclos
	  if ((dtype(1:2) .ne. '/X')  .and.
	1      (dtype(1:3) .ne. '/GW') .and.
	1      (dtype(1:4) .ne. '/CGW'))
	1      call sys_rd_tmo(char(0), buf, l)
	  call sys_getenv('PG_SWITCH', buf)
	  if (buf .ne. ' ') call sys_cmd(buf)
	  gradev=0
	endif
	end


	subroutine gra_close
	call gra_close_it
	end


	subroutine gra_print_exit
	call gra_print
	end



	subroutine gra_mult(px, py, perr, mk, n, nstyl, styl, ncol)

C 24.2.95 M. Zolliker

	implicit none

	real px(*), py(*), perr(*)
	integer mk(*), nstyl, styl(*)
	integer n,ncol

	integer i,sty,set,maxlin
	logical pen

	real nan_value

	real gra_trf

	call pgbbuf

	call gra_get_nan_value(nan_value)

	maxlin=0
	do i=1,n
	  if (ncol .gt. 0) then
	    call pgsci(mod(mk(i)-1, ncol)+1)
	  endif
	  set=mk(i)
	  sty=set
	  if (sty .gt. nstyl) sty=mod(sty-1,nstyl)+1
	  sty=styl(sty)
	  if (abs(sty) .ge. 10 .and. set .gt. maxlin) maxlin=set
	  if (sty .eq. -8) then	! draw point
	    call pgpt1(px(i), gra_trf(py(i), mk(i)), -1)
	  else
	    if (sty .gt. 0) then				! draw error bars
	      call pgmove(px(i), gra_trf(py(i)+perr(i), mk(i)))
	      call pgdraw(px(i), gra_trf(py(i)-perr(i), mk(i)))
	    endif
	    sty=mod(abs(sty),10)
	    if (sty .lt. 8 .and. sty .ne. 0) then ! draw marker (do not draw for point or none)
	      call pgpt1(px(i), gra_trf(py(i), mk(i)), sty)
	    endif
	  endif
	enddo
	do set=1,maxlin
	  if (ncol .gt. 0) then
	    call pgsci(mod(set-1, ncol)+1)
	  endif
	  sty=set
	  if (sty .gt. nstyl) sty=mod(sty-1,nstyl)+1
	  if (abs(styl(sty)) .ge. 10) then
	    pen=.false.
	    do i=1,n
	      if (mk(i) .eq. set) then
	        if (py(i) .eq. nan_value) then
	          pen=.false.
	        elseif (pen) then
	          call pgdraw(px(i), gra_trf(py(i), mk(i)))
	        else
	          call pgmove(px(i), gra_trf(py(i), mk(i)))
	          pen=.true.
	        endif
	      endif
	    enddo
	  endif
	enddo

	if (ncol .gt. 0) then
	  call pgsci(1)
	endif
	call pgebuf
	end


	subroutine gra_get_nan_value(val)

	real val

	val=-0.5**31
	end


!!
	subroutine gra_auto_fmt(ndig, vx, erx, units
	1                      , n, inc, val, ince, err, pdig) !!
!!
!! determine number of digits needed after decimal point for a table of values and errors
!!
	implicit none

	real vx, erx            !! (out) graphical width of largest value/error
	integer units           !! pglen units for vx/erx
	real val(*), err(*)	!! values, errors
	integer ndig            !! (out) digits after decimal point
	integer n, inc, ince    !! number of values, index increment (of value, of error)
	integer pdig            !! number of digits to add if error is 0.0

	real vmax, emax, emin, ex, ey
	integer i, j, k, n1, n3, numl, l, lmax, mdig
	character*24 num

	vmax=0.0
	emin=10.0
	emax=0.0
	j=1
	k=1
	do i=1,n
	  if (abs(val(j)) .gt. vmax) then
	    vmax=abs(val(j))			! highest value
	  endif
	  if (err(k) .ne. 0) then
	    if (abs(err(k)) .gt. emax) emax=abs(err(k))	! highest error
	    if (abs(err(k)) .lt. emin) emin=abs(err(k))	! lowest error
	  elseif (abs(val(j)) .lt. emin .and. val(j) .ne. 0) then
	    emin=abs(val(j))			! if error=0 then check for value
	  endif
	  j=j+inc
	  k=k+ince
	enddo

	if (vmax .eq. 0.0) then
	  n3=0
	else
	  n3=int(3.3-log10(vmax))		! largest number min. 3.3 digits
	  if (emax .eq. 0.0) n3=n3+pdig		! add more digits if error undefined
	endif
	n1=int(1-log10(emin))			! smallest error min. 1 digit
	ndig=max(n1,n3)
	if (emax .gt. 0) then
	  n1=int(2-log10(emax))
	  if (n1 .gt. ndig) ndig=n1		! largest error min. 2 digits
	endif
	if (n3+5 .lt. ndig) ndig=n3+5		! largest number max. 8 digits
	if (ndig .gt. 7) ndig=7			! max. 7 digits after decimal point

! determine max. width of columns
	j=1
	k=1
	vx=0.0
	erx=0.0
	mdig=ndig
	do i=1,n
	  call cvt_real_str(num, numl, val(j), 1, ndig, 1, 0)
	  l=numl
	  lmax=numl-mdig
	  do while (l .gt. lmax .and. num(l:l) .eq. '0')
	    l=l-1
	  enddo
	  if (num(l:l) .eq. '.') l=l-1
	  mdig=numl-l
	  call pglen(units, num(1:numl), ex, ey)
	  if (ex .gt. vx) vx=ex
	  if (err(k) .ne. 0.0) then
	    call cvt_real_str(num, numl, err(k), 1, ndig, 1, 0)
	    call pglen(units, num(1:numl), ex, ey)
	    if (ex .gt. erx) erx=ex
	  endif
	  j=j+inc
	  k=k+ince
	enddo
	if (mdig .gt. 0 .and. erx .eq. 0) then
	  ndig=ndig-mdig
	  if (ndig .lt. 0) ndig=0
	  call pglen(units, num(numl-mdig+1:numl), ex, ey)
	  vx=vx-ex
	endif

	end

!!
	subroutine gra_show_par(gx, gy, columns, txt	 !!
	1                       , n, inc, val, ince, err) !!
!!
!! output parameter(s) 
!!
	implicit none

	real gx, gy		!! (in) start coordinates, (out) end coordinates (gx only)
	real val(*), err(*)	!! values, errors
	integer columns         !! 0: name above value(s), else: name to the left, number of columns
	integer n, inc, ince    !! number of values, index increment (of value, of error)
	character*(*) txt	!! name

	real ex, ey, top, vx, vx0, ve, erx, px, dx
	integer i, j, k, ndig, numl
	character*24 num
	real x1, x2, y1, y2

	call pgqvp(0, x1, x2, y1, y2)	! save viewport
	call pgsvp(0.0, 1.0, 0.0, 1.0)	! set viewport to max.

	call gra_auto_fmt(ndig, vx, erx, 0, n, inc, val, ince, err, 1)

	top=gy-1.0
	if (txt .ne. ' ') then ! show parameter name
	  call pgmtxt('T', top, gx, 0.0, txt)
	  if (columns .eq. 0) top=top-1.0
	endif

	call pglen(0, char(177), px, ey) ! length of +- character
	if (columns .ne. 0) then
	  dx=vx+erx+px*1.5
	  call pglen(0, txt, ex, ey)
	  vx=vx+ex+px
	else
	  vx=vx+px*0.1
	endif
	vx=gx+vx
	vx0=vx
	ve=vx+erx
	if (erx .gt. 0) ve=ve+px*1.5
	j=1
	k=1
	do i=1,n
	  call cvt_real_str(num, numl, val(j), 1, ndig, 1, 2)
	  call pgmtxt('T', top, vx, 1.0, num(1:numl))
	  if (err(k) .ne. 0) then
	    call pgmtxt('T', top, vx+px*0.75, 0.5, char(177))	! show +/- 
	    call cvt_real_str(num, numl, err(k), 1, ndig, 1, 2)
	    call pgmtxt('T', top, ve, 1.0, num(1:numl))
	  endif
	  if (columns .eq. 0) then
	    top=top-1.0
	  else
	    gx=ve+px
	    if (mod(i,columns) .eq. 0) then
	      vx=vx0
	      ve=vx+erx+px*1.5
	      top=top-1.0
	    else
	      vx=vx+dx
	      ve=ve+dx
	    endif
	  endif
	  j=j+inc
	  k=k+ince
	enddo
	call pgsvp(x1, x2, y1, y2)
	end


	subroutine gra_trfmode(modus, shift)

	integer modus
	real shift

	integer mode
	real step
	common /gra_com/mode, step

	data mode/0/
	data step/0.0/

	mode=max(0,min(1,modus))
	step=shift
	end


	real function gra_trf(y, n)

	real y
	integer n

	integer mode
	real step
	common /gra_com/mode, step

	if (mode .eq. 0) then
	  gra_trf=y+(n-1)*step
	elseif (mode .eq. 1) then ! logarithmic
	  gra_trf=log10(max(1.e-30,y))+(n-1)*step
	else
	  stop 'unimplemented transfer function'
	endif
	end


	real function gra_itrf(y)

	real y

	integer mode
	real step
	common /gra_com/mode, step

	if (mode .eq. 0) then
	  gra_itrf=y
	elseif (mode .eq. 1) then ! logarithmic
	  gra_itrf=10**(max(-30.,min(30.,y)))
	else
	  stop 'unimplemented transfer function'
	endif
	end