*************************************************************************
*			CALCOMP TO HIPLOT INTERFACE			*
*									*
* Developed by   Dr. Michael Peterson					*
*		     Remote Sensing Applications Lab			*
*		     University of Nebraska at Omaha			*
*		     Omaha, Nebraska   68182				*
*									*
*			(402) 554-2726					*
*									*
* 		This set of subroutines can be substituted for          *
*	the supplied Calcomp subroutines for use with HIPLOT   		*
*	compatible devices. ( Houston Instruments DMP series )		*
*									*
*	Known differences:						*
*									*
*		1) The subroutine OUTPUT must be used to initialize     *
*		   the output unit in the package.			*
*									*
*		2) In subroutines SYMBOL and NUMBER, the position       *
*		   999,999 in the parameter list do not work as         *
*		   the Calcomp version. (Since the location of the      *
*		   plotters pen is not applicable to HIPLOT devices)	*
*									*
*	To compile these subroutines with source code:     		*
*									*
*		   $ FORTRAN filename+HILIB				*
*		   $ LINK filename					*
*		   $ run filename					*
*			or						*
*		   $ FORTRAN filename					*
*		   $ LINK filename,HILIB				*
*		   $ RUN filename					*
*									*
      subroutine squares
      dimension x(4),y(4),xd(4),yd(4)
c..
c .. initialize first square
      data x/4.,4.,-4.,-4./
      data y/4.,-4.,-4.,4./
c ..
c .. initialize plot and go to center 
      call output (ifile)
      call plots (0,0,ifile)
      call plot(5.0,5.0,-3)
c .. 
c .. set mu value and produce 20 squares
      smu=0.1
      rmu=1.0-smu
c
      do 3 i=1,21
        call plot(x(4),y(4),3)
c ..
c .. draw the square and calculate the coordinates of the next square
        do 1 j=1,4
          call plot(x(j),y(j),2)
          nj=mod(j,4)+1
          xd(j)=rmu*x(j)+smu*x(nj)
          yd(j)=rmu*y(j)+smu*y(nj)
  1     continue
c
c .. reset square coordinates
        do 2 j=1,4
          x(j)=xd(j)
          y(j)=yd(j)
  2     continue
  3   continue
c
      call efplot
      return
      end
      subroutine spiro 
      dimension x(100), y(100)
      call output (ifile)
      write(6,*) '  Enter the number of vertices:'
      read(5,*) n
      call plots(0,0,ifile)
c ..
c .. center plot
c ..
      thinc=6.283185307/float(n)
      theta=0.0
      do 1 i=1,n
        theta=theta+thinc
        x(i)=6.0*cos(theta)
        y(i)=6.0*sin(theta)
   1  continue
      call xyscale (x,y,n,0.0,9.0,0.0,9.0)
c ..
c .. join point i to point j for all i.lt.j and i.lt.n
      nm1=n-1
      do 3 i=1,nm1
        ip1=i+1
        do 2 j=ip1,n
          call plot(x(i),y(i),3)
          call plot(x(j),y(j),2)
   2    continue
   3  continue
      call efplot
      return
      end
      subroutine output (ifile)
      write(6,890)
  890 format(//5x,'Output unit selection:')
c     2 /8x,'0 = no graphics initialization, output to screen',
c     3 /8x,'5 = input unit, not definable as output unit',
c     4 /8x,'6 = graphics output to screen',
c     5 /8x,'All other numbers will open an output file called PLOT.OUT',
c     6 //5x,'Specify unit number:',/8x)
      read(5,*) ifile
      if(ifile .eq. 5) ifile=6
      return
      end
      subroutine init
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      write(iunit,*) ' ;: A'
      return
      end
      subroutine error(istring,ilen)
      call symbol(.1,.30,.0625,'*** ERROR ***',0.0,13)
      call symbol(.1,0.1,.0625,istring,0.0,ilen)
      return
      end
      subroutine screen
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      write(iunit,900) 
  900 format(' H')
      return
      end
      subroutine plots (x,y,iplot)
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      data fact,xorigin,yorigin,iunit/1.0,0.0,0.0,6/
      if (iplot .eq. 999) goto 999
      if (iplot .eq. 0) goto 200
      if (iplot .ne. 6) open (unit=iplot,file='plot.out',status='new')
      iunit=iplot
c ..
c .. initiate graphics
      call init
  200 call screen
      return
  999 call efplot
      return
      end
      subroutine plot (x1,y1,iplot)
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      x=x1
      y=y1
c ..
c .. relate x and y to current origin
      x=(x*fact)+xorigin
      y=(y*fact)+yorigin
c ..
c .. save the current x and y value
      ysave=y
      xsave=x
c ..
c .. convert to plot coordinates
      ix=nint(x*200.)
      iy=nint(y*200.)
c
c ..
c .. check if ix and iy within visible window
      if (iy.lt.0.or.iy.gt.2200.or.ix.lt.0.or.ix.gt.3400)
     2  call error ('plot: x or y exceeds screen coordinates',39)      
c ..
c .. switch-off based on iplot parameter
      if (iplot .eq. 2 .or. iplot .eq. 3) goto 200
      if (iplot .eq.-2 .or. iplot .eq.-3) goto 500
      call error('plot: iplot does not have an acceptable value',45)
      return
  200 if (iplot .eq. 3) goto 300
c ..
c .. plot line to x,y position
      write(iunit,900) ix,iy
  900 format(' D ',i4,',',i4)
      return
c ..
c .. move to x,y position
  300 write(iunit,910) ix,iy
  910 format(' U ',i4,',',i4)
      return
c ..
c .. reset origin with a plotted line
  500 if (iplot .eq. -3) goto 600
      write(iunit,900) ix,iy
      xorigin=x
      yorigin=y
      return
c ..
c .. reset origin with no line plotted
  600 write(iunit,910) ix,iy
      xorigin=x
      yorigin=y
      return
      end
      subroutine symbol (x,y,hgt,istr,ang,nchar)
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      dimension istr(20)
c ..
c .. check angle
      if (ang .ge. 0.0 .and. ang .le. 360.0) goto 50
      call error ('symbol: angle not within 0 to +360 range',40)
      return
c ..
c .. classify height into one of the 4 allowable hiplot heights
   50 call clashgt (hgt, iht, is, ang, iang)
c ..
c .. move to position for lettering
      call plot (x,y,3)
c ..
c .. compute length of character string
      write(iunit,900) iang,iht,(istr(kk),kk=1,ilen),char(94)
  900 format(' S',i2,i2,20a4,a1)
      return
      end
      subroutine number (x,y,hgt,fpn,ang,ndec)
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      character*12 istring
c ..
c .. convert floating point number to string
      encode (12,900,istring) fpn
  900 format(f12.<ndec>)
c ..
c .. check angle
      if (ang .ge. 0.0 .and. ang .le. 360.0) go to 50
      call error ('number: angle not within 0 to +360 range',40)
      return
c ..
c .. classify height into one of the 17 allowable gigi heights
  50  call clashgt (hgt, iht, is, ang, iang)
c ..
c .. move to position for symbol plotting
      call plot (x,y,3)
c ..
c .. write number using gigi format
      write(iunit,910) iang,iht,istring,char(94)
  910 format(' S',i2,i2,a12,a1)
      return
      end
      subroutine clashgt (hgt, iht, is, ang, iang)
      dimension heights(5)
      data heights /.07,.14,.28,.56,1.12/
      if (hgt .lt. heights(1)) iht=0
      if (hgt .ge. heights(5)) iht=5
      if (hgt .lt. heights(1) .or. hgt .gt. heights(5)) goto 400
      do 100 i=1,4
        hmid=(heights(i)+heights(i+1))/2.0
        if (hgt .lt. heights(i+1) .and. hgt .ge. hmid) goto300
        if (hgt .lt. hmid .and. hgt .ge. heights(i)) goto200
  100 continue
  300 iht = i + 1
      goto 400
  200 iht = i
c ..
c .. compute angle as a number between 1 and 4
  400 if (ang. gt. 315.0 .and. ang .le. 45.0) iang=1
      if (ang .gt. 45.0 .and. ang .le. 135.0) iang=2
      if (ang .gt. 135.0 .and. ang .le. 225.0) iang=3
      if (ang .gt. 225 .and. ang .le. 315.0) iang=4
      return
      end
      subroutine efplot
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      write(iunit,*) '@@@@'
      return
      end
      subroutine newpen (icolor)
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      if (icolor .lt. 0 .or. icolor .gt. 7) icolor=7
      write(iunit,900) icolor
  900 format(' !w(i',i1,')')
      return
      end
      subroutine maxmin (v,n,vmin,vmax)
      dimension v(n)
      vmax=-1000000.0
      vmin=1.0e15
      do 100 i=1,n
        vmax = amax1 (v(i), vmax)
        vmin = amin1 (v(i), vmin)
  100 continue
      return
      end
      subroutine xyscale (x,y,n,xlow,xhigh,ylow,yhigh)
      dimension x(n), y(n)
      call maxmin (x,n,xmin,xmax)
      call maxmin (y,n,ymin,ymax)
c ..
c .. determine delta-x and delta-y in old and new coordinate system
      dxold = xmax - xmin
      dyold = ymax - ymin
c     write(6,*) 'dxold, xmax,xmin', dxold,xmax,xmin
      dxnew = xhigh - xlow
      dynew = yhigh - ylow
c     write(6,*) 'dxnew,dynew,xhigh,xlow,yhigh,ylow'
c     write(6,*) dxnew,dynew,xhigh,xlow,yhigh,ylow
c ..
c .. check for division by zero
      if (dxold .eq. 0.0 .or. dyold .eq. 0.0) return
      factx = dxnew / dxold
      facty = dynew / dyold
c     write(6,*) 'factx,facty,dxnew,dxold,dynew,dyold'
c     write(6,*) factx,facty,dxnew,dxold,dynew,dyold
c ..
c .. scale to new coordinate system
      do 100 i=1,n
        x(i) = abs (x(i) - xmin) * factx + xlow
        y(i) = abs (y(i) - ymin) * facty + ylow
c       write(6,*) x(i),y(i)
  100 continue
      return
      end
      subroutine xyline (x,y,n)
      dimension x(n), y(n)
      call plot (x(1),y(1),3)
      do 100 i = 2,n
        call plot (x(i),y(i),2)
  100 continue
      return
      end
      subroutine rotate (x,y,n,ang)
c **************************************************************
c
c .. ang = counterclockwise angle that brings the original
c          points into the desired rotation.
c
c **************************************************************
      dimension x(n), y(n)
      data radconv /57.29578/
c ..
c .. if ang out of range, no rotation is performed
      if (ang .le. 0.0 .or. ang .ge. 360) return
c ..
c .. convert ang to a radian measure
      radang = ang / radconv
c ..
c .. rotate coordinates by ang
      u = cos (radang)
      v = sin (radang)
      do 100 i = 1, n
        tempx = x(i) * u + y(i) * v
        y(i) = -x(i) * v + y(i) * u
        x(i) = tempx
  100 continue
      return
      end
      subroutine dash (idash)
      common /plotvar/fact,xorigin,yorigin,iunit,xsave,ysave
      if (idash .lt. 0 .or. idash .gt. 9) return
      write(iunit,900) idash
  900 format (' L',i2)
      return
      end
      subroutine xscale (x,y,n,xlow,xhigh)
      dimension x(n), y(n)
      call maxmin (x,n,xmin,xmax)
      call maxmin (y,n,ymin,ymax)
c ..
c .. determine delta-x in old and new coordinate system
      dxold = xmax - xmin
      dxnew = xhigh - xlow
c ..
c .. check for division by zero
      if (dxold .eq. 0.0 ) return
      factx = dxnew / dxold
c ..
c .. for this application, facty = factx
      facty = factx
c ..
c .. scale to new coordinate system
      do 100 i=1,n
        x(i) = abs (x(i) - xmin) * factx + xlow
        y(i) = abs (y(i) - ymin) * facty + xlow
  100 continue
      return
      end
      subroutine yscale (x,y,n,ylow,yhigh)
      dimension x(n), y(n)
      call maxmin (x,n,xmin,xmax)
      call maxmin (y,n,ymin,ymax)
c ..
c .. determine delta-y in old and new coordinate system
      dyold = ymax - ymin
      dynew = yhigh - ylow
c ..
c .. check for division by zero
      if (dyold .eq. 0.0 ) return
      facty = dynew / dyold
c ..
c .. for yscale, factx=facty
      factx=facty
c ..
c .. scale to new coordinate system
      do 100 i=1,n
        x(i) = abs (x(i) - xmin) * factx + ylow
        y(i) = abs (y(i) - ymin) * facty + ylow
  100 continue
      return
      end