Program ELFread
c Example program to read ELF data files
c This version setup to read 2-channel data only

c local variables
      integer*4 numrecs, iw
	
c disk image name 
      character*80 imagefile
c header common
      integer*4 BlockNum(4096), NumChan(4096), FreqNum(4096)
      integer*4 NBlocks(4096), tapeblock
      integer*4 address(4096),time1(4096),time2(4096),ymdhms(6,4096)
      common /TheHeader/ BlockNum, NumChan, FreqNum, NBlocks, address,
     *   time1, time2, ymdhms, tapeblock

      real gains(6), tgains(6)
      integer*2 idata(256000)  
	dimension fdata(6,256000)
      real center(14), vcogain(14), freqs(14)
     
c  center frequencies of the VCOs at various sampling frequencies
      data center /3906.25, 7812.5, 15625., 31250.,7812.5, 15625., 
     * 31250.,7812.5, 15625., 31250.,1953.125, 3906.25, 7812.5,15625./
        
c get the file name
      write(*,*) ' enter file name'
      read (*,'(a)') imagefile
      open (12, file=imagefile, form='unformatted', 
     *              access='direct', recl=512, readonly)

c read the header
      iw = 1
      call GetHeader(iw, numrecs)
      
      write(*,*) ' Enter records to read'
      read (*,*) nnbb1, nnbb2
c      write(*,*) ' Enter gains for ',NumChan(1),' channels (Volts/physical unit)'
c      read (*,*) (gains(i), i=1,NumChan(1))
c hardwire the gains for 2 channel AC electric fields with a 10m dipole
      gains(1) = 9000000.
      gains(2) = 9000000.
      	
c set up a frequency lookup table for the ELFs
      do i = 1,14
	  freqs(i) = 512./float(2**i)
	end do
c find the frequency for these data
	freq = freqs(FreqNum(nnbb1))

c set up the VCO gains for a +- 5V swing
      do i = 1,14
        vcogain(i) = center(i)/5.
	end do
c pick off the center frequency for these data
      cent = center(FreqNum(nnbb1)) - 32768.     
c set up the gains
      nc = NumChan(nnbb1)
      do i = 1, nc
        tgains(i) = 1./gains(i)/vcogain(FreqNum(nnbb1))
	end do

c initialize all the counters
      nall = 245760
      ndata = nall/nc 
	nrec = ndata/512
		
	open(13,file='elfdata.ch1',form='unformatted',
     *     access='direct',recl=512*4)
      open(14,file='elfdata.ch2',form='unformatted',
     *     access='direct',recl=512*4)

      irec = 0
      do k = nnbb1,nnbb2
c get the data
        write(*,*) ' getting record ', k
        call GetDatai(k, idata, ndat)
c  demultiplex the data 
        do ic = 1, nc
	    j = 1
          do i = ic, nall, nc
	      fdata(ic,j) = idata(i)
	      j = j+1
	    end do
	  end do
	  do j = 0, nrec-1
	    irec = irec + 1
          write(13,rec=irec) (fdata(1,i),i=512*j+1,512*j+512)
          write(14,rec=irec) (fdata(2,i),i=512*j+1,512*j+512)
	  end do
      end do
	
      close(12)
      close(13)
      stop
      end
      
c----------------------------------------------------------------------      
      subroutine GetDatai(nrec, idata, ndat)
c----------------------------------------------------------------------      

c header common
      integer*4 BlockNum(4096), NumChan(4096), FreqNum(4096)
      integer*4 NBlocks(4096), tapeblock
      integer*4 address(4096), time1(4096), time2(4096), ymdhms(6,4096)
      common /TheHeader/ BlockNum, NumChan, FreqNum, NBlocks, address,
     *   time1, time2, ymdhms, tapeblock
c arguments
      integer*4 nrec
      integer*2 idata(256000)
c local variables
      integer*4  NumBlocks, StartBlock
c binary data buffer
      character*1 buffer(512000)

c reading from disk image file

        StartBlock = address(nrec)
        NumBlocks = NBlocks(nrec)
        do 20 l = 1,NumBlocks
          iblock = StartBlock+l
          read(12,rec=iblock+2,err=19)(buffer(k),k=512*(l-1)+1,512*l)
20      continue
              
c byteswap data and offset by 256*127 to turn unsigned into signed integer*2
        ndata = 491520
        i = 0
        do 800 index = 0, ndata-2, 2
          i = i+1
800       idata(i) = ichar(buffer(index+1)) +
     *           256*(ichar(buffer(index+2))-127)

        ndat = i
        return
19      continue
        write(*,*) ' read error from image file'
        stop
        end
c----------------------------------------------------------------------      
      subroutine GetHeader (iw, numrecs)
c----------------------------------------------------------------------      
c reads an ELF header from a disk image file opened on unit 12 
c if iw .eq. 0 then the header info is returned in the common block 
c without writing any diagnostics
c if iw .eq. 1 then the entire header is written to screen.
c if iw .eq. 2 then the first and last entries are written to screen.
c
c on output, numrecs is the number of disk writes (records) for this
c deployment
c
c as of Dec 2000, the image files are offset 2 blocks from the raw
c disk data to accomodate an additional header structure in the 
c second block containing information about deployment location, etc. 

c header common
      integer*4 BlockNum(4096), NumChan(4096), FreqNum(4096)
      integer*4 NBlocks(4096), tapeblock
      integer*4 address(4096), time1(4096), time2(4096), ymdhms(6,4096)
      common /TheHeader/ BlockNum, NumChan, FreqNum, NBlocks, address,
     *   time1, time2, ymdhms, tapeblock
c argument list
      integer*4 numrecs, iw
c read buffer
      character*1 buffer(131072)
c local variables
      integer*4 ibuff(512)

      if (iw. ne. 0) write(*,*) ' Reading header'

      do l = 1,256
         read (12,rec=l+2,err=21)(buffer(k),k=512*(l-1)+1,512*l)
      end do
c
c decompose the header
c
      lastblock = 0      
      do 500 k = 1, 4096
c
c pull of numerical data from the buffer
c
        index = (k-1)*32
        do 90 i = 1,32
  90      ibuff(i) = ichar(buffer(index+i))
c
c get the record parameters into the right arrays
c
        i0 = 1
        BlockNum(k) = ibuff(i0+0) + 256*ibuff(i0+1)
        do 111 i=1,6        
	    ihigh = iand(int2(jishftc(ibuff(i0+1+i), -4, 32)),$000F)
	    ilow =  iand(int2(ibuff(i0+1+i)),$000F)
111       ymdhms(i,k) = 10*ihigh + ilow
        NumChan(k) = ibuff(i0+8)
        FreqNum(k) = ibuff(i0+9)
        NBlocks(k) = ibuff(i0+10) + 256*ibuff(i0+11)
        address(k) = ibuff(i0+14) + 256*(ibuff(i0+15) + 
     *       256*(ibuff(i0+12) + 256*ibuff(i0+13)))
        time1(k) = ibuff(i0+18) + 256*(ibuff(i0+19) +
     *       256*(ibuff(i0+16) + 256*ibuff(i0+17)))
        time2(k) = ibuff(i0+22) + 256*(ibuff(i0+23) + 
     *       256*(ibuff(i0+20) + 256*ibuff(i0+21)))
c
c test to see if the data still make sense
c
        if (BlockNum(k) .ne. lastblock+1) then
          numrecs = k-1
          goto 600
        end if
        lastblock = BlockNum(k)
500   continue
600   continue      

c write it out
      if (iw .ge. 1) then
        inc = 1
        if (iw .eq. 2) then
          inc = numrecs-1
          write(*,*) 'first and last records:'
        end if
        do 700 k = 1, numrecs, inc
700          write(*,'(i5,1x,5(i2,'':''),i2,2i3,i4,i8,2i12)') 
     *      BlockNum(k), (ymdhms(j,k),j=1,6),
     *      NumChan(k), FreqNum(k), NBlocks(k), 
     *      address(k), time1(k), time2(k)
        write(*,*) 'Number of good records =',numrecs 
      end if
            
      return
21    continue
      write(*,*) ' read error from image file'
      stop
      end