c rainbow.for
c
c creates red, green, and blue density slices from a standard grid
c
c Jeff Phillips 8 Aug 91
      dimension x(3000),r(3000),g(3000),b(3000),an(6)
      dimension id(14),pgm(2)
      character*50 ifile
      data ddval/1.e37/
      print *,'Enter input grid name: '
      read(5,100) ifile
  100 format(a)
      open(10,file=ifile,form='unformatted',status='old')
      read(10) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      am=0.
      as=0.
      n=0
      an(1)=1.e38
      an(6)=-1.e38
      do 10 j=1,nrow
      read(10) dlt,(x(i),i=1,ncol)
      do 10 i=1,ncol
      if(x(i).gt.ddval) go to 10
      am=am+x(i)
c      as=as+x(i)**2
      n=n+1
      if(x(i).gt.an(6)) an(6)=x(i)
      if(x(i).lt.an(1)) an(1)=x(i)
   10 continue
      am=am/n
      rewind(10)
      read(10) id
      do 11 j=1,nrow
      read(10) dlt,(x(i),i=1,ncol)
      do 11 i=1,ncol
      if(x(i).gt.ddval) go to 11
      as=as+(x(i)-am)**2
   11 continue
      rewind(10)
      as=sqrt(as/(n-1))
      print *,'data range = ',an(1),an(6)
      print *,'mean = ',am
      print *,'s.d. = ',as
      an(6)=min(an(6),am+3.*as)
      an(1)=max(an(1),am-3.*as)
      print *,'new data range = ',an(1),an(6)
      do 20 i=2,5
      an(i)=an(1)+(i-1)*(an(6)-an(1))/5.
   20 continue
      open(11,file='red.grd',form='unformatted',status='unknown')
      open(12,file='green.grd',form='unformatted',status='unknown')
      open(13,file='blue.grd',form='unformatted',status='unknown')
      read(10) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      write(11) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      write(12) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      write(13) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      do 40 j=1,nrow
      read(10) dlt,(x(i),i=1,ncol)
      do 30 i=1,ncol
      if(x(i).lt.an(1)) x(i)=an(1)
      if(x(i).lt.an(2)) then
           r(i)=0.
           g(i)=(x(i)-an(1))/(an(2)-an(1))
           b(i)=1.
      else if (x(i).lt.an(3)) then
           r(i)=0.
           g(i)=1.
           b(i)=(an(3)-x(i))/(an(3)-an(2))
      else if (x(i).lt.an(4)) then
           r(i)=(x(i)-an(3))/(an(4)-an(3))
           g(i)=1.
           b(i)=0.
      else if (x(i).lt.an(5)) then
           r(i)=1.
           g(i)=(an(5)-x(i))/(an(5)-an(4))
           b(i)=0.
      else if (x(i).lt.ddval) then
           if(x(i).gt.an(6)) x(i)=an(6)
           r(i)=1.
           g(i)=0.
           b(i)=(x(i)-an(5))/(an(6)-an(5))
      else
           r(i)=x(i)
           g(i)=x(i)
           b(i)=x(i)
      endif
   30 continue
      write(11) dlt,(r(i),i=1,ncol)
      write(12) dlt,(g(i),i=1,ncol)
      write(13) dlt,(b(i),i=1,ncol)
   40 continue
      stop
      end
