c Copyright (c) 1990 by James A. D. Connolly, Institute for Mineralogy and
c Petrography, Swiss Federal Insitute of Technology, CH-8092 Zurich,
c SWITZERLAND. All rights reserved.
 
c PSVDRAW - a program to generate postscript code for x-y plots and chemographic
c phase diagrams.  The input data format is consistent with the files generated
c by the programs VERTEX and FRENDLY.
 
c Please do not distribute any part of this source.
 
 
      PROGRAM PVDRAW
 
      integer h9
 
      parameter (k1=1000,h9=18)
 
      character*80 fname, yes*1, names*8, sname*10
 
      common / cst8 /names(k1)/ csta7 /sname(h9)
      common / xyrat /xfac
 
      n4 = 24
 
      write (*,1020) 
      read (*,1000) fname
 
      open (n4,err=90,file=fname,status='old')
      icopt = 1
 
      call psopen (fname,ier)
      goto (999), ier
  
      write (*,1030) 
      read (*,1010) yes
      iop0 = 0
      if (yes.eq.'y'.or.yes.eq.'Y') iop0 = 1

      xfac = 1.0
      call psxypl (n4,iop0,icopt)
  
      goto 99
90    write (*,1040) fname
      goto 999
 
99    call psclos
 
      close (99)
 
1000  format (a80)
1010  format (a1)
1020  format (/,' Enter plot file name: ',/)
1030  format (/,' Modify the default plot (y/n)? ',/)
1040  format (/,' no such file as: ',a14,/)
1050  format (/,a14,' is not formatted correctly.',/)
 
999   end
c---------------------------------------------------------------------
      subroutine psxypl (n4,iop0,icopt)
 
c psxypl - subroutine to output x-y plot.
 
      character*72 fname(4),vname*30, names*8, sname*10
 
      integer h9, k1
 
      parameter (k1=1000,h9=18,k9=2000)
 
      common / cst8 /names(k1)/ csta7 /sname(h9)/ phase /ikp(k1)
      common/ vs    /vmin(5),vmax(5),iv/ vsa /vname(5)
      common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3)
      common / xypts /xpt(k9),ypt(k9),isym(k9),ipts
 
      iv = 2
      ipts = 1
      vmin(1) = 1.d99
      vmax(1) = -1.d99
      vmin(2) = 1.d99
      vmax(2) = -1.d99
      vname(1) = 'x axis'
      vname(2) = 'y axis'
      fname(1) = 'title'

10    read (n4,*,end=30) isym(ipts),xpt(ipts),ypt(ipts)
      if (xpt(ipts).lt.vmin(1)) vmin(1) = xpt(ipts)
      if (xpt(ipts).gt.vmax(1)) vmax(1) = xpt(ipts)
      if (ypt(ipts).gt.vmax(2)) vmax(2) = ypt(ipts)
      if (ypt(ipts).lt.vmin(2)) vmin(2) = ypt(ipts)
      ipts = ipts + 1
      goto 10
c                                 get some options and
c                                 set up transformations
30    call psaxop (iop0,iop1,jop0,jwidth,icopt)

      ipts = ipts - 1

      call psvdrw (n4,fname,iop0,iop1,jwidth,jop0)
 
      call psaxes (jwidth,jop0)
 
1000  format (a72)
1010  format (10a8)
1020  format (2x,a30)
1030  format (8a10)

      end
c----------------------------------------------------------------------
      subroutine psaxop (iop0,iop1,jop0,jwidth,icopt)

c psaxop - subroutine to make graphics transformation and get some options

      character yes*1, vname*30

      common / wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen
      common/ vs    /vmin(5),vmax(5),iv/ vsa /vname(5)
      common / xyrat /xfac


      dsx = 0.18
      dsy = 0.18
      dtx = 160.
      dty = 220.
      drot = 0.
      jop0 = 0
      jwidth = 1 

      if (icopt.ne.3.and.iop0.eq.1) then
c                                     modify default drafting
c                                     options?          
         write (*,1090) 
         read (*,1030) yes
         if (yes.eq.'y'.or.yes.eq.'Y') jop0 = 1

       else if (icopt.eq.3) then

          jop0 = iop0

       end if 

       if (jop0.eq.1) then
 
         if (icopt.eq.3) goto 10 

         write (*,1070) 
         read (*,1030) yes
         iop1 = 0
 
         if (yes.eq.'y'.or.yes.eq.'Y') then
            write (*,1040) vname(1),vmin(1),vmax(1)
            read (*,*) vmin(1),vmax(1)
            write (*,1040) vname(2),vmin(2),vmax(2)
            read (*,*) vmin(2),vmax(2)
            iop1 = 1
            write (*,1080) 
         end if

10       write (*,1100)
         read (*,1030) yes

         if (yes.eq.'y'.or.yes.eq.'Y') then

            write (*,1110) dsx,dsy,dtx,dty,drot
            read (*,*) dsx,dsy,dtx,dty,drot

         end if 

         write (*,1130)
         read (*,1030) yes

         if (yes.eq.'y'.or.yes.eq.'Y') then

            write (*,1140) 
            read (*,*) xfac

         end if 

         if (icopt.eq.3) goto 20 

         write (*,1120) 
         read (*,1030) yes
         if (yes.eq.'y'.or.yes.eq.'Y') jwidth = 0
 
      end if

20    xmax = vmax(1)
      xmin = vmin(1)
      ymax = vmax(2)
      ymin = vmin(2)
      xlen = xmax - xmin
      ylen = ymax - ymin
c                                     set default char size, based
c                                     on font 7 (helv 14) the axes
c                                     length is about 85 characters
      dcx = xlen / 85.
      dcy = ylen / 85.
 
      call psssc2 (xmin,xmax,ymin,ymax)

c      call pssscm (xfac, 1.0)

      dsx = dsx * xfac
 
      call psstrn (dsx,dsy,dtx,dty,drot)

1030  format (a1)
1040  format (/,' Enter new min and max for ',a8,' old values ',
     *        ' were: ',2(g11.5,1x),/)
1070  format (/,' Modify x-y limits (y/n)? ',/)
1080  format (/,' OK, but this may look sloppy. ',/)
1090  format (/,' Modify default drafting options (y/n)?',/)
1100  format (/,' Modify default picture transformation (y/n)?',/)
1110  format (/,' The default values for x-y scaling, x-y',
     *          ' translation, and rotation',/,' are, respectively:',
     *        /,(5(4x,g8.3)),/,' Enter the new values:',/)
1120  format (/,' Use minimum width lines (y/n)?',/)
1130  format (/,' Modify ratio of x to y axis length (y/n)',/)
1140  format (/,' Enter fractional length of x axis relative',
     *          ' to the y axis:',/)

      end 

c----------------------------------------------------------------------
      subroutine psaxes (iwidth,jop0)
 
c psaxes - subroutine to output (sloppy) axes.
 
      character*30 vname, record, yes*1
 
      common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen
      common/ vs    /vmin(5),vmax(5),iv/ vsa /vname(5)
      common / xyrat /xfac
 
      ihalf = 1
      itenf = 0
      igrid = 0
      rline = 1.

      csn = 1.
      csa = 1.

      x0 = xmin
      y0 = ymin
      dx = xlen / 5.
      dy = ylen / 5.

      xtic  = xlen/45.
      xtic1 = xtic*.67
      xtic2 = xtic1*.67
 
      ytic  = ylen/45.
      ytic1 = ytic*.67
      ytic2 = ytic1*.67 

      if (jop0.eq.1) then
         write (*,1010) 
         read (*,1020) yes
         if (yes.ne.'y'.and.yes.ne.'Y') goto 10

         write (*,1030) 'X', x0, dx
         read (*,*) x0, dx
         write (*,1030) 'Y', y0, dy
         read (*,*) y0, dy

         write (*,1040)
         read (*,1020) yes
         if (yes.eq.'y'.or.yes.eq.'Y') ihalf = 0

         write (*,1050) 
         read (*,1020) yes
         if (yes.eq.'y'.or.yes.eq.'Y') itenf = 1

         write (*,1070) 
         read (*,1020) yes
         if (yes.eq.'y'.or.yes.eq.'Y') then 
            write (*,1060) 
            read (*,*) csn, csa
         end if

         write (*,1080) 
         read (*,1020) yes 
         if (yes.eq.'y'.or.yes.eq.'Y') igrid = 1
     
      end if 
c                                       draw axes box
10    call psrect (xmin,xmax,ymin,ymax,1.,iwidth,0)
c                                       draw left vertical axis tics
      call psytic (xmin, y0, dy, xtic, xtic1, xtic2, rline, 
     *             ihalf, itenf, iwidth)
c                                       draw right vertical axis
      call psytic (xmax, y0, dy, -xtic, -xtic1, -xtic2, rline, 
     *             ihalf, itenf, iwidth)
c                                       draw bottom horizontal axis
      call psxtic (ymin, x0, dx, ytic, ytic1, ytic2, rline,
     *             ihalf, itenf, iwidth)
c                                       draw top horizontal axis
      call psxtic (ymax, x0, dx, -ytic, -ytic1, -ytic2, rline, 
     *             ihalf, itenf, iwidth)
 
      call pssctr (7,csn*1.2,csn*1.2,0.)

      cx = dcx * csn
      cy = dcy * csn 
c                                       numeric axis labels:
      call psylbl (xmin, xmax, y0, ymax, dy, cx, cy, tmin, igrid) 
 
      call psxlbl (x0, dx, cx, 2.3*cy, igrid)

      cx = dcx * csa
      cy = dcy * csa
c                                       x-axis name
      call pssctr (7,csa*2.,csa*2.,0.0)
      t = xmin + 0.5 * xlen - (7. * cx) / xfac
      call pstext (t, ymin - 0.07 * ylen, vname(1),0)
c                                       y-axis name
      call pssctr (7,csa*2.,csa*2., 90.)
      t = tmin - 5. * cx / xfac 
      call pstext (t, ymin + 0.43 * ylen,vname(2),0)
c                                       sectioning constraints
      if (iv.gt.2) then
 
         call pssctr (8,csa*1.2,csa*1.2,0.)
         y = ymax + 0.12 * ylen
         do 1 i = 3, iv
            write (record,1000) vname(i),vmin(i)
            call pstext (xmin,y,record,0)
1           y = y - 0.037 * ylen
 
      end if
 
1000  format (a8,'=',g9.3)
1010  format (/,' Modify default axes (y/n)?',/)
1020  format (a1)
1030  format (/,' Enter the starting value and interval for',
     *          ' major tick marks on',/,' the ',a1,'-axis (the'
     *          ' current values are:',2(1x,g9.3),')',/, 
     *          ' Enter the new values:',/)
1040  format (/,' Cancel half interval tick marks (y/n)?',/)
1050  format (/,' Draw tenth interval tick marks (y/n)?',/)
1060  format (/,' Enter scale factors for numeric and text',
     *          ' axes labels:') 
1070  format (/,' Rescale axes text (y/n)?',/)
1080  format (/,' Turn gridding on (y/n)?',/)
      end
c------------------------------------------------------------------
      subroutine psytic (x0, y0, dy, tic, tic1, tic2, rline, ihalf,
     *                   itenf, iwidth)

      common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen

c psytic - subroutine to draw y-axis ticks.

      y = y0 
      call psmove (x0,y0)

      if (ihalf.eq.0.and.itenf.eq.0) then 

10       if (y.ge.ymax) goto 99
         call psrlin (tic,0.,rline,iwidth)
         call psrmov (-tic,dy)
         y = y + dy 
         goto 10

      else if (itenf.eq.0) then

         dy2 = dy / 2.

20       if (y.ge.ymax) goto 30
         call psrlin (tic,0.,rline,iwidth)
         call psrmov (-tic,dy2)
         y = y + dy2
         if (y.ge.ymax) goto 30
         call psrlin (tic1, 0., rline, iwidth)
         call psrmov (-tic1,dy2)
         y = y + dy2
         goto 20

30       if (y0-dy2.gt.ymin)
     *      call psline (x0,y0-dy2,x0+tic1,y0-dy2,rline,iwidth)
         goto 99

      else

         dy2 = dy / 10.

40       if (y.ge.ymax) goto 50
         call psrlin (tic, 0., rline, iwidth)
         call psrmov (-tic, dy2)
         y = y + dy2

         do 60 i = 1, 4
            if (y.ge.ymax) goto 50
            call psrlin (tic2, 0., rline, iwidth)
            call psrmov (-tic2, dy2)
60          y = y + dy2

         if (y.ge.ymax) goto 50
         call psrlin (tic1, 0., rline, iwidth)
         call psrmov (-tic1, dy2)
         y = y + dy2

         do 70 i = 1, 4
            if (y.ge.ymax) goto 50
            call psrlin (tic2, 0., rline, iwidth)
            call psrmov (-tic2, dy2)
70          y = y + dy2

         goto 40

50       if (y0-dy2.lt.ymin) goto 99

         y = y0 - dy2
         call psmove (x0,y)

         do 80 i = 1, 4
            if (y.le.ymin) goto 99
            call psrlin (tic2, 0., rline, iwidth)
            call psrmov (-tic2, -dy2)
80          y = y - dy2

         if (y.le.ymin) goto 99
         call psrlin (tic1, 0., rline, iwidth)
         call psrmov (-tic1, -dy2)
         y = y - dy2 

         do 90 i = 1, 4
            if (y.le.ymin) goto 99
            call psrlin (tic2, 0., rline, iwidth)
            call psrmov (-tic2, -dy2)
90          y = y - dy2

      end if 

99    end

c------------------------------------------------------------------
      subroutine psxtic (y0, x0, dx, tic, tic1, tic2, rline, ihalf,
     *                   itenf, iwidth)

      common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen

c psxtic - subroutine to draw x-axis ticks.

      x = x0 
      call psmove (x0,y0)

      if (ihalf.eq.0.and.itenf.eq.0) then 

10       if (x.ge.xmax) goto 99
         call psrlin (0., tic, rline, iwidth)
         call psrmov (dx, -tic)
         x = x + dx 
         goto 10

      else if (itenf.eq.0) then

         dx2 = dx / 2.

20       if (x.ge.xmax) goto 30
         call psrlin (0., tic, rline, iwidth)
         call psrmov (dx2, -tic)
         x = x + dx2
         if (x.ge.xmax) goto 30
         call psrlin (0., tic1, rline, iwidth)
         call psrmov (dx2, -tic1)
         x = x + dx2
         goto 20

30       if (x0-dx2.gt.xmin)
     *      call psline (x0-dx2, y0, x0+tic1, y0, rline, iwidth)
         goto 99

      else

         dx2 = dx / 10.

40       if (x.ge.xmax) goto 50
         call psrlin (0., tic, rline, iwidth)
         call psrmov (dx2, -tic)
         x = x + dx2

         do 60 i = 1, 4
            if (x.ge.xmax) goto 50
            call psrlin (0., tic2, rline, iwidth)
            call psrmov (dx2, -tic2)
60          x = x + dx2

         if (x.ge.xmax) goto 50
         call psrlin (0., tic1, rline, iwidth)
         call psrmov (dx2, -tic1)
         x = x + dx2

         do 70 i = 1, 4
            if (x.ge.xmax) goto 50
            call psrlin (0., tic2, rline, iwidth)
            call psrmov (dx2, -tic2)
70          x = x + dx2

         goto 40

50       if (x0-dx2.lt.xmin) goto 99

         x = x0 - dx2

         call psmove (x, y0)

         do 80 i = 1, 4
            if (x.le.xmin) goto 99
            call psrlin (0., tic2, rline, iwidth)
            call psrmov (-dx2, -tic2)
80          x = x - dx2

         if (x.le.xmin) goto 99
         call psrlin (0., tic1, rline, iwidth)
         call psrmov (dx2, -tic1) 
         x = x - dx2 

         do 90 i = 1, 4
            if (x.le.xmin) goto 99
            call psrlin (0., tic2, rline, iwidth)
            call psrmov (-dx2, -tic2)
90          x = x - dx2

      end if 

99    end
 
      subroutine psylbl (xmin, xmax, y0, ymax, dy, 
     *                   dcx, dcy, tmin, igrid)
 
c psylbl - subroutine to put on y-axis labels.
 
      character*12 numbs(40)
 
      integer ic(40)

      common / xyrat /xfac

      xdc = 1.4 * dcx
      ydc = .8 * dcy
      tmin = 1.d30
 
      call psnum (y0,ymax,dy,ic,iy1,numbs)
 
      y = y0
 
      do 10 i = 1, iy1
 
         t = xmin - (float(ic(i)+1) * xdc) / xfac
         if (t.lt.tmin) tmin = t
         call pstext (t, y + ydc, numbs(i), ic(i))
         if (igrid.eq.1) call psline (xmin, y, xmax, y, 10., 0)
 
10    y = y + dy
 
      end
 
      subroutine psnum (rmin,rmax,dr,ic,i1,numbs)
 
      character*1 text(12), numbs(40)*12, next(12)
 
      integer ic(40), k

      i1 = 1 + ifix ((rmax - rmin) / dr)
 
      if (rmax.gt. 9999.9 .and.rmax.lt. 99999.9 .and.rmin.gt.-rmax) then
         icase = 1
      else if (rmax.gt. 999.9 .and.rmax.le. 9999.9
     *                        .and.rmin.gt.-rmax) then
         icase = 2
      else if (rmax.gt. 99.9 .and.rmax.le. 999.9
     *                       .and.rmin.gt.-rmax) then
         icase = 3
      else
         icase = 4
      end if
 
      r = rmin
 
      do 10 i = 1, i1
         goto (1,2,3,4) icase
1        write (numbs(i),1001) ifix(r)
         goto 20
2        write (numbs(i),1002) ifix(r)
         goto 20
3        write (numbs(i),1003) ifix(r)
         goto 20
4        write (numbs(i),1004) r
20       read (numbs(i),1010) text
         k = 0
         do 30 j = 1, 12
            if (text(j).eq.' ') goto 30
            k = k + 1
            next(k) = text(j)
30       continue
         ic(i) = k
         write (numbs(i),1010) (next(j),j=1,k)
10    r = r + dr
 
1001  format (i5)
1002  format (i4)
1003  format (i3)
1004  format (g9.3)
1010  format (12a1)
 
      end
 
      subroutine psxlbl (x0, dx, dc, dyl, igrid)
 
c psxlbl - subroutine to put on x-axis labels.
 
      character*12 numbs(40)

      common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen
 
      integer ic(40)

      common / xyrat /xfac
c
      y = ymin - dyl
 
      x = x0
 
      xdc = .75 * dc
 
      call psnum (x0,xmax,dx,ic,ix1,numbs)
 
      do 10 i = 1, ix1
         if (x.eq.xmin.or.x.eq.xmax) goto 10
         t = x - (float(ic(i)) * xdc) / xfac
         call pstext (t, y, numbs(i), ic(i))
         if (igrid.eq.1) call psline (x, ymin, x, ymax, 10., 0)
10    x = x + dx
 
      end
 
      subroutine psvdrw (n4,pname,iop0,iop1,jwidth,jop0)
C----------------------------------------------------------------
 
c psvdrw - subroutine to output schreinemakers diagrams.
 
      character*72 pname(4), yes*1, prompt*14

      common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen
      common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3)
c---------------------------------------------------------------
 
      iop2 = 0
      iop3 = 0
      iop5 = 0
      iop6 = 0
      iop7 = 0
      iop8 = 1
 
      fac1 = .6
      fac2 = .2
      fac3 = .05
  
      call pssctr (8,1.2,1.2,0.)
 
      call psipts 
 
      call pssctr (8,1.6,1.6,0.)

      x = xmin - .15 * xlen

      call pstext (x, ymax+.32*ylen, pname(1), 0)
      call pssctr (8,1.,1.,0.)

      y = ymax + 0.27 * ylen

      do 30 i = 2, 4
         call pstext (x,y,pname(i),0)
30       y = y - 0.03 * ylen
 
 
1000  format (a1)
1030  format (/,' Enter minimum fraction of the axes length that a',/
     *        ' curve must be to receive a text label (0-1): ',/)
1070  format (/,' Change default labelling of curve segments (y/n)?',/)
1080  format (/,' Reset default size for text labels (y/n)? ',/)
1090  format (/,' Scale default text size by a factor of: ',/)
1110  format (/,' Enter minimum fraction of the axes length that a',/
     *        ' curve must be to receive a numeric label (0-',f5.3,
     *        '): ',/)
1120  format (/,' Show only without phases (y/n)? ',/)
1130  format (/,' Show only with phases (y/n)? ',/)
1140  format (/,' Show only with assemblage (y/n)? ',/)
1150  format (/,' Restrict phase fields (y/n)? ',/)
1160  format (/,' Suppress high variance equilibria (y/n)? ',/)
1170  format (/,' Suppress pseudounivariant curves (y/n)? ',/)
1180  format (/,' Suppress pseudoinvariant points (y/n)? ',/)
1190  format (/,' Suppress point labels (y/n)? ',/)
1200  format (/,1x,g13.6,' is an invalid value.',/)
1210  format (/,' Fit curves with B-splines (y/n)? ',/)
1220  format (/,' Suppress curve labels (y/n)? ',/)
1230  format (/,' Modify default equilibrium labelling (y/n)?',/)
 
      end 




c---------------------------------------------------------------
      subroutine psipts 
 
c psipts - subprogram to output invariant points.

      parameter (k9=2000)
 
      character*10 xnams
 
      integer kex(9), k

      common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3)
      common/ excl4 /xnams(50,3)
      common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen
      common / xypts /xpt(k9),ypt(k9),isym(k9),ipts
  
      do 20 i = 1, ipts

         x = xpt(i)
         y = ypt(i)
         r = 0.78

 
         if (isym(i).lt.4) then
            if (isym(i).eq.0) then

               call pselip (x,y,r*dcx,r*dcy,1.,0,7)

            else if (isym(i).eq.1) then
 
               call psrect (x-r*dcx,x+r*dcx,y-r*dcy,y+r*dcy,1.,0,7)

            else if (isym(i).eq.2) then

               call pselip (x,y,r*dcx,r*dcy,1.,0,1)

            else if (isym(i).eq.3) then

               call psrect (x-r*dcx,x+r*dcx,y-r*dcy,y+r*dcy,1.,0,1)

            end if 

         else 
 
            r = 0.38

            if (isym(i).eq.4) then

               call pselip (x,y,r*dcx,r*dcy,1.,0,7)

            else if (isym(i).eq.5) then
 
               call psrect (x-r*dcx,x+r*dcx,y-r*dcy,y+r*dcy,1.,0,7)

            else if (isym(i).eq.6) then

               call pselip (x,y,r*dcx,r*dcy,1.,0,1)

            else if (isym(i).eq.7) then

               call psrect (x-r*dcx,x+r*dcx,y-r*dcy,y+r*dcy,1.,0,1)

            end if 

         end if 

 20    continue
 
99    end