      program ask

C                  histpar is the number of bins in the [Fe/H] histogram
      integer      histpar
      parameter    (histpar = 28)
      real         metalref(histpar),ageref(histpar),shist(histpar)
      character*80 text


      print*, ' '
      print*, ' '
      print*, 'ASK returns a flux, colors, SBF mags, and indices'
      print*, 'for combinations of single-burst populations. '
      print*, 'This the Padova isochrone version. '
      print*, ' '
      print*, 'You can enter up to ',histpar,
     &        ' separate populations, along '
      print*, 'with weights (which the program will normalize so that '
      print*, 'the total population mass is a million Mo). '
      print*, 'Metallicities are restricted to -1.7 < [Fe/H] < +0.4,'
      print*, 'and ages to between 0.0039 and 19.9 Gyr.'
      print*, ' '
      print*, 'Enter trios of numbers separated by spaces: '
      print*, '( age  [Fe/H]  weight ). '
      print*, 'Enter a blank to exit the loop '
      print*, ' ' 

      do i=1,histpar
 11      print*, 'Enter (age [Fe/H] weight) no. ',i
         read(*,'(a80)') text
         if (text.eq.' ') then
            npop =  i-1
            goto 22
         end if
         read(text,*,err=99)  ageref(i),metalref(i),shist(i)
         goto 12
 99      print*, 'I had trouble reading that. Try again.'
         goto 11
 12   end do
      npop = 28

C     normalize mass histogram to integral one
 22   xsum = 0.0
      print*, ' npop = ',npop
      do j=1,npop
         xsum = xsum + shist(j)
      end do
      do j=1,npop
         shist(j) = shist(j)/xsum
      end do

      call asksub(metalref,ageref,shist,npop)

      stop 
      end

      subroutine asksub(metalref,ageref,shist,npop)
C                  histpar is the number of bins in the [Fe/H] histogram
C                  nind is the number of indices
C                  npop is the number of populations to combine
      integer      nind,histpar,npop
      parameter    (nind = 25, histpar = 28)
      real         metalref(histpar),ageref(histpar),shist(histpar)

C     primary passed variables
      real         hlam(1221),wave(1221), bolo, z, y, age
      real         mags(11),Lsq(11),Lavg(11)
      real         widths(nind),Fcont(nind),Fline(nind)
      integer      indexunits(nind)
      character*10 labels(nind)
      real         mass
C     other passed variables
      real         imfpar, zsolar, ysolar, yprimo, imflow, imfupp
      real         hbmorph, pagb, eta

C    
      real         zin, agein
      character*25 optdir
C     variables for storing things.
      real         ahlam(1221,histpar), abolo(histpar)
      real         amags(11,histpar)
      real         aLsq(11,histpar),aLavg(11,histpar)
      real         aFcont(nind,histpar),aFline(nind,histpar)
      real         amass(histpar),aind(nind,histpar)
C     variables for output
      real         index(nind), Lbar(11)

C     other variables
      character*24   ctime,t
      integer        time,status

C     variables for Spinrad's UV break indices
      real         sp(2,4), spavg(4)
      data      sp / 2600.0, 2630.0,    2645.0, 2675.0,
     z               2855.0, 2885.0,    2915.0, 2945.0   /

C     get and store the interpolated population models
      do i=1,npop
         zin = 0.0169*10**(metalref(i))
         agein = ageref(i)
c        optdir is 25 characters long, gives subdirectory for model grid
         optdir = 'padmods+00'

         call getmodel(zin,agein,optdir,
     &          hlam,wave,bolo,z,y,age,mags,Lsq,Lavg,
     &       widths,Fcont,Fline,indexunits,labels,mass,imfpar,zsolar,
     &       ysolar,yprimo,imflow,imfupp,hbmorph,pagb,eta,status)


C        store important results
         do j=1,1221
            ahlam(j,i) = hlam(j)
         end do
         abolo(i) = bolo
         amass(i) = mass
c         print*, ' bolo, mass = ',bolo,mass
         do j=1,11
            amags(j,i) = mags(j)
            aLsq(j,i) = Lsq(j)
            aLavg(j,i) = Lavg(j)
c            print*, ' mag, Lsq, Lavg = ',mags(j),Lsq(j),Lavg(j)
         end do
         do j=1,nind
            aFcont(j,i) = Fcont(j)
            aFline(j,i) = Fline(j)
            if (indexunits(j).eq. 1) then
               aind(j,i) = widths(j)*(1.0-Fline(j)/Fcont(j))
            else
               aind(j,i) = -2.5*alog10(Fline(j)/Fcont(j))
            end if
c            print*, 'index, Fline, Fcont = ',aind(j,i),
c     &                      Fline(j),Fcont(j)
         end do
C     end of loop over ingredient populations
      end do


C     prepare to sum by clearing arrays
      do j=1,1221
         hlam(j) = 0.0
      end do
      bolo = 0.0
      mass = 0.0
      do j=1,11
         mags(j) = 0
         Lsq(j) = 0.0
         Lavg(j) = 0.0
      end do
      do j=1,nind
         Fcont(j) = 0.0
         Fline(j) = 0.0
      end do

C     loop over histogram bins and sum
      do iz = 1,npop
C        if there are few stars in a bin, skip it.
         if (shist(iz).lt.0.0001) goto 55
         do j=1,1221
            hlam(j) = hlam(j) + shist(iz)*ahlam(j,iz)
         end do
C        note that we go to non-log units
c         print*, iz,bolo, shist(iz),abolo(iz)
         bolo = bolo + shist(iz)*10**(abolo(iz))
         mass = mass + shist(iz)*amass(iz)
         do j=1,11
            mags(j) = mags(j) +shist(iz)*10**(-0.4*amags(j,iz))
            Lsq(j) = Lsq(j) + shist(iz)*aLsq(j,iz)
            Lavg(j) = Lavg(j) + shist(iz)*aLavg(j,iz)
         end do
         do j=1,nind
            Fcont(j) = Fcont(j) + shist(iz)*aFcont(j,iz)
            Fline(j) = Fline(j) + shist(iz)*aFline(j,iz)
         end do
 55      continue
      end do

C     form Spinrad's UV break amplitudes.
C     B(2640) = ave Fnu flux (2645 to 2675) / (2630 to 2600)
C     B(2900) = ave Fnu flux (2915 to 2945) / (2855 to 2885)  (angstroms)
C     30-Angstrom bandpasses. Do 30 points of interpolation.
      do i=1,4
c         print*, 'bandpass number ',i
         delta = (sp(2,i) - sp(1,i))/30.0
         spavg(i) = 0.0
         do j=1,30
            xwav = 0.10*(sp(1,i) + (j-1)*delta + 0.5*delta)
            call locate(wave,1221,xwav,jw)
            frac = (xwav - wave(jw))/(wave(jw+1) - wave(jw))
            xx = (1.0-frac)*hlam(jw)*xwav**2 + frac*hlam(jw+1)*xwav**2
            spavg(i) = spavg(i) + xx
c            print*, j, xwav, xx/xwav**2, xx, spavg(i)
         end do
      end do
      b2640 = spavg(2)/spavg(1)
      b2900 = spavg(4)/spavg(3)
      write(*,'(a,f7.3)') ' Break B(2640) = ',b2640
      write(*,'(a,f7.3)') ' Break B(2900) = ',b2900


C     open some output files
      open(unit=45,file='ask_most.out',status='unknown')
      open(unit=48,file='ask_sed.out',status='unknown')
      t = ctime(time())
      write(45,'(a,a,a24)') 'Worthey machinery',
     &         ' operating on Padova isochrones ',t
      write(45,'(x)')
      write(45,'(a)') 'Composite population comprised of: '
      write(45,'(a)') ' Age (Gyr)    [Fe/H]    weight'
      write(48,'(a,a,a24)') 'Worthey machinery',
     &         ' operating on Padova isochrones ',t
      write(48,'(x)')
      write(48,'(a)') 'Composite population comprised of: '
      write(48,'(a)') ' Age (Gyr)    [Fe/H]    weight'
      do i=1,npop
         write(45,'(2f10.3,f10.5)') ageref(i), metalref(i), shist(i)
         write(48,'(2f10.3,f10.5)') ageref(i), metalref(i), shist(i)
      end do
      write(45,'(x)')
      write(45,'(x)')
      write(48,'(x)')
      write(48,'(x)')

C     see if it is legal to output spectral indices
      ibad = 0
      do i=1,npop
         if (ageref(i).lt.0.4) ibad = -1
      end do
      if (ibad.eq.-1) then
         write(45,'(a,a)') 'No spectral indices output: ',
     &    ' indices good only for ages older than 0.4 Gyr '
      else
C        form indices and output them to a file
         do k=1,nind
            if (indexunits(k).eq. 1) then
               index(k) = widths(k)*(1.0-Fline(k)/Fcont(k))
            else
               index(k) = -2.5*alog10(Fline(k)/Fcont(k))
            end if
         end do
         write(45,'(a)') 'Spectral Indices'
         write(45,'(7(a10,1x))') (labels(k),k=1,7)
         write(45,'(7(2x,f8.4,1x))') (index(k),k=1,7)
         write(45,'(x)')
         write(45,'(7(a10,1x))') (labels(k),k=8,14)
         write(45,'(7(2x,f8.4,1x))') (index(k),k=8,14)
         write(45,'(x)')
         write(45,'(7(a10,1x))') (labels(k),k=15,21)
         write(45,'(7(2x,f8.4,1x))') (index(k),k=15,21)
         write(45,'(x)')
         write(45,'(7(a10,1x))') (labels(k),k=22,25)
         write(45,'(7(2x,f8.4,1x))') (index(k),k=22,25)
      end if
      write(45,'(x)')
      write(45,'(x)')
C     form Lbar and output

      do k=1,11
         Lbar(k) = -2.5 * alog10( Lsq(k)/Lavg(k) )
      end do
      write(45,'(a)') ' SBF magnitudes '
      write(45,'(a,a,a)') 
     &      '   Ubar   Bbar   Vbar   Rbar   Ibar',
     &      '   Jbar   Hbar   Kbar   Lbar   Lpbar',
     &       '  Mbar'
      write(45,'(11f7.3)')  
     z         (Lbar(k),k=1,11)
      write(45,'(x)')
      write(45,'(x)')

C     output magnitudes and colors
      do k=1,11
         mags(k) = -2.5*alog10(mags(k))
      end do
      write(45,'(a)') ' Ordinary magnitudes '
      write(45,'(a,a,a)') 
     &      '       U       B       V       Rc      I',
     &      'c      J       H       K'
      write(45,'(11f8.3)') (mags(k),k=1,8)
      write(45,'(x)')
      write(45,'(a,a)') '       L       Lp','      M'
      write(45,'(11f8.3)') (mags(k),k=9,11)
      write(45,'(x)')
      write(45,'(x)')
      write(45,'(a)') 'Some colors'
      write(45,'(a,a)') 
     &      '     U-V     B-V     V-Rc    V-Ic    V-K',
     &      '     J-K     H-K'
      write(45,'(7f8.3)') mags(1)-mags(3),mags(2)-mags(3),
     &  mags(3)-mags(4),mags(3)-mags(5),mags(3)-mags(8),
     &  mags(6)-mags(8),mags(7)-mags(8)
      if (mags(2)-mags(3).gt.0.6) then
         write(45,'(x)')
         write(45,'(a,f6.3,a)') 'A B-V of ',mags(2)-mags(3)-0.05,
     &      ' is probably more accurate (see Worthey 94).'
         write(45,'(x)')
      end if
      
      write(45,'(x)')
      write(45,'(x)')
      write(45,'(a,f9.0)') 'Total population mass: ',mass
      write(45,'(a,f7.3)') 'log Lbol/Lbol_sun    : ',alog10(bolo)
      write(45,'(a,f7.3)') 'IMF power-law slope  : ',imfpar
      write(45,'(a,f7.3)') 'IMF lower mass cut   : ',imflow
      write(45,'(a,f7.3)') 'IMF upper mass cut   : ',imfupp
      write(45,'(x)')


C     output SED
      write(48,'(a)')'Spectral energy distribution in wavelength units.'
      write(48,'(x)')
      write(48,'(a)')'Normalized so that the integral over wavelength'
      write(48,'(a)')'in nm is equal to the bolometric luminosity.'
      write(48,'(a)')'Interpolation error can cause up to 1% mismatch.'
      write(48,'(x)')
      write(48,'(a)')'Column 1: wavelength in Angstroms '
      write(48,'(a)')'Column 2: log wavelength (in Angstroms) '
      write(48,'(a)')'Column 3: Flux in wavelength units'
      write(48,'(a)')'Column 4: log of column 3'
      write(48,'(x)')
      write(48,'(x)')
      do j=1,1221
         if ( hlam(j).gt.0.0 ) then
            x1 = alog10(hlam(j))
         else
            x1 = 0.0
         end if
         write(48,'(4e14.5)') 10.0*wave(j),alog10(wave(j))+1.0,
     &        hlam(j),x1
      end do

C     let us do some checks, shall we? Does integral of flux = Lbol?
c      xsum = 0.0
c      do j=2,1221
c         xsum = xsum + 0.5*(hlam(j)+hlam(j-1))*(wave(j)-wave(j-1))
c      end do
c      print*, ' integral = ',xsum,'   bolo =',bolo


      print*, 'Indices, colors, SBF mags are in "ask_most.out"'
      print*, 'SED is in "ask_sed.out"'
      print*, ' Bye!'

      close(45)
      close(48)

      return
      end
