      program xdat

!  2D data read added into pixelled array, eps printed
      implicit real*8 (a-h,o-z)
      parameter (nspec=5,ncd=10) ! elements
      parameter (ndim=1000,nzbin=1000,natt=3100)
      common /colordefs/ rgbs(nspec,ncd,4),isp(nspec),irgbs(nspec)
      character*150 title
      dimension xv(3),yv(3),zv(3),zz(2),ba(3,3),bn(3,3)
      dimension vx(3),vy(3),vz(3),rr(3),ri(3),rs(3)
      dimension rgb(3)
      dimension r(3),xx(3),xa(3),w(3),ww(3),xw(3),cce(3)
      dimension d(0:ndim,0:ndim)
      dimension r0(3,natt),ra(3,natt),dri(3),ial(natt),nma(natt)
      dimension cn(4,2),iex(3)
      dimension zh(0:nzbin),dw(3,natt),dw2(3,natt)
      dimension iv(3,1000)
      dimension u1(3),u2(3),u3(3),u4(3),ev(2,3,natt),sc(2,1:2000)
      dimension xma(natt),xmag(natt),wfac(natt)
      integer isis(1000,natt)
      integer count,nsym
      logical xb,yb,zb
      character*1 hco(2),hw(6),hws(60)
!! colormap stuff; in rgbs, "4" recolength is 3 for rgb, 4th for 
!! intensity limit
      dimension trgb(3)
      dimension ss(192,3,4)
      common/stra/ss

      psofx=0
      psofy=0
      half=0.499999
      cn(1,1)=-half
      cn(1,2)=-half
      cn(2,1)=-half
      cn(2,2)=half
      cn(3,1)=half
      cn(3,2)=-half
      cn(4,1)=half
      cn(4,2)=half
      do ko=1,3
         ww(ko)=w(ko)
      enddo
      nmax=400
      open(2,file='xperp.inp')
      read(2,*)psw
      read(2,*)xmin,xmax,ymin,ymax
      dx=xmax-xmin
      dy=ymax-ymin
      if(dx>dy) then
         pswx=psw
         pswy=psw*dy/dx
      else
         pswy=psw
         pswx=psw*dx/dy
      endif
      read(2,*)sthr
      read(2,*)satur
!      satur=1.0
      do
         read(*,*,end=3339)x1,y1
         w(2)=(x1-xmin)/dx
         w(1)=(y1-ymin)/dy
!         write(*,*)'a ',w
!         w(1)=w(1)+0.5
!         w(2)=w(2)+0.5
         w1=w(1)*nmax
         w2=w(2)*nmax
         call rround(w1,m1,dx1)
         call rround(w2,m2,dx2)
!        write(*,1151)nmax,m1,m2,dx1,dx2
 1151    format('hgt ',3i8,2f8.4)
         tarea=0d0
         do jo=1,4
            w1=w(1)*nmax+cn(jo,1)
            w2=w(2)*nmax+cn(jo,2)
            call rround(w1,k1,dk1)
            call rround(w2,k2,dk2)
!            write(*,1152)nmax,k1,k2,dk1,dk2
 1152       format('hgt2 ',3i8,2f8.4)
            dm1=(k1-m1+dk1)
            dm2=(k2-m2+dk2)
            if(dm1.lt.0)dm1=-dm1
            if(dm2.lt.0)dm2=-dm2
            area=(1-dm1)*(1-dm2)
            if(k1.lt.0.or.k1.gt.nmax)goto 9174
            if(k2.lt.0.or.k2.gt.nmax)goto 9174
 610        format('aaa ',2f9.3,2i12,2f9.3,f9.4)
            d(k1,k2)=d(k1,k2)+area
            tarea=tarea+area
            if(d(k1,k2).gt.dmx)dmx=d(k1,k2)
            mm=mm+1
         enddo
 9174    continue
      enddo
 3339 continue
!      write(*,*)'tarea',tarea
      i0=int(psofx)+1
      i1=int(psofy)+1
      iwi=int(pswx)+1
      ihe=int(pswy)+1
!     offsy(1)=psofy
!     offsy(2)=psofy+pswy+ygap
      xx0=psofx                 !leave space for ticks
      yy0=psofy                 !leave space for ticks
      xx1=psofx+pswx
      yy1=psofy+pswy
      open(3,file='perp.eps')
      write(3,1600)
      write(3,500)xx0,yy0,xx1,yy1 !bbox
      write(3,1601)
      write(3,1602)i0,i1        ! translate
      write(3,1603)iwi,ihe
      write(3,1604)
      write(3,1605)
      write(3,1606)nmax
      write(3,1607)nmax
      write(3,1608)
      write(3,1609)
      write(3,1610)nmax,nmax
      write(3,1611)
      write(3,1612)
      write(3,1613)
 500  format('%%BoundingBox: ',4f6.1)
 1600 format('%!PS-Adobe-2.0')
 1601 format('/DeviceRGB setcolorspace')
 1602 format(2i6,' translate')
 1603 format(2i6,' scale')
 1604 format('<<')
 1605 format('/ImageType 1')
 1606 format('/Width',i6)
 1607 format('/Height',i6)
 1608 format('/BitsPerComponent 8')
 1609 format('/Decode [0 1 0 1 0 1]')
 1610 format('/ImageMatrix [',i4,' 0 0 ',i4,' 0 0]')
 1611 format('/DataSource currentfile /ASCIIHexDecode filter')
 1612 format('/Interpolate true')
 1613 format('>>image')
 1614 format('>')
      zero=0d0
      dxx=(1.-sthr)/4.
      x1=dxx
      x2=2*dxx
      x3=3*dxx
      x4=4*dxx
      zero=1d-4
      ih=0
      nnw=0
      ispec=0
      do ie=1,nmax
         do ik=1,nmax
            if(ispec.gt.0) then
               rho=0
               do io=1,3
                  rgb(io)=0d0
               enddo
               denom=d(ie,ik)
               if(denom.lt.1d-6) then
                  do io=1,3
                     rgb(io)=1d0
                  enddo
                  goto 499
               endif
               ipp=0
               if(dmx.lt.1d-6)goto 19445
               xxx=d(ie,ik)/dmx/satur
               if(xxx.gt.1d0)xxx=1d0 ! saturation correction
!               fac=xxx/d(ie,ik)
               call ucolormap(xxx,io,trgb(1),ipo)   ! need to have trgb() !!!
               if(ipo.gt.0) then
                  ipp=1
                  do jo=1,3
                     rgb(jo)=rgb(jo)+trgb(jo)*d(ie,ik)/denom
!     rgb(jo)=rgb(jo)+trgb(jo)
                  enddo
                  xg=d(ie,ik)
!     write(*,*)'ffff ',io,xxx,xg,denom,(rgb(jo),jo=1,3)
               endif
19445          continue
               if(ipp.eq.0) then
                  do io=1,3
                     rgb(io)=1d0
                  enddo
               endif
            else
               denom=d(ie,ik)
               if(denom.lt.1d-6) then
                  do io=1,3
                     rgb(io)=1d0
                  enddo
                  goto 499
               endif
               rho=d(ie,ik)/dmx/satur
               if(rho.gt.1d0)rho=1d0
               xps=1.d0-rho
               if(xps.lt.zero)xps=zero
               call colormap(xps,dxx,x1,x2,x3,x4,rgb)
            endif
 499        continue
            call rgbtohex(rgb,hw)
            do io=1,6
               hws(ih+io)=hw(io)
            enddo
            ih=ih+6
            if(mod(ih,60).eq.0) then
               write(3,1248)hws
               ih=0
               nnw=nnw+1
            endif
 1248       format(60a)
            goto 57
 57         continue
         enddo
      enddo
      if(ih.gt.0) then
         write(3,1248)(hws(io),io=1,ih)
!         write(*,*)'last line has ',ih/6,' records'
      endif
      write(3,1614)
      close(3)
      stop
      end

      subroutine rround(x,i,dx)
      IMPLICIT double precision (A-H,O-Z)
      integer i
      i=x
      dx=x-float(i)
      if(dx.gt.0.5)then 
         i=i+1
         dx=x-i
      endif
      if(dx.lt.-0.5)then
         i=i-1
         dx=x-i
      endif
!      write(*,222)x,i,dx
 222  format('inround ',f10.3,i10,f10.3)
      return
      end

      subroutine colormap(q,dx,x1,x2,x3,x4,rgb)
      implicit real*8 (a-h,o-z)
      dimension rgb(*)
      if(q.lt.x1)then           ! start from red, add green
         rgb(1)=1.
         sl=q/dx
         rgb(2)=sl
         rgb(3)=0.
         goto 56
      endif
      if(q.lt.x2)then           !subtract red
         sl=(q-x1)/dx
         rgb(1)=1.-sl
         rgb(2)=1.
         rgb(3)=0.
         goto 56
      endif
      if(q.lt.x3)then           ! add blue
         rgb(1)=0.
         rgb(2)=1.
         sl=(q-x2)/dx
         rgb(3)=sl
         goto 56
      endif
      if(q.lt.x4)then           ! subtract green
         rgb(1)=0.
         sl=(q-x3)/dx
         rgb(2)=1.-sl
         rgb(3)=1.
         goto 56
      endif
      xxx=(q-x4)/(1-x4)
      rgb(1)=xxx
      rgb(2)=xxx
      rgb(3)=xxx
 56   continue
      return
      end

      subroutine ucolormap(q,is,rgb,ii)
! ii tells which interval was used
      implicit real*8 (a-h,o-z)
      parameter (nspec=5,ncd=10) ! elements
      common /colordefs/ rgbs(nspec,ncd,4),isp(nspec),irgbs(nspec)
      dimension w0(3),w1(3),rgb(3)
      rgb(1)=1d0
      rgb(2)=1d0
      rgb(3)=1d0
      eps=1d-6
      ii=0
!! there must be at least two definitions !!
      do io=2,irgbs(is)
         r0=rgbs(is,io-1,4)+eps
         r1=rgbs(is,io,4)+eps
         do jo=1,3
            w0(jo)=rgbs(is,io-1,jo)
            w1(jo)=rgbs(is,io,jo)
         enddo
         if(q.gt.r0.and.q.lt.r1.and.r1-r0.gt.eps) then
            ii=io
            dd=(q-r0)/(r1-r0)
            do jo=1,3
               rgb(jo)=w0(jo)+dd*(w1(jo)-w0(jo))
            enddo
         endif
      enddo
!      if(ii.gt.0)write(*,*)'ucol ',ii,(rgb(jo),jo=1,3)
      return
      end

      subroutine rgbtohex(rgb,hw)
      implicit real*8 (a-h,o-z)
      character*1 hw(6),hco(2)
      dimension rgb(3)
      ih=0
      do io=1,3
         icval=int(rgb(io)*255)
         if(rgb(io)*255-icval.gt.0.5)icval=icval+1
!     write(*,*)'icval ',icval
         call dtoh(icval,hco)
         ih=ih+1
         hw(ih)=hco(1)
         ih=ih+1
         hw(ih)=hco(2)
!         write(*,*)io,rgb(io),' ',icval,' ',hco(1),hco(2)
      enddo
      return
      end

!     ******************************************
!     ROUTINE TO CONVERT INTEGER VALUES TO HEX
!     WRITTEN BY:  ERIK PETER OOSTERWAL, SEP 12, 1988
      SUBROUTINE DTOH(INUM,ARRY)
      CHARACTER*1 ARRY(2),HEX(16)
      integer inum
      DATA HEX/'0','1','2','3','4','5','6','7','8',
     *     '9','a','b','c','d','e','f'/
      DO I = 2,1,-1
         ARRY(I)=HEX(MOD(INUM,16)+1)
         INUM = INUM/16
      ENDDO 
      RETURN
      END

