program rodplot
c
c  plots sets of Rodrigues triples in 7 triangles
c  output= postscript file
c
	integer type
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
	parameter (m=16)
	parameter (n=10000)
	real energy(m,m,m),correct(m,m,m)
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
c  ix 99 - not relevant here, but keep in case we want to weight each point
	real sigma(n,3)
c  SIGMA contains the sigma value of each boundary, if any
	real ap(n,3,3),bp(n,3,3)
c  AP and BP contain the A and B boundary plane normals
	integer mresolve,ndata
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
c
	character*1 ten,one,ctype,cdummy
	character*16 fname
c
	real   tx(10),ty(10)
c
      pi=4.*atan(1.)
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
c
	tx(1)=1.
	ty(1)=1.
	tx(2)=1.
	ty(2)=3.1
	tx(3)=1.
	ty(3)=5.2
	tx(4)=1.
	ty(4)=7.3
	tx(5)=3.5
	ty(5)=3.1
	tx(6)=3.5
	ty(6)=5.2
	tx(7)=3.5
	ty(7)=7.3
c  7 triangles, ll corners
c
      pi=4.*atan(1.)
c
	type=0			!  points only
	write(*,*) 'Name of input file? [enclose in quotes; try rodhemi.out]'
	read(*,*) fname
	open(unit=8,file=fname,status='old')
c
c	read (8,*) 
	do 1130, i=1,n
	do 1130, j=1,3
	read (8,*,end=1131) 
     &	idum,idum
     &	,sigma(i,j)
     &	,(disorient(i,j,mn),mn=1,4)
     &	,(ap(i,j,ia),ia=1,3)
     &	,(bp(i,j,ia),ia=1,3)
1124	format(2(x,i4),12f9.3)
	write(*,*) 'misorient, Rodr., sigma '
     &	,(disorient(i,j,mn),mn=1,4),sigma(i,j)
	ratio(i,j,1)=1.
	ratio(i,j,2)=1.
c uniform weighting for now
c
1130	continue
1131	write(*,*) 'End of input file, ',i-1,' boundaries read in.'
	ndata=i-1
c
	write(*,*) 'Do you want sigma values identified [0=NO]?'
	read(*,*) iqsigma
c
c	ten=char(48+(mresolve/10))
c	one=char(48+mod(mresolve,10))
	ctype=char(48+type)
c	write(*,*) 'Name of output file? [enclose in quotes]'
c	read(*,*) fname
	write(*,*) 'Do you want colour [=1] or B&W [=0]?'
	read(*,*) iqcolour
c
c	do 2010, ijk=1,16
c	write(*,*) fname(ijk:ijk)
c	if(fname(ijk:ijk).eq.' ') then
c		iflag=ijk
c	endif
2010	continue
	open(unit=9,file='rfp'//'.'//ctype//'.ps',status='unknown')
c
	write(9,4010) 
4010	format('%!PS-Adobe-3.0')
	write(9,4020) 
4020	format('%%BoundingBox: 36 36 599 792')
c  box for .5" margins on a 8.5x11 page
	write(9,4030) 
4030	format('%%Creator:RFplot.f')
	write(9,4040) 
4040	format('%%EndComments')
c
	write(9,5010) 
5010	format('/rod2pix {347.646753 mul} def')
c  transforms RF units (0.41 = 45 degr about 100) to 2" or 144 pixels in PS space
	write(9,5020) 
5020	format('/inch {72 mul} def')
c  transforms inches to 72 pixels in PS space
	tri_size=2.0
c  sets size of triangle
	write(9,5030) tri_size,tri_size,-1*tri_size,-1*tri_size
5030	format('/triangle    %stack: x y == ll corner'/
     &	'{ newpath moveto'/
     &	' ',f8.3,' inch 0 rlineto'/
     &	' 0 ',f8.3,' inch rlineto'/
     &	' ',f8.3,' inch ',f8.3,' inch rlineto'/
     &	' closepath } def')
c  defines a triangle of size 2 inches, for plotting within
	box_size=0.2
	write(9,5040) box_size
5040	format('/bsz ',f8.3,' def')
	write(9,5050) 
5050	format('/box    %stack: x y == corner'/
     &	'{ newpath moveto'/
     &	' bsz 2 div inch 0 rmoveto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' bsz -1 mul inch 0 rlineto'/
     &	' 0 bsz -1 mul inch rlineto'/
     &	' bsz inch 0 rlineto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' closepath } def')
c  defines a box of size 0.2 inches, for use as a symbol
	write(9,5060) 
5060	format('/rbox    %draws a box wherever we are!'/
     &	'{ '/
     &	' bsz 2 div inch 0 rmoveto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' bsz -1 mul inch 0 rlineto'/
     &	' 0 bsz -1 mul inch rlineto'/
     &	' bsz inch 0 rlineto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' closepath } def')
c  defines a box of size 0.2 inches, for use as a symbol
	WRITE(9,4050) 
4050	format('/Times-Roman findfont 18 scalefont setfont')
	write(9,4090) 
4090	format('%%EndProlog')
	rfincr=1/3./7.
c  increment of R3 between sections
	call pltpts(disorient,
     &	ratio,sigma,ndata,rfincr,tx,ty,tri_size,iqcolour,iqsigma)
	do 100, i=1,7
	rflabel=float(i-1)*rfincr
c  current R3 value
	rtmp=tri_size*rfincr
	write(9,7010) tx(i),ty(i)
7010	format(f8.3,' inch ',f8.3,' inch triangle '/
     &	'2 setlinewidth stroke ')
	if(i.eq.1) write(9,7020) tx(i)+0.2,ty(i)+1.6,rflabel
	if(i.gt.1) write(9,7021) tx(i)+0.2,ty(i)+1.6,rflabel
7020	format(f8.3,' inch ',f8.3,' inch moveto (R3 =',f6.3,') show')
7021	format(f8.3,' inch ',f8.3,' inch moveto (       ',f6.3,') show')
c  should put the label above and to the right of the ll corner
	write(9,7026) tx(i)+0.1,ty(i)+0.6,char(96+i)
7026	format(f8.3,' inch ',f8.3,' inch moveto ((',a1,')) show')
c
	write(9,7030) tx(i),ty(i)
7030	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	write(9,7035) rflabel,rflabel
7035	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point on the diagonal
	write(9,7040) tend-rflabel
7040	format(f8.3,' rod2pix 0 rlineto',
     &	     ' [] 0 setdash 2 setlinewidth stroke')
c  draw dashed line to the edge
	write(9,7050) tx(i),ty(i)
7050	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	write(9,7060) rflabel+rfincr,rflabel+rfincr
7060	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point on the diagonal
	write(9,7070) tend-rflabel-rfincr
7070	format(f8.3,' rod2pix 0 rlineto',
     &	     ' [2] 1 setdash 2 setlinewidth stroke [] 0 setdash')
c  draw solid line to the edge
c
	if(rflabel.gt.0.1715) then
c  we need the cut at the top right corner
	write(9,7110) tx(i),ty(i)
7110	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	rod2=1.-rflabel-tend
c   based on the cutoff plane being defined by X+Y+Z=1
	write(9,7115) tend,rod2
7115	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point at edge of triangle
	rod2p=(1.-rflabel)/2.
	write(9,7120) (rod2p-tend),(rod2p-rod2)
7120	format(f8.3,' rod2pix ',f8.3,' rod2pix rlineto',
     &	     ' [] 0 setdash 2 setlinewidth stroke')
c  draw dashed line to the diagonal
	write(9,7130) tx(i),ty(i)
7130	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	rod2=1.-(rflabel+rfincr)-tend
	write(9,7135) tend,rod2
7135	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point at edge of triangle
	rod2p=(1.-(rflabel+rfincr))/2.
	write(9,7140) rod2p-tend,rod2p-rod2
7140	format(f8.3,' rod2pix ',f8.3,' rod2pix rlineto',
     &	    ' [2] 1 setdash 2 setlinewidth stroke')
c  draw solid line to the diagonal
	endif
c    draws lines in upper right corner of triangle
c	  
100	continue
c
	write(9,8010) tx(1)+2.7,ty(1)+.5,type
8010	format(f8.3,' inch ',f8.3,' inch moveto (plot type = ',i3,') show')
	write(9,8020) tx(1)+2.7,ty(1)+.9,mresolve
8020	format(f8.3,' inch ',f8.3,' inch moveto (resol. = ',i4,') show')
	write(9,8030) tx(1)+2.7,ty(1)+1.3,error
8030	format(f8.3,' inch ',f8.3,' inch moveto (error = ',f10.6,') show')
	write(9,9010) 
9010	format('showpage  '/'%%PageEnd'/'%%EOF')
	close(9)
	write(*,*) 'Look for rfp.0.ps'
	call exit
	end
c
c ______________________________
c
	subroutine pltpts(angle,ratio,sigma,
     &	ndata,rfincr,tx,ty,tri_size,iqcolour,iqsigma)
	parameter (n=10000)
	real tx(10),ty(10),box_size
	real angle(n,3,5),ratio(n,3,2),sigma(n,3),rfincr,tri_size
	integer ndata,iqcolour,iqsigma
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      pi=4.*atan(1.)
      pi2=0.5*pi
      degrad=180./pi
c
	write(*,*) 'Enter the box size for plotting [0.05]'
	read(*,*) box_size
	write(9,*) '/bsz ',box_size,' def'
c  sets boxsize
c
	scale=0.1
	do 100, i=1,ndata
	do 100, j=1,3
	ntriang=1+int(0.5+angle(i,j,4)/rfincr)
	if(ntriang.lt.1) ntriang=1
	if(ntriang.gt.7) ntriang=7
c	write(*,*) 'i,j,ntriang  ',i,j,ntriang
	write(9,6010) tx(ntriang),ty(ntriang)
6010	format('newpath ',f8.3,' inch ',f8.3,' inch moveto ')
c  puts us at the ll corner of the triangle
c
	if(sigma(i,j).le.0..or.iqsigma.eq.0) then
c  plot colored box if not a CSL
c
	write(9,6020) angle(i,j,2),angle(i,j,3)
6020	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto rbox')
	tmp=alog10(sqrt(ratio(i,j,1)*ratio(i,j,2)))
c  for now, assume that we are interested in a scale of 0.5 to 2 in ratios
	if(tmp.gt.scale) tmp=scale
	if(tmp.lt.-1*scale) tmp=-1*scale
	thue=scale-tmp
	tsat=thue
c  colour (hue) ranges from about red to blue
	tbright=4.*abs(tmp)/scale
	if(tbright.gt.1.0) tbright=1.0
	if(iqcolour.eq.0) then
		write(9,6028) tsat
6028	format(' gsave ',f8.3,' setgray fill grestore ')
c  outline a square symbol, set color and grey level, fill
c  outline a square symbol,  fill with gray
c
	elseif(iqcolour.eq.1) then
		write(9,6030) thue,1.0,tbright
	endif
c  arrange for bright colors unless close to ratio=1
6030	format(' gsave ',3f8.3,' sethsbcolor fill grestore ')
c  outline a square symbol, set color and grey level, fill
c
	else
	write(9,6120) angle(i,j,2),angle(i,j,3)
6120	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto -0.1 0 rmoveto')
c	write(*,*) 'writing sigma number for ',i,j,'   sigma =',sigma(i,j)
	if(iqcolour.eq.0) then
c  set gray level by the twist angle, with a max grey of 0.8 (to stay away from white)
	write(9,6415) tsat, int(sigma(i,j))
 6415	format(f8.3,' setgray (',i2,') show ')
c 6415	format('  .01388888889 .01388888889 scale '/
c     &  f8.3,' setgray (',i2,') show '/' 72 72 scale')
	elseif(iqcolour.eq.1) then
	write(9,6420)  thue,1.0,1.0, int(sigma(i,j))
6420	format(3f8.3,' sethsbcolor (',i2,') show ')
c6420	format(' .01388888889 .01388888889 scale'/ 
c     &   3f8.3,' sethsbcolor (',i2,') show '/' 72 72 scale')
	   endif
	 endif
c
100	continue
c
	write(9,*) '/bsz 0.5 def'
c  sets boxsize to 0.5" for legend
c
	write(*,*) 'Do you want the legend? [0=NO]'
	read(*,*) iqleg
	if(iqleg.ne.0) then
	do 200, i=-4,4
	rtmp=(i+4)*0.5+0.25
	tmp=scale*float(i)/4.
	if(mod(i,2).eq.0) write(9,8010) 
     &	    tx(4)+rtmp-0.25,ty(4)+2.15,exp(tmp*2.3026)
8010	format(f8.3,' inch ',f8.3,' inch moveto (',f6.2,') show')
	write(9,8020) tx(4)+rtmp,ty(4)+2.6
8020	format(f8.3,' inch ',f8.3,' inch box')
	thue=scale-tmp
c  colour (hue) ranges from about red to blue
	tbright=4.*abs(tmp)/scale
	if(tbright.gt.1.0) tbright=1.0
	write(9,6030) thue,1.0,tbright
c  arrange for bright colors unless close to ratio=1
c  outline a square symbol, set color and grey level, fill

200	continue
	endif
	return
	end