program pickround C written by Phil Massey, Lowell Observatory, June2003 C finds the median value in "tempround" and writes it into "temproundans" implicit double precision(a-h,o-z) dimension data(200000) open(unit=1,file='tempround',status='old') open(unit=2,file='temproundans',status='new') ix=0 1 read(1,*,end=99,err=1) xval ix=ix+1 data(ix)=xval go to 1 99 continue call median(data,ix,amed) write(6,101) amed write(2,102) amed 101 format('Median roundness is:',f6.2) 102 format(f6.2) stop end subroutine median(data,n,amed) implicit double precision(a-h,o-z) dimension data(*),indx(200000) 71 format('Calling sort',i5,f8.3) call sort(n,data,indx) if( ((n/2)*2).eq.n) then amed=(data(indx(n/2))+data(indx(n/2+1)))/2. else amed=data(indx(n/2+1)) end if return end subroutine sort(n,arrin,indx) implicit double precision(a-h,o-z) dimension arrin(*),indx(*) do 11 j=1,n indx(j)=j 11 continue if(n.eq.1) return l=n/2+1 ir=n 10 continue if(l.gt.1) then l=l-1 indxt=indx(l) q=arrin(indxt) ELSE indxt=indx(ir) q=arrin(indxt) indx(ir)=indx(1) ir=ir-1 if(ir.eq.1) then indx(1)=indxt return endif endif i=l j=l+l 20 if(j.le.ir) then if(j.lt.ir) then if(arrin(indx(j)).lt.arrin(indx(j+1)))j=j+1 endif if(q.lt.arrin(indx(j)))then indx(i)=indx(j) i=j j=j+j else j=ir+1 endif go to 20 endif indx(i)=indxt go to 10 end