C 
C Subroutines to linearly interpolate discrete colour maps in RGB space
C
C Version: 5.50
C Date: 2000, November 8
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C                                                    
C Description of the data file specifying the colour map:
C (1) N0,N1,N2,...,NK,/
C     Colour map is specified on the (K+1)-dimensional rectangular grid
C     of N0*N1*N2*...*Nk points situated within (K+1)-dimensional unit
C     cube.
C (2) (K+1) times (2.1), for J=0,1,...,K
C (2.1) CI(1),CI(2),...,CI(NJ)
C     Grid coordinates.  Must be 0.LE.CI(1).LE.CI(2).LE. ... .LE.CI(NI).
C (3) R(I0,I1,...,IK),G(I0,I1,...,IK),B(I0,I1,...,IK) for I0=1,2,...,N0;
C     for I1=1,2,...,N1; ...; for IK=1,2,...,NK
C     RGB components of grid values of the colour map.  The colours are
C     linearly (bilinearly, trilinearly, ...) interpolated between the
C     gridpoints and  are constant between the sides of the unit cube
C     and the nearest gridpoint.
C (4) (K+1) times (2.1), for J=0,1,...,K
C (4.1) REPL,REPR,CREF,VREFNORM
C     REPL,REPR... Colours of interval (0,1) are cyclically repeated
C             along the interval (-REPL,1+REPR).  Colours left to -REPL
C             are constant, colours right to 1+REPR are constant.
C     CREF... Default reference colour.  It may be changed by the input
C             SEP file using parameters CREF, CREF1, CREF2, ...
C     VREFNORM... Determines the default reference value.  The reference
C             value may be changed by the input SEP file using
C             parameters VREF, VREF1, VREF2, ...
C For an example refer to colour map hsv.dat
C
C                                                     
C Data file 'SEP' has the form of the SEP (Stanford Exploration Project)
C parameter file:
C     All the data are specified in the form of PARAMETER=VALUE, e.g.
C     N1=50, with PARAMETER directly preceding = without intervening
C     spaces and with VALUE directly following = without intervening
C     spaces.  The PARAMETER=VALUE couple must be delimited by a space
C     or comma from both sides.
C     The PARAMETER string is not case-sensitive.
C     PARAMETER= followed by a space resets the default parameter value.
C     All other text in the input files is ignored.  The file thus may
C     contain unused data or comments without leading comment character.
C     Everything between comment character # and the end of the
C     respective line is ignored, too.
C     The PARAMETER=VALUE couples may be specified in any order.
C     The last appearance takes precedence.
C Data specifying the colour scale:
C     VADD=real... Controls the default value of VPER.
C             Default: VADD=0.
C     VMUL=real... Controls the default value of VPER.
C             Default: VMUL=1.
C     VPER=real... Period of values corresponding to one period in the
C             colour map.
C             Default: VPER=VMUL*(GMAX-GMIN+VADD),
C             where GMIN and GMAX are the minimum and maximum values to
C             be displayed in colours.
C     VREF=real... Reference value.  It will be displayed in the
C             reference colour.
C             Default: VREF=GMIN+(GMAX-GMIN)*VREFNORM,
C             where VREFNORM is taken from the colour map file.
C             For COLORS='hsv.dat', default VREF=GMIN.
C     CREF=real... Reference colour.
C             Default value is taken from the colour map file.
C             For COLORS='hsv.dat', default CREF=0.666667 (blue).
C     VADD1=real, VADD2=real, ..., VMUL1=real, VMUL2=real, ...,
C     VPER1=real, VPER2=real, ..., VREF1=real, VREF2=real, ...,
C     CREF1=real, CREF2=real, etc... Analogous to VADD, VMUL, VPER, VREF
C             and CREF for multidimensional colour maps.
C             Default values are also analogous.
C             For COLORS='hsv.dat' (three-dimensional colour map),
C             defaults are determined by VREFNORM1=1., VREFNORM2=1.,
C             CREF1=1. (maximum saturation), CREF2=1. (maximum
C             brightness).
C
C=======================================================================
C
      SUBROUTINE COLOR1(LU,MRAM,IRAM,RAM,NVALUE,VALMIN,VALMAX)
C
      INTEGER LU,MRAM,IRAM(0:MRAM-1),NVALUE
      REAL RAM(0:MRAM-1),VALMIN(NVALUE),VALMAX(NVALUE)
C
C-----------------------------------------------------------------------
C
C     External functions and subroutines:
      EXTERNAL RSEP3T,RSEP3R,ERROR
C
C     Storage in array (I)RAM:
C     IRAM(0)=K... Dimensionality of the colour map.
C     IRAM(1:K)=(L1,L2,...,LK)... Last indices of the grid coordinates
C             of the colour map.  L0=K, L(i)=L(i-1)+N(i) where
C             N1,N2,...,NK are the numbers of gridpoints of the colour
C             map.
C     RAM(K+1:LK)... Grid coordinates.
C     RAM(LK+1:LK+3*N1*N2*...*NK)... Gridpoint RGB colours.
C     RAM(LK+3*N1*N2*...*NK+1:LK+3*N1*N2*...*NK+NREF*K)... For each grid
C             coordinate: the left and right colour repetitions, the
C             reference colour, the reference normalized value switched
C             into the reference value, the period of values, auxiliary
C             storage locations.
C
C.......................................................................
C
      CHARACTER*80 FILE
      CHARACTER*5  TEXT
      REAL VADD,VMUL,VPER,VREF,CREF
      INTEGER NGRID,NREF,I1,I2,I
      PARAMETER (NREF=8)
C
C     NGRID=LK+3*N1*N2*...*NK... Location of the last grid colour in
C             array RAM.
C     NREF... Number of reference values for each grid coordinate.
C
C.......................................................................
C
      IF(10.GT.MRAM-1) THEN
C       COLORS-01
        CALL ERROR('COLORS-01: Too small array RAM')
      END IF
      CALL RSEP3T('COLORS',FILE,'hsv.dat' )
      OPEN(LU,FILE=FILE)
      DO 11 I=1,10
        IRAM(I)=0
   11 CONTINUE
      READ(LU,*) (IRAM(I),I=1,10)
      NGRID=1
      DO 12 I=1,10
        IF(IRAM(I).LE.0) THEN
          IRAM(0)=I-1
          GO TO 13
        END IF
        NGRID=NGRID*IRAM(I)
   12 CONTINUE
C       COLORS-02
        CALL ERROR('COLORS-02: More than 9 colour coordinates')
   13 CONTINUE
      DO 14 I=1,IRAM(0)
        IRAM(I)=IRAM(I-1)+IRAM(I)
   14 CONTINUE
      NGRID=IRAM(IRAM(0))+3*NGRID
      IF(NGRID+NREF*IRAM(0).GT.MRAM-1) THEN
C       COLORS-03
        CALL ERROR('COLORS-03: Too small array RAM')
      END IF
C     Reading grid coordinates
      DO 15 I2=1,IRAM(0)
        READ(LU,*) (RAM(I1),I1=IRAM(I2-1)+1,IRAM(I2))
   15 CONTINUE
C     Reading grid RGB values
      READ(LU,*) (RAM(I1),I1=IRAM(IRAM(0))+1,NGRID)
C     Reading repetitions, reference colours and reference values
      DO 17 I2=NGRID,NGRID+NREF*(IRAM(0)-1),NREF
        READ(LU,*) (RAM(I1),I1=I2+1,I2+4)
   17 CONTINUE
C
      TEXT=' '
      DO 21 I2=1,MIN0(IRAM(0),NVALUE)
        I=NGRID+NREF*(I2-1)
        IF(I2.GT.1) THEN
          TEXT(5:5)=CHAR(ICHAR('0')+I2-1)
        END IF
        TEXT(1:4)='CREF'
        CREF=RAM(I+3)
        CALL RSEP3R(TEXT,RAM(I+3),CREF)
        TEXT(1:4)='VREF'
        VREF=VALMIN(I2)+(VALMAX(I2)-VALMIN(I2))*RAM(I+4)
        CALL RSEP3R(TEXT,RAM(I+4),VREF)
        TEXT(1:4)='VADD'
        CALL RSEP3R(TEXT,VADD,0.)
        TEXT(1:4)='VMUL'
        CALL RSEP3R(TEXT,VMUL,1.)
        TEXT(1:4)='VPER'
        VPER=(VALMAX(I2)-VALMIN(I2)+VADD)*VMUL
        CALL RSEP3R(TEXT,RAM(I+5),VPER)
   21 CONTINUE
      DO 22 I2=MIN0(IRAM(0),NVALUE)+1,IRAM(0)
        I=NGRID+NREF*(I2-1)
        IF(I2.GT.1) THEN
          TEXT(5:5)=CHAR(ICHAR('0')+I2-1)
        END IF
        TEXT(1:4)='CREF'
        CREF=RAM(I+3)
        CALL RSEP3R(TEXT,RAM(I+3),CREF)
   22 CONTINUE
      CLOSE(LU)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE COLOR2(MRAM,IRAM,RAM,NVALUE,VALUE,R,G,B)
C
      INTEGER MRAM,IRAM(0:MRAM-1),NVALUE
      REAL RAM(0:MRAM-1),VALUE(NVALUE),R,G,B
C
C-----------------------------------------------------------------------
C
      REAL VPER,VREF,CREF,COLOR,W
      INTEGER NREF,NGRID,IGRID,I1,I2,I,J,N
      PARAMETER (NREF=8)
C
C     NREF... Number of reference values for each grid coordinate.
C     NGRID...Location of the last grid colour in array RAM.
C
C.......................................................................
C
C     Index of the grid origin = initial index of the cube origin
      IGRID=IRAM(IRAM(0))
C     Index of the last grid value
      NGRID=1
      DO 11 I=1,IRAM(0)
        NGRID=NGRID*(IRAM(I)-IRAM(I-1))
   11 CONTINUE
      NGRID=IGRID+3*NGRID
C
      N=3
C     N is the shift between indices of neighbouring gridpoints
      DO 39 I2=1,IRAM(0)
        I=NGRID+NREF*(I2-1)
        CREF=RAM(I+3)
        VREF=RAM(I+4)
        VPER=RAM(I+5)
        IF(I2.LE.NVALUE) THEN
          COLOR=(VALUE(I2)-VREF)/VPER+CREF
        ELSE
          COLOR=CREF
        END IF
        COLOR=AMAX1(COLOR,  -RAM(I+1))
        COLOR=AMIN1(COLOR,1.+RAM(I+2))
        IF(COLOR.LT.0.) THEN
          COLOR=COLOR-AINT(COLOR)+1.
        ELSE IF(COLOR.GT.1.) THEN
          COLOR=COLOR-AINT(COLOR)
          IF(COLOR.LE.0.) THEN
            COLOR=1.
          END IF
        END IF
        DO 31 I1=IRAM(I2-1)+1,IRAM(I2)
          IF(COLOR.LE.RAM(I1)) THEN
C           Colour is located left to the grid coordinate RAM(I1)
            IF(I1.LE.IRAM(I2-1)+1) THEN
              RAM(I+7)=1.
              RAM(I+8)=0.
              J=0
            ELSE
              J=I1-IRAM(I2-1)-2
              RAM(I+8)=(COLOR-RAM(I1-1))/(RAM(I1)-RAM(I1-1))
              RAM(I+7)=1.-RAM(I+8)
            END IF
            GO TO 32
          END IF
   31   CONTINUE
        J=IRAM(I2)-IRAM(I2-1)-2
        RAM(I+7)=0.
        RAM(I+8)=1.
   32   CONTINUE
        IGRID=IGRID+N*J
        IRAM(I+6)=N
        N=N*(IRAM(I2)-IRAM(I2-1))
   39 CONTINUE
C
      R=0.
      G=0.
      B=0.
C     Loop over the vertices of the IRAM(0)-dimensional cube:
      DO 49 I2=0,2**IRAM(0)-1
        N=IGRID
        W=1.
        J=I2
        DO 45 I1=0,IRAM(0)-1
          I=MOD(J,2)
          J=J/2
          N=N+I*IRAM(NGRID+NREF*I1+6)
          W=W*RAM(NGRID+NREF*I1+7+I)
   45   CONTINUE
        R=R+W*RAM(N+1)
        G=G+W*RAM(N+2)
        B=B+W*RAM(N+3)
   49 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE COLOR3(MRAM,IRAM,RAM,NVALUE,IREF,IRGB)
C
      INTEGER MRAM,IRAM(0:MRAM-1),NVALUE,IREF,IRGB
      REAL RAM(0:MRAM-1)
C
C     IREF... Index of CREF in array RAM:
C             RAM(IREF  )=CREF,  RAM(IREF+1)=VREF,  RAM(IREF+2)=VPER,
C             RAM(IREF+8)=CREF1, RAM(IREF+9)=VREF1, RAM(IREF+10)=VPER1,
C             etc.
C     IRGB... Index of the first output grid value in array RAM:
C             RAM(IREF  )=CREF,  RAM(IREF+1)=VREF,  RAM(IREF+2)=VPER,
C             RAM(IREF+8)=CREF1, RAM(IREF+9)=VREF1, RAM(IREF+10)=VPER1,
C             etc.
C
C-----------------------------------------------------------------------
C
      INTEGER NREF,I1,I2,I,N
      PARAMETER (NREF=8)
C
C     NREF... Number of reference values for each grid coordinate.
C
C.......................................................................
C
C     Index of the first reference value:
      IREF=1
      DO 11 I1=1,IRAM(0)
        IREF=IREF*(IRAM(I1)-IRAM(I1-1))
   11 CONTINUE
      IREF=IRAM(IRAM(0))+3*IREF+3
C
C     Index of the first output grid value:
      IRGB=IREF-2+NREF*IRAM(0)
C
      N=1
      DO 21 I1=1,NVALUE
        N=N*(IRAM(I1)-IRAM(I1-1))
   21 CONTINUE
      IF(IRGB+3*N+NVALUE.GT.MRAM) THEN
C       COLORS-04
        CALL ERROR('COLORS-04: Too small array RAM')
      END IF
      DO 29 I2=0,N-1
        I=I2
        DO 28 I1=0,NVALUE-1
          RAM(IRGB+3*N+I1)=(RAM(IRAM(I1)+MOD(I,IRAM(I1+1)-IRAM(I1))+1)
     *                -RAM(IREF+8*I1))*RAM(IREF+8*I1+2)+RAM(IREF+8*I1+1)
          I=I/(IRAM(I1+1)-IRAM(I1))
   28   CONTINUE
        CALL COLOR2(MRAM,IRAM(0),RAM(0),NVALUE,RAM(IRGB+3*N),
     *                 RAM(IRGB+3*I2),RAM(IRGB+3*I2+1),RAM(IRGB+3*I2+2))
   29 CONTINUE
      RETURN
      END
C
C=======================================================================
C