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 implicit real (a-g,o-z),integer (h-n) parameter (k1=5000,h9=30) character*80 fname, yes*1, names*8, sname*10, chars(80)*1, xname common / cst8 /names(k1)/ csta7 /sname(h9) common / basic /n4,ifont,icopt,jwidth,iop0 common / xyrat /xfac n4 = 24 c default font ifont = 7 c default line width jwidth = 1 write (*,1020) read (*,1000) fname xname = fname open (n4,err=90,file=fname,status='old') read (n4,*,err=91) icopt 20 call psopen (fname,ier) goto (999), ier c if icopt = 5, open c polygon file: ipoly = 0 if (icopt.eq.5) then read (xname,1060) chars write (xname,1070) (chars(i),i=2,80) open (n4+1,err=92,file=xname,status='old') ipoly = 1 end if 30 if (icopt.ne.0) then write (*,1030) read (*,1010) yes iop0 = 0 if (yes.eq.'y'.or.yes.eq.'Y') iop0 = 1 end if if (icopt.eq.1.or.icopt.ge.5) then xfac = 1.0 call psxypl (ipoly) else if (icopt.eq.0) then call pschem else if (icopt.eq.3) then xfac = 0.5 call psmixd else goto 91 end if goto 99 90 write (*,1040) fname goto 999 91 write (*,1050) fname goto 999 92 write (*,1080) xname read (*,1010) yes if (yes.eq.'y'.or.yes.eq.'Y') goto 30 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.',/) 1060 format (80a1) 1070 format ('b',79a1) 1080 format (/,' No file named ',a14,' run vertex to create it.',/, * ' Continue without it (y/n)? ',/) 999 end c--------------------------------------------------------------------- subroutine psxypl (ipoly) c psxypl - subroutine to output x-y plot. implicit real (a-g,o-z),integer (h-n) character*162 fname(4),vname*30, names*8, sname*10 parameter (k1=5000,k5=12,h9=30) common/ cst8 /names(k1)/ csta7 /sname(h9)/ phase /ikp(k1) * / vs /vmin(5),vmax(5),iv/ vsa /vname(5) * / excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ bulk /ipid,ipvar,amode(k5),vmode(k5),wtp(k5),rho, * ipds(k5),ivol,incts,pmn,dp common/ basic /n4,ifont,icopt,jwidth,iop0 n5 = n4 + 1 if (icopt.eq.6) read (n4,*) incts, ivol read (n4,*) iphct, isoct if (iphct.ne.0) then read (n4,1010) (names(i),i=1,iphct) read (n4,*) (ikp(i),i=1,iphct) end if if (isoct.ne.0) read (n4,1030) (sname(i),i=1,isoct) read (n4,1000) fname read (n4,*) iv read (n4,*) (vmax(i),vmin(i),i=1,iv) read (n4,1020) (vname(i),i=1,iv) if (ipoly.eq.1) then read (n5,*) icopt read (n5,*) incts, ivol read (n5,*) iphct, isoct if (iphct.ne.0) then read (n5,1010) (names(i),i=1,iphct) read (n5,*) (ikp(i),i=1,iphct) end if if (isoct.ne.0) read (n5,1030) (sname(i),i=1,isoct) read (n5,1000) fname read (n5,*) iv read (n5,*) (vmax(i),vmin(i),i=1,iv) read (n5,1020) (vname(i),i=1,iv) end if c get some options and c set up transformations call psaxop (iop1,jop0) if (icopt.eq.1) then call psvdrw (fname,iop1,jop0) else if (icopt.eq.5) then call psbplt (fname,ipoly,iop1,jop0) else if (icopt.eq.6) then call pscplt (fname,iop1,jop0) end if call psaxes (jop0) 1000 format (a162) 1010 format (10a8) 1020 format (2x,a30) 1030 format (8a10) end c---------------------------------------------------------------------- subroutine psaxop (iop1,jop0) c psaxop - subroutine to make graphics transformation and get some options implicit real (a-g,o-z),integer (h-n) 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/ basic /n4,ifont,icopt,jwidth,iop0 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) 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 pschem c pschem - subroutine to output ternary chemographies. implicit real (a-g,o-z),integer (h-n) parameter (k1=5000,k2=40000,k5=12,h9=30) character names*8, vname*30, title*162, record*72, yes*1, * sname*10, xname(k5)*12 real x3(3), y3(3), xx(1000), yy(1000) integer idss(10), iperm(2,3) common/ vs /vmin(5),vmax(5),iv/ vsa /vname(5) common/ basic /n4,ifont,icopt,jwidth,iop0 common/ cst8 /names(k1)/ csta7 /sname(h9)/ phase /ikp(k1) common/ asmbl /idf(3,k2),ib,iasmbl(k2),ivchk(k1),x(2,k1) save x3,y3,dyt,yt,xt,iperm data x3,y3,dyt,yt,xt,iperm/0.0,0.5,1.0,0.0,0.866025,0.0, * 300.,980.,100.,1,2,1,3,2,3/ c ----------------------------------------------- c read header information read (n4,*) icp,istct,iphct,ipoint,ifyn,isyn,ipot,isoct if (icp.ne.3) then write (*,*) ' Sorry, PSVDRAW only works', * ' for ternary composition diagrams.' stop end if if (iphct.gt.k1.or.iv.gt.5 * .or.ib.gt.k2.or.ipot.gt.5) then write (*,*) ' dimensioning error in pschem' goto 99 end if read (n4,1000) (vname(i),i=1,ipot) read (n4,1010) title icp1 = icp - 1 c names read (n4,1020) (names(i), i = 1, iphct) c phase coordinates read (n4,*) ((x(j,i),j = 1, 2), i=istct, iphct) read (n4,*) (ikp(i), i = 1, iphct) c solution names if (isoct.ne.0) read (n4,1090) (sname(i),i = 1, isoct) c write graphics code names for x variables: read (n4,1080) (xname(i), i = 1, icp) c ----------------------------------------------- c convert to equilateral coordinates: do i = istct, iphct x(1,i) = x(1,i) + 0.5 * x(2,i) x(2,i) = x(2,i) * 0.866025 end do c draw tie lines? write (*,1060) read (*,1070) yes iop1 = 0 c iop1 = 1 = yes if (yes.eq.'y'.or.yes.eq.'Y') iop1 = 1 iflag = 0 c read simple variables 10 read (n4,*,end=99) (vmin(i), i = 1, ipot) read (n4,*) ib c read phase configurations read (n4,*) ((idf(j,i),j = 1, icp), i = 1, ib) do i = 1, iphct ivchk(i) = 0 end do c read assemblage flags read (n4,*) (iasmbl(i), i = 1, ib) c read saturated phase id's if (isyn.ne.1) then read (n4,*) isat read (n4,*) (idss(i), i = 1, isat) end if iflag = (iflag-1)**2 if (iflag.eq.0) then xt = 375. else yt = yt - dyt xt = 75. end if call psssc2 (0.,1.,0.,1.) call psstrn (0.08,0.08,xt,yt,0.) if (iop1.eq.0) then c don't draw tielines c first divariant do i = 1, ib if (iasmbl(i).eq.3) then call merger (i,iperm,ivert,ipoint,xx,yy) call pspygn (xx,yy,ivert,0.,0,7) end if end do c second miscibility do i = 1, ib if (iasmbl(i).eq.4) then call merger (i,iperm,ivert,ipoint,xx,yy) call pspygn (xx,yy,ivert,1.,0,0) end if end do c third univariant do i = 1, ib if (iasmbl(i).eq.2) then call merger (i,iperm,ivert,ipoint,xx,yy) call pspygn (xx,yy,ivert,1.,0,3) end if end do c last invariant do i = 1, ib if (iasmbl(i).eq.1) then do j = 1, 3 id = idf(j,i) ivchk(id) = 1 xx(j) = x(1,id) yy(j) = x(2,id) end do call pspygn (xx,yy,3,1.,0,0) end if end do else c old mode, draw tielines: do i = 1, ib ias = iasmbl(i) do j = 1, 3 id = idf(j,i) if (ias.eq.1.or.id.le.ipoint) ivchk(id) = 1 xx(j) = x(1,id) yy(j) = x(2,id) end do rias = float(ias) c draw tie lines: if (ias.eq.1) then c invariant call pspygn (xx,yy,3,rias,0,0) else if (ias.eq.2) then c univariant call pspygn (xx,yy,3,rias,0,0) c draw a thin line c between pseudocompounds of the c same phase do j = 1, 3 id = idf(iperm(1,j),i) jd = idf(iperm(2,j),i) c draw the line: if (ikp(jd).eq.ikp(id)) call psline * (xx(iperm(1,j)),yy(iperm(1,j)), * xx(iperm(2,j)),yy(iperm(2,j)),1.,0) end do else if (ias.eq.3) then c divariant call pspygn (xx,yy,3,rias,0,7) else if (ias.eq.4) then c immiscibility call pspygn (xx,yy,3,rias,0,0) end if end do end if do i = istct, iphct if (ivchk(i).ne.0) then call pssctr (8,1.2,1.2,0.) x1 = x(1,i) + 0.015 y1 = x(2,i) call pstext (x1,y1,names(i),8) call pselip (x1-0.015,y1,0.007,0.007,0.,0,7) end if end do c draw bounding triangle call pspygn (x3,y3,3,1.,0,0) c sectioning constraints call pssctr (8,1.2,1.2,0.) y = 1. do i = 1, ipot write (record,1040) vname(i),vmin(i) nchar = 0 call psublk (record,nchar) call pstext (0.75,y,record,nchar) y = y - 0.06 end do if (ifyn.eq.0) then call pstext (0.75,y,'(fluid saturated)',17) y = y - 0.06 end if if (isyn.eq.0) then write (record,1050) (names(idss(i)),i=1,isat) nchar = 72 call psublk (record,nchar) call pstext (0.75,y,record,nchar) end if goto 10 1000 format (2x,a30) 1010 format (a162) 1020 format (10a8) 1030 format (bn,i80) 1040 format (a8,'=',g9.3) 1050 format ('+ ',6(a8,' ')) 1060 format (' Draw tielines (y/n[recommended])?') 1070 format (a1) 1080 format (5a18) 1090 format (8a10) 99 end subroutine merger (i,iperm,ivert,ipoint,xx,yy) implicit real (a-g,o-z),integer (h-n) parameter (k1=5000,k2=40000,k5=12,h9=30) real xx(1000), yy(1000) integer iperm(2,3), idv(1000), idp(3), jdv(3) common/ asmbl /idf(3,k2),ib,iasmbl(k2),ivchk(k1),x(2,k1) * / phase /ikp(k1) c merge high variance fields c load first part into polygon: ias = iasmbl(i) ivert = 3 do j = 1, 3 id = idf(j,i) idv(j) = id c identify the phases if (ikp(id).ne.0) then idp(j) = -ikp(id) c tag used endmember compounds: if (id.le.ipoint) then ivchk(id) = 1 end if else idp(j) = id c tag used compounds: ivchk(id) = 1 end if end do c now search for the common simplexes: 50 do 20 j = 1, ib jas = iasmbl(j) c no match, reject: if (jas.ne.ias.or.j.eq.i) goto 20 do 30 k = 1, 3 id = idf(k,j) jdv(k) = id c check if the phases match: if (ikp(id).ne.0) then if (id.lt.ipoint) ivchk(id) = 1 do l = 1, 3 if (-ikp(id).eq.idp(l)) goto 30 end do c no match, reject: goto 20 else do l = 1, 3 if (id.eq.idp(l)) goto 30 end do c no match, reject: goto 20 end if 30 continue c ok, the simplex has the same assemblage c now match vertices: jvert = ivert do k = 1, ivert if (k.lt.ivert) then kp = k + 1 else kp = 1 end if do l = 1, 3 if (jdv(iperm(1,l)).eq.idv(k).and. * jdv(iperm(2,l)).eq.idv(kp).or. * jdv(iperm(2,l)).eq.idv(k).and. * jdv(iperm(1,l)).eq.idv(kp)) then if (kp.eq.1) then c tag the common join on the end idv(ivert+1) = jdv(4-l) else c the simplex might share the next c segment as well: if (ias.gt.2) then if (kp.lt.ivert) then kp2 = kp + 1 else kp2 = 1 end if do l2 = 1, 3 if (jdv(iperm(1,l2)).eq.idv(kp).and. * jdv(iperm(2,l2)).eq.idv(kp2).or. * jdv(iperm(2,l2)).eq.idv(kp).and. * jdv(iperm(1,l2)).eq.idv(kp2)) then c the simplex matches two segments c remove the matched point do l3 = kp+1, ivert idv(l3-1) = idv(l3) end do iasmbl(j) = -1 jvert = jvert - 1 goto 40 end if end do end if c found a common join, displace the old vertices do l1 = ivert + 1, kp, -1 idv(l1) = idv(l1-1) end do c insert the new vertex idv(kp) = jdv(4-l) end if c reset the assemblage flag: iasmbl(j) = -1 jvert = ivert + 1 goto 40 end if end do end do 20 continue c if got to here then there were no c no matches, save the polygon do j = 1, ivert c the polygon now has ivert vertices c load the vertices xx(j) = x(1,idv(j)) yy(j) = x(2,idv(j)) end do goto 99 c if here found a match, check the c list for more 40 ivert = jvert goto 50 99 end c---------------------------------------------------------------------- subroutine psaxes (jop0) c psaxes - subroutine to output (sloppy) axes. implicit real (a-g,o-z),integer (h-n) 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/ basic /n4,ifont,icopt,iwidth,iop0 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) c draw right vertical axis call psytic (xmax, y0, dy, -xtic, -xtic1, -xtic2, rline, * ihalf, itenf) c draw bottom horizontal axis call psxtic (ymin, x0, dx, ytic, ytic1, ytic2, rline, * ihalf, itenf) c draw top horizontal axis call psxtic (ymax, x0, dx, -ytic, -ytic1, -ytic2, rline, * ihalf, itenf) 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 subroutine psytic (x0, y0, dy, tic, tic1, tic2, rline, ihalf, * itenf) implicit real (a-g,o-z),integer (h-n) common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen common/ basic /n4,ifont,icopt,iwidth,iop0 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) implicit real (a-g,o-z),integer (h-n) common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen common/ basic /n4,ifont,icopt,iwidth,iop0 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-dx2, y0+tic1, 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. implicit real (a-g,o-z),integer (h-n) 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) implicit real (a-g,o-z),integer (h-n) 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. implicit real (a-g,o-z),integer (h-n) 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 c----------------------------------------------------------------- subroutine psalbl (x,y,ipr,ivar,ird,ityp,rline,iop9,fac1) c psalbl - subroutine to output reaction(ityp=0) or ip(ityp=1) c name label. implicit real (a-g,o-z),integer (h-n) common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen common/ basic /n4,ifont,icopt,jwidth,iop0 real x(ipr),y(ipr) character*4 lnms, pnms*6, string*400 if (ivar.eq.1) then fac = 1.0 * fac1 else fac = 0.75 * fac1 end if if (ityp.eq.0) then c get curve midpoint dx = (x(ipr) - x(1)) dy = (y(ipr) - y(1)) delta = 1d30 if (abs(dx)/xlen.gt.abs(dy)/ylen) then xmid = x(1) + dx / 2. do 1 i = 2, ipr if (abs(x(i)-xmid).lt.delta) then delta = abs(x(i)-xmid) imid = i end if 1 continue else ymid = y(1) + dy / 2. do 2 i = 2, ipr if (abs(y(i)-ymid).lt.delta) then delta = abs(y(i)-ymid) imid = i end if 2 continue end if if (iop9.eq.0) then c numeric label: write (lnms,1000) ird call pssctr (ifont,fac,fac,0.0) if (ird.lt.10) then call pselip (x(imid),y(imid), * fac*dcx,fac*1.2*dcy,rline,0,1) call pstext (x(imid)-2.45*dcx*fac, * y(imid)+0.8*dcy*fac,lnms,4) else if (ird.lt.100) then call pselip (x(imid),y(imid),fac*1.74*dcx,fac*1.6*dcy, * rline,0,1) call pstext (x(imid)-2.71*dcx*fac, * y(imid)+0.8*dcy*fac,lnms,4) else if (ird.lt.1000) then call pselip (x(imid),y(imid), * fac*2.75*dcx,fac*1.83*dcy,rline,0,1) call pstext (x(imid)-3.04*dcx*fac, * y(imid)+0.8*dcy*fac,lnms,4) else call pselip (x(imid),y(imid), * fac*3.75*dcx,fac*2.*dcy,rline,0,1) call pstext (x(imid)-4.*dcx*fac, * y(imid)+0.8*dcy*fac,lnms,4) end if else if (x(imid).eq.x(imid-1)) then theta = 1.5708 else if (y(imid).eq.y(imid-1)) then theta = 0.0 else c text label, get angle: theta = atan ( ((y(imid)-y(imid-1)) / ylen) * / ((x(imid)-x(imid-1)) / xlen) ) end if call rxntxt (string,iend) c normalized x text length: clx = iend * dcx * fac / xlen / 2. c normalized y text displacement: cly = .5 * dcy / ylen c get displacements from midpoint: dy = (sin(theta) * clx + cos(theta) * cly) * ylen dx = (cos(theta) * clx - sin(theta) * cly) * xlen theta = 57.29578 * theta call pssctr (ifont,fac,fac,theta) call pstext (x(imid)-dx,y(imid)-dy,string,iend) end if else c invariant point labels: if (ird.gt.999) then write (pnms,1040) ird else if (ird.gt.99) then write (pnms,1010) ird else if (ird.gt.9) then write (pnms,1020) ird else write (pnms,1030) ird end if call pssctr (ifont,1.2*fac,1.2*fac,0.0) call pstext (x(1)+.7*dcx,y(1)+.5*dcy,pnms,5) end if 1000 format (i4) 1010 format ('(',i3,') ') 1020 format ('(',i2,') ') 1030 format ('(',i1,') ') 1040 format ('(',i4,')') end subroutine psvdrw (pname,iop1,jop0) c---------------------------------------------------------------- c psvdrw - subroutine to output schreinemakers diagrams. implicit real (a-g,o-z),integer (h-n) character*162 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) common/ basic /n4,ifont,icopt,jwidth,iop0 c--------------------------------------------------------------- iop2 = 0 iop3 = 0 iop5 = 0 iop6 = 0 iop7 = 0 iop8 = 1 fac1 = .6 fac2 = .2 fac3 = .05 if (jop0.eq.1) then write (*,1210) read (*,1000) yes if (yes.ne.'y'.and.yes.ne.'Y') iop8 = 0 write (*,1230) read (*,1000) yes if (yes.ne.'y'.and.yes.ne.'Y') goto 40 write (*,1220) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') iop2 = 1 if (iop2.eq.0) then write (*,1070) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then 10 write (*,1030) read (*,*) fac2 if (fac2.lt.0.0.or.fac2.gt.1.0) then write (*,1200) fac2 goto 10 end if 20 write (*,1110) fac2 read (*,*) fac3 if (fac3.lt.0.0.or.fac3.gt.fac2) then write (*,1200) fac3 goto 20 end if end if write (*,1080) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then 25 write (*,1090) read (*,*) fac1 if (fac1.le.0.) then write (*,1200) fac1 goto 25 end if end if end if write (*,1190) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') iop3 = 1 end if 40 if (iop0.eq.1) then write (*,1150) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then write (*,1160) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then write (*,1170) read (*,1000) yes iop4c = 0 if (yes.eq.'y'.or.yes.eq.'Y') iop4c = 1 write (*,1180) read (*,1000) yes iop4p = 0 if (yes.eq.'y'.or.yes.eq.'Y') iop4p = 1 end if write (*,1140) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop5 = 1 prompt='present in the' call rname (1,prompt) end if write (*,1120) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop6 = 1 prompt=' absent in all' call rname (2,prompt) end if write (*,1130) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop7 = 1 prompt='present in all' call rname (3,prompt) end if end if end if call pssctr (8,1.2,1.2,0.) call pscurv (iop1,iop2,iop4c,iop5,iop6,iop7,iop8, * fac1,fac2,fac3,iend) if (iend.eq.1) goto 100 call psipts (iop1,iop3,iop4p,iop5,iop6,iop7,fac1) 100 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 rname (iw,prompt) implicit real (a-g,o-z),integer (h-n) character*10 xnams, xnam, prompt*14 common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ excl4 /xnams(50,3) jxct = 1 110 write (*,1040) prompt read (*,1020) xnam if (xnam.eq.' ') goto 210 call matchi (xnam,iex(jxct,iw),jex(jxct,iw)) if (iex(jxct,iw).eq.-1) then write (*,1100) xnam goto 110 end if xnams(jxct,iw) = xnam jxct = jxct + 1 goto 110 210 ixct(iw) = jxct - 1 1020 format (a10) 1040 format (/,' Enter the name of a phase ',a14,' fields ', * /,' (left justified, to finish): ',/) 1100 format (/,' No such entity as ',a10,', try again: ',/) end c--------------------------------------------------------------- subroutine pscurv (iop1,iop2,iop4,iop5,iop6,iop7, * iop8,fac1,fac2,fac3,iend) c pscurv - subroutine to output curves. implicit real (a-g,o-z),integer (h-n) parameter (k1=5000,k2=6000,k5=12,k7=k5+1,k0=4000) integer k, kex(k7), idr, lchk(k2) real x(k0), y(k0) character*10 names*8, xnams common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen common/ rxn /vnu(k7),idr(k7),ivct,iplus(k7),iminus(k7) common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ excl4 /xnams(50,3)/ phase /ikp(k1)/ cst8 /names(k1) common/ basic /n4,ifont,icopt,jwidth,iop0 data lchk/k2*0/ iend = 0 do 11 i = 1, 3 11 ict(i) = 0 10 read (n4,*,end=100) ipt,ird,ivar,ivct,(idr(i),i=1,ivct), * (vnu(i),i=1,ivct) if (ipt.eq.1) then goto 99 else if (ipt.eq.0) then goto 10 end if ipr = ipt/2 if (ird.gt.k2) then write (*,*) ' you gotta problem, increase k2 in pscurv' stop end if if (ipr.gt.k0*2) then write (*,*) ' you gotta problem, increase k0 in pscurv' write (*,*) ' a curve with ',ipr,' nodes?' stop end if read (n4,*) (x(i),y(i),i=1,ipr) if (iop4.eq.1.and.ivar.gt.1) goto 10 if (iop5.eq.1.or.iop6.eq.1.or.iop7.eq.1) then ifail = 1 if (iop5.eq.1.and.ivct.ge.ixct(1)) then c show only with assemblage: itic = 0 do 2 i = 1, ivct c if you get a match then look at c look for the next phase call checki (1,idr(i),kex(i)) if (kex(i).ge.0) then if (isoct.gt.0.and.kex(i).ne.0.and. * i.gt.1.and.itic.gt.0) then do 1 k = 1, i - 1 1 if (kex(k).eq.kex(i)) goto 2 end if itic = itic + 1 end if 2 continue c itic must = ixct if all phases match if (itic.lt.ixct(1)) goto 50 c if here the assemblage contains c all the phases requested: ict(1) = ict(1) + 1 ifail = 0 end if 50 if (iop6.eq.1.and.ivct.ne.0) then do 5 j = 1, ivct call checki (2,idr(j),itis) if (itis.ge.0) then c name matches a phase in the field ict(2) = ict(2) + 1 goto 60 end if 5 continue c no match ifail = 0 end if 60 if (iop7.eq.1.and.ivct.ne.0) then do 3 j = 1, ivct call checki (3,idr(j),itis) if (itis.ge.0) then c name matches a phase in the field ict(3) = ict(3) + 1 ifail = 0 goto 70 end if 3 continue c no match end if 70 if (ifail.eq.1) goto 10 end if if (iop1.eq.1) then jpr = 0 do 20 i = 1, ipr if (x(i).gt.xmax.or.x(i).lt.xmin.or. * y(i).gt.ymax.or.y(i).lt.ymin) goto 20 jpr = jpr + 1 x(jpr) = x(i) y(jpr) = y(i) 20 continue if (jpr.eq.0) goto 10 ipr = jpr end if if (ivar.eq.1) then iwidth = 1 else if (ivar.eq.98) then ivar = 1 iwidth = 2 else if (ivar.eq.99) then ivar = 1 iwidth = 0 else iwidth = 0 end if rline = float(ivar) iwidth = iwidth * jwidth if (iop8.eq.0) then call pspyln (x,y,ipr,rline,iwidth,0) else call psbspl (x,y,ipr,rline,iwidth,0) end if if (ird.eq.0) goto 10 if (iop2.eq.0.and.lchk(ird).eq.0) then yfrac = abs(y(ipr)-y(1)) / ylen xfrac = abs(x(ipr)-x(1)) / xlen if (xfrac.lt.fac3.and.yfrac.lt.fac3) then c curve is too short goto 10 else if (xfrac.lt.fac2.and.yfrac.lt.fac2) then c use a numeric label call psalbl (x,y,ipr,ivar,ird,0,rline,0,fac1) else c use a text label if (ivct.ne.0) call psalbl (x,y,ipr,ivar,ird,0,rline,1,fac1) end if if (ird.gt.k2) then write (*,1000) k2 stop end if lchk(ird) = 1 end if goto 10 100 write (*,1010) iend = 1 99 if (iop5.eq.1) write (6,*) ict(1), * ' curves have the assemblage: ', * (xnams(k,1),' ',k = 1, ixct(1)) if (iop6.eq.1) write (6,*) ict(2), * ' curves have none of the phases: ', * (xnams(k,2),' ',k = 1, ixct(2)) if (iop7.eq.1) write (6,*) ict(3), * ' curves have one of the phases: ', * (xnams(k,3),' ',k = 1, ixct(3)) 1000 format (' Too many reactions, increase parameter k2 (', * i4,') in routine pscurv',/) 1010 format (' End-of-file! if this is a PeRpLeX file ', * 'then the PeRpLeX',/, * ' program terminated incorrectly',/) end c--------------------------------------------------------------- subroutine plumin (ip,im) c subprogram to order a reactants according to sign) implicit real (a-g,o-z),integer (h-n) parameter (k5=12,k7=k5+1) real vp(k7), vm(k7) common/ rxn /vnu(k7),idr(k7),ivct,iplus(k7),iminus(k7) ip = 0 im = 0 c load id's of phases with + or - coeffs do 10 i = 1, ivct if (vnu(i).gt.0.0) then ip = ip + 1 vp(ip) = vnu(i) iplus(ip) = idr(i) else im = im + 1 vm(im) = vnu(i) iminus(im) = idr(i) end if 10 continue do 20 i = 1, im vnu(i) = vm(i) 20 idr(i) = iminus(i) do 30 i = 1, ip vnu(i+im) = vp(i) 30 idr(i+im) = iplus(i) end c--------------------------------------------------------------- subroutine rxntxt (string,iend) c subprogram to write a text label for a reaction implicit real (a-g,o-z),integer (h-n) parameter (k1=5000,k5=12,k7=k5+1,h9=30) character*8 names, text(400)*1, fname*10, string*400 common/ csta7 /fname(h9)/ cst8 /names(k1)/ phase /ikp(k1) common/ rxn /vnu(k7),idr(k7),ivct,iplus(k7),iminus(k7) c load id's of phases with + or - coeffs call plumin (ip,im) c now dump the names into the array text ist = 1 is = 1 80 do 20 i = is, im id = idr(i) if (ikp(id).eq.0) then c simple compound: read (names(id),1000) (text(j), j = ist, ist + 7) text(ist+8) = ' ' ist = ist + 9 else c solution phases: read (fname(ikp(id)),1000) (text(j), j = ist, ist + 9) text(ist+10) = '(' read (names(id),1000) (text(j), j = ist + 11, ist + 18) text(ist+19) = ')' text(ist+20) = ' ' ist = ist + 21 end if 20 continue if (is.eq.1) then text(ist) = '=' text(ist+1) = ' ' ist = ist + 2 is = im + 1 im = ivct goto 80 else text(ist) = ' ' end if c now filter out blanks iend = 1 do 40 i = 2, ist if (text(i).eq.' '.and.text(i+1).eq.' ') then goto 40 else if (text(i).eq.' '.and.text(i+1).eq.')') then goto 40 else if (text(i).eq.' '.and.text(i+1).eq.'(') then goto 40 end if iend = iend + 1 text(iend) = text(i) 40 continue if (iend.gt.400) iend = 400 write (string,1000) (text(i),i=1,400) 1000 format (400a1) end c--------------------------------------------------------------- subroutine psipts (iop1,iop3,iop4,iop5,iop6,iop7,fac1) c psipts - subprogram to output invariant points. implicit real (a-g,o-z),integer (h-n) parameter (k5=12,k8=k5+2) character*10 xnams real x(1), y(1) integer ipds(k8), kex(k8), k common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ basic /n4,ifont,icopt,jwidth,iop0/ excl4 /xnams(50,3) common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen data iop9/0/ do 11 i = 1, 3 11 ict(i) = 0 read (n4,*) ipct,incts if (ipct.eq.0) goto 99 do 20 i = 1, ipct read (n4,*) ipid,ipvar,(ipds(j), j=1,incts),x,y if (iop4.eq.1.and.ipvar.ne.1) goto 20 if (iop1.eq.1) then if (x(1).gt.xmax.or.x(1).lt.xmin.or. * y(1).gt.ymax.or.y(1).lt.ymin) goto 20 end if if (iop5.eq.1.or.iop6.eq.1.or.iop7.eq.1) then ifail = 1 if (iop5.eq.1.and.incts.ge.ixct(1)) then itic = 0 do 2 j = 1, incts call checki (1,ipds(j),kex(j)) if (kex(j).ge.0) then if (isoct.gt.0.and.kex(j).ne.0 * .and.j.gt.1.and.itic.gt.0) then do 1 k = 1, j-1 1 if (kex(k).eq.kex(j)) goto 2 end if itic = itic + 1 end if 2 continue if (itic.lt.ixct(1)) goto 50 c if here the assemblage contains c all the phases requested: ifail = 0 ict(1) = ict(1) + 1 end if 50 if (iop6.eq.1.and.incts.ne.0) then do 5 j = 1, incts call checki (2,ipds(j),itis) if (itis.ge.0) then c name matches a phase in the field ict(2) = ict(2) + 1 goto 60 end if 5 continue c no match ifail = 0 end if 60 if (iop7.eq.1.and.incts.ne.0) then do 3 j = 1, incts call checki (3,ipds(j),itis) if (itis.ge.0) then c name matches a phase in the field ict(3) = ict(3) + 1 ifail = 0 goto 70 end if 3 continue end if 70 if (ifail.eq.1) goto 20 end if if (ipvar.eq.1) then r = .78 else r = .39 end if rline = ipvar call pselip (x(1),y(1),r*dcx,r*dcy,0.,0,7) if (iop3.eq.0) call psalbl (x,y,1,ipvar,ipid,1, * rline,iop9,fac1) 20 continue if (iop5.eq.1) write (6,*) ict(1), * ' points have the assemblage: ', * (xnams(k,1),' ',k = 1, ixct(1)) if (iop6.eq.1) write (6,*) ict(2), * ' points do not have any of the phases: ', * (xnams(k,2),' ',k = 1, ixct(2)) if (iop7.eq.1) write (6,*) ict(3), * ' points have one of the phases: ', * (xnams(k,3),' ',k = 1, ixct(3)) 99 end subroutine matchi (unnown,itis,icpd) c---------------------------------------------------------------- c matchi - subroutine to determine if the string unnown is a valid c solution or compound name. c itis = 0 if compound, = ikp if solution, and =-1 if invalid c icpd = the id of the compound if itis = 0 implicit real (a-g,o-z),integer (h-n) character*10 sname, names*8, unnown parameter (k1=5000,h9=30) common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ cst8 /names(k1)/ csta7 /sname(h9)/ phase /ikp(k1) if (isoct.ne.0) then do 10 i = 1, isoct if (unnown.eq.sname(i)) then itis = i goto 99 end if 10 continue end if do 20 i = 1, iphct if (unnown.eq.names(i)) then itis = 0 icpd = i goto 99 end if 20 continue itis = -1 99 end subroutine checki (iw,iun,itis) c---------------------------------------------------------------- c checki - subroutine to determine if the cpd iun is in the c list iex or jex solution or compound name. c itis = -1 if not, itis>=0 if it is. implicit real (a-g,o-z),integer (h-n) parameter (k1=5000) common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ phase /ikp(k1) if (ikp(iun).ne.0) then id = ikp(iun) do 10 i = 1, ixct(iw) if (id.eq.iex(i,iw)) then itis = iex(i,iw) goto 99 end if 10 continue end if do 20 i = 1, ixct(iw) if (iun.eq.jex(i,iw)) then itis = 0 goto 99 end if 20 continue itis = -1 99 end subroutine psmixd implicit real (a-g,o-z),integer (h-n) c psmixd - subroutine to draw binary mixed variable diagrams parameter (k1=5000,k5=12,k7=k5+1,h9=30) character names*8, vname*30, title*162, sname*10, string*400, y*1 character tname(5)*8,xname(k5)*30 integer idss(10), idf(3), jphi(k1) real v(5) common/ rxn /vnu(k7),idr(k7),ivct,iplus(k7),iminus(k7) common/ vs /vmin(5),vmax(5),iv/ vsa /vname(5) common/ cst8 /names(k1)/ csta7 /sname(h9)/ phase /ikp(k1) common/ tx /x(k1),iphi(k1),ivph(k1),iphct,ib,istct common / wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen common/ basic /n4,ifont,icopt,jwidth,iop0 / xyrat /xfac data idf,igo/4*0/ c start-of-header c ------------------------------ c read simple variables read (n4,*) icp,istct,iphct,ipoint,ifyn,isyn,ipot,isoct if (icp.ne.2) then write (*,*) ' Sorry, PSVDRAW only works', * ' for binary mixed variable diagrams.' goto 99 end if if (iv.gt.5.or.ib.gt.k1.or.ipot.gt.5) then write (*,*) ' dimensioning error in psmixd ' goto 99 end if if (isoct.ne.0) write (*,1060) read (n4,1000) (tname(i), i = 1, ipot) read (n4,1010) title c names read (n4,1020) (names(i),i=1,iphct) c phase coordinates read (n4,*) (x(i), i = istct, iphct) c read solution identifier flags read (n4,*) (ikp(i),i=1,iphct) c read solution names if (isoct.ne.0) read (n4,1050) (sname(i), i=1, isoct) read (n4,1030) (xname(i),i = 1, icp) c end-of-header c ------------------------------ read (n4,*) (v(i), i = 1, ipot) c read stable phase ids read (n4,*) ib read (n4,*) (iphi(i), i = 1, ib) c read saturated phase id's goto (10),isyn read (n4,*) isat read (n4,*) (idss(i),i=1,isat) c read index of y variable 10 read (tname(1),1020) vname(2) vname(1) = xname(2) read (n4,*) vmin(2), vmax(2) vmin(1) = 0.0 vmax(1) = 1.0 call psaxop (iop1,jop0) c get some other options iop1 = 1 iop2 = 0 iop3 = 0 fac = 1. rx = 0.5 * dcx/2./xfac ry = 0.5 * dcy/2. if (iop0.eq.1) then write (*,2050) read (*,2010) y if (y.eq.'y'.or.y.eq.'Y') iop2 = 1 write (*,2000) read (*,2010) y if (y.eq.'y'.or.y.eq.'Y') then write (*,2020) read (*,*) fac end if write (*,2030) read (*,2010) y if (y.eq.'y'.or.y.eq.'Y') then iop1 = 0 else write (*,2040) read (*,2010) y if (y.eq.'y'.or.y.eq.'Y') iop1 = 2 end if if (isoct.gt.0) then write (*,2060) read (*,2010) y if (y.eq.'y'.or.y.eq.'Y') iop3 = 1 end if end if t0 = vmin(2) tlast = t0 - t0 c get the variance of each compound: do 30 i = 1, ib 30 call getvar (i,iop3) 130 read (n4,*) ird, ivct, ivar, t1, (idr(i), i = 1, ivct), * (vnu(i), i = 1, ivct) if (ivct.eq.1) goto 90 do 135 i = 1, ib 135 jphi(i) = iphi(i) jb = ib call plumin (jplus,jminus) do 20 i = 2, ib - 1 rline = ivph(i) 20 if (rline.ne.0.0) call psline * (x(iphi(i)),t0,x(iphi(i)),t1,rline,0) isum = 0 idif = 0 do 50 i = 1, ivct if (ikp(idr(i)).ne.0) then isum = isum + 1 do 60 j = 1, idif 60 if (idf(j).eq.ikp(idr(i))) goto 50 idif = idif + 1 idf(idif) = ikp(idr(i)) end if 50 continue c consider different cases: if (ivct.eq.2) then c get index (id) of reacted phase: do 40 i = 1, ib if (iphi(i).eq.idr(1)) then i1 = i - 1 if (i1.lt.1) i1 = 1 id1 = iphi(i1) it = i i2 = i + 1 if (i2.gt.ib) i2 = ib id2 = iphi(i2) goto 55 end if 40 continue write (*,*) 'nope' write (*,*) (iphi(i),i=1,ib) write (*,*) idr(1), ird,t1 goto 99 c check if only a compositional c change in a solution: 55 if (isum.eq.2.and.ikp(idr(1)).eq.ikp(idr(2))) then c draw an open dot, but don't c do anything else: call pselip (x(idr(1)),t1,rx,ry,1.,0,0) iphi(it) = idr(2) goto 130 end if c draw dot call pselip (x(idr(1)),t1,rx,ry,1.,0,7) c draw line call psline (x(id1),t1,x(id2),t1,1.,0) c replace reacted phase: iphi(it) = idr(2) call getvar (it,iop3) if (it.gt.1) call getvar (it-1,iop3) if (it.lt.ib) call getvar (it+1,iop3) else c classify reaction, jminus = 2 c eutectoid, else peritectoid. if (jminus.eq.2) then i1 = idr(1) i2 = idr(2) i0 = idr(3) c find indices of old phases: it = 0 do 70 i = 1, ib if (iphi(i).eq.i1) id1 = i 70 if (iphi(i).eq.i2) id2 = i else i0 = idr(1) i1 = idr(2) i2 = idr(3) c find indices of old phases: do 80 i = 1, ib if (iphi(i).eq.i1) id1 = i 80 if (iphi(i).eq.i2) id2 = i end if if (id1.gt.id2) then it = id1 jt = i1 i1 = i2 i2 = jt id1 = id2 id2 = it if (jminus.eq.2) then idr(1) = i1 idr(2) = i2 else idr(2) = i1 idr(3) = i2 end if end if if (isum.lt.2.or.idif.eq.3) then c all phases different c draw line call psline (x(i1),t1,x(i2),t1,1.,0) call pselip (x(i0),t1,rx,ry,1.,0,7) else if (isum.eq.2) then if (idif.eq.2) then c all phases different call psline (x(i1),t1,x(i2),t1,1.,0) call pselip (x(i0),t1,rx,ry,1.,0,7) else c one solution + one compound if (ikp(i0).eq.0) then c compound in the middle c change variance of reaction ivar = 1 call pselip (x(i0),t1,rx,ry,1.,0,7) call psline (x(i1), t1, x(i2), t1, 1., 0) else c pseudoinvariant if (ikp(i0).eq.0) then call psline (x(i1),t1,x(i2),t1,1.,0) else if (ikp(i1).eq.0) then call psline (x(i0),t1,x(i2),t1,1.,0) else call psline (x(i0),t1,x(i1),t1,1.,0) end if end if end if else if (isum.eq.3.and.idif.eq.1.and.igo.eq.0) then c reaction on solvus c check if on left line itot = 0 call miscib (x(i1),x(i0),ikp(i0),imis,igo) if (imis.eq.0) call psline (x(i0),t1,x(i1),t1,1.,0) itot = itot + imis call miscib (x(i0),x(i2),ikp(i0),imis,igo) itot = itot + imis if (imis.eq.0) call psline (x(i0),t1,x(i2),t1,1.,0) if (itot.eq.0) then ivar = 1 call pselip (x(i0),t1,rx,ry,1.,0,7) end if else c two solutions present. c solvus possible, find c the pseudocompounds of c the same solution: if (ikp(i0).ne.ikp(i1).and.ikp(i0).ne.ikp(i2)) then c i1 and i2 are the same: call miscib (x(i1),x(i2),ikp(i1),imis,igo) else if (ikp(i0).eq.ikp(i1)) then c i0 and i1 are the same: call miscib (x(i1),x(i0),ikp(i0),imis,igo) else c i0 and i2 are the same: call miscib (x(i0),x(i2),ikp(i0),imis,igo) end if if (imis.eq.1) then ivar = 1 c eutectic/peritectic with solvus call psline (x(i1),t1,x(i2),t1,1.,0) call pselip (x(i0),t1,rx,ry,1.,0,7) else if (ikp(i1).eq.ikp(i0)) then c pseudounivariant, left call psline (x(i1),t1,x(i0),t1,1.,0) else if (ikp(i2).eq.ikp(i0)) then c pseudounivariant, right call psline (x(i2),t1,x(i0),t1,1.,0) else c else center c compound in the middle c change variance: itoc = 1 do 210 j = 1, ib if (iphi(j).eq.i0) goto 210 if (ikp(iphi(j)).eq.ikp(i0)) itoc = 0 210 continue if (itoc.eq.1) ivar = 1 if (ivar.eq.1) call pselip * (x(idr(1)),t1,rx,ry,1.,0,7) call psline ( x(iphi(id1)), t1, * x(iphi(id2)), t1, 1., 0) end if end if c get new chemography: if (jminus.eq.2) then c eutectoid, number of phases c increases. do 100 i = ib, id2, -1 iphi(i+1) = iphi(i) 100 ivph(i+1) = ivph(i) ib = ib + 1 iphi(id1+1) = i0 do 110 i = id1, id1 + 2 110 call getvar (i,iop3) else c peritectoid, number of phases c decreases. do 120 i = id2, ib iphi(i-1) = iphi(i) 120 ivph(i-1) = ivph(i) ib = ib - 1 do 140 i = id1, id1 + 1 140 call getvar (i,iop3) end if end if if (iop2.eq.1) then do 170 i = 1, ib - 1 if (ikp(iphi(i)).eq.0) goto 170 i1 = i + 1 if (ikp(iphi(i)).ne.ikp(iphi(i1))) goto 170 do 160 j = 1, jb - 1 if (x(jphi(j)).gt.x(iphi(i))) then goto 170 else if (x(jphi(j)).eq.x(iphi(i))) then j1 = j + 1 if (ikp(jphi(j)).eq.ikp(iphi(i)).and. * ikp(jphi(j)).eq.ikp(jphi(j1))) then c draw a rectangle from c i-i1-j-j1 call psrect (x(iphi(i)),x(iphi(i1)), * t0,t1,0.,iwidth,2) end if end if 160 continue 170 continue end if t0 = t1 if (iop1.eq.0) then goto 130 else if (iop1.eq.1) then if (ivar.ne.1) goto 130 call rxntxt (string,iend) else call rxntxt (string,iend) end if call pssctr (7,fac,fac,0.) ttext = t1 + .9 * dcy if (ttext.lt.tlast+2.1*dcy) ttext = tlast + 2.1*dcy call pstext (xmax+2.*dcx,ttext,string,iend) tlast = ttext goto 130 90 t1 = vmax(2) c fill in last part of one phase regions if (iop2.eq.1) then do 180 i = 1, ib - 1 if (ikp(iphi(i)).eq.0) goto 180 i1 = i + 1 if (ikp(iphi(i)).ne.ikp(iphi(i1))) goto 180 c draw a rectangle from c i-i1-j-j1 call psrect (x(iphi(i)),x(iphi(i1)), * t0,t1,0.,iwidth,2) 180 continue end if do 150 i = 2, ib - 1 rline = ivph(i) if (rline.eq.0.0) goto 150 call psline (x(iphi(i)),t0,x(iphi(i)),t1,rline,0) 150 continue 99 call psaxes (jop0) 1000 format (2x,a8) 1010 format (a162) 1020 format (10a8) 1030 format (5a18) 1050 format (8a10) 1060 format (/,' WARNING: Psvdraw may not draw t-x diagrams', * ' correctly if immiscibility',/,' occurs in a', * ' projected ternary or higher order solution, to', * ' avoid',/,' problems turn the miscibility test', * ' off (modify default options).',/) 2000 format (/,' Modify default text scaling (y/n)?',/) 2010 format (a1) 2020 format (/,' Enter the scaling factor:',/) 2030 format (/,' Suppress reaction labels (y/n)?',/) 2040 format (/,' Show pseudounivariant reaction labels (y/n)', * /) 2050 format (/,' Try to fill one phase fields (y/n)?',/) 2060 format (/,' Skip tests for immiscibility (y/n)?',/) end c------------------------------------------------------------------- subroutine getvar (i,igo) c getvar - subroutine to get variance of compound i in a binary row c ivphi set 1 if true compound c ivphi set 0 if pseudocompound c ivphi set 2 if solvus limb and igo = 0. c if igo = 1 skip miscibility test implicit real (a-g,o-z),integer (h-n) parameter (k1=5000) common/ phase /ikp(k1)/tx/x(k1),iphi(k1),ivph(k1),iphct,ib,istct if (i.eq.1.or.i.eq.ib) then ivph(i) = 1 goto 30 end if i1 = i - 1 i2 = i + 1 id1 = iphi(i1) id = iphi(i) id2 = iphi(i2) if (ikp(id).eq.0) then c true compound ivph(i) = 1 else if (ikp(id1).eq.ikp(id).and.ikp(id2).eq.ikp(id)) then c only one solution: ivph(i) = 0 call miscib (x(id1), x(id), ikp(id), imis, igo) if (imis.eq.1) ivph(i) = 1 call miscib (x(id), x(id2), ikp(id), imis, igo) if (imis.eq.1) ivph(i) = 1 else ivph(i) = 1 end if 30 end c------------------------------------------------------------------- subroutine miscib (x1, x2, ids, imis, igo) c miscib - subroutine check for immiscibility in a binary phase c imis = 1 => solvus between solution(ids) compositions x1 and x2 implicit real (a-g,o-z),integer (h-n) parameter (k1=5000) character*8 names common/ cst8 /names(k1) common/ phase /ikp(k1) common/ tx /x(k1),iphi(k1),ivph(k1),iphct,ib,istct imis = 0 c disabled cause of projections goto (99), igo do 10 i = istct, iphct if (ikp(i).ne.ids) goto 10 if (x(i).gt.x1.and.x(i).lt.x2) then imis = 1 goto 99 end if 10 continue 99 end subroutine psbplt (pname,ipoly,iop1,jop0) c---------------------------------------------------------------- c psbplt - subroutine to fixed bulk composition sections implicit real (a-g,o-z),integer (h-n) parameter (k5=12) integer ilop(9) character*162 pname(4), yes*1, prompt*14, propty(3)*60 common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen * / excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) * / basic /n4,ifont,icopt,jwidth,iop0 common/ bulk /ipid,ipvar,amode(k5),vmode(k5),wtp(k5),rho, * ipds(k5),ivol,incts,pmn,dp save propty data propty/'1 - Density (g/cm3) in solid assemblage', * '2 - Wt % of volatile 1 in solid assemblage', * '3 - Wt % of volatile 2 in solid assemblage'/ c--------------------------------------------------------------- iop5 = 0 iop6 = 0 iop7 = 0 ivmn = 0 ivmx = 100 do i = 1, 9 ilop(i) = 0 end do fac1 = .6 if (jop0.eq.1) then c drafting options c polygon fill mode: write (*,1010) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then write (*,1030) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then c contour by property ilop(8) = 2 c choose property if (ivol.eq.0) then ilop(9) = 0 else write (*,1050) do i = 1, 1 + ivol write (*,1060) propty(i) end do read (*,*) iprop ilop(9) = iprop - 1 end if c get range of property write (*,1070) read (*,*) pmn, dp dp = dp - pmn else c contour by variance write (*,1040) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(8) = 1 end if end if c text label size write (*,1080) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then 25 write (*,1090) read (*,*) fac1 if (fac1.le.0.) then write (*,1200) fac1 goto 25 end if end if c modify label: write (*,1180) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then write (*,1220) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(1) = 1 write (*,1230) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(2) = 1 write (*,1240) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(3) = 1 write (*,1250) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(4) = 1 write (*,1260) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(5) = 1 write (*,1270) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(6) = 1 if (ilop(5).eq.1.and.ilop(6).eq.0) then write (*,1280) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(7) = 1 end if end if end if 40 if (iop0.eq.1) then c phase field restrictions write (*,1150) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then c variance restrictions write (*,1160) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then write (*,1170) read (*,*) ivmn, ivmx end if c restrict by assemblage write (*,1140) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop5 = 1 prompt='present in the' call rname (1,prompt) end if c restrict by phase presence write (*,1120) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop6 = 1 prompt=' absent in all' call rname (2,prompt) end if c restrict by phase absence write (*,1130) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop7 = 1 prompt='present in all' call rname (3,prompt) end if end if end if call pssctr (8,1.2,1.2,0.) call psbpts (ipoly,iop1,iop5,iop6,iop7,ilop,ivmn,ivmx,fac1) 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 i = 2, 4 call pstext (x,y,pname(i),0) y = y - 0.03 * ylen end do 1000 format (a1) 1010 format ( ' Change polygon fill mode (y/n)? ',/) 1020 format (/,' Suppress polygon fill (y/n)? ',/) 1030 format (/,' Fill according to the value of a property (y/n)? ',/) 1040 format (/,' Fill according to variance (y/n)? ',/) 1050 format (/,' Fill according to the value of which property? ',/) 1060 format (' ',a60) 1070 format (/,' Enter the range (min->max) of the property to be', * ' contoured by false color:',/) 1080 format (/,' Reset size for text labels (y/n)? ',/) 1090 format (/,' Scale default text size by a factor of: ',/) 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 (/,' Restrict phase fields by variance (y/n)? ',/) 1170 format (/,' Enter minimum (2) and maximum (c+2) variance :',/) 1180 format (/,' Modify assemblage labeling (y/n):',/) 1200 format (/,1x,g13.6,' is an invalid value.',/) 1220 format (' Suppress numerical label (y/n)? ') 1230 format (' Suppress variance label (y/n)? ') 1240 format (' Suppress density (y/n)? ') 1250 format (' Suppress volatile concentrations (y/n)? ') 1260 format (' Suppress modes (y/n)? ') 1270 format (' Suppress phase names (y/n)? ') 1280 format (' Show only solution names (i.e. suppress compositions' * ,', y/n)? ') end subroutine pscplt (pname,iop1,jop0) c---------------------------------------------------------------- c pscplt - subroutine to plot bulk compositions implicit real (a-g,o-z),integer (h-n) parameter (k5=12) integer ilop(9) character*162 pname(4), yes*1, prompt*14, propty(3)*60 common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen * / excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) * / basic /n4,ifont,icopt,jwidth,iop0 common/ bulk /ipid,ipvar,amode(k5),vmode(k5),wtp(k5),rho, * ipds(k5),ivol,incts,pmn,dp save propty data propty/'1 - Density (g/cm3) in solid assemblage', * '2 - Wt % of volatile 1 in solid assemblage', * '3 - Wt % of volatile 2 in solid assemblage'/ c--------------------------------------------------------------- iop5 = 0 iop6 = 0 iop7 = 0 ivmn = 0 ivmx = 100 do i = 1, 9 ilop(i) = 0 end do fac1 = .6 if (jop0.eq.1) then c drafting options c text label size write (*,1080) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then 25 write (*,1090) read (*,*) fac1 if (fac1.le.0.) then write (*,1200) fac1 goto 25 end if end if c modify label: write (*,1180) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then write (*,1220) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(1) = 1 write (*,1230) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(2) = 1 write (*,1240) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(3) = 1 write (*,1250) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(4) = 1 write (*,1260) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(5) = 1 write (*,1270) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(6) = 1 if (ilop(5).eq.1.and.ilop(6).eq.0) then write (*,1280) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') ilop(7) = 1 end if end if end if 40 if (iop0.eq.1) then c phase field restrictions write (*,1150) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then c variance restrictions write (*,1160) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then write (*,1170) read (*,*) ivmn, ivmx end if c restrict by assemblage write (*,1140) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop5 = 1 prompt='present in the' call rname (1,prompt) end if c restrict by phase presence write (*,1120) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop6 = 1 prompt=' absent in all' call rname (2,prompt) end if c restrict by phase absence write (*,1130) read (*,1000) yes if (yes.eq.'y'.or.yes.eq.'Y') then iop7 = 1 prompt='present in all' call rname (3,prompt) end if end if end if call pssctr (8,1.2,1.2,0.) call psclbl (iop1,iop5,iop6,iop7,ilop,ivmn,ivmx,fac1) 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 i = 2, 4 call pstext (x,y,pname(i),0) y = y - 0.03 * ylen end do 1000 format (a1) 1010 format ( ' Change polygon fill mode (y/n)? ',/) 1020 format (/,' Suppress polygon fill (y/n)? ',/) 1030 format (/,' Fill according to the value of a property (y/n)? ',/) 1040 format (/,' Fill according to variance (y/n)? ',/) 1050 format (/,' Fill according to the value of which property? ',/) 1060 format (' ',a60) 1070 format (/,' Enter the range (min->max) of the property to be', * ' contoured by false color:',/) 1080 format (/,' Reset size for text labels (y/n)? ',/) 1090 format (/,' Scale default text size by a factor of: ',/) 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 (/,' Restrict phase fields by variance (y/n)? ',/) 1170 format (/,' Enter minimum (2) and maximum (c+2) variance :',/) 1180 format (/,' Modify assemblage labeling (y/n):',/) 1200 format (/,1x,g13.6,' is an invalid value.',/) 1220 format (' Suppress numerical label (y/n)? ') 1230 format (' Suppress variance label (y/n)? ') 1240 format (' Suppress density (y/n)? ') 1250 format (' Suppress volatile concentrations (y/n)? ') 1260 format (' Suppress modes (y/n)? ') 1270 format (' Suppress phase names (y/n)? ') 1280 format (' Show only solution names (i.e. suppress compositions' * ,', y/n)? ') end c--------------------------------------------------------------- subroutine psbpts (ipoly,iop1,iop5,iop6,iop7,ilop,ivmn,ivmx,fac1) c psbpts - subprogram to output bulk composition data. implicit real (a-g,o-z),integer (h-n) parameter (k0=4000,k2=6000,k5=12,k7=k5+1,k8=k5+2) real x(k0), y(k0) character*10 xnams integer k, kex(k8), ilop(9) common/ rxn /vnu(k7),idr(k7),ivct,iplus(k7),iminus(k7) common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ basic /n4,ifont,icopt,jwidth,iop0/ excl4 /xnams(50,3) common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen common/ bulk /ipid,ipvar,amode(k5),vmode(k5),wtp(k5),rho, * ipds(k5),ivol,incts,pmn,dp do i = 1, 3 ict(i) = 0 end do n5 = n4 + 1 rline = 2. iwidth = 0 propmn = 1e30 propmx = -1e30 20 read (n4,*,end=90) ipt,ird,ivar,ivct,(idr(i),i=1,ivct), * (vnu(i),i=1,ivct) if (ipt.eq.1) then goto 99 else if (ipt.eq.0) then goto 20 end if ipr = ipt/2 if (ird.gt.k2) then write (*,*) ' you gotta problem, increase k2 in pscurv' stop end if if (ipr.gt.k0*2) then write (*,*) ' you gotta problem, increase k0 in pscurv' write (*,*) ' a curve with ',ipr,' nodes?' stop end if read (n4,*) (x(i),y(i),i=1,ipr) if (ipoly.eq.0) then call pspygn (x,y,ipr,rline,iwidth,3) goto 20 end if read (n5,*,end=900) ipid,ipvar,(ipds(i), i = 1, incts) read (n5,*) (vmode(i), i = 1, incts) read (n5,*) (amode(i), i = 1, incts) read (n5,*) xx, yy read (n5,*) rho if (ivol.gt.0) read (n5,*) (wtp(i),i=1,ivol) c variance restriction if (ipvar+2.gt.ivmx.or.ipvar+2.lt.ivmn) goto 20 if (iop1.eq.1.and.(xx.gt.xmax.or.xx.lt.xmin.or. * yy.gt.ymax.or.yy.lt.ymin)) goto 20 if (iop5.eq.1.or.iop6.eq.1.or.iop7.eq.1) then ifail = 1 if (iop5.eq.1.and.incts.ge.ixct(1)) then itic = 0 do 2 j = 1, incts call checki (1,ipds(j),kex(j)) if (kex(j).ge.0) then if (isoct.gt.0.and.kex(j).ne.0. * and.j.gt.1.and.itic.gt.0) then do 1 k = 1, j-1 1 if (kex(k).eq.kex(j)) goto 2 end if itic = itic + 1 end if 2 continue if (itic.lt.ixct(1)) goto 50 c if here the assemblage contains c all the phases requested: ifail = 0 ict(1) = ict(1) + 1 end if 50 if (iop6.eq.1.and.incts.ne.0) then do 5 j = 1, incts call checki (2,ipds(j),itis) if (itis.ge.0) then c name matches a phase in the field ict(2) = ict(2) + 1 goto 60 end if 5 continue c no match ifail = 0 end if 60 if (iop7.eq.1.and.incts.ne.0) then do 3 j = 1, incts call checki (3,ipds(j),itis) if (itis.ge.0) then c name matches a phase in the field ict(3) = ict(3) + 1 ifail = 0 goto 70 end if 3 continue end if 70 if (ifail.eq.1) goto 20 end if r = .78 if (ilop(8).eq.0) then c no fill call pspyln (x,y,ipr,1.,1,0) else if (ilop(8).eq.1) then c variance based fill if (ipvar.eq.0) then call pspygn (x,y,ipr,rline,iwidth,0) else if (ipvar.le.5) then c for di-septa-variant c use gray-scale rfill = 1. - float (ipvar) * 0.2 call pspygr (x,y,ipr,rline,iwidth,rfill) else c use pattern fills call pspygn (x,y,ipr,rline,iwidth,3+ipvar) end if else if (ilop(8).eq.2) then if (ilop(9).eq.0) then prop = rho else prop = wtp(ilop(9)) end if if (prop.lt.propmn) propmn = prop if (prop.gt.propmx) propmx = prop rfill = (prop-pmn) / dp if (rfill.lt.0.) then rfill = 0. else if (rfill.gt.1) then rfill = 1. end if rfill = 1. - rfill call pspygr (x,y,ipr,rline,iwidth,rfill) end if goto 20 90 if (ilop(8).eq.2) write (*,1010) propmn,propmx,pmn,pmn+dp if (iop5.eq.1) write (6,*) ict(1), * ' fields have the assemblage: ', * (xnams(k,1),' ',k = 1, ixct(1)) if (iop6.eq.1) write (6,*) ict(2), * ' fields do not have any of the phases: ', * (xnams(k,2),' ',k = 1, ixct(2)) if (iop7.eq.1) write (6,*) ict(3), * ' fields do have one of the phases: ', * (xnams(k,3),' ',k = 1, ixct(3)) goto 99 900 write (*,1000) c stick the labels on to the c points after drawing the polygons: 99 if (ipoly.ne.0) call psblbl (iop1,iop5,iop6,iop7,ilop, * ivmn,ivmx,fac1) 1000 format (/,' your polygon file is inconsistent with ', * ' your assemblage file, i am quitting.',/) 1010 format (/,' The property you chose for false color ranges from', * g12.6,'->',g12.6,/,' The range you specified for scaling', * ' was ',g12.6,'->',g12.6,/) end subroutine ftext (text,ist,iend) c----------------------------------------------------------------------- c subprogram to filter blanks from text c----------------------------------------------------------------------- implicit real (a-g,o-z),integer (h-n) character*1 text(400) itic = ist - 1 igot = 0 jend = iend do 40 i = ist, iend-1 if (text(i).eq.' '.and.text(i + 1).eq.' '.or. * text(i).eq.' '.and.text(i + 1).eq.')'.or. * text(i).eq.' '.and.text(i + 1).eq.'('.or. * igot.eq.0.and.text(i).eq.' ') goto 40 if (i.gt.ist) then if (text(i-1).eq.'-'.and.text(i).eq.' ') goto 40 end if itic = itic + 1 igot = 1 text(itic) = text(i) 40 continue if (text(iend).ne.' ') then itic = itic + 1 text(itic) = text(iend) end if iend = itic + 1 do i = iend, jend text(i) = ' ' end do end subroutine psbtxt (string,iend,ilop) c----------------------------------------------------------------------- c subprogram to write a text labels for bulk composition output c----------------------------------------------------------------------- implicit real (a-g,o-z),integer (h-n) parameter (k1=5000,k5=12,h9=30) character*8 names, text(400)*1, fname*10, string*400, char*1 integer ilop(9), iwas(k5) common/ cst8 /names(k1)/ csta7 /fname(h9)/ phase /ikp(k1) common/ bulk /ipid,ipvar,amode(k5),vmode(k5),wtp(k5),rho, * ipds(k5),ivol,incts,pmn,dp 1000 format (400a1) 1020 format (1x,f6.2,1x,a8) 1030 format (a10) 1040 format (a8) 2000 format (i5) 2010 format (f6.3) 2020 format (f6.2,1x,f6.2) 2030 format (a1,i2) ist = 1 iend = 0 do i = 1, 400 text(i) = ' ' end do if (ilop(1).eq.0) then write (string,2000) ipid iend = ist + 4 read (string,1000) (text(j),j=ist,iend) call ftext (text,ist,iend) ist = iend end if if (ilop(2).eq.0) then char = ' ' if (ilop(1).eq.0) char = '-' write (string,2030) char, ipvar+2 iend = ist + 2 read (string,1000) (text(j), j = ist, iend) call ftext (text,ist,iend) ist = iend + 1 else if (ilop(1).eq.0) then ist = ist + 1 end if if (ilop(3).eq.0.or.ivol.ne.0.and.ilop(4).eq.0) then text(ist) = '(' ist = ist + 1 end if if (ilop(3).eq.0) then write (string,2010) rho iend = ist + 5 read (string,1000) (text(j), j= ist, iend) call ftext (text,ist,iend) ist = iend + 1 end if if (ivol.ne.0.and.ilop(4).eq.0) then write (string,2020) (wtp(i), i = 1, ivol) iend = ist + 5 + (ivol-1)*7 read (string,1000) (text(j),j = ist, iend) call ftext (text,ist,iend) end if if (ilop(3).eq.0.or.ivol.ne.0.and.ilop(4).eq.0) then text(iend) = ')' text(iend+1) = ' ' iend = iend + 1 end if if (ilop(6).ne.1) then if (ilop(4).ne.1) then do i = 1, incts write (string,1020) vmode(i),names(ipds(i)) ist = iend + 1 iend = ist + 15 read (string,1000) (text(j),j=ist,iend) call ftext (text,ist,iend) end do else if (ilop(7).eq.1) then ict = 0 do 10 i = 1, incts idsol = ikp(ipds(i)) if (idsol.ne.0) then do j = 1, ict if (iwas(j).eq.idsol) goto 10 end do ict = ict + 1 iwas(ict) = idsol write (string,1030) fname(idsol) else write (string,1030) names(ipds(i)) end if ist = iend + 1 iend = ist + 11 read (string,1000) (text(j),j=ist,iend) call ftext (text,ist,iend) 10 continue else do i = 1, incts write (string,1040) names(ipds(i)) ist = iend + 1 iend = ist + 9 read (string,1000) (text(j),j=ist,iend) call ftext (text,ist,iend) end do end if end if write (string,1000) (text(j),j=1,iend) end c--------------------------------------------------------------- subroutine psblbl (iop1,iop5,iop6,iop7,ilop,ivmn,ivmx,fac1) c psblbl - subprogram to output bulk composition label data. implicit real (a-g,o-z),integer (h-n) parameter (k0=4000,k2=6000,k5=12,k7=k5+1,k8=k5+2) real rnum character*10 xnams, text*400, char*10 integer k, kex(k8), ilop(9), int common/ rxn /vnu(k7),idr(k7),ivct,iplus(k7),iminus(k7) common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ basic /n4,ifont,icopt,jwidth,iop0/ excl4 /xnams(50,3) common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen common/ bulk /ipid,ipvar,amode(k5),vmode(k5),wtp(k5),rho, * ipds(k5),ivol,incts,pmn,dp n5 = n4 + 1 rewind (n5) c reread header read (n5,*) icopt, incts, ivol, iphct, isoct if (iphct.ne.0) then read (n5,1010) (char, i = 1, iphct) read (n5,*) (int, i = 1, iphct) end if if (isoct.ne.0) read (n5,1030) (char, i = 1, isoct) read (n5,1000) (char, i = 1, 4) read (n5,*) iv read (n5,*) (rnum, rnum, i = 1, iv) read (n5,1000) (char, i = 1 , iv) 1000 format (a10) 1010 format (10a8) 1030 format (8a10) c start reading assemblages rline = 2. iwidth = 0 20 read (n5,*,end=99) ipid,ipvar,(ipds(i), i = 1, incts) read (n5,*) (vmode(i), i = 1, incts) read (n5,*) (amode(i), i = 1, incts) read (n5,*) xx, yy read (n5,*) rho if (ivol.gt.0) read (n5,*) (wtp(i),i=1,ivol) c variance restriction if (ipvar+2.gt.ivmx.or.ipvar+2.lt.ivmn) goto 20 if (iop1.eq.1.and.(xx.gt.xmax.or.xx.lt.xmin.or. * yy.gt.ymax.or.yy.lt.ymin)) goto 20 if (iop5.eq.1.or.iop6.eq.1.or.iop7.eq.1) then ifail = 1 if (iop5.eq.1.and.incts.ge.ixct(1)) then itic = 0 do 2 j = 1, incts call checki (1,ipds(j),kex(j)) if (kex(j).ge.0) then if (isoct.gt.0.and.kex(j).ne.0. * and.j.gt.1.and.itic.gt.0) then do 1 k = 1, j-1 1 if (kex(k).eq.kex(j)) goto 2 end if itic = itic + 1 end if 2 continue if (itic.lt.ixct(1)) goto 50 c if here the assemblage contains c all the phases requested: ifail = 0 end if 50 if (iop6.eq.1.and.incts.ne.0) then do 5 j = 1, incts call checki (2,ipds(j),itis) if (itis.ge.0) then c name matches a phase in the field goto 60 end if 5 continue c no match ifail = 0 end if 60 if (iop7.eq.1.and.incts.ne.0) then do 3 j = 1, incts call checki (3,ipds(j),itis) if (itis.ge.0) then c name matches a phase in the field ifail = 0 goto 70 end if 3 continue end if 70 if (ifail.eq.1) goto 20 end if r = .5 call pselip (xx,yy,r*dcx,r*dcy,0.,0,7) call psbtxt (text,iend,ilop) call pssctr (ifont,1.2*fac1,1.2*fac1,0.0) call pstext (xx+.7*dcx,yy+.5*dcy,text,iend) goto 20 99 end c--------------------------------------------------------------- subroutine psclbl (iop1,iop5,iop6,iop7,ilop,ivmn,ivmx,fac1) c psclbl - subprogram to output bulk composition label data, without c polygon data. implicit real (a-g,o-z),integer (h-n) parameter (k0=4000,k2=6000,k5=12,k7=k5+1,k8=k5+2) real rnum character*10 xnams, text*400, char*10 integer k, kex(k8), ilop(9), int common/ rxn /vnu(k7),idr(k7),ivct,iplus(k7),iminus(k7) common/ excl1 /ixct(3),iphct,isoct,iex(50,3),jex(50,3),ict(3) common/ basic /n4,ifont,icopt,jwidth,iop0/ excl4 /xnams(50,3) common/ wsize /xmin,xmax,ymin,ymax,dcx,dcy,xlen,ylen common/ bulk /ipid,ipvar,amode(k5),vmode(k5),wtp(k5),rho, * ipds(k5),ivol,incts,pmn,dp do i = 1, 3 ict(i) = 0 end do c start reading assemblages 20 read (n4,*,end=99) ipid,ipvar,(ipds(i), i = 1, incts) read (n4,*) (vmode(i), i = 1, incts) read (n4,*) (amode(i), i = 1, incts) read (n4,*) xx, yy read (n4,*) rho if (ivol.gt.0) read (n4,*) (wtp(i),i=1,ivol) c variance restriction if (ipvar+2.gt.ivmx.or.ipvar+2.lt.ivmn) goto 20 if (iop1.eq.1.and.(xx.gt.xmax.or.xx.lt.xmin.or. * yy.gt.ymax.or.yy.lt.ymin)) goto 20 if (iop5.eq.1.or.iop6.eq.1.or.iop7.eq.1) then ifail = 1 if (iop5.eq.1.and.incts.ge.ixct(1)) then itic = 0 do 2 j = 1, incts call checki (1,ipds(j),kex(j)) if (kex(j).ge.0) then if (isoct.gt.0.and.kex(j).ne.0. * and.j.gt.1.and.itic.gt.0) then do 1 k = 1, j-1 1 if (kex(k).eq.kex(j)) goto 2 end if itic = itic + 1 end if 2 continue if (itic.lt.ixct(1)) goto 50 c if here the assemblage contains c all the phases requested: ifail = 0 ict(1) = ict(1) + 1 end if 50 if (iop6.eq.1.and.incts.ne.0) then do 5 j = 1, incts call checki (2,ipds(j),itis) if (itis.ge.0) then c name matches a phase in the field ict(2) = ict(2) + 1 goto 60 end if 5 continue c no match ifail = 0 end if 60 if (iop7.eq.1.and.incts.ne.0) then do 3 j = 1, incts call checki (3,ipds(j),itis) if (itis.ge.0) then c name matches a phase in the field ict(3) = ict(3) + 1 ifail = 0 goto 70 end if 3 continue end if 70 if (ifail.eq.1) goto 20 end if call pselip (xx,yy,0.5*dcx,0.5*dcy,0.,0,7) call psbtxt (text,iend,ilop) call pssctr (ifont,1.2*fac1,1.2*fac1,0.0) call pstext (xx+.7*dcx,yy+.5*dcy,text,iend) goto 20 99 if (iop5.eq.1) write (6,*) ict(1), * ' fields have the assemblage: ', * (xnams(k,1),' ',k = 1, ixct(1)) if (iop6.eq.1) write (6,*) ict(2), * ' fields do not have any of the phases: ', * (xnams(k,2),' ',k = 1, ixct(2)) if (iop7.eq.1) write (6,*) ict(3), * ' fields do have one of the phases: ', * (xnams(k,3),' ',k = 1, ixct(3)) end