c   Reads the HEAO-1 A2 output files and creates summed spectra
c   for the specified detector. Creates LFOV, SFOV and LFOV-SFOV
c   spectra.
c   Uses all observations listed in an input file.

	integer*4 MXCHAN
	parameter (MXCHAN=64)

	real*4 input(MXCHAN), errors(MXCHAN)
	real*4 stime, ltime
	real*4 small(MXCHAN-1), small2(MXCHAN-1)
	real*4 large(MXCHAN-1), large2(MXCHAN-1)

	integer*4 ierr, ilun, jlun, lenact, idet
	integer*4 i, j, iobs, iparse, nret, iflag, idelim

	character*72 outfil, lstfil, cfile, infil
	character*3 cdet
	character*512 string

	logical qend, qexist, qfound

	do i = 1, MXCHAN-1
	    small(i) = 0.
	    small2(i) = 0.
	    large(i) = 0.
	    large2(i) = 0.
	end do
	stime = 0.
	ltime = 0.

	call xinird('List file, detector (HD1,HD2,MED,HD3), output file : ',
     &	 string, iparse)
	if (ierr .ne. 0) call exit(1)
	call xgtstr(string, iparse, 1, 'List file', 1, lstfil, nret,
     &	  iflag, idelim)
	if (nret .eq. 0) call exit(1)
	call xgtstr(string, iparse, 1, 'detector (HD1,HD2,MED,HD3)', 1, cdet,
     &	  nret, iflag, idelim)
	if (nret .eq. 0) call exit(1)
	call xgtstr(string, iparse, 1, 'Output file', 1, outfil, nret,
     &	  iflag, idelim)
	if (nret .eq. 0) call exit(1)

	call getlun(jlun)
	call openwr(jlun,lstfil,'old',' ',' ',0,0,ierr)
	if (ierr .ne. 0) then
	    write (*,*) 'IOSTAT= ', ierr
	    write (*,*) lstfil
	    call exit(2)
	end if

c  loop over all observations

	read (jlun,*) iobs
	qend = .false.
	call getlun(ilun)

	do while(.not.qend)

	    infil = char(92)//'/ros2/HEAO1A2/SPECTRA/OBS/obs.'

c  set up filename and check whether the file exists

	    if(iobs.lt.10) then
		infil(32:33)='00'
	        write(infil(34:34),'(i1)') iobs
	    else if(iobs.lt.100) then
		infil(32:32)='0'
	        write(infil(33:34),'(i2)') iobs
	    else
	        write(infil(32:34),'(i3)') iobs
	    end if

	    inquire(file=infil(2:34), exist=qexist)
	    if (qexist) then
		call openwr(ilun,infil,'old',' ',' ',0,0,ierr)
	    else
		ierr = 1
		write (*,*) 'Failed to find '//infil(2:lenact(infil))
	    end if

	    if (ierr .ne. 0) then
		write (*,*) 'Failed to open '//infil(:lenact(infil))
	    end if

	    if (ierr .eq. 0) then

c  if file exists read it and find the data for the selected detector

		read (ilun,'(a)') string
		qfound = .false.

		do while (.not.qfound)

		    if(string(2:4) .eq. cdet) then
			qfound = .true.
			do i = 1, 8
			    read(ilun,*) (input((i-1)*8+j),j=1,8)
			end do
		
c  include in the running sums

cd
	write (*,*) iobs, input(1)
cd

			ltime = ltime + input(1)
			do i = 2, MXCHAN-1
			    large(i) = large(i) + input(i)
			    large2(i) = large2(i) + input(i)*input(i)
			end do

c  jump over background data

			read (ilun,*)
			do i = 1, 8
			    read(ilun,*)
			end do

c  now the SFOV

			read(ilun,*)
			do i = 1, 8
			    read(ilun,*) (input((i-1)*8+j),j=1,8)
			end do

cd
	write (*,*) iobs, input(1)
cd
			stime = stime + input(1)
			do i = 2, MXCHAN-1
			    small(i) = small(i) + input(i)
			    small2(i) = small2(i) + input(i)*input(i)
			end do

			read(ilun,*)		
			do i = 1, 8
			    read(ilun,*)
			end do

c  else if not selected data jump over

		    else
			do i = 1, 17
			    read(ilun,*)
			end do
		    end if

c  read header for next section of the input file

		    read(ilun,'(a)',iostat=ierr) string
		    if (ierr .ne. 0) qfound = .true.

		end do

	    end if

c  read next file to be used

	    read(jlun,*,iostat=ierr) iobs
	    if (ierr.ne.0) qend=.true.

	end do

c now create the data files. Calculate a spectrum from LFOV-SFOV
c simply counting rate errors.

	if (cdet .eq. 'HD1') then
	    idet = 1
	else if (cdet .eq. 'HD2') then
	    idet = 4
	else if (cdet .eq. 'MED') then
	    idet = 7
	else if (cdet .eq. 'HD3') then
	    idet = 10
	end if

c  first the small FOV

	do i = 1, MXCHAN-1
	    input(i) = small(i)
	    errors(i) = sqrt(small(i))
	end do

	do i = 1, MXCHAN-1
	    if (errors(i)**2.ge.input(i)) then
		errors(i) = sqrt(errors(i)**2-input(i))/input(i)
	    else
		errors(i) = 0.
	    end if
	end do

	call locase(cdet)
	cfile=outfil(:lenact(outfil))//cdet//'s'
	call xtend(cfile,'.pha')
	call write_pha (input, errors, stime, cfile, 'none',
     &	  'none', idet, MXCHAN-1)

c  then the large FOV

	do i = 1, MXCHAN-1
	    input(i) = large(i)
	    errors(i) = sqrt(large(i))
	end do

	do i = 1, MXCHAN-1
	    if (errors(i)**2.ge.input(i)) then
		errors(i) = sqrt(errors(i)**2-input(i))/input(i)
	    else
		errors(i) = 0.
	    end if
	end do

	call locase(cdet)
	cfile=outfil(:lenact(outfil))//cdet//'l'
	call xtend(cfile,'.pha')
	call write_pha (input, errors, ltime, cfile, 'none',
     &	  'none', idet, MXCHAN-1)

c and finally the LFOV-SFOV

	do i = 1, MXCHAN-1
	    input(i) = large(i) - (ltime/stime)*small(i)
	    errors(i) = sqrt(large(i) + (ltime/stime)**2*small(i))
	end do

	do i = 1, MXCHAN-1
	    if (errors(i)**2.ge.input(i)) then
		errors(i) = sqrt(errors(i)**2-input(i))/input(i)
	    else
		errors(i) = 0.
	    end if
	end do

	call locase(cdet)
	cfile=outfil(:lenact(outfil))//cdet
	call xtend(cfile,'.pha')
	call write_pha (input, errors, ltime, cfile, 'none',
     &	  'none', idet, MXCHAN-1)

	call exit(0)
	end

