FFLAGS = -O -Zf -Zg 
LIBS = -lrcgos -lmlib -lutil -lgrafixf -lvdif

OBJECT = arcdump.o colmanw.o initialize.o makemap.o selection.o

colenhanc: $(OBJECT)
	f77 $(DEBUG) $(FFLAGS) -o colenhanc $(OBJECT) $(LIBS)
	mv colenhanc ../.

clean:
	rm -f *.o made

.f.o:
	f77 $(DEBUG) $(FFLAGS) -c $*.f
      subroutine arcdump(fullscreen)
      integer*2 colormax,colormin,datamax,datamin,idim,jdim
      integer*2 localmax,localmin,numpix,totcolor
      integer*2 cmap(0:1023,3),pixels(1024)
      integer*2 mandelset(1:2049,1:2049)
      integer*4 mandelcolor,mapfunc,paint,whole
      integer*4 i,lstpix(1024),green,temp,y
      integer*4 colman,decimal,max,min,xwinlen,xwinorg,ywinlen,ywinorg
      logical fullscreen,gray,interp
      real datascale,power
      common/colormap/ datascale,decimal,paint,gray,interp,mapfunc,
     .                 power,whole
      common/extreme/ colormax,colormin,datamax,datamin,localmax,
     .                localmin,max,min
      common/image/ mandelset
      common/window/ colman,idim,jdim,xwinlen,xwinorg,ywinlen,ywinorg
      call inigra
      call opegra('test.ras')
      call rasdim(1024,768)
      call rasori(1,1)
      call bacras(0.0,0.0,0.0)
      totcolor=1023
      do 10 i=0,totcolor
         call getmco(i,cmap(i,1),cmap(i,2),cmap(i,3))
   10 continue
      do 20 y=0,767
         if(fullscreen) then
            call cmov2i(0,y)
            numpix=readpi(1024,pixels)
         endif
         do 30 i=1,1024
            if(fullscreen) then
               mandelcolor=pixels(i)
            else
               mandelcolor=nint((mandelset(i,y+1)-min)/datascale)+32
               if(mandelcolor.lt.32) mandelcolor=32
               if(mandelcolor.gt.1023) mandelcolor=1023
            endif
c red
            temp=cmap(mandelcolor,1)
            temp=ishft(temp,16)
c green
            green=cmap(mandelcolor,2)
            green=ishft(green,8)
            temp=ior(temp,green)
c blue
            temp=ior(temp,cmap(mandelcolor,3))
            lstpix(i)=temp
   30    continue
         call raslin(1,y+1,1024,lstpix)
   20 continue
      call fragra
      call clogra
      call endgra
      return
      end
      program colmanw
      character file*80,extrlabel*12,extrlist(2)*12,funclabel(5)*12
      character label*132,yes*1
      integer*1 onebyte,mandelset1(1:2049)
      integer*2 ax,ay,bx,by,colormax,colormin,datamax,datamin
      integer*2 diff,eatnum,i,idim,ioff,itemp,j,jdim,joff
      integer*2 localmax,localmin,maxband,oldskip,one,skip,start
      integer*2 tempi,tempj,val
      integer*2 mandelset(1:2049,1:2049)
      integer*2 numpix,totcolor,yi
      integer*4 bottom,colman,decimal,dev,gcolor(0:1,0:1),index
      integer*4 left,mandelcolor,mapfunc,paint,oldpaint,whole
      integer*4 max,min,right,top
      integer*4 xwinlen,xwinorg,ywinlen,ywinorg
      logical byte,dimentioninfo,eat,gray,gtoggle,interp
      logical scaleinfo,test,temp,topbottom
      real colorscale,datascale,power
      real scale,scalex,scaley
      real xorigin,xscale,y,yorigin,yscale
      real spectrum(32,3)
      common/colorband/ spectrum,maxband
      common/colormap/ datascale,decimal,paint,gray,interp,mapfunc,
     .                 power,whole
      common/extreme/ colormax,colormin,datamax,datamin,localmax,
     .                localmin,max,min
      common/image/ mandelset
      common/labels/ extrlabel,extrlist,funclabel
      common/plot/ bottom,left,right,skip,top,xorigin,xscale,
     .             yorigin,yscale
      common/window/ colman,idim,jdim,xwinlen,xwinorg,ywinlen,ywinorg
$ include /usr/include/fgl.h
$ include /usr/include/fdevice.h
c
c initialize some parameters
c
      data extrlist/'Interpolated','Truncated'/
      data funclabel/'  Log(V)^P','  V^(1/P)','    V^P','Exp(V*2^P)-1',
     .               'Cos(V)^P+1'/
      colorscale=741.0/991.0
      decimal=0
      extrlabel=extrlist(1)
      paint=0
      oldpaint=paint
      gray=.false.
      interp=.true.
      maxband=8
      mapfunc=32
      oldskip=1
      one=1
      power=whole+real(decimal)/10.0
      skip=oldskip
c initial color bands
c white
      spectrum(1,1)=255
      spectrum(1,2)=255
      spectrum(1,3)=255
c magenta
      spectrum(2,1)=255
      spectrum(2,2)=0
      spectrum(2,3)=255
c red
      spectrum(3,1)=255
      spectrum(3,2)=0
      spectrum(3,3)=0
c yellow
      spectrum(4,1)=255
      spectrum(4,2)=255
      spectrum(4,3)=0
c green
      spectrum(5,1)=0
      spectrum(5,2)=255
      spectrum(5,3)=0
c cyan
      spectrum(6,1)=0
      spectrum(6,2)=255
      spectrum(6,3)=255
c blue
      spectrum(7,1)=0
      spectrum(7,2)=0
      spectrum(7,3)=255
c black
      spectrum(8,1)=0
      spectrum(8,2)=0
      spectrum(8,3)=0
      whole=1
      write(*,"('Enter file name:',$)")
      read(*,'(a80)') file
      open(2,file=file,form='binary',status='old')
      rewind 2
      write(*,"('Is there dimention information? ',$)")
      read(*,'(a1)') yes
      dimentioninfo=.true.
      if(yes.eq.'n'.or.yes.eq.'N') dimentioninfo=.false.
      write(*,"('Is there scale information? ',$)")
      read(*,'(a1)') yes
      scaleinfo=.true.
      if(yes.eq.'n'.or.yes.eq.'N') scaleinfo=.false.
      write(*,"('Is this byte information? ',$)")
      read(*,'(a1)') yes
      byte=.true.
      if(yes.eq.'n'.or.yes.eq.'N') byte=.false.
      write(*,"('Eat bytes? ',$)")
      read(*,'(a1)') yes
      eat=.false.
      if(yes.eq.'y'.or.yes.eq.'Y') eat=.true.
      if(eat) then
         write(*,"('Eat how many bytes of information? ',$)")
         read(*,*) eatnum
      endif
      write(*,"('Reverse top and bottom? ',$)")
      read(*,'(a1)') yes
      topbottom=.false.
      if(yes.eq.'y'.or.yes.eq.'Y') topbottom=.true.
      if(dimentioninfo) then
         if(byte) then
            read(2) idim1,jdim1
            idim=idim1
            jdim=jdim1
         else
            read(2) idim,jdim
         endif
         write(*,"('From input file: idim=',i5,' jdim=',i5)") idim,jdim
      else
         write(*,"('Enter idim,jdim:',$)")
         read(*,*) idim,jdim
      endif
      write(*,"('Enter left,right,bottom,top:',$)")
      read(*,*) left,right,bottom,top
      if(eat) read(2)(onebyte,i=1,eatnum)
      if(scaleinfo) then
         read(2) xorigin,yorigin,xscale,yscale
      else
         xorigin=0.0
         yorigin=0.0
         xscale=1.0
         yscale=1.0
      endif
      if(byte) then
         if(topbottom) then
            do 10 j=jdim,1,-1
               read(2) (mandelset1(i),i=1,idim)
               do 10 i=1,idim
                  if(mandelset1(i).lt.0) then
                     mandelset(i,j)=mandelset1(i)+256
                  else
                     mandelset(i,j)=mandelset1(i)
                  endif
   10       continue
         else
            do 20 j=1,jdim
               read(2) (mandelset1(i),i=1,idim)
               do 20 i=1,idim
                  if(mandelset1(i).lt.0) then
                     mandelset(i,j)=mandelset1(i)+256
                  else
                     mandelset(i,j)=mandelset1(i)
                  endif
   20       continue
         endif
      else
         if(topbottom) then
            read(2) ((mandelset(i,j),i=1,idim),j=jdim,1,-1)
         else
            read(2) ((mandelset(i,j),i=1,idim),j=1,jdim)
         endif
      endif
      call ringbe
      do 30 j=1,jdim
         do 30 i=1,idim
            itemp=mandelset(i,j)
            if(i.eq.1.and.j.eq.1) then
               datamax=itemp
               datamin=itemp
            else if(itemp.gt.datamax) then
               datamax=itemp
            else if(itemp.lt.datamin) then
               datamin=itemp
            endif
   30 continue
      localmax=datamax
      localmin=datamin
      max=datamax
      min=datamin
      call ringbe
      close(2)
      call initialize
      call setval(mousex,xwinorg+xwinlen/2,0,1023)
      call setval(mousey,ywinorg+ywinlen/2,0,767)
      call selection
      oldskip=skip
      oldpaint=paint
   50 continue
      call cursof
      call qreset
      call getori(xwinorg,ywinorg)
      call getsiz(xwinlen,ywinlen)
      call viewpo(0,xwinlen-1,0,ywinlen-1)
      call ortho2(0.0,1023.0,0.0,767.0)
      call lookat(0.0,0.0,1.0,0.0,0.0,0.0,0.0)
C     call color(black)
C     call clear
c
c Color scale
c
      do 60 index=32,1023
         y=anint(real(index-32)*colorscale+14.0)
         call color(index)
         call pmv2(988.0,y)
         call pdr2(1023.0,y)
         call pdr2(1023.0,y+colorscale)
         call pdr2(988.0,y+colorscale)
         call pclos
   60 continue
      call color(white)
      write(label,1010) colormin
      call cmov2i(988,2)
      call charst(label,4)
      write(label,1010) colormax
      call cmov2i(988,756)
      call charst(label,4)
c
c Draw data
c
      if(real(jdim)/real(idim).ge.0.75) then
         call viewpo(0,nint(0.75*real(idim)/real(jdim)*988*xwinlen/
     .               1024)-1,28*ywinlen/768,ywinlen-1)
      else
         call viewpo(0,988*xwinlen/1024-1,28*ywinlen/768,
     .               nint(real(jdim)/real(idim)/0.75*ywinlen)-1)
      endif
      call ortho2(real(left),real(right),real(bottom),real(top))
      call lookat(0.0,0.0,1.0,0.0,0.0,0.0,0.0)
      test=.true.
      if(paint.eq.0) then
         call color(black)
         call clear
         do 70 j=bottom,top,skip
            do 80 i=left,right,skip
               itemp=mandelset(i,j)
               if(test) then
                  localmax=itemp
                  localmin=itemp
                  test=.false.
               else if(itemp.gt.localmax) then
                  localmax=itemp
               else if(itemp.lt.localmin) then
                  localmin=itemp
               endif
               mandelcolor=nint((itemp-min)/datascale)+32
               if(mandelcolor.lt.32) mandelcolor=32
               if(mandelcolor.gt.1023) mandelcolor=1023
               call color(mandelcolor)
               call pnt2s(i,j)
   80       continue
   70    continue
      else if(paint.eq.1) then
         do 90 j=bottom,top,skip
            tempj=j+one*skip
            do 100 i=left,right,skip
               tempi=i+one*skip
               itemp=mandelset(i,j)
               if(test) then
                  localmax=itemp
                  localmin=itemp
                  test=.false.
               else if(itemp.gt.localmax) then
                  localmax=itemp
               else if(itemp.lt.localmin) then
                  localmin=itemp
               endif
               mandelcolor=nint((itemp-min)/datascale)+32
               if(mandelcolor.lt.32) mandelcolor=32
               if(mandelcolor.gt.1023) mandelcolor=1023
               call color(mandelcolor)
c              call pnt2s(i,j)
               call pmv2s(i,j)
               call pdr2s(tempi,j)
               call pdr2s(tempi,tempj)
               call pdr2s(i,tempj)
               call pclos
  100       continue
   90    continue
      else if(paint.eq.2) then
         do 110 j=bottom,top-skip,skip
            tempj=j+one*skip
            do 120 i=left,right-skip,skip
               tempi=i+one*skip
c find local min/max
               itemp=mandelset(i,j)
               if(test) then
                  localmax=itemp
                  localmin=itemp
                  test=.false.
               else if(itemp.gt.localmax) then
                  localmax=itemp
               else if(itemp.lt.localmin) then
                  localmin=itemp
               endif
c use last colors calculated
               if(i.eq.left) then
                  start=0
               else
                  start=1
                  gcolor(0,0)=gcolor(1,0)
                  gcolor(0,1)=gcolor(1,1)
               endif
c calculate colors
               do 130 ioff=start,1
                  do 140 joff=0,1
                     itemp=mandelset(i+ioff*skip,j+joff*skip)
                     mandelcolor=nint((itemp-min)/datascale)+32
                     if(mandelcolor.lt.32) mandelcolor=32
                     if(mandelcolor.gt.1023) mandelcolor=1023
                     gcolor(ioff,joff)=mandelcolor
  140             continue
  130          continue
               call setsha(gcolor(0,0))
               call pmv2s(i,j)
               call setsha(gcolor(1,0))
               call pdr2s(tempi,j)
               call setsha(gcolor(1,1))
               call pdr2s(tempi,tempj)
               call setsha(gcolor(0,1))
               call pdr2s(i,tempj)
               call spclos
  120       continue
  110    continue
      endif
      call viewpo(0,xwinlen-1,0,ywinlen-1)
      call ortho2(0.0,1023.0,0.0,767.0)
      call lookat(0.0,0.0,1.0,0.0,0.0,0.0,0.0)
      call color(white)
      write(label,1000) real(left-1)*xscale+xorigin,
     .                  real(bottom-1)*yscale+yorigin,
     .                  real(right-1)*xscale+xorigin,
     .                  real(top-1)*yscale+yorigin,
     .                  funclabel(mapfunc-29),power,extrlabel,
     .                  localmax,localmin
      call cmov2i(0,7*ywinlen/768)
      call charst(label,109)
  150 if(qtest().ne.0) then
         dev=qread(val)
         if(dev.eq.qkey) then
            goto 900
         else if(dev.eq.ckey) then
            call curson
            call setcur(1,1,1023)
         else if(dev.eq.dkey.and.val.gt.0) then
	    call ringbe
	    call arcdump(.true.)
	    call ringbe
         else if(dev.eq.rkey.and.val.gt.0) then
	    call ringbe
	    call arcdump(.false.)
	    call ringbe
         else if(dev.eq.rightm.and.val.gt.0) then
            call selection
            if(oldskip.ne.skip) then
               oldskip=skip
               goto 50
            endif
            if(oldpaint.ne.paint) then
               oldpaint=paint
               goto 50
            endif
         else if(dev.eq.leftmo) then
            call curson
            if(val.gt.0) then
               dev=qread(bx)
               bx=bx-xwinorg
               dev=qread(by)
               by=by-ywinorg-nint(28.0*ywinlen/768.0)
               call setcur(2,1,1023)
            else
               dev=qread(ax)
               ax=ax-xwinorg
               dev=qread(ay)
               ay=ay-ywinorg-nint(28.0*ywinlen/768.0)
               if(ax.lt.bx) then
                  itemp=ax
                  ax=bx
                  bx=itemp
               endif
               if(ay.lt.by) then
                  itemp=ay
                  ay=by
                  by=itemp
               endif
               scalex=real(right-left)/(988.0*xwinlen/1024.0-1.0)
               scaley=real(top-bottom)/(740.0*ywinlen/768.0-1.0)
               if((ax-bx)*.75*scalex.gt.(ay-by)*scaley) then
                  diff=ax-bx
                  scale=scalex
               else
                  diff=nint((ay-by)/0.75)
                  scale=scaley
               endif
	       if(diff.eq.0) diff=1
               left=left+nint(scalex*real(bx))
               if(left.lt.1) left=1
               bottom=bottom+nint(scaley*real(by))
               if(bottom.lt.1) bottom=1
               right=left+nint(scale*real(diff))
               top=bottom+nint(scale*real(diff)*0.75)
               if(right.gt.idim) right=idim
               if(top.gt.jdim) top=jdim
               call qreset
               goto 50
            endif
         else if(dev.eq.middle) then
            call curson
            if(val.gt.0) then
               dev=qread(bx)
               bx=bx-xwinorg
               dev=qread(by)
               by=by-ywinorg-nint(28.0*ywinlen/768.0)
               call setcur(2,1,1023)
            else
               dev=qread(ax)
               ax=ax-xwinorg
               dev=qread(ay)
               ay=ay-ywinorg-nint(28.0*ywinlen/768.0)
               if(ax.lt.bx) then
                  itemp=ax
                  ax=bx
                  bx=itemp
               endif
               if(ay.lt.by) then
                  itemp=ay
                  ay=by
                  by=itemp
               endif
               scalex=real(right-left)/(988.0*xwinlen/1024.0-1.0)
               scaley=real(top-bottom)/(740.0*ywinlen/768.0-1.0)
               if((ax-bx)*.75*scalex.gt.(ay-by)*scaley) then
                  diff=ax-bx
		  if(diff.eq.0) diff=1
                  scale=real(right-left)/real(diff)
               else
                  diff=nint((ay-by)/0.75)
		  if(diff.eq.0) diff=1
                  scale=real(top-bottom)/real(diff)
               endif
               left=left-nint(scale*real(bx))
               if(left.lt.1) left=1
               bottom=bottom-nint(scale*real(by))
               if(bottom.lt.1) bottom=1
               right=left+nint(scale*(988.0*xwinlen/1024.0-1.0))
               top=bottom+nint(scale*(ywinlen-28.0*ywinlen/768.0))
               if(right.gt.idim) right=idim
               if(top.gt.jdim) top=jdim
               call qreset
               goto 50
            endif
         endif
      endif
      goto 150
  900 call color(black)
      call clear
      call setcur(0,1,1023)
      print *,real(left-1)*xscale+xorigin,real(bottom-1)*yscale+yorigin,
     .        real(right-1)*xscale+xorigin,real(top-1)*yscale+yorigin
      print *,xwinorg,xwinlen,ywinorg,ywinlen
      stop
 1000 format('l-b=(',e13.7,',',e13.7,') r-t=(',e13.7,',',e13.7,') ',
     .       a12,' P=',f4.1,1x,a5,' sx=',i4,' sn=',i4)
 1010 format(i4)
      end
      subroutine initialize
      integer*2 curs(16,2),idim,jdim
      integer*4 colour,decimals,extr,func,funcextr,funcmap,menu,skipnum
      integer*4 numbers,range,painting
      integer*4 colman,colorbar,oldwin,xwinlen,xwinorg,ywinlen,ywinorg
      common/menus/ colour,colorbar,decimals,extr,func,menu,numbers,
     .              range,skipnum,painting
      common/window/ colman,idim,jdim,xwinlen,xwinorg,ywinlen,ywinorg
$ include /usr/include/fgl.h
$ include /usr/include/fdevice.h
      curs(1,1)=65535
      curs(2,1)=64512
      curs(3,1)=63488
      curs(4,1)=63488
      curs(5,1)=64512
      curs(6,1)=56832
      curs(7,1)=36608
      curs(8,1)=34688
      curs(9,1)=33728
      curs(10,1)=33248
      curs(11,1)=33008
      curs(12,1)=32888
      curs(13,1)=32828
      curs(14,1)=32798
      curs(15,1)=32782
      curs(16,1)=32772
      curs(1,2)=8193
      curs(2,2)=28673
      curs(3,2)=30721
      curs(4,2)=15361
      curs(5,2)=7681
      curs(6,2)=3841
      curs(7,2)=1921
      curs(8,2)=961
      curs(9,2)=481
      curs(10,2)=241
      curs(11,2)=123
      curs(12,2)=63
      curs(13,2)=31
      curs(14,2)=31
      curs(15,2)=63
      curs(16,2)=65535
      call foregr
      call keepas(1024,768)
      colman=winope('colmanw',7)
      oldwin=winatt()
      call double
      call gconfi
      call getori(xwinorg,ywinorg)
      call getsiz(xwinlen,ywinlen)
      call frontb(.true.)
      call defcur(1,curs(1,1))
      call curori(1,0,0)
      call defcur(2,curs(1,2))
      call curori(2,16,16)
      call setcur(0,1,1023)
      call qdevic(ckey)
      call qdevic(dkey)
      call qdevic(qkey)
      call qdevic(rkey)
      call qdevic(leftmo)
      call qdevic(middle)
      call qdevic(rightm)
      call tie(leftmo,mousex,mousey)
      call tie(middle,mousex,mousey)
      colorbar=newpup()
      call addtop(colorbar,"White %x50",10,1)
      call addtop(colorbar,"Magenta %x51",12,1)
      call addtop(colorbar,"Red %x52",8,1)
      call addtop(colorbar,"Yellow %x53",11,1)
      call addtop(colorbar,"Green %x54",10,1)
      call addtop(colorbar,"Cyan %x55",9,1)
      call addtop(colorbar,"Blue %x56",9,1)
      call addtop(colorbar,"Black %x57",10,1)
      call addtop(colorbar,"Exit %x69",9,1)
      painting=newpup()
      call addtop(painting,"Painting %t",11,1)
      call addtop(painting,"Points %x70",11,1)
      call addtop(painting,"Solid %x71",10,1)
      call addtop(painting,"Gouraud %x72",12,1)
      colour=newpup()
      call addtop(colour,"Color Map %t",12,1)
      call addtop(colour,"Colors %x10",11,1)
      call addtop(colour,"Gray Scale %x11",15,1)
      extr=newpup()
      call addtop(extr,"Handling Extremes %t",20,1)
      call addtop(extr,"Interpolation %x20",18,1)
      call addtop(extr,"Truncation %x21",15,1)
      func=newpup()
      call addtop(func,"Color Map Functions %t",22,1)
      call addtop(func,"Log(value)^P %x30",17,1)
      call addtop(func,"value^(1/P) %x31",16,1)
      call addtop(func,"value^P %x32",12,1)
      call addtop(func,"Exp(value*2^P)-1 %x33",21,1)
      call addtop(func,"Cos(value)^P+1 %x34",19,1)
      decimals=newpup()
      call addtop(decimals,"-9%x-109",8,1)
      call addtop(decimals,"-8%x-108",8,1)
      call addtop(decimals,"-7%x-107",8,1)
      call addtop(decimals,"-6%x-106",8,1)
      call addtop(decimals,"-5%x-105",8,1)
      call addtop(decimals,"-4%x-104",8,1)
      call addtop(decimals,"-3%x-103",8,1)
      call addtop(decimals,"-2%x-102",8,1)
      call addtop(decimals,"-1%x-101",8,1)
      call addtop(decimals,"0%x100",6,1)
      call addtop(decimals," 1%x101",7,1)
      call addtop(decimals," 2%x102",7,1)
      call addtop(decimals," 3%x103",7,1)
      call addtop(decimals," 4%x104",7,1)
      call addtop(decimals," 5%x105",7,1)
      call addtop(decimals," 6%x106",7,1)
      call addtop(decimals," 7%x107",7,1)
      call addtop(decimals," 8%x108",7,1)
      call addtop(decimals," 9%x109",7,1)
      numbers=newpup()
      call addtop(numbers,"-9%x-209",8,1)
      call addtop(numbers,"-8%x-208",8,1)
      call addtop(numbers,"-7%x-207",8,1)
      call addtop(numbers,"-6%x-206",8,1)
      call addtop(numbers,"-5%x-205",8,1)
      call addtop(numbers,"-4%x-204",8,1)
      call addtop(numbers,"-3%x-203",8,1)
      call addtop(numbers,"-2%x-202",8,1)
      call addtop(numbers,"-1%x-201",8,1)
      call addtop(numbers,"0%x200",6,1)
      call addtop(numbers," 1%x201",7,1)
      call addtop(numbers," 2%x202",7,1)
      call addtop(numbers," 3%x203",7,1)
      call addtop(numbers," 4%x204",7,1)
      call addtop(numbers," 5%x205",7,1)
      call addtop(numbers," 6%x206",7,1)
      call addtop(numbers," 7%x207",7,1)
      call addtop(numbers," 8%x208",7,1)
      call addtop(numbers," 9%x209",7,1)
      skipnum=newpup()
      call addtop(skipnum," 1%x301",7,1)
      call addtop(skipnum," 2%x302",7,1)
      call addtop(skipnum," 3%x303",7,1)
      call addtop(skipnum," 4%x304",7,1)
      call addtop(skipnum," 5%x305",7,1)
      call addtop(skipnum," 6%x306",7,1)
      call addtop(skipnum," 7%x307",7,1)
      call addtop(skipnum," 8%x308",7,1)
      call addtop(skipnum," 9%x309",7,1)
      call addtop(skipnum,"10%x310",7,1)
      call addtop(skipnum,"11%x311",7,1)
      call addtop(skipnum,"12%x312",7,1)
      call addtop(skipnum,"13%x313",7,1)
      call addtop(skipnum,"14%x314",7,1)
      call addtop(skipnum,"15%x315",7,1)
      call addtop(skipnum,"16%x316",7,1)
      call addtop(skipnum,"17%x317",7,1)
      call addtop(skipnum,"18%x318",7,1)
      call addtop(skipnum,"19%x319",7,1)
      call addtop(skipnum,"20%x320",7,1)
      range=newpup()
      call addtop(range,"Set Color Map Extremes %t",25,1)
      call addtop(range,"Input Extremes %x40",19,1)
      call addtop(range,"Use Local Extremes %x41",23,1)
      call addtop(range,"Use File Extremes %x42",22,1)
      menu=newpup()
      call addtop(menu,"Main Menu %t",12,1)
      call addtop(menu,"Exit %x999",10,1)
      call addtop(menu,"Number %m",9,numbers)
      call addtop(menu,"Decimal %m",10,decimals)
      call addtop(menu,"Functions %m",12,func)
      call addtop(menu,"Extremes %m",11,extr)
      call addtop(menu,"Range %m",8,range)
      call addtop(menu,"Color/Grayscale %m",18,colour)
      call addtop(menu,"Color Bar %m",12,colorbar)
      call addtop(menu,"Shading %m",10,painting)
      call addtop(menu,"Skip %m",7,skipnum)
      call color(black)
      call clear
      return
      end
      subroutine makemap
      double precision pi,temp2
      integer*2 band,maxband
      integer*2 colormax,colormin,datamax,datamin,localmax,localmin
      integer*4 decimal,i,mapfunc,max,min,whole,paint,pred,pgreen,pblue
      logical gray,interp
      real datascale,deltamax,deltamin,deltastep,funcmax,funcmin
      real funcvalue,power,slope,spectrum(32,3),step(0:31),temp
      common/colorband/ spectrum,maxband
      common/colormap/ datascale,decimal,paint,gray,interp,mapfunc,
     .                 power,whole
      common/extreme/ colormax,colormin,datamax,datamin,localmax,
     .                localmin,max,min
      data pi/3.14159265358979323d0/
      dataminmax=datamax-datamin
      if(min.lt.0) min=0
      datascale=(max-min)/991.0
      if(mapfunc.eq.30) then
         step(0)=log(real(min+2))**power
         step(maxband-1)=log(real(max+2))**power
      else if(mapfunc.eq.31) then
         temp=2.0**(-abs(power))
         if(power.lt.0) then
            temp=-temp
         endif
         step(0)=real(min+1)**temp
         step(maxband-1)=real(max+1)**temp
      else if(mapfunc.eq.32) then
         temp=1.3**abs(power)
         if(power.lt.0) then
            temp=-temp
         endif
         step(0)=real(min+1)**temp
         step(maxband-1)=real(max+1)**temp
      else if(mapfunc.eq.33) then
         step(0)=exp(dble(min-datamin)/dble(dataminmax)*
     .               (2.0d0**power))-1.0d0
         step(maxband-1)=exp(dble(max-datamin)/dble(dataminmax)*
     .               (2.0d0**power))-1.0d0
      else if(mapfunc.eq.34) then
         temp2=cos(dble(min-datamin)/dble(dataminmax)*pi)
         temp=abs(temp2)**(1.3d0**power)
         step(0)=1.0-sign(temp,real(temp2))
         temp2=cos(dble(max-datamin)/dble(dataminmax)*pi)
         temp=abs(temp2)**(1.3d0**power)
         step(maxband-1)=1.0-sign(temp,real(temp2))
      endif
      if(step(maxband-1).lt.step(0)) then
         temp=step(maxband-1)
         step(maxband-1)=step(0)
         step(0)=temp
      endif
      if(gray) then
         deltastep=step(maxband-1)-step(0)
         colormin=min
         colormax=max
      else
         deltastep=(step(maxband-1)-step(0))/real(maxband-1)
         do 10 i=1,maxband-2
            step(i)=step(0)+real(i)*deltastep
   10    continue
         if(interp) then
            colormin=localmin
            colormax=localmax
            if(mapfunc.eq.30) then
               funcmin=log(real(localmin+2))**power
               funcmax=log(real(localmax+2))**power
            else if(mapfunc.eq.31) then
               temp=2.0**(-abs(power))
               if(power.lt.0) then
                  temp=-temp
               endif
               funcmin=real(localmin+1)**temp
               funcmax=real(localmax+1)**temp
            else if(mapfunc.eq.32) then
               temp=1.3**abs(power)
               if(power.lt.0) then
                  temp=-temp
               endif
               funcmin=real(localmin+1)**temp
               funcmax=real(localmax+1)**temp
            else if(mapfunc.eq.33) then
               funcmin=exp(dble(localmin-datamin)/dble(dataminmax)*
     .                     (2.0d0**power))-1.0d0
               funcmax=exp(dble(localmax-datamin)/dble(dataminmax)*
     .                     (2.0d0**power))-1.0d0
            else if(mapfunc.eq.34) then
               temp2=cos(dble(localmin-datamin)/dble(dataminmax)*pi)
               temp=abs(temp2)**(1.3d0**power)
               funcmin=1.0-sign(temp,real(temp2))
               temp2=cos(dble(localmax-datamin)/dble(dataminmax)*pi)
               temp=abs(temp2)**(1.3d0**power)
               funcmax=1.0-sign(temp,real(temp2))
            endif
         else
            colormin=min
            colormax=max
            funcmin=step(0)
            funcmax=step(maxband-1)
         endif
         if(funcmin.gt.funcmax) then
           temp=funcmin
           funcmin=funcmax
           funcmax=temp
         endif
         deltamin=step(1)-funcmin
         deltamax=funcmax-step(maxband-2)
      endif
      datascale=real(colormax-colormin)/991.0
      do 30 i=32,1023
         funcvalue=real(i-32)*datascale+(colormin+1)
         if(mapfunc.eq.30) then
            funcvalue=log(funcvalue+1)**power
         else if(mapfunc.eq.31) then
            temp=2.0**(-abs(power))
            if(power.lt.0) then
               temp=-temp
            endif
            funcvalue=funcvalue**temp
         else if(mapfunc.eq.32) then
            temp=1.3**abs(power)
            if(power.lt.0) then
               temp=-temp
            endif
            funcvalue=funcvalue**temp
         else if(mapfunc.eq.33) then
            funcvalue=exp((funcvalue-1.0d0-datamin)/dble(dataminmax)*
     .                    (2.0d0**power))-1.0d0
         else if(mapfunc.eq.34) then
            temp2=cos(dble(funcvalue-1.0d0-datamin)/dble(dataminmax)*pi)
            temp=abs(temp2)**(1.3d0**power)
            funcmin=1.0-sign(temp,real(temp2))
         endif
         if(gray) then
            pred=nint(255.0*(funcvalue-step(0))/deltastep)
            pgreen=pred
            pblue=pred
         else
            do 40 band=1,maxband-2
               if(funcvalue.lt.step(band)) goto 50
   40       continue
   50       if(band.eq.1) then
               slope=(funcvalue-funcmin)/deltamin
               pred=nint(spectrum(1,1)+
     .              (spectrum(2,1)-spectrum(1,1))*slope)
               pgreen=nint(spectrum(1,2)+
     .                (spectrum(2,2)-spectrum(1,2))*slope)
               pblue=nint(spectrum(1,3)+
     .               (spectrum(2,3)-spectrum(1,3))*slope)
            else if(band.eq.maxband-1) then
               slope=(funcvalue-step(maxband-2))/deltamax
               pred=nint(spectrum(maxband-1,1)+
     .              (spectrum(maxband,1)-spectrum(maxband-1,1))*slope)
               pgreen=nint(spectrum(maxband-1,2)+
     .                (spectrum(maxband,2)-spectrum(maxband-1,2))*slope)
               pblue=nint(spectrum(maxband-1,3)+
     .               (spectrum(maxband,3)-spectrum(maxband-1,3))*slope)
            else
               slope=(funcvalue-step(band-1))/deltastep
               pred=nint(spectrum(band,1)+
     .              (spectrum(band+1,1)-spectrum(band,1))*slope)
               pgreen=nint(spectrum(band,2)+
     .                (spectrum(band+1,2)-spectrum(band,2))*slope)
               pblue=nint(spectrum(band,3)+
     .               (spectrum(band+1,3)-spectrum(band,3))*slope)
            endif
         endif
         if(pred.lt.0) then
            pred=0
         else if(pred.gt.255) then
            pred=255
         endif
         if(pgreen.lt.0) then
            pgreen=0
         else if(pgreen.gt.255) then
            pgreen=255
         endif
         if(pblue.lt.0) then
            pblue=0
         else if(pblue.gt.255) then
            pblue=255
         endif
         call mapcol(i,pred,pgreen,pblue)
   30 continue
      return
      end
      subroutine selection
      character extrlabel*12,extrlist(2)*12,funclabel(5)*12,label*132
      integer*2 colormax,colormin,datamax,datamin,idim,jdim
      integer*2 localmax,localmin,maxband,skip
      integer*4 bottom,left,mapfunc,max,min,right,top
      integer*4 decimal,dev,whole
      integer*4 colour,colorbar,decimals,extr,func,funcextr,funcmap
      integer*4 menu,numbers,range,skipnum,painting
      integer*4 colman,nullwin,oldwin,paint
      integer*4 xwinlen,xwinorg,ywinlen,ywinorg
      logical gray,interp
      real colorscale,datascale
      real power,xorigin,xscale,yorigin,yscale
      real spectrum(32,3)
      common/colorband/ spectrum,maxband
      common/colormap/ datascale,decimal,paint,gray,interp,mapfunc,
     .                 power,whole
      common/extreme/ colormax,colormin,datamax,datamin,localmax,
     .                localmin,max,min
      common/labels/ extrlabel,extrlist,funclabel
      common/menus/ colour,colorbar,decimals,extr,func,menu,numbers,
     .              range,skipnum,painting
      common/plot/ bottom,left,right,skip,top,xorigin,xscale,
     .             yorigin,yscale
      common/window/ colman,idim,jdim,xwinlen,xwinorg,ywinlen,ywinorg
$ include /usr/include/fgl.h
$ include /usr/include/fdevice.h
      colorscale=741.0/991.0
   10 dev=dopup(menu)
      if(abs(dev).ge.100.and.abs(dev).le.109) then
c
c decimals
c
         dev=sign(abs(dev)-100,dev)
   30    if((mapfunc.eq.30.and.whole.eq.0.and.dev.eq.0).or.
     .      (mapfunc.eq.33.and.(2.**(whole+real(dev)/10.)).gt.88.722))
     .                                                              then
            call ringbe
            dev=dopup(decimals)
            dev=sign(abs(dev)-100,dev)
            goto 30
         endif
         decimal=dev
         power=real(whole)+real(decimal)/10.0
      else if(abs(dev).ge.200.and.abs(dev).le.209) then
c
c numbers
c
         dev=sign(abs(dev)-200,dev)
   40    if((mapfunc.eq.30.and.dev.eq.0.and.decimal.eq.0).or.
     .      (mapfunc.eq.33.and.(2.**(dev+real(decimal)/10.)).gt.88.722))
     .                                                              then
            call ringbe
            dev=dopup(numbers)
            dev=sign(abs(dev)-200,dev)
            goto 40
         endif
         whole=dev
         power=real(whole)+real(decimal)/10.0
      else if(dev.ge.301.and.dev.lt.400) then
c
c skiping step
c
         skip=dev-300
c
c color/grayscale
c
      else if(dev.eq.10) then
c colors
         gray=.false.
      else if(dev.eq.11) then
c gray scale
         gray=.true.
      else if(dev.eq.20) then
c
c interpolation to file extremes
c
         interp=.true.
         extrlabel=extrlist(1)
      else if(dev.eq.21) then
c
c truncation to specified limits
c
         interp=.false.
         extrlabel=extrlist(2)
      else if(dev.ge.30.and.dev.le.39) then
c
c distribution functions
c
         mapfunc=dev
         if((mapfunc.eq.30.or.mapfunc.eq.31.or.mapfunc.eq.32).and.
     .      power.eq.0) then
            decimal=0
            whole=1
            power=1.0
         endif
         if(mapfunc.eq.33.and.(2.0**power).gt.88.722) then
            decimal=4
            whole=6
            power=6.4
         endif
c
c set colormap extremes
c
      else if(dev.eq.40) then
c
c Input Extremes
c
         call noport
         nullwin=winope('null',4)
         oldwin=winatt()
         call winclo(nullwin)
         call ringbe
         read(*,*) max,min
         call winset(colman)
         oldwin=winatt()
         call qreset
      else if(dev.eq.41) then
c
c Use Local Extremes
c
         max=localmax
         min=localmin
      else if(dev.eq.42) then
c
c Use File Extremes
c
         max=datamax
         min=datamin
      else if(dev.ge.50.and.dev.le.69) then
c
c choose color bar colors
c
         if(dev.ne.69) maxband=0
   50    if(dev.eq.50) then
c White
            maxband=maxband+1
            spectrum(maxband,1)=255
            spectrum(maxband,2)=255
            spectrum(maxband,3)=255
         else if(dev.eq.51) then
c Magenta
            maxband=maxband+1
            spectrum(maxband,1)=255
            spectrum(maxband,2)=0
            spectrum(maxband,3)=255
         else if(dev.eq.52) then
c Red
            maxband=maxband+1
            spectrum(maxband,1)=255
            spectrum(maxband,2)=0
            spectrum(maxband,3)=0
         else if(dev.eq.53) then
c Yellow
            maxband=maxband+1
            spectrum(maxband,1)=255
            spectrum(maxband,2)=255
            spectrum(maxband,3)=0
         else if(dev.eq.54) then
c Green
            maxband=maxband+1
            spectrum(maxband,1)=0
            spectrum(maxband,2)=255
            spectrum(maxband,3)=0
         else if(dev.eq.55) then
c Cyan
            maxband=maxband+1
            spectrum(maxband,1)=0
            spectrum(maxband,2)=255
            spectrum(maxband,3)=255
         else if(dev.eq.56) then
c Blue
            maxband=maxband+1
            spectrum(maxband,1)=0
            spectrum(maxband,2)=0
            spectrum(maxband,3)=255
         else if(dev.eq.57) then
c Black
            maxband=maxband+1
            spectrum(maxband,1)=0
            spectrum(maxband,2)=0
            spectrum(maxband,3)=0
         else if(dev.eq.69.and.maxband.ge.2) then
            goto 10
         endif
         if(maxband.lt.32) then
            dev=dopup(colorbar)
            goto 50
         endif
c Painting type
      else if(dev.eq.70) then
         paint=0
      else if(dev.eq.71) then
         paint=1
      else if(dev.eq.72) then
         paint=2
      else if(dev.eq.999) then
c Exit
         call makemap
         call getori(xwinorg,ywinorg)
         call getsiz(xwinlen,ywinlen)
         call viewpo(0,xwinlen-1,0,ywinlen-1)
         call ortho2(0.0,1023.0,0.0,767.0)
         call lookat(0.0,0.0,1.0,0.0,0.0,0.0,0.0)
         call color(black)
         call rectf(988.0,0.0,1023.0,767.0)
         call rectf(0.0,0.0,1023.0,27.0)
         do 60 index=32,1023
            y=anint(real(index-32)*colorscale+14.0)
            call color(index)
            call pmv2(988.0,y)
            call pdr2(1023.0,y)
            call pdr2(1023.0,y+colorscale)
            call pdr2(988.0,y+colorscale)
            call pclos
   60    continue
         call color(white)
         write(label,1000) real(left-1)*xscale+xorigin,
     .                     real(bottom-1)*yscale+yorigin,
     .                     real(right-1)*xscale+xorigin,
     .                     real(top-1)*yscale+yorigin,
     .                     funclabel(mapfunc-29),power,extrlabel,
     .                     localmax,localmin
         call cmov2i(0,7*ywinlen/768)
         call charst(label,109)
         write(label,1010) min
         call cmov2i(988,2)
         call charst(label,4)
         write(label,1010) max
         call cmov2i(988,756)
         call charst(label,4)
         return
      endif
      goto 10
 1000 format('l-b=(',e13.7,',',e13.7,') r-t=(',e13.7,',',e13.7,') ',
     .       a12,' P=',f4.1,1x,a5,' sx=',i4,' sn=',i4)
 1010 format(i4)
      end