      program mash3
C     Inquires of the user stellar parameters and returns colors, with errors

      integer      nind,nvk,nfe,ng
      parameter    (nind=9,nvk=86,nfe=7,ng=13)
C     The 8th nind is the temperature. 
C     There are actually only 75 vk points, but don't change the nvk parameter.
C     The 9th nind is BC_v
      real         a(nvk,nfe,ng,nind)
      real         grav,feh,teff,theta,clrs(nind),cerr(nind)
      real         ete(8),eterr(8),cpres(9),bcefloor(8)
C     out-of-bounds flag. Nominal is zero.
      integer      iflag

C     theta and assumed percentage error arrays:
      data ete /0.1008,0.252,0.504,0.84,1.26,1.44,1.68,2.52 /
c         Teff= 50000, 20000,10000,6000,4000,3500,3000,2000
      data eterr / 4.0, 2.5, 1.0, 0.5, 0.5, 1.0, 1.5, 4.0 /
C     The cool errors probably ought to be a bit bigger for dwarfs, but that's
C       not implemented...

C     color (or BC) intrinsic precisions (entry 8, 9 are dummies)
C      data cpres / .014, .007, .007, .007, .007, .007, .007, .0, .05 /
      data cpres / .071, .017, .010, .010, .011, .004, .002, .0, .05 /
C     make bolometric correction floor error depend on theta
      data bcefloor / 0.2, 0.1, 0.07, 0.05, 0.05, 0.07, 0.10, 0.2 /

C     set extrapolation flag to nominal value
      iflag = 0

C     load data table
      call readtable(a)

C     inquire as to input parameters
      print*, 'enter log g '
      read*, grav
      print*, 'enter [Fe/H]'
      read*, feh
      print*, 'enter Teff '
      read*, teff
      print*, ' '
C     switch Teff to THETA
      theta = 5040./teff

C     get colors for these parameters
      call teffinterp(a,grav,feh,theta,clrs,iflag)
C     iflag warns of extrapolations.
C     Nominally zero, it becomes -N if too hot or N if too cool,
C     where N is the number of point spacings extrapolated.
C     Obviously, abs(iflag) > 1 is of concern . . .
      if ( iflag .ne. 0 ) print*, 'Extrapolation warning! Iflag =',iflag
C     and colors for 1% lower temperature to compute errors
      call teffinterp(a,grav,feh,1.01*theta,cerr,iflag)
      call linear(ete,eterr,8,theta,xerr,iOK)
      if (iOK.le.-2) print*, 'Panic! Interpolation error!'
      do i=1,7
         cerr(i) = sqrt(cpres(i)**2 + (xerr**2)*(cerr(i)-clrs(i))**2)
      end do
C     BCv error
      call linear(ete,bcefloor,8,theta,xbcerr,iOK)
      cerr(9) = sqrt(xbcerr**2 + (xerr**2)*(cerr(9)-clrs(9))**2)

      write(*,'(a6,f7.3,a5,f5.3)') 'U-B = ',clrs(1),' +/- ',cerr(1)
      write(*,'(a6,f7.3,a5,f5.3)') 'B-V = ',clrs(2),' +/- ',cerr(2)
      write(*,'(a6,f7.3,a5,f5.3)') 'V-R = ',clrs(3),' +/- ',cerr(3)
      write(*,'(a6,f7.3,a5,f5.3)') 'V-I = ',clrs(4),' +/- ',cerr(4)
      write(*,'(a6,f7.3,a5,f5.3)') 'J-K = ',clrs(5),' +/- ',cerr(5)
      write(*,'(a6,f7.3,a5,f5.3)') 'H-K = ',clrs(6),' +/- ',cerr(6)
      write(*,'(a6,f7.3,a5,f5.3)') 'V-K = ',clrs(7),' +/- ',cerr(7)
      write(*,'(a6,f7.3,a5,f5.3)') 'BCv = ',clrs(9),' +/- ',cerr(9)

      stop
      end

C--------------------------------------------------------------------------
      subroutine readtable(a)

      integer      nind,nvk,nfe,ng,i
      parameter    (nind=9,nvk=86,nfe=7,ng=13)
C     the 8th nind is the temperature
C     the 9th nind is BC_v
      real         a(nvk,nfe,ng,nind),g(ng),fe(nfe)

      data g / -0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5 /
      data fe /  -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5  /

      open(unit=19,file='mash2.out',status='old')
      do ife = 1,nfe
         do ig = 1,ng
            do ivk=1,75
               read(19,*) x1,x2,xteff,x3,
     z                 (a(ivk,ife,ig,k),k=1,7),a(ivk,ife,ig,9)
               if ( abs(x1 - fe(ife)).gt. 0.01 ) print*, 'bad fe'
               if ( abs(x2 - g(ig)).gt. 0.01 ) print*, 'bad g'
C              transform to THETA rather than Teff . . . 
               a(ivk,ife,ig,8) = 5040./xteff
            end do
         end do
      end do

      return
      end
C----------------------------------------------------------------------------

C     subroutine teffinterp
      subroutine teffinterp(a,grav,feh,theta,clrs,iflag)
C     given a THETA (=5040/Teff), return colors
      integer      nind,nvk,nfe,ng,iflag
      parameter    (nind=9,nvk=86,nfe=7,ng=13)
      real         grav,feh,theta,clrs(nind)
      real         a(nvk,nfe,ng,nind),g(ng),fe(nfe),vk(nvk)
C     the 8th nind is the temperature, 9th BC_v
C     local variables
      integer      jg,jfe
      real         c(nvk,nind)

      data g / -0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5 /
      data fe /  -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5  /
C     clrs are  1 U-B  2 B-V  3 V-R  4 V-I  5 J-K  6 H-K  7 V-K  (8 Teff)
C     a() is organised: a(nvk,nfe,ng,nind)                       (9 BCv)

C     find jfe and jg interpolation corners
      call locate(g,ng,grav,jg)
      if (jg.eq.0) jg = 1
      if (jg.eq.ng) jg = ng-1
      call locate(fe,nfe,feh,jfe)
      if (jfe.eq.0) jfe=1
      if (jfe.eq.nfe) jfe=nfe-1
C     fill c array with bilinear-interp results
      ffe = ( feh - fe(jfe) )/( fe(jfe+1)-fe(jfe)  )
      gg  = ( grav - g(jg) )/( g(jg+1) - g(jg)  )
      do ivk=1,75
         do ind=1,nind
            c(ivk,ind) = (1.0-ffe)*(1.0-gg)*a(ivk,jfe,jg,ind)
     z                 +  ffe     *(1.0-gg)*a(ivk,jfe+1,jg,ind)
     z                 +  ffe     *gg      *a(ivk,jfe+1,jg+1,ind)
     z                 + (1.0-ffe)*gg      *a(ivk,jfe,jg+1,ind)
         end do
      end do

C     find temperature (it's really THETA) and interpolate colors
      call locate(c(1,8),75,theta,jt)
      jt0 = jt
      if (jt.eq.0) jt = 1
      if (jt.gt.70) jt = 71
      do i=1,7
         call polint(c(jt,8),c(jt,i),5,theta,clrs(i),dy)
      end do
      call polint(c(jt,8),c(jt,9),5,theta,clrs(9),dy)
C     if requested temperature is out-of-bounds, use linear interpolation to
C     extrapolate. Return iflag = int(number of segments beyond the tabulated)
      if ( jt0 .eq. 0 ) then
         frac = (theta - c(1,8))/(c(2,8)-c(1,8))
         iflag = -1 + int(frac)
         do i=1,7
            clrs(i) = (1.0-frac)*c(1,i) + frac*c(2,i)
         end do
            clrs(9) = (1.0-frac)*c(1,9) + frac*c(2,9)
      end if
      if ( jt0 .eq. 75 ) then
         frac = (theta - c(74,8))/(c(75,8)-c(74,8))
         iflag = 1 + int(frac)
         do i=1,7
            clrs(i) = (1.0-frac)*c(74,i) + frac*c(75,i)
         end do
            clrs(9) = (1.0-frac)*c(74,9) + frac*c(75,9)
      end if

      return
      end
C-----------------------------------------------------------------------

C     subroutine vkinterp
      subroutine vkinterp(a,grav,feh,theta,clrs)
C     given a V-K (in clrs(7)), return colors and THETA = 5040/Teff
      integer      nind,nvk,nfe,ng
      parameter    (nind=9,nvk=86,nfe=7,ng=13)
      real         grav,feh,theta,clrs(nind)
      real         a(nvk,nfe,ng,nind),g(ng),fe(nfe),vk(nvk)
C     the 8th nind is the temperature, 9th BC_v
C     local variables
      integer      jg,jfe
      real         c(nvk,nind)

      data g / -0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5 /
      data fe /  -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5  /
C     clrs are  1 U-B  2 B-V  3 V-R  4 V-I  5 J-K  6 H-K  7 V-K  (8 Teff)
C     a() is organised: a(nvk,nfe,ng,nind)                       (9 BCv)

C     find jfe and jg interpolation corners
      call locate(g,ng,grav,jg)
      if (jg.eq.0) jg = 1
      if (jg.eq.ng) jg = ng-1
      call locate(fe,nfe,feh,jfe)
      if (jfe.eq.0) jfe=1
      if (jfe.eq.nfe) jfe=nfe-1
C     fill c array with bilinear-interp results
      ffe = ( feh - fe(jfe) )/( fe(jfe+1)-fe(jfe)  )
      gg  = ( grav - g(jg) )/( g(jg+1) - g(jg)  )
      do ivk=1,75
         do ind=1,nind
            c(ivk,ind) = (1.0-ffe)*(1.0-gg)*a(ivk,jfe,jg,ind)
     z                 +  ffe     *(1.0-gg)*a(ivk,jfe+1,jg,ind)
     z                 +  ffe     *gg      *a(ivk,jfe+1,jg+1,ind)
     z                 + (1.0-ffe)*gg      *a(ivk,jfe,jg+1,ind)
         end do
      end do

C     find V-K and interpolate THETA and colors
      call locate(c(1,7),75,clrs(7),jt)
      if (jt.eq.0) jt = 1
      if (jt.gt.70) jt = 71
      do i=1,6
         call polint(c(jt,7),c(jt,i),5,clrs(7),clrs(i),dy)
      end do
      call polint(c(jt,7),c(jt,8),5,clrs(7),clrs(8),dy)
      call polint(c(jt,7),c(jt,9),5,clrs(7),clrs(9),dy)
      theta = clrs(8)

      return
      end
C-----------------------------------------------------------------------
C     subroutine viinterp
      subroutine viinterp(a,grav,feh,theta,clrs)
C     given a V-I (in clrs(4)), return colors and THETA = 5040/Teff
      integer      nind,nvk,nfe,ng
      parameter    (nind=9,nvk=86,nfe=7,ng=13)
      real         grav,feh,theta,clrs(nind)
      real         a(nvk,nfe,ng,nind),g(ng),fe(nfe),vk(nvk)
C     the 8th nind is the temperature, 9th BCv
C     local variables
      integer      jg,jfe
      real         c(nvk,nind)

      data g / -0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5 /
      data fe /  -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5  /
C     clrs are  1 U-B  2 B-V  3 V-R  4 V-I  5 J-K  6 H-K  7 V-K  (8 Teff)
C     a() is organised: a(nvk,nfe,ng,nind)

C     find jfe and jg interpolation corners
      call locate(g,ng,grav,jg)
      if (jg.eq.0) jg = 1
      if (jg.eq.ng) jg = ng-1
      call locate(fe,nfe,feh,jfe)
      if (jfe.eq.0) jfe=1
      if (jfe.eq.nfe) jfe=nfe-1
C     fill c array with bilinear-interp results
      ffe = ( feh - fe(jfe) )/( fe(jfe+1)-fe(jfe)  )
      gg  = ( grav - g(jg) )/( g(jg+1) - g(jg)  )
      do ivk=1,75
         do ind=1,nind
            c(ivk,ind) = (1.0-ffe)*(1.0-gg)*a(ivk,jfe,jg,ind)
     z                 +  ffe     *(1.0-gg)*a(ivk,jfe+1,jg,ind)
     z                 +  ffe     *gg      *a(ivk,jfe+1,jg+1,ind)
     z                 + (1.0-ffe)*gg      *a(ivk,jfe,jg+1,ind)
         end do
      end do

C     find V-I and interpolate THETA and colors
      call locate(c(1,4),75,clrs(4),jt)
      if (jt.eq.0) jt = 1
      if (jt.gt.70) jt = 71
      do i=1,3
         call polint(c(jt,4),c(jt,i),5,clrs(4),clrs(i),dy)
      end do
      do i=5,9
         call polint(c(jt,4),c(jt,i),5,clrs(4),clrs(i),dy)
      end do
      theta = clrs(8)

      return
      end
C-----------------------------------------------------------------------
C------------------------------------------------------------------------
C     subroutine linear
C     quick linear interpolation
      subroutine linear(x,y,nxy,xin,yout,iOK)
      integer    nxy,iOK
      real       x(nxy),y(nxy),xin,yout
C     x and y are nxy long
C     xin is the nontabulated input x value
C     yout is the interpolated y guess.
C     iOK is -3 if X is not in ascending order, -2 if input x is
C     out-of-bounds by more than 1 xpoint spacing, -1 if out-of-bounds
C     by less than 1 xpoint spacing, 0 if all is OK

C     local variables
      integer    j,jl,ju,jm
      real       frac

      iOK = 0
      if ( x(2).lt.x(1) ) then
         print*, 
     z   'Error. Sub LINEAR. Input X array must be in ascending order.'
         yout = 0.0
         iOK = -3
         return
      end if

C     locate correct array element by bisection
      jl=0
      ju=nxy+1
10    if (ju-jl.gt.1) then
         jm = (ju+jl)/2
         if ((x(nxy).gt.x(1)).eqv.(xin.gt.x(jm))) then
            jl=jm
         else
            ju=jm
         end if
         goto 10
      end if
      j = jl
C     j is 0 or nxy if xin is off the grid

C     if off-grid, reset j and set output flag iOK
      if ( j.eq.0) then
         if ( xin.lt.(x(1)-(x(2)-x(1))) ) then
            iOK = -2
         else
            iOK = -1
         end if
         j=1
      endif
      if ( j.eq.nxy) then
         if ( xin.gt.(x(nxy)+(x(nxy)-x(nxy-1))) ) then
            iOK = -2
         else
            iOK = -1
         end if
         j = nxy-1
      end if

C     now interpolate/extrapolate
      frac = (xin - x(j))/(x(j+1)-x(j))
      yout = (1.0-frac)*y(j) + frac*y(j+1)


      return
      end
C     end subroutine linear ---------------------------------------------

C     -----NUMERICAL RECIPES routines: locate and polint
      SUBROUTINE LOCATE(XX,N,X,J)
      DIMENSION XX(N)
      JL=0
      JU=N+1
10    IF(JU-JL.GT.1)THEN
        JM=(JU+JL)/2
        IF((XX(N).GT.XX(1)).EQV.(X.GT.XX(JM)))THEN
          JL=JM
        ELSE
          JU=JM
        ENDIF
      GO TO 10
      ENDIF
      J=JL
      RETURN
      END

      SUBROUTINE POLINT(XA,YA,N,X,Y,DY)
      PARAMETER (NMAX=10) 
      DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
      NS=1
      DIF=ABS(X-XA(1))
      DO 11 I=1,N 
        DIFT=ABS(X-XA(I))
        IF (DIFT.LT.DIF) THEN
          NS=I
          DIF=DIFT
        ENDIF
        C(I)=YA(I)
        D(I)=YA(I)
11    CONTINUE
      Y=YA(NS)
      NS=NS-1
      DO 13 M=1,N-1
        DO 12 I=1,N-M
          HO=XA(I)-X
          HP=XA(I+M)-X
          W=C(I+1)-D(I)
          DEN=HO-HP
          IF(DEN.EQ.0.)PAUSE
          DEN=W/DEN
          D(I)=HP*DEN
          C(I)=HO*DEN
12      CONTINUE
        IF (2*NS.LT.N-M)THEN
          DY=C(NS+1)
        ELSE
          DY=D(NS)
          NS=NS-1
        ENDIF
        Y=Y+DY
13    CONTINUE
      RETURN
      END
