program sod4con

c	This program reads in .SOD and .COD files and plot these
c       in Bunge-convention ODF space.  At current stage, sections
c       are made every 5 degrees (phi2 for .COD, and phi1 for .SOD).

c	SODCON written by	Paul S. Lee,	November 1998
c  much edited by a d rollett,  i 02
c  this version for plotting COD or SOD with no sample symmetry
c  edited again xii 04 to get a plot on one page

	real fint(73,20,20), vodf(73,20,20)
	real scint(0:6),scleg(0:6),scleg2(0:6)
	integer nvint(73,20,20)
	character afile*40, file*19, out1*22, outfile*30

	real array(73,19),cval(0:10),xold,yold,xdiff,ydiff,xorig,yorig
	integer nval
	real colrval(3,0:10)
c  colors are in COLRVAL
	character*5 nomen,seclab,label
	real delth,rm,delom,pm
	integer iw,jw,iper(3),iavg,ngr
	logical rescale
	integer iycount	! counter for y direct.
	character page*1	! page number, 0-5
	integer ll,ml,nl
	real data_min,data_max,data_delta,auto_val(6)
	real y_crd_fl,y_crd_sc  !  y-coords for legends
	logical l_com_line,l_ticks,l_lines
	integer q_contours
	character arg_line*20,ticks*1,lines*1

      common/conpar/ispec,ioffpp,spvall,ilegg,ilabb,nhii,ndeccn,nlbll,
     +mscall,ldsh,hgtlab

	l_com_line = .false.
	CALL GetArg(1,afile)
c	print*,'debug: afile = ',afile
	if(afile.ne.'') l_com_line = .true.
c  detect an argument for the input file name; look for user input if not present
	CALL GetArg(2,arg_line)
c	print*,'debug: arg_line = ',arg_line
	read(arg_line,"(i4)") q_contours
c	print*,'debug: q_contours = ',q_contours
	CALL GetArg(3,ticks)
c	print*,'debug: ticks = ',ticks
	l_ticks = .false.
	if(ticks.eq.'1') l_ticks = .true.
	CALL GetArg(4,lines)
c	print*,'debug: lines = ',lines
	l_lines = .false.
	if(lines.eq.'1') l_lines = .true.

	ispec=1
	ilegg=0
	ilabb=0
	nhii=-1
c turns off legend, labels, hi/lo labels
	data_min=9999999.
	data_max=0.
	if(.not.l_com_line) then
	   print *,' '
	   print *,' '
	   print *,'     ***********************************************'
	   print *,'     *                                             *'
	   print *,'     *     .COD and .SOD plotting in ODF space     *'
	   print *,'     *            by Paul Seungyong Lee            *'
	   print *,'     *            and A.D. (Tony) Rollett          *'
	   print *,'     *                                             *'
	   print *,'     *         Carnegie Mellon University          *'
	   print *,'     *                                             *'
	   print *,'     ***********************************************'
	   print *,' '
	   print *,'     for plotting triclinic sample symmetry COD '
	   
c	print *,'Name of input file (.SOD, .SMH or .COD) ?'
	   print *,'Name of input file (.CMH or .COD) ?'
	   read(*,1011) afile
 1011	   format(a)
	   print *,' '
c	else
c	   continue
c	   CALL GetArg(1,afile)
	endif

c	print *,'Choose the scaling method.'
c	print *,'    Linear ------------ 1 '
c	print *,'    Log --------------- 2 '
c	read(*,*)nscale
	nscale = 1  !  stick to linear!

	kk=40
	do 154, i=40,1,-1
c  find the right-most occurrence of a "."
	   if(afile(i:i).eq.'.') then
	   	kk=i
	   	goto 155
	   endif
 154	continue
155	continue
	   file=afile(kk+1:kk+3)
c	   if(file.eq.'sod'.or.file.eq.'SOD'.or.file.eq.'smh') then
c	      ntype=1
c	      goto 888
c	   endif
	   if(file.eq.'cod'.or.file.eq.'COD'.or.file.eq.'cmh') then
	      ntype=2
	      goto888
	   endif
c	   if(file.eq.'son'.or.file.eq.'SON'.or.file.eq.'snh') then
c	      ntype=3
c	      goto 888
c	   endif
c	   if(file.eq.'con'.or.file.eq.'CON'.or.file.eq.'cnh') then
c	      ntype=4
c	      goto888
c	   endif
	   write(*,*)'This program digests only SOD,COD,SON and CON files.'
	   goto 932

 888	open(11,file=afile,status='old')
c
c	write(*,*) 'Enter a name for the OUTPUT file, e.g. sample.cmh.ps'
c	read(*,1011) afile
c	outfile=afile(1:(kk+3))//'.ps'
c	write(*,*) 'Using ',outfile,' as output filename'
c	call newdev(outfile,3)
c	call newdev(afile,3)

	if(ntype.eq.2.or.ntype.eq.4) then
	   ff1=1
	   ff2=2
	   goto 887
	endif
	   ff1=2
	   ff2=1
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	read in the input data
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 887	continue
c
c  see below for dealing with header lines
c
c 	do 21, i=1,2
c	read(11,*)
c      read(11,"(a5)") nomen
c
c 21	continue
c
c     read the first header lines
c
	do 725,k=1,19
c
	read(11,*)
	read(11,1332) nomen,delth,rm,delom,pm,iw,jw
     &,(iper(i),i=1,3),iavg,ngr,seclab,label
 1332 format(a5,4f5.1,5i2,2i5,2a5)

      if(pm.ne.360.) then
	 write(*,*) 'The range of phi1 is not 360, but =',pm
	 write(*,*) 'therefore this is not suitable for SOD4CON'
	 print*,'If phi1 goes to 90 degrees, try SODCON or SODCONX'
	 stop
      endif

c  on the first section, detect the type of Euler angles in use
      if(k.eq.1) then
	 angle_type=1
	 if(nomen(5:5).eq.'B'.or.nomen(4:4).eq.'B') then
	    write(*,*) 'detected Bunge angles'
	 endif
c  assume Bunge angles
	 if(nomen(5:5).eq.'K'.or.nomen(4:4).eq.'K') then
	    write(*,*) 'detected Kocks angles'
	    angle_type=2
	 endif
c  Kocks angles
      endif

      scalef=100./float(iavg)
      rescale=(abs(scalef-1.0).gt.0.05)
      if(rescale) write(*,*) 'Rescaling intensities for k=',k,scalef
c  re-scale all the intensities in the section by the ratio
c  of the nominal average (100) divided by IAVG in the file
c  but require more than a 5% change to do it

      do 724, j=1,19
	 read(11,31,err=932)(nvint(i,j,k),i=1,54)
	 read(11,339,err=932)(nvint(i,j,k),i=55,73)
 31	 format(1x,18i4)
 339	 format(1x,19i4)
	 if(rescale) then
	    do 723, i=1,73
	       nvint(i,j,k)=int(float(nvint(i,j,k))*scalef)
 723	    continue
	 endif

	 do i=1,73
	    if(nvint(i,j,k).lt.data_min) data_min = nvint(i,j,k)
	    if(nvint(i,j,k).gt.data_max) data_max = nvint(i,j,k)
	 enddo

 724  continue
      read(11,*,end=725)
c  blank line after each section

 725  continue

	 print*,' minimum value = ',data_min,' max value = ',data_max
	 data_delta=(data_max-data_min)/6.
	 print*,'difference/6. = ',data_delta
	 print*,' automatic levels = '
      do i=1,6
c	 auto_val(i)=(data_min+(float(i-1)*data_delta)
c     1   	 +(data_delta*0.5))*0.0001
	 auto_val(i)=(data_min+(float(i-1)*data_delta)
     1   	 +(data_delta*0.5))*0.01
	 print*,' auto contour level ',i,' = ',auto_val(i)
      enddo

      if(.not.l_com_line) then
	 print*,'Now we deal with contour levels'
	 print*,'NB: values given here are multiples of random (uniform)'
	 print*,'  commonly written as MRD units'
	 print*,'NB: we assume popLA format files, such that the'
	 print*,'intensities are stored as i4 integers, scaled such'
	 print*,'100 == 1MRD, unless IAVG is not equal to 100'
	 print*,'in which case, inspect the code for meaning of'
	 print*,'the re-scaling operation'
	 print*
      endif
      nvmx=nvint(1,1,1)
      do 722, k=1,19
	 do 722, j=1,19
	    do 729, i=1,73
	       if(nvint(i,j,k).le.0.and.nscale.eq.2) nvint(i,j,k)=1
c  reset intensity to 1 if we are using a log scale
	       if(nvint(i,j,k).gt.nvmx) nvmx=nvint(i,j,k)
 729	    continue
 722	 continue


c	write(*,492)nvmx
c 492	format('Maximum intensity =',i7,/,
c     &         'Put the maximum intensity to be plotted.')
c	read(*,*)npmx
	npmx = nvmx + 1
	write(*,*)
	pmx=float(npmx)

c    making scale bar

	if(l_com_line) then
	   iq0 = q_contours
	else
	   write(*,*) 'Automatic choice of levels? 0=yes'
	   write(*,*) 'Enter 1 for 0.5 / 1 / 2 / 4 / 8 / 16 '
	   write(*,*) 'Enter 2 for -0.5 / -0.25 / 0 / 0.25 / 0.50 / 1.0'
	   write(*,*) 'Enter 3 for 0.04 / 0.1 / 0.3 / 0.5 / 0.7 / 0.9 '
	   write(*,*) 'Enter 4 for 0.009 /0.013 / 0.04 / 0.2 / 0.6 /0.9'
	   print*,    'Enter 5 for auto levels above, based on min-max'
	   read(*,*) iq0
	endif
	if(iq0.eq.0) then
	   do 385,i=0,5
	      xii=float(i)/5.
	      scint(i)=xii
	      scleg(i)=(1.-xii)*pmx
	      scleg2(i)=exp((1.-xii)*log(pmx))
c	write(*,*)scint(i),scleg(i),scleg2(i)
	      if(nscale.eq.1) then 
		 cval(i+1)=scleg(i)
	      else
		 cval(i+1)=scleg2(i)
	      endif
 385	   continue

	elseif(iq0.eq.1) then
	   scleg(0) = 0.5
	   cval(1) = 0.5
	   scleg(1) = 1.
	   cval(2) = 1.
	   scleg(2) = 2.
	   cval(3) = 2.
	   scleg(3) = 4.
	   cval(4) = 4.
	   scleg(4) = 8.
	   cval(5) = 8.
	   scleg(5) = 16.
	   cval(6) = 16.
	   nbelow0=1
	   nscale=1
	elseif(iq0.eq.2) then
	   scleg(0)=-0.50
	   cval(1)=-0.50
	   scleg(1)=-0.25
	   cval(2)=-0.25
	   scleg(2)=0.
	   cval(3)=0.
	   scleg(3)=0.25
	   cval(4)=0.25
	   scleg(4)=0.50
	   cval(5)=0.5
	   scleg(5)=1.
	   cval(6)=1.
	   nbelow0=2
	   nscale=1
	elseif(iq0.eq.3) then
	   scleg(0)=0.0400
	   cval(1)=0.0400
	   scleg(1)=0.1000
	   cval(2)=0.1000
	   scleg(2)=0.3000
	   cval(3)=0.3000
	   scleg(3)=0.5000
	   cval(4)=0.5000
	   scleg(4)=0.7000
	   cval(5)=0.7000
	   scleg(5)=0.9000
	   cval(6)=0.9000
	   nbelow0=1
	   nscale=1
        elseif(iq0.eq.4) then
	   scleg(0)=0.0090
	   cval(1)=0.0090
	   scleg(1)=0.0130
	   cval(2)=0.0130
	   scleg(2)=0.0400
	   cval(3)=0.0400
	   scleg(3)=0.2000
	   cval(4)=0.2000
	   scleg(4)=0.6000
	   cval(5)=0.6000
	   scleg(5)=0.9000
	   cval(6)=0.9000
	   nbelow0=1
	   nscale=1
	elseif(iq0.eq.5) then
	   do i=1,6
	      scleg(i-1)=auto_val(i)
	      cval(i)=auto_val(i)
	   enddo
	   nbelow0=1
	   nscale=1
	   
	else
	   nbelow0=0
	   nscale=1
	   write(*,*) 'Enter 6 contour values [1. = 1 MRD] '
	   do 387, i=0,5
	      write(*,*) 'Enter the value for the ',i,'th value'
	      read(*,*) scleg(i)
	      cval(i+1)=scleg(i)
	      if(scleg(i).le.1.0) nbelow0 = nbelow0+1
c     was counting contours below 100;  ADR  dec 01
 387	   continue
	endif

c    end of scale bar information

	colrval(1,0)=0.0
	colrval(2,0)=0.0
	colrval(3,0)=0.0	! black
	colrval(1,1)=0.0
	colrval(2,1)=0.0
	colrval(3,1)=1.0	! blue
	colrval(1,2)=0.0
	colrval(2,2)=1.0
	colrval(3,2)=1.0	! cyan
	colrval(1,3)=0.0
	colrval(2,3)=1.0
	colrval(3,3)=0.0	!  green
	colrval(1,4)=1.0
	colrval(2,4)=1.0
	colrval(3,4)=0.0	! yellow
	colrval(1,5)=1.0
	colrval(2,5)=0.35
	colrval(3,5)=0.0	!  orange
	colrval(1,6)=1.0
	colrval(2,6)=0.0
	colrval(3,6)=0.0	!  red
	colrval(1,7)=1.0
	colrval(2,7)=0.0
	colrval(3,7)=1.0	! magenta
c  add more values if you use more than 6 contour values!

	if(.not.l_com_line) then
	   l_ticks = .false.
	   print*,'Tick marks?  1 = YES '
	   read(*,"(a)") ticks
	   if(ticks.eq.'1') l_ticks = .true.
	   print*,'Check: l_ticks = ',l_ticks
	endif

	if(.not.l_com_line) then
	   l_ticks = .false.
	   print*,'Line contours (vs solid)?  1 = YES '
	   read(*,"(a)") lines
	   if(lines.eq.'1') l_lines = .true.
	endif

c  changed the scaling to have 1.0 = 1 MRD
	do 149,nk=1,19
	   do nj=1,19
	      do ni=1,73
		 if(nvint(ni,nj,nk).gt.npmx) then
		   nvint(ni,nj,nk)=npmx
		 endif
		 if(nscale.eq.2) then
		    fint(ni,nj,nk)=log(float(nvint(ni,nj,nk))*0.01)
c	  write(*,*)fint(ni,nj,nk)
		 else
		    fint(ni,nj,nk)=float(nvint(ni,nj,nk))*0.01
c	  write(*,*) fint(ni,nj,nk)
		 endif
	      enddo
	   enddo
 149    continue

	write(*,32) afile
 32	format('Status : Data read in from ',a)

c       end of data input

c  new loop for plotting ALL the sections, but first just
c  a one-page selection

	do 4000, kmj=1,5

	   page=char(48+kmj-1)
	   outfile=afile(1:(kk+3))//'.p'//page//'.ps'
	   write(*,*) 'Using ',outfile,' as output filename'
	   call newdev(outfile,(kk+3+2+1+3))

	   call psinit(.true.)

c     initiating post script plotting

	   write(*,239)
 239	   format('Status : Making labels')

	   red=0.
	   green=0.
	   blue=0.
	   call setcolr(red,green,blue)
	   call setlw(.035)

	   call arrow(0.1,10.0,1.1,10.0,.08,20.,1)
	   call arrow(0.1,10.0,0.1,9.0,.08,20.,1)

c     call arrow(5.7,2.4,6.4,2.4,.08,20.,1)
c     call arrow(5.7,2.4,5.7,1.7,.08,20.,1)
	   if(angle_type.eq.1) then 
	      call grksym(0.3,9.3,.15,21,0.,1,1)
	      call grksym(0.55,9.75,.15,45,0.,1,1)
	      call keknum(0.65,9.65,.09,ff1,0.,-1,1)
	   endif
	   if(angle_type.eq.2) then
	      call grksym(0.2,9.3,.15,8,0.,1,1)
	      if(ntype.eq.2) call grksym(0.35,9.3,.15,23,0.,1,1)
	      if(ntype.eq.1) call grksym(0.35,9.3,.15,45,0.,1,1)
	   endif
	   call keknum(0.15,8.7,.15,90.,0.,-1,1)
	   call keknum(1.5,9.85,.15,360.,0.,-1,1)

	   y_crd_fl = 0.15
	   call setfnt(31)
	   call keksym(0.9,y_crd_fl,.14,10hInput file,0.,10,0)
c     call keksym(0.7,0.20,.10,19hMask size (microns),0.,19,0)
c     call keksym(0.7,0.05,.10,20hPixel size (degrees),0.,20,0)
c     call keksym(0.7,-0.1,.10,15hSample symmetry,0.,15,0)
c     call keksym(0.7,-0.1,.1,21hOrientation tolerance,0.,21,0)

	   call keksym(1.7,y_crd_fl,.14,1h:,0.,1,0)
c     call keksym(2.4,0.20,.10,1h=,0.,1,0)
c     call keksym(2.4,0.05,.10,1h=,0.,1,0)
c     call keksym(2.4,-0.1,.10,1h=,0.,1,0)
c     call keksym(2.4,-0.25,.10,1h=,0.,1,0)

c     call keknum(2.7,0.10,.10,dia,0.,-1,0)
c     call keknum(2.7,0.05,.10,resol,0.,-1,0)
c     call keknum(2.7,-0.1,.10,fsymm,0.,-1,0)
c     call keknum(2.7,-0.1,.1,atol,0.,-1,0)
	   call keksym(2.5,y_crd_fl,.14,afile,0.,20,0)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     
c     making scale list
c     
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	   y_crd_sc = 0.45
	   call keksym(0.9,y_crd_sc,.14,12hContours at ,0.,12,0)
	   red=0.
	   green=0.
	   blue=1.
	   call setcolr(red,green,blue)
	   do 386, i=0,5
	      if(i.gt.1) then
		 red=0.
		 green=0.
		 blue=0.
		 call setcolr(red,green,blue)

	      endif
	      if(i.ge.4) then
		 red=1.
		 green=0.
		 blue=0.
		 call setcolr(red,green,blue)
	      endif
	      scx = 2.5+float(i)/10.*6.5
c     call rectfilg(scx,.31,scx+.7,.31,.3,scint(i))
	      if(nscale.eq.1) then
		 call keknum(scx+.1,y_crd_sc,.14,scleg(i),0.,3,0)
	      endif
	      if(nscale.eq.2) then
		 call keknum(scx+.1,y_crd_sc,.14,scleg2(i),0.,3,0)
	      endif
 386	   continue
	   red=0.
	   green=0.
	   blue=0.
	   call setcolr(red,green,blue)

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     
c     start of 3 orientation loops (phi1,phi,phi2)
c     phi2 is multiples of 5.
c     phi1 and phi are defined by the resolution of VODF space.
c     
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	   dx=2.5/90.*1.5
c  size of the boxes
	   write(*,39)
 39	   format('Status : Plugging in data ')

	   xdiff=0.
	   ydiff=0.
c     do 26,k=1,19
	   icount=0
c  now we set the increments and limits
	   if(kmj.eq.1) then  !  =1 means first page
	      ll = 1
	      ml = 19
	      nl = 1
	      xsize = 0.70
	      ysize = xsize
	      xold = 0.
	      yold = 0.
	   else
	      ll=(kmj-2)*5+1
	      ml=(kmj-2)*5+5
	      if(ml.gt.19) ml=19
	      nl=1
	      xsize = 1.4
	      ysize = xsize
c	      xold = -1.
c	      yold = -1.
	      xold = 0.
	      yold = 0.
	   endif
	   do 26, k=ll,ml,nl  !  each section 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     
c     vodfmx is the maximum vodf.  vodf will be normalized by vodfmx
c     because gray scale should be in the range of 0 and 1.
c     
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	      nsec=k-1

c     nsec is the number of ODF section from top left to bottom right.
c     Its range is from 0 to 18.

c     nys=5-nsec/4
c     nxs=nsec+1
c     if(nxs.gt.16) nxs=nxs-16
c     if(nxs.gt.12) nxs=nxs-12
c     if(nxs.gt.8) nxs=nxs-8
c     if(nxs.gt.4) nxs=nxs-4

	      if(kmj.eq.1) then	!  =1 means first page
		 nxs = nsec/10 + 1
		 nys = 10 - mod(nsec,10)  ! 10 sections per column
		 xcen = 1.2 + (nxs-1)*5.7*xsize
		 ycen = 0.6 + float(nys)*(ysize+0.05)
	      else  !  the other pages
		 nxs=1
		 nys=5-icount
		 icount=icount+1
		 xcen = 1.2 + nxs*xsize
		 ycen = 0.05 + (nys)*(ysize+0.4)
	      endif
	      print*,'nsec,nxs, nys = ',nxs,nys

c     nxs is the x-section number from left.
c     nys is the y-section number from top.

c	      xsize=1.2
c	      ysize=1.2  !  see above - vary by page
c	      xcen=0.2+nxs*xsize
c	      ycen = (nys)*(ysize+0.4)
	      xsta = xcen-(xsize/2.)
c	      ysta=ycen+(ysize/2.)
	      ysta = ycen - (ysize/2.)
	      xdiff = xsta-xold
c	      ydiff=(ysta-ysize)-yold
	      ydiff = ysta-yold
	      xold = xsta
c	      yold=ysta-ysize
	      yold = ysta
	      write(*,*)'x- y-centre= ',xcen,ycen
	      print*,'   x-, y-starting points = ',xsta,ysta
	      print*,'        xdiff, ydiff = ',xdiff,ydiff
	      call plot(xdiff,ydiff,-3)

c     xcen, and ycen are the coordinates of centers of ODF sections.
c     These are the same locations as the centers in ODF frames.
c     xsta, and ysta are the coordinates of starting point of each
c     ODF section.

	      red=0.
	      green=0.
	      blue=0.
	      call setcolr(red,green,blue)
	      call setlw(.025)

	      xss = -0.9
	      yss = 0.2
	      if(angle_type.eq.1) then 
		 call grksym(xss,yss,.15,45,0.,1,1)
		 call keknum(xss+.12,yss-.07,.09,ff2,0.,-1,1)
	      elseif(angle_type.eq.2) then
		 if(ntype.eq.2) call grksym(xss,yss,.15,45,0.,1,1)
		 if(ntype.eq.1) call grksym(xss,yss,.15,23,0.,1,1)
	      endif
	      xsss = xss + 0.23
	      call keksym(xsss,yss,.15,3h  = ,0.,3,1)
c	      ky = 6-j
c     fnum=float((i+(ky-1)*4)*5-5)
c	      fnum=float((ll-1+((5-j))*nl)*5)
	      fnum = 5.*float(k-1)
	      call keknum(xsss+.4,yss-0.02,.15,fnum,0.,-1,1)

c     do 1020, ijk=1,19
	      do 1020, ijk=1,73
		 do  jkl=1,19
		    array(ijk,jkl)=fint(ijk,(20-jkl),k)
c     have to go from bottom to top to obtain the correct order
		 enddo
 1020	      continue
c     copy section into array for contouring

	      CALL SETLW (.005)
	     if(l_ticks) then
		call border((4.*xsize),ysize,1111,1111,4,9,1,9)
	     else
		call border((4.*xsize),ysize,0000,1111,4,9,1,9)
	     endif
c     if(nbelow0.gt.0) then

	   if(l_lines) then

c	      ldsh=6
	      ldsh = 0  ! solid lines for now
	      red=0.
	      green=0.
	      blue = 0.8
	      call setcolr(red,green,blue)
	      CALL SETLW (.015)
	      call conrec(array,73,73,19,(4.*xsize),ysize,cval(1),2)
c     endif  !  nbelow0 above
	      ldsh=0
	      red=0.
	      green=0.
	      blue=0.
	      call setcolr(red,green,blue)
	      CALL SETLW (.02)
	      call conrec(array,73,73,19,(4.*xsize),ysize,cval(3),2)

	      ldsh=0
	      red=1.
	      green=0.
	      blue=0.
	      call setcolr(red,green,blue)
	      CALL SETLW (.03)
	      call conrec(array,73,73,19,(4.*xsize),ysize,cval(5),2)

	   else
	      call concolr(array,73,73,19,(4.*xsize),ysize,cval(0)
     +           ,colrval,7,0,0.)
	      red=0.
	      green=0.
	      blue=0.
	      call setcolr(red,green,blue)
	      CALL SETLW (.01)
c	      call conrec(array,73,73,19,(4.*xsize),ysize,cval(0),6)
	   endif

 26	   continue

c	   xorig = -1.*xsta
c	   yorig = -2.*(ysta-ysize)
c	   call plot(xorig,yorig,-3)
c	   print*, 'resetting origin to ',xorig,yorig

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     
c     frame setting
c     
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	   call plotnd

 4000	continue

	print*,'plotting completed without errors'
	call exit

 932	write(*,*)'Input file is not in the right format.'
	stop
	end

	include '/code/paul.lee.codes/psplot.txt'