      program mean_med ! Bernd Berg Sep 08 2008.
c Mean and median for exp[xr*ln(10000)], x uniform in [0,1).
      include '../../ForLib/implicit.sta'
      include '../../ForLib/constants.par'
      parameter(iuo=6,iud=10,miuo=-iuo,iseed1=1,iseed2=0,ndat=10**5+1)
      dimension a(ndat)
      if(2*(ndat/2).eq.ndat) stop "mean_med: chose ndat odd."
c
      fact1=10000*one
      fact2=log(fact1)
      write(iuo,'(/," ndat =",1I9,/)') ndat
      call rmaset(miuo,iud,iseed1,iseed2,'ranmar.d')
      do idat=1,ndat
        call ranmar(xr)
        a(idat)=fact1*exp(xr*fact2)
      enddo
      call heapsort(ndat,a)
      i025=nint(0.025d00*ndat)
      i15=nint(0.15d00*ndat)
      i50=1+ndat/2
      i85=nint(0.85d00*ndat)
      i975=nint(0.975d00*ndat)
      write(iuo,'(" 95% and 70% confidence ranges:",/)') 
      write(iuo,'(6X,"0.025%",7X,"15%",7X,"Median",8X,"85%",8X,"97.5%",
     &  /,/,5G14.6)') a(i025),a(i15),a(i50),a(i85),a(i975)
      write(iuo,'(" at array locations:",/,5I12)') 
     &                         i025,i15,i50,i85,i975
      amean=stmean(ndat,a)
      write(iuo,'(/," Mean value =",G14.6,/)') amean
      call qdf_gnu(iud,ndat,a)
c
      stop
      end

      include '../../ForLib/heapsort.f'
      include '../../ForLib/qdf_gnu.f'
      include '../../ForLib/rmaset.f'
      include '../../ForLib/ranmar.f'
      include '../../ForLib/rmasave.f'
      include '../../ForLib/stmean.f'
      include '../../ForLib/write_progress.f'
