program do_fe

*     calculates spectral indices, fe and errors for segue spectra

*     error treatment uses segue errors for each pixel and works
*     out first bit analytically, then monte-carlos bits that need
*     approximation (division to produce index and index-to-fe)
*     also uses trimsig for more robust estimates now

*     measures spectral indices for spectra from various  runs for 
*     luminosity and metallicity calibration

*     includes fixes to CaI and UV CN indices to make right 5/00
*     includes proper treatment of the fact that the iraf wavelength 
*     is the start of the pixel, not the middle, 9/05

      implicit none

      real c,fe_kp,fe_mg,fe_ca,fe_kp2,fe_ca2,fe_mg2,sig,lamvac
      real lam, val,lam2,lam1,sumred1,sumblue1,sumfeat1
      real sumred,sumblue,lamtest1,lamtest2,frac
      real sumfeat6,sumfeat12,sumfeat18,sumfeat
      real kprime6err,kprime12err,kprime18err,kerr
      real kprime(6), kprime12(2), kprime18(2),uvcn(4)
      real ca1(6),  mg(6),vel,k_monte,ca_monte,mg_monte
      real caindex,k,mgbh,dellam,gr0,gr0err,newgr0,fe,g,ra,dec,s
      real value(4000), error(4000), lambda(4000), gasdev, ran
      real kplam(500),calam(500),mglam(500)
      real feksn(1000),fecasn(1000),femgsn(1000),mean,sigk,sigmg,sigca
      real crval1,cd1_1,mgerr,ca1err,featerr,trimmed(900),order(1000)
      real vhelio, verr,inst_col,vmag,uvcnindex,cnerr

      integer i,j,idum,nn,i1,i2,i3,ii,ncols,nullval,group,nbuffer
      integer status,unitin,unitout,readwrite,blocksize,naxes(2),nfound
      integer ntrim

      logical anynull

      character*23 infile
      character*42 comment,errfile,fitsfile

      external gasdev

      data c / 3.0e5 /

      data uvcn / 3846.,3881.5, 3884.5, 3912. /
      data kprime / 3908., 3918., 3930.7, 3936.7, 4010., 4025. /
      data kprime12 / 3927.7, 3939.7 /
      data kprime18 / 3924.7,3942.7 /
      data ca1 / 4147., 4164., 4221.7, 4231.7, 4240., 4247. /
      data mg / 4935., 4975., 5130., 5200., 5217., 5258/

* input file has location in all difft run dirs plus color and velocity

      open(11,file='indices.in',status='old')      

      open(22,file='dofe.out',status='unknown')

      idum = -81

      write(22,10)
 10   format('star  g-r_0  UVCN  err  kprime   err   [Fe/H]K    err  '
     :   '    CaI   err   ' 
     :    '  [Fe/H]Ca  err   Mg index  err  [Fe/H]Mg  err  ')


      do j=1,500
         read(11,25,end=99) infile,vel, gr0,g,ra,dec
 25      format(a13,f14.3,f13.6,f13.4,f13.3,f13.4)

         vel = vel / c
         if(j.eq.1) then
            print*,gr0
         end if
*     read in data

*     it reads data into the array "value" and the errors into array "error"

         open(3,file=infile,status='old')

         do i=1,8
            read(3,*)
         end do
 
         do i=1,5000

* read in lambda , value, error and convert lambda to wavelengths in air

           read(3,*,end=98) lamvac, value(i),error(i)
           lambda(i) = lamvac/(1.0+2.735182E-4+131.4182/lamvac**2 +
     :                               2.76249E8 / lamvac**4)

*   AIR = VAC / (1.0 + 2.735182E-4 + 131.4182 / VAC^2 + 2.76249E8 / VAC^4)

         end do

 98      close(3)
         ncols = i - 1

*     transform to rest lambda. 

         do i=1,ncols
            lambda(i) = lambda(i)*(1-vel)
         end do

*     then call each subroutine to calculate index for this chunk        

* violet CN:

         call accumcn(lambda,value,error,ncols,uvcn,uvcnindex,idum,cnerr)

* K prime:
         call accumkp(lambda,value,error,ncols,kprime,kprime12,kprime18,
     :                            k,idum,kerr)

*********NB**********
******* I have not at this point changed the subroutines that calc [Fe/H]
******* to deal with g-r rather than m-t2
******* need to do this when photometry calibration is done
*********NB**********

         call kprime_interp(gr0, k, fe_kp)
         
*     CaI:
         call accumca(lambda,value,error,ncols,ca1,caindex,idum,ca1err)
         
         call ca1_interp(gr0, caindex, fe_ca)
         
*     Mg:
         
         call accummg(lambda,value,error,ncols,mg,mgbh,idum,mgerr)
         
         call mg2_interp(gr0, mgbh, fe_mg)      
         
         if(j.eq.1) then
            print*,gr0
         end if
*     monte carlo calculations for errors
         
         do i=1,1000
            
            ran = gasdev(idum)*kerr
            k_monte = k + ran
            newgr0 = gr0 + gasdev(idum)*gr0err
            call kprime_interp(newgr0, k_monte, fe_kp2)
            feksn(i) = fe_kp2

            ran = gasdev(idum)* ca1err
            ca_monte = caindex + ran
            newgr0 = gr0 + gasdev(idum)*gr0err
            call ca1_interp(newgr0, ca_monte, fe_ca2)
            fecasn(i) = fe_ca2
            
            ran = gasdev(idum)* mgerr
            mg_monte = mgbh + ran
            newgr0 = gr0 + gasdev(idum)*gr0err
            call mg2_interp(newgr0, mg_monte, fe_mg2)
            femgsn(i) = fe_mg2
         end do

         if(j.eq.1) then
            print*,gr0
         end if

* work out trimmed sigma so rare flyers dont bias it too much

         ntrim=50

         call sort(1000,feksn,order)
         do i=1,900
            trimmed(i) = order(i+ntrim)
         end do

         call stdev(trimmed,900,mean,sig)
         sigk = sig/0.7894


         call sort(1000,fecasn,order)
         do i=1,900
            trimmed(i) = order(i+ntrim)
         end do

         call stdev(trimmed,900,mean,sig)
         sigca = sig/0.7894


         call sort(1000,femgsn,order)
         do i=1,900
            trimmed(i) = order(i+ntrim)
         end do

         call stdev(trimmed,900,mean,sig)
         sigmg = sig/0.7894



*         write(22,15) infile(6:13), gr0, uvcnindex,cnerr, k, kerr,fe_kp,sigk, 
*     :           caindex, ca1err, fe_ca,sigca,
*     :            mgbh, mgerr, fe_mg, sigmg, g,ra,dec,s
*         
* 15      format(a8, 15f8.3, 5f9.3)

         write(22,15) infile(6:13), gr0, uvcnindex,cnerr, k, kerr, 
     :           caindex, ca1err, 
     :            mgbh, mgerr,g,ra,dec
         
 15      format(a8, 13f8.3, 5f9.3)

 19      continue
      end do

 99   close(22)
      close(11)

      end

      subroutine accumcn(lambda,value,error,nsize,uvcn,uvcnindex,idum,cnerr)

      implicit none

      real lam, val,lam2,lam1,del,bluerr,rederr,top,bottom,mid,alpha
      real lamtest1,lamtest2,frac,conterr,ans,ca1err
      real sumred,sumblue,sumrederr,sumbluerr,err,sumfeaterr
      real sumfeat,featerr
      real uvcn(4),uvcnindex,cnerr
      real dellam
      real value(4000), lambda(4000),error(4000)
      real cont,feat,index(100),mean
      real gasdev

      integer nsize, i,idum

      external gasdev

      sumred=0.
      sumrederr=0.
      sumfeat=0.
      sumfeaterr=0.


      do i=1,nsize
         val = value(i)
         lam = lambda(i)
         err = error(i)
         if(i.gt.1) then
            del = lambda(i)-lambda(i-1)
         else
            del = 0.873
         end if
         
* area under curve for violet CN bands

         lamtest1 = lam + del/2.
         lamtest2 = lam - del/2.

         if((lamtest1).gt.uvcn(1)) then
            if((lamtest2).lt.uvcn(2)) then
               call fraction(lamtest1,lamtest2,uvcn(1),uvcn(2),del,frac)
               sumfeat = sumfeat + frac*val*del
               sumfeaterr = sumfeaterr + (err*del)**2
            else if((lamtest1).gt.uvcn(3).and.(lamtest2).lt.uvcn(4)) then
               call fraction(lamtest1,lamtest2,uvcn(3),uvcn(4),del,frac)
               sumred = sumred + frac*val*del
               sumrederr = sumrederr + (err*del)**2
            end if
         end if

      end do

* calculate error in continuum level

      conterr = sqrt(sumrederr) 

* now for Cn band:

      featerr = sqrt(sumfeaterr)

*calculate index value

      uvcnindex = -2.5 * alog10(sumfeat/sumred)

*     now calc error on index with a monte carlo thingy

      do i=1,100
         cont = sumred + gasdev(idum)*conterr
         feat = sumfeat + gasdev(idum)*featerr
         index(i) = -2.5 * alog10(abs(feat/cont))
      end do

      call stdev(index,100,mean,cnerr)
      
      return
      end


      subroutine accumkp(lambda,value,error,nsize,kprime,kprime12,kprime18,
     :              k,idum,kerr)

      implicit none

      real lam, val,lam2,lam1,del,bluerr,rederr,top,bottom,mid,alpha
      real lamtest1,lamtest2,frac,conterr,ans,kprime6err,kprime12err,kprime18err
      real sumred,sumblue,sumrederr,sumbluerr,err,sumfeat6err
      real sumfeat12err,sumfeat18err
      real sumfeat6,sumfeat12,sumfeat18
      real kprime(6), kprime12(2), kprime18(2)
      real k,dellam,mt2,fe,mt2err
      real value(4000), lambda(4000),error(4000)
      real cont,feat,index1(100),index2(100),index3(100),mean,featerr
      real gasdev,kerr
      real k6,k12,k18

      integer nsize, i,idum

      external gasdev

      sumblue=0.
      sumbluerr=0.
      sumred=0.
      sumrederr=0.
      sumfeat6=0.
      sumfeat6err=0.
      sumfeat12=0.
      sumfeat12err=0.
      sumfeat18=0.
      sumfeat18err=0.

      do i=1,nsize
         val = value(i)
         lam = lambda(i)
         err = error(i)
         if(i.gt.1) then
            del = lambda(i)-lambda(i-1)
         else
            del = 0.873
         end if

         lamtest1 = lam + del/2.
         lamtest2 = lam - del/2.

         if((lamtest1).gt.kprime(1)) then
            if((lamtest2).lt.kprime(2)) then
               call fraction(lamtest1,lamtest2,kprime(1),kprime(2),del,frac)
               sumblue = sumblue + frac*val*del/(kprime(2)-kprime(1))
               sumbluerr = sumbluerr + (err*del/(kprime(2)-kprime(1)))**2
            else 
               if((lamtest1).gt.kprime(3).and.(lamtest2).lt.kprime(4)) then
                  call fraction(lamtest1,lamtest2,kprime(3),kprime(4),del,frac)
                  sumfeat6 = sumfeat6 + frac*val*del/(kprime(4)-kprime(3))
                  sumfeat6err = sumfeat6err + (err*del/(kprime(4)-kprime(3)))**2
               end if
               if((lamtest1).gt.kprime12(1).and.(lamtest2).lt.kprime12(2)) then
                  call fraction(lamtest1,lamtest2,kprime12(1),kprime12(2),del,frac)
                  sumfeat12 = sumfeat12 + frac*val*del/(kprime12(2)-kprime12(1))
                  sumfeat12err = sumfeat12err + (err*del/(kprime12(2)-kprime12(1)))**2
               end if
               if((lamtest1).gt.kprime18(1).and.(lamtest2).lt.kprime18(2)) then
                  call fraction(lamtest1,lamtest2,kprime18(1),kprime18(2),del,frac)
                  sumfeat18 = sumfeat18 + frac*val*del/(kprime18(2)-kprime18(1))
                  sumfeat18err = sumfeat18err + (err*del/(kprime18(2)-kprime18(1)))**2
               end if
            end if
            if((lamtest1).gt.kprime(5).and.(lamtest2).lt.kprime(6)) then
               call fraction(lamtest1,lamtest2,kprime(5),kprime(6),del,frac)
               sumred = sumred + frac*val*del/(kprime(6)-kprime(5))
               sumrederr = sumrederr + (err*del/(kprime(6)-kprime(5)))**2
            end if
         end if
         
      end do

*     calculate error in continuum level

      bluerr = sqrt(sumbluerr)
      rederr = sqrt(sumrederr)
      top = (kprime(5) + kprime(6))/2.
      bottom = (kprime(1) + kprime(2))/2.
      mid = (kprime(3)+kprime(4))/2.
      alpha = (mid-bottom) / (top-bottom)

      conterr = sqrt( (bluerr*(1-alpha))**2 + (rederr*alpha)**2)

*     now calc EW 

      call interp(kprime,sumblue,sumred,ans)

      k6  =  6*(1- sumfeat6/ans)
      k12 = 12*(1-sumfeat12 / ans)
      k18 = 18*(1-sumfeat18 / ans)

* now monte carlo to get estimate of error in EW

      if(k6.le.2) then
         k = k6
         do i=1,100
            cont = ans + gasdev(idum)*conterr
            feat = sumfeat6 + gasdev(idum)*sqrt(sumfeat6err)
            index3(i) = 6.*(1-feat/cont)
         end do
         call stdev(index3,100,mean,kerr)

      else if(k12.le.5.) then
         k = k12
         do i=1,100
            cont = ans + gasdev(idum)*conterr
            feat = sumfeat12 + gasdev(idum)*sqrt(sumfeat12err)
            index2(i) = 12.*(1-feat/cont)
         end do
         call stdev(index2,100,mean,kerr)

      else
         k = k18
         do i=1,100
            cont = ans + gasdev(idum)*conterr
            feat = sumfeat18 + gasdev(idum)*sqrt(sumfeat18err)
            index1(i) = 18.*(1-feat/cont)
         end do
         call stdev(index1,100,mean,kerr)
      end if

      return
      end


      subroutine accumca(lambda,value,error,nsize,ca1,caindex,idum,
     :                       ca1err)

      implicit none

      real lam, val,lam2,lam1,del,bluerr,rederr,top,bottom,mid,alpha
      real lamtest1,lamtest2,frac,conterr,ans,ca1err
      real sumred,sumblue,sumrederr,sumbluerr,err,sumfeaterr
      real sumfeat,featerr
      real ca1(6),caindex
      real dellam
      real value(4000), lambda(4000),error(4000)
      real cont,feat,index(100),mean
      real gasdev

      integer nsize, i,idum

      external gasdev


      sumred=0.
      sumrederr=0.
      sumfeat=0.
      sumfeaterr=0.

      do i=1,nsize
         val = value(i)
         lam = lambda(i)
         err = error(i)
         if(i.gt.1) then
            del = lambda(i)-lambda(i-1)
         else
            del = 0.873
         end if

         lamtest1 = lam + del/2.
         lamtest2 = lam - del/2.

         if((lamtest1).gt.ca1(1)) then
            if((lamtest1).gt.ca1(3).and.(lamtest2).lt.ca1(4)) then
               call fraction(lamtest1,lamtest2,ca1(3),ca1(4),del,frac)
               sumfeat = sumfeat + frac*val*del
               sumfeaterr = sumfeaterr + (err*del)**2
            else if((lamtest1).gt.ca1(5).and.(lamtest2).lt.ca1(6)) then
               call fraction(lamtest1,lamtest2,ca1(5),ca1(6),del,frac)
               sumred = sumred + frac*val*del
               sumrederr = sumrederr + (err*del)**2
            end if
         end if

      end do

* calculate error in continuum level

      conterr = sqrt(sumrederr) 

* now for Ca band:

      featerr = sqrt(sumfeaterr)

*calculate index value

      caindex = -2.5 * alog10(sumfeat/sumred)

*     now calc error on index with a monte carlo thingy

      do i=1,100
         cont = sumred + gasdev(idum)*conterr
         feat = sumfeat + gasdev(idum)*featerr
         index(i) = -2.5 * alog10(abs(feat/cont))
      end do

      call stdev(index,100,mean,ca1err)
      
      return
      end


      subroutine accummg(lambda,value,error,nsize,mg,mgbh,idum,sig)

      implicit none

      real lam, val,lam2,lam1,del,bluerr,rederr,top,bottom,mid,alpha
      real lamtest1,lamtest2,frac,conterr,ans,featerr
      real sumred,sumblue,sumrederr,sumbluerr,err,sumfeaterr
      real sumfeat,differr,diffpercen
      real mg(6),mgbh
      real dellam,mt2,fe,mt2err
      real value(4000), lambda(4000),error(4000)
      real cont,feat,index(100),mean,sig,realmgbh
      real gasdev

      integer nsize, i,idum

      external gasdev

      sumblue=0.
      sumbluerr=0.
      sumred=0.
      sumrederr=0.
      sumfeat=0.
      sumfeaterr=0.

      do i=1,nsize
         val = value(i)
         lam = lambda(i)
         err = error(i)
         if(i.gt.1) then
            del = lambda(i)-lambda(i-1)
         else
            del = 0.873
         end if
         
         lamtest1 = lam + del/2.
         lamtest2 = lam - del/2.

         if(lamtest1.gt.mg(1)) then
            if((lamtest2).lt.mg(2)) then
               call fraction(lamtest1,lamtest2,mg(1),mg(2),del,frac)
               sumblue = sumblue + frac*val*del/(mg(2)-mg(1))
               sumbluerr = sumbluerr + (err*del/(mg(2)-mg(1)))**2
            else if((lamtest1).gt.mg(3).and.(lamtest2).lt.mg(4)) then
               call fraction(lamtest1,lamtest2,mg(3),mg(4),del,frac)
               sumfeat = sumfeat + frac*val*del/(mg(4)-mg(3))
               sumfeaterr = sumfeaterr + (err*del/(mg(4)-mg(3)))**2
            else if((lamtest1).gt.mg(5).and.(lamtest2).lt.mg(6)) then
               call fraction(lamtest1,lamtest2,mg(5),mg(6),del,frac)
               sumred = sumred + frac*val*del/(mg(6)-mg(5))
               sumrederr = sumrederr + (err*del/(mg(6)-mg(5)))**2
            end if
         end if
      end do
            
* calculate error in continuum level

      bluerr = sqrt(sumbluerr)
      rederr = sqrt(sumrederr)
      top = (mg(5) + mg(6))/2.
      bottom = (mg(1) + mg(2))/2.
      mid = (mg(3)+mg(4))/2.
      alpha = (mid-bottom) / (top-bottom)

      conterr = sqrt( (bluerr*(1-alpha))**2 + (rederr*alpha)**2)
      featerr = sqrt(sumfeaterr)

* calc index

      call interp(mg,sumblue,sumred,ans)
      mgbh  =  (mg(4)-mg(3))*(1 - sumfeat/ans)


*     now calc error on EW with a monte carlo thingy

      do i=1,100
         cont = ans + gasdev(idum)*conterr
         feat = sumfeat + gasdev(idum)*featerr
         index(i) = (mg(4)-mg(3))*(1-feat/cont)
      end do
      call stdev(index,100,mean,sig)

      return
      end


      subroutine fraction(lamtest1,lamtest2,lhs,rhs,del,frac)

      real lamtest1,lamtest2,lhs,rhs,del,frac

      if(lamtest1.gt.lhs.and.lamtest2.lt.rhs) then
         frac = 1.0
      end if

      if(lamtest1.gt.lhs.and.lamtest2.lt.lhs) then
         frac = (lamtest1-lhs)/del
      end if

      if(lamtest1.gt.rhs.and.lamtest2.lt.rhs) then
         frac = (rhs - lamtest2)/del
      end if

      return
      end



      subroutine interp(band,sumblue,sumred,ans)
            
      real band(6),top,bottom,mid,ans

      top = (band(5) + band(6))/2.
      bottom = (band(1) + band(2))/2.
      mid = (band(3)+band(4))/2.

      ans = sumblue + (mid - bottom)* (sumred-sumblue)/(top-bottom)
      
      return
      end



      subroutine kprime_interp(mt2,index,fe)

c     interpolate in globular cluster grid to find metallicity
c     with Kprime index

      implicit none
      
      integer npts,i

      parameter (npts = 10)

      real hix(npts),hiy(npts),medx(npts),medy(npts),lowx(npts),lowy(npts)
      real vlowx(npts),vlowy(npts)
      real midx(npts),midy(npts),femid, lowishx(npts),lowishy(npts),felowish
      real mt2, indexlow, indexmed, indexhi, index, gradient, fe
      real fehi, femed, felow, fevlow, indexvlow,indexmid,indexlowish

*     Lucy's new values from WHT run May 2005:


*     (-0.7)
      data hix / 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8 /
      data hiy / 9.9, 10.4, 10.9, 11.4, 11.7, 11.9, 12.0, 12.1, 12.3, 12.5  /
      data fehi / -0.7 /

*     (-1.4)
      data midx / 0.9, 1.0, 1.1, 1.2, 1.3,  1.4, 1.5, 1.6, 1.7, 1.8 /
      data midy / 8.7, 9.4, 10., 10.6,11.2, 11.4,11.45,11.5,11.55,11.6  /
      data femid / -1.4 /

*     (-2.0)
      data lowx / 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8 /
      data lowy / 5.1, 5.6, 6.1, 6.8, 7.7, 6.6, 9.2, 9.5, 10.1,10.7 /
      data felow / -2.0 /

*     synthetic -3.0 isochrone 
      data vlowx / 0.9, 1.0, 1.108, 1.282,1.411,1.46,1.536, 1.6, 1.7, 1.8 /
      data vlowy / 3.,3.85, 4.929, 5.996, 6.704, 6.9,7.156, 7.6, 8.0, 8.4 /
      data fevlow / -3.0 /

*     47tuc (-0.7):  old values
*     data hix / 1.1,  1.2,  1.31,  1.44, 1.6 /
*     data hiy / 10.5, 11.0, 11.6, 12.1, 12.6 /
*     data fehi / -0.7 /

*     N1851(-1.08)  old values
*     data medx / 1.1, 1.21, 1.31, 1.6, 1.6 /
*     data medy / 9.6, 10.6, 11.1, 12.0, 12.0/
*     data femed / -1.08 /

*     M2/M3(-1.35)  old values
*     data midx / 1.1, 1.2, 1.26, 1.45, 1.6 /
*     data midy / 8.9,  9.6,9.89, 10.76,10.5 /
*     data femid / -1.35 /

*     N6397(-1.82)  old values
*     data lowishx / 1.1, 1.2, 1.31, 1.4, 1.6 /
*     data lowishy / 6.9, 7.55, 8.34, 8.63, 9.3 /
*     data felowish / -1.82 /

*     N4590(-1.99)  old values
*     data lowx / 1.1, 1.13, 1.43, 1.6, 1.6 /
*     data lowy / 6.46, 6.6, 8.46, 9.56, 9.0 /
*     data felow / -1.99 /

*     synthetic -3.0 isochrone  old values
*     data vlowx / 1.108, 1.282, 1.411, 1.536, 1.6 /
*     data vlowy / 2.529, 3.596, 4.304, 4.756, 5.2 /
*     data fevlow / -3.0 /

c     find index value for five clusters at this mt2 value
*     extrapolate outside m=t2 range

      if(mt2.lt.vlowx(1)) then
         gradient = (vlowy(2)-vlowy(1))/(vlowx(2)-vlowx(1)) 
         indexvlow = vlowy(1) + gradient*(mt2-vlowx(1))
      else
         do i=2,npts
            if(mt2.lt.vlowx(i)) then
               gradient = (vlowy(i)-vlowy(i-1))/(vlowx(i)-vlowx(i-1)) 
               indexvlow = vlowy(i-1)+gradient*(mt2-vlowx(i-1))
               go to 10
            end if
         end do
         gradient = (vlowy(npts)-vlowy(npts-1))/(vlowx(npts)-vlowx(npts-1))
         indexvlow = vlowy(npts)+gradient*(mt2-vlowx(npts))          
      end if
      
 10   if(mt2.lt.lowx(1)) then
         gradient = (lowy(2)-lowy(1))/(lowx(2)-lowx(1)) 
         indexlow = lowy(1) + gradient*(mt2-lowx(1))
      else
         do i=2,npts
            if(mt2.lt.lowx(i)) then
               gradient = (lowy(i)-lowy(i-1))/(lowx(i)-lowx(i-1)) 
               indexlow = lowy(i-1)+gradient*(mt2-lowx(i-1))
               go to 15
            end if
         end do
         gradient = (lowy(npts)-lowy(npts-1))/(lowx(npts)-lowx(npts-1))
         indexvlow = lowy(npts)+gradient*(mt2-lowx(npts))          
      end if
      
 15   if(mt2.lt.midx(1)) then
         gradient = (midy(2)-midy(1))/(midx(2)-midx(1)) 
         indexlow = midy(1) + gradient*(mt2-midx(1))
      else
         do i=2,npts
            if(mt2.lt.midx(i)) then
               gradient = (midy(i)-midy(i-1))/(midx(i)-midx(i-1)) 
               indexmid = midy(i-1)+gradient*(mt2-midx(i-1))
               go to 25
            end if
         end do
         gradient = (midy(npts)-midy(npts-1))/(midx(npts)-midx(npts-1))
         indexvlow = midy(npts)+gradient*(mt2-midx(npts))          
      end if

 25   if(mt2.lt.hix(1)) then
         gradient = (hiy(2)-hiy(1))/(hix(2)-hix(1)) 
         indexlow = hiy(1) + gradient*(mt2-hix(1))
      else
         do i=2,npts
            if(mt2.lt.hix(i)) then
               gradient = (hiy(i)-hiy(i-1))/(hix(i)-hix(i-1)) 
               indexhi = hiy(i-1)+gradient*(mt2-hix(i-1))
               go to 40
            end if
         end do
         gradient = (hiy(npts)-hiy(npts-1))/(hix(npts)-hix(npts-1))
         indexvlow = hiy(npts)+gradient*(mt2-hix(npts))          
      end if


*     now interpolate to get metallicity estimate

 40      if(index.le.indexlow) then
            gradient = (felow - fevlow)/(indexlow-indexvlow)
            fe = fevlow + gradient * (index- indexvlow)
         else if(index.le.indexmid) then
            gradient = (femid - felow)/(indexmid-indexlow)
            fe = felow + gradient * (index- indexlow)
         else
            gradient = (fehi - femid)/(indexhi-indexmid)
            fe = femid + gradient * (index- indexmid)
         end if

         return
         end

      subroutine ca1_interp(mt2,index,fe)

c     interpolate in globular cluster grid to find metallicity
c     with CaI index

      implicit none

      integer npts,i

      parameter (npts = 6)

      real hix(npts),hiy(npts),medx(npts),medy(npts)
      real lowx(npts),lowy(npts),vlowx(npts),vlowy(npts)
      real midx(npts),midy(npts),femid
      real mt2, indexlow, indexmed, indexhi, index, gradient, fe
      real fehi, femed, felow, fevlow, indexvlow,indexmid


*     (-0.7)
      data hix / 0.9, 1.0, 1.2, 1.4, 1.6, 1.8 /
      data hiy / -0.23, -0.21, -0.185, -0.14, -0.1, -0.07 /
      data fehi / -0.7 /

*     (-1.2)
      data medx / 0.9, 1.0, 1.2, 1.4, 1.6, 1.8 /
      data medy / -0.33, -0.32, -0.29, -0.25, -0.2, -0.15 /
      data femed / -1.2 /

*     (-2.0)
      data lowx / 0.9, 1.0, 1.2, 1.4, 1.6, 1.8 /
      data lowy / -0.36, -0.36, -0.36, -0.35, -0.31, -0.25 /
      data felow / -2.0 /

*     synthetic -3.0 isochrone 
      data vlowx / 0.9,  1.08, 1.282, 1.411, 1.536,  1.8 /
      data vlowy / -0.365,  -0.36, -0.355, -0.35, -0.345, -0.34 /
      data fevlow / -3.0 /

*     47tuc (-0.7): old values
*     data hix / 1.0, 1.1,  1.2,  1.31,  1.44, 1.6 /
*     data hiy / -0.3,-0.3, -0.25, -0.18, -0.14, -0.09 /
*     data fehi / -0.7 /

*     N1851(-1.08) old values
*     data medx / 1.0,  1.1,   1.21,   1.31,  1.6,  1.6 /
*     data medy / -0.33,-0.33, -0.32, -0.29, -0.22, -0.22 /
*     data femed / -1.08 /

*     M2/M3(-1.35) (not used at present; CaI changes slowly here old values
*     data midx / 1.0,  1.1,  1.2, 1.26, 1.45, 1.6 /
*     data midy / -.34,-.34,-.36,-.32, -.30, -.26 /
*     data femid / -1.35 /

*     N4590(-1.99) old values
*     data lowx / 1.0,  1.1,  1.13,   1.43,   1.6,   1.6 /
*     data lowy / -0.35,-0.35, -0.35, -0.323, -0.323, -0.323 /
*     data felow / -1.99 /

*     synthetic -3.0 isochrone old values
*     data vlowx / 1.0,  1.108, 1.282, 1.411, 1.536, 1.6 /
*     data vlowy / -0.37, -0.37, -0.365, -0.36, -0.355, -0.355 /
*     data fevlow / -3.0 /

      if(mt2.lt.vlowx(1)) then
         gradient = (vlowy(2)-vlowy(1))/(vlowx(2)-vlowx(1)) 
         indexvlow = vlowy(1) + gradient*(mt2-vlowx(1))
      else
         do i=2,npts
            if(mt2.lt.vlowx(i)) then
               gradient = (vlowy(i)-vlowy(i-1))/(vlowx(i)-vlowx(i-1)) 
               indexvlow = vlowy(i-1)+gradient*(mt2-vlowx(i-1))
               go to 10
            end if
         end do
         gradient = (vlowy(npts)-vlowy(npts-1))/(vlowx(npts)-vlowx(npts-1))
         indexvlow = vlowy(npts)+gradient*(mt2-vlowx(npts))          
      end if


 10   if(mt2.lt.lowx(1)) then
         gradient = (lowy(2)-lowy(1))/(lowx(2)-lowx(1)) 
         indexlow = lowy(1) + gradient*(mt2-lowx(1))
      else
         do i=2,npts
            if(mt2.lt.lowx(i)) then
               gradient = (lowy(i)-lowy(i-1))/(lowx(i)-lowx(i-1)) 
               indexlow = lowy(i-1)+gradient*(mt2-lowx(i-1))
               go to 20
            end if
         end do
         gradient = (lowy(npts)-lowy(npts-1))/(lowx(npts)-lowx(npts-1))
         indexlow = lowy(npts)+gradient*(mt2-lowx(npts))          
      end if
         
 20   if(mt2.lt.medx(1)) then
         gradient = (medy(2)-medy(1))/(medx(2)-medx(1)) 
         indexmed = medy(1) + gradient*(mt2-medx(1))
      else
         do i=2,npts
            if(mt2.lt.medx(i)) then
               gradient = (medy(i)-medy(i-1))/(medx(i)-medx(i-1)) 
               indexmed = medy(i-1)+gradient*(mt2-medx(i-1))
               go to 30 
            end if
         end do
         gradient = (medy(npts)-medy(npts-1))/(medx(npts)-medx(npts-1))
         indexmed = medy(npts)+gradient*(mt2-medx(npts))          
      end if

 30   if(mt2.lt.hix(1)) then
         gradient = (hiy(2)-hiy(1))/(hix(2)-hix(1)) 
         indexhi = hiy(1) + gradient*(mt2-hix(1))
         do i=2,6
            if(mt2.lt.hix(i)) then
               gradient = (hiy(i)-hiy(i-1))/(hix(i)-hix(i-1)) 
               indexhi = hiy(i-1)+gradient*(mt2-hix(i-1))
               go to 40
            end if
         end do
         gradient = (hiy(npts)-hiy(npts-1))/(hix(npts)-hix(npts-1))
         indexhi = hiy(npts)+gradient*(mt2-hix(npts))          
      end if


*     now interpolate to get metallicity estimate

 40      if(index.le.indexlow) then
            gradient = (felow - fevlow)/(indexlow-indexvlow)
            fe = fevlow + gradient * (index- indexvlow)
         else if(index.le.indexmed) then
            gradient = (femed - felow)/(indexmed-indexlow)
            fe = felow + gradient * (index- indexlow)
         else
            gradient = (fehi - femed)/(indexhi-indexmed)
            fe = femed + gradient * (index- indexmed)
         end if

         return
         end

      subroutine mg2_interp(mt2,index,fe)

c     interpolate in globular cluster grid to find metallicity
c     with new Mg index; no synthetic data yet

      implicit none

      integer npts,i

      parameter (npts = 8)

      real hix(npts),hiy(npts),medx(npts),medy(npts)
      real lowx(npts),lowy(npts),vlowx(npts),vlowy(npts)
      real midx(npts),midy(npts),femid
      real mt2, indexlow, indexmed, indexhi, index, gradient, fe
      real fehi, femed, felow, fevlow, indexvlow,indexmid


*     (-0.75)
      data hix / 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.6, 1.8  / 
      data hiy / 4.5, 4.8, 5.6, 6.5, 6.7, 7.0, 8.5, 12.5  /
      data fehi / -0.75 /

*     (-1.4)
      data medx / 0.9, 1.1, 1.2,  1.3, 1.4, 1.5, 1.6, 1.8 /
      data medy / 3.0, 3.1, 3.2, 3.4, 3.45, 3.65, 4.2, 5.8 /
      data femed / -1.4 /

*     (-1.9)
      data lowx / 0.9, 1.0, 1.1, 1.2, 1.4, 1.5, 1.6, 1.8 /
      data lowy / 2.6, 2.7, 2.6, 2.5, 2.6, 2.8, 3.0, 3.4 /
      data felow / -1.9 /

*     new point for 47 Tuc at Mt2=1.6 from M71 wht data

*     47tuc (-0.7)/6171(-0.9; red end only): old values
*     data hix / 1.0, 1.07,  1.2,  1.31, 1.45, 1.6 / 
*     data hiy / 3.1, 3.1,  5.8,  6.5,  6.55, 7.5 /
*     data fehi / -0.8 /

*     N1851/M2/M3(-1.3) old values
*     data medx / 1.0, 1.04, 1.19, 1.28, 1.4, 1.75 /
*     data medy / 1.84, 1.84, 2.5,  3.,  3.44, 5.52 /
*     data femed / -1.3 /

*     M2/M3(-1.35) old values
*     data midx / 1.0, 1.1,  1.2, 1.26, 1.45, 1.6 /
*     data midy / 1.9, 1.9, 2.7,  2.6,  3.3, 4.2 /
*     data femid / -1.35 /

*     N4590/6397(-1.9) old values
*     data lowx / 1.0, 1.13, 1.20, 1.31, 1.40, 1.6 /
*     data lowy / 1.2, 1.3,  1.45, 1.6, 1.7, 3.0 /
*     data felow / -1.9 /

*     synthetic -3.0 isochrone old values
*     data vlowx / 1.108, 1.282, 1.411, 1.536, 1.6 /
*     data vlowy / -0.37, -0.365, -0.36, -0.355, -0.355 /
*     data fevlow / -3.0 /

c     find index value for three clusters at this mt2 value
*     extrapolate outside m=t2 = 1.1, 1.6; almost never a real problem

*     currently dont have synthetic vals for Mg
*     do i=2,5
*     if(mt2.lt.vlowx(i)) then
*     gradient = (vlowy(i)-vlowy(i-1))/(vlowx(i)-vlowx(i-1)) 
*     indexvlow = vlowy(i-1)+gradient*(mt2-vlowx(i-1))
*     go to 10
*     end if
*     end do
*     gradient = (vlowy(5)-vlowy(4))/(vlowx(5)-vlowx(4))
*     indexvlow = vlowy(4)+gradient*(mt2-vlowx(4))          

 10   if(mt2.lt.lowx(1)) then
         gradient = (lowy(2)-lowy(1))/(lowx(2)-lowx(1)) 
         indexlow = lowy(1) + gradient*(mt2-lowx(1))
      else
         do i=2,npts
            if(mt2.lt.lowx(i)) then
               gradient = (lowy(i)-lowy(i-1))/(lowx(i)-lowx(i-1)) 
               indexlow = lowy(i-1)+gradient*(mt2-lowx(i-1))
               go to 20
            end if
         end do
         gradient = (lowy(npts)-lowy(npts-1))/(lowx(npts)-lowx(npts-1))
         indexlow = lowy(npts)+gradient*(mt2-lowx(npts))          
      end if
         

 20   if(mt2.lt.medx(1)) then
         gradient = (medy(2)-medy(1))/(medx(2)-medx(1)) 
         indexmed = medy(1) + gradient*(mt2-medx(1))
      else
         do i=2,npts
            if(mt2.lt.medx(i)) then
               gradient = (medy(i)-medy(i-1))/(medx(i)-medx(i-1)) 
               indexmed = medy(i-1)+gradient*(mt2-medx(i-1))
               go to 30 
            end if
         end do
         gradient = (medy(npts)-medy(npts-1))/(medx(npts)-medx(npts-1))
         indexmed = medy(npts)+gradient*(mt2-medx(npts))          
      end if


 30   if(mt2.lt.hix(1)) then
         gradient = (hiy(2)-hiy(1))/(hix(2)-hix(1)) 
         indexmed = hiy(1) + gradient*(mt2-hix(1))
      else
         do i=2,npts
            if(mt2.lt.hix(i)) then
               gradient = (hiy(i)-hiy(i-1))/(hix(i)-hix(i-1)) 
               indexhi = hiy(i-1)+gradient*(mt2-hix(i-1))
               go to 40
            end if
         end do
         gradient = (hiy(npts)-hiy(npts-1))/(hix(npts)-hix(npts-1))
         indexhi = hiy(npts)+gradient*(mt2-hix(npts))          
      end if

*     now interpolate to get metallicity estimate

*     40   if(index.le.indexlow) then
*     gradient = (felow - fevlow)/(indexlow-indexvlow)
*     fe = fevlow + gradient * (index- indexvlow)
 40      if(index.le.indexmed) then
            gradient = (femed - felow)/(indexmed-indexlow)
            fe = felow + gradient * (index- indexlow)
         else
            gradient = (fehi - femed)/(indexhi-indexmed)
            fe = femed + gradient * (index- indexmed)
         end if

         return
         end


      subroutine stdev(xx,n,mean,sig)

      real*8  variance
      real*8 sumx,sumxsq,x(1000)
      real xx(1000),mean,sig
      integer n

      sumx = 0.
      sumxsq = 0.

      do i=1,n
         x(i) = dble(xx(i))
      end do

      do i=1,n        
         sumx = sumx + x(i)
         sumxsq = sumxsq + x(i)**2
      end do

      mean = sngl(sumx/float(n))
      variance = (sumxsq - ((sumx)**2)/float(n))/float(n-1)
      sig = sngl(sqrt(variance))


      return
      end

        Subroutine sort(n,inarray,outarray)

* from Blair Phillips, much better than Kavan's

        implicit none
        integer n
        real*4 inarray(n),outarray(n)
        integer i,sorted,ismall
        real*4 small        

        do i = 1,n
          outarray(i) = inarray(i)
        enddo
*
*  outarray(1..sorted-1) is sorted
*
        do sorted = 1,n-1
          ismall = sorted
          small = outarray(ismall)
          do i = sorted+1,n
            if (outarray(i).lt.small) then
              ismall = i
              small = outarray(ismall)
            endif
          enddo
          outarray(ismall) = outarray(sorted)
          outarray(sorted) = small
        enddo
        return
        end



      FUNCTION gasdev(idum)
      INTEGER idum
      REAL gasdev
CU    USES ran1
      INTEGER iset
      REAL fac,gset,rsq,v1,v2,ran3
      SAVE iset,gset
      DATA iset/0/
      if (idum.lt.0) iset=0
      if (iset.eq.0) then
1       v1=2.*ran3(idum)-1.
        v2=2.*ran3(idum)-1.
        rsq=v1**2+v2**2
        if(rsq.ge.1..or.rsq.eq.0.)goto 1
        fac=sqrt(-2.*log(rsq)/rsq)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      else
        gasdev=gset
        iset=0
      endif
      return
      END

      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END

      FUNCTION ran2(idum)
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software #725.

      FUNCTION ran3(idum)
      INTEGER idum
      INTEGER MBIG,MSEED,MZ
C     REAL MBIG,MSEED,MZ
      REAL ran3,FAC
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
C     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
      INTEGER i,iff,ii,inext,inextp,k
      INTEGER mj,mk,ma(55)
C     REAL mj,mk,ma(55)
      SAVE iff,inext,inextp,ma
      DATA iff /0/
      if(idum.lt.0.or.iff.eq.0)then
        iff=1
        mj=MSEED-iabs(idum)
        mj=mod(mj,MBIG)
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.MZ)mk=mk+MBIG
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
12        continue
13      continue
        inext=0
        inextp=31
        idum=1
      endif
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.MZ)mj=mj+MBIG
      ma(inext)=mj
      ran3=mj*FAC
      return
      END
C     (C) Copr. 1986-92 Numerical Recipes Software #725.