program get_EUV
c
c          Read the Hinteregger data file SC21OBSC.EUVS into internal
c       arrays. This program can be made into a subroutine that returns
c       the AE-E EUVS data.
c
      character*128 buffer
      integer*4 YYDDD, UT
      integer*4 lenstr
      logical st_data /.false./
      logical reference /.false./
      real*8 year_j, YYDDD2year, waveg
      parameter ( nptmax=4096, ngroups=15 )
      dimension waveg(ngroups)
      dimension flux(nptmax, ngroups), year(nptmax, ngroups)
      dimension ireg(2), reg(ngroups), ngroup(ngroups), UT(nptmax)
      dimension YYDDD(nptmax), npgrp(ngroups), FJUL76(ngroups)
c
c          WAVEG is the central wavelength of group i.
c
      data waveg/ 179.d0, 198.d0, 230.5d0, 277.5d0, 304.d0, 545.d0,
     $            584.d0, 625.d0, 1026.d0,  335.d0, 284.d0, 202.d0,
     $           180.5d0, 171.d0, 1216.d0/
c
c          Open the combined AE-E EUVS data file.
c
      open (unit=11, file='sc21obsc.euvs', status='old', readonly )
c
c          Each wavelength group has a different number of points so
c       keep track of each group with a separate counter.
c
      npts = 0
      do i=1,ngroups
         npgrp(i) = 0
      enddo
      do while ( npts .le. nptmax )
         read(11,1000,end=900) buffer
         if( .not. reference .and. buffer(1:5) .eq. ' WLG#' ) then
c
c          Get the reference spectrum values for the assorted wavelength
c       groups.
c
            reference = .true.
            read(11,1000) buffer
            read(11,1100) (ngroup(i), FJUL76(i),i=1,ngroups)
            write(6,1100) (ngroup(i), FJUL76(i),i=1,ngroups)
         elseif( reference .and. .not. st_data .and.
     $                           buffer(1:6) .eq. '  DATE' ) then
c
c          Past the reference spectrum and hitting the keyword '  DATE',
c       the rest of the file is ratios of irradiance in the different 
c       groups.
c
            st_data = .true.
         elseif( reference .and. st_data .and.
     $                           lenstr(buffer) .gt. 0 ) then
c
c          The LENSTR guarantees no blank lines in the data listing
c
            read(buffer, 2100) ireg(1), ireg(2), (reg(i),i=1,ngroups)
            npts = npts + 1
            YYDDD(npts) = ireg(1)
            UT(npts) = ireg(2)
c
c          YEAR_J is the (year - 1900)
c
            year_j = YYDDD2year( YYDDD(npts) ) +
     $                               float(UT(npts))/3.15576d7
c
c          Cycle through the groups and keep only the points that have
c       data.
c
            do j=1, ngroups
               if( reg(j) .gt. 0 ) then
                  npgrp(j) = npgrp(j)+1
                  flux(npgrp(j),j) = reg(j)*FJUL76(j)
                  year(npgrp(j),j) = year_j
               endif
            enddo
         endif
      enddo
c
c          End-of-File or maximum number of points stored
c
 900  continue
      write(6,2000) npts
      close(11)
c
c          In this section of the program you can manipulate the data
c
c          As an example of data manipulation, write out a 2-column file
c       with the Lyman-alpha data as a function of time.
c
      open (unit=15, file='ae-e.lya', status='new' )
      write(15,4000) (year(i,15), flux(i,15), i=1,npgrp(15))
      stop
c
 1000 format(a)
 1100 format(5(i5,2x,f7.3))
 2100 format(1x,i5,1x,i5,1x,15(2x,f5.2))
 2000 format(1x,'Read ',i5,' points from datafile')
 4000 format(1x,0pf9.6,1x,1pe12.4)
      end
      real*8 function YYDDD2year( idate )
c
c          Returns the year + fraction for the year and day number
c       passed in YYDDD. The first two digits are the year - 1900 and
c       the last three digits give the day number.
c
      iyear = idate/1000
      iday = idate - 1000*iyear
      if( iyear .lt. 100 ) iyear = iyear + 1900
      if( mod(iyear,100) .eq. 0 ) then
c
c          Century leap year must be divisible by 400
c
         if( mod(iyear,400) .eq. 0 ) then
            ndays = 366
         else
            ndays = 365
         endif
      else
c
c          Normal leap year must be divisible by 4
c
         if( mod(iyear,4) .eq. 0 ) then
            ndays = 366
         else
            ndays = 365
         endif
      endif
      YYDDD2year = float(iyear-1900) + float(iday-1)/float(ndays)
      return
      end
      function lenstr( string )
      character*(*) string
      character*1 blank, null
      parameter ( blank=' ', null=char(0) )
c
c          Returns the position of the right-most non-blank 
c       character in string.
c
      ls = len( string )
      if( ls .eq. 0 ) then
         lreg = ls
      else
         lreg = ls
         do 10 i=ls,1,-1
            if( string(i:i) .ne. blank .and. 
     $          string(i:i) .ne. null ) goto 20
            lreg = i-1
  10     continue
         lreg = 0
      endif
  20  continue
      lenstr = lreg
      return
      end