C order_bdm.f May 2004 S.G. C C This program orders bdm catalogs by mass CHARACTER*120 name1,name2,f_tmp CHARACTER*120 line CHARACTER*90 trash PARAMETER (nmax = 120000) DIMENSION r_bdm(3,nmax),r_prim(3,3), + v_fof(nmax), + v_bdm3(3,nmax),v_bdm(nmax), + total_mass_bdm(nmax),rh_bdm(nmax), Vrms_bdm(nmax), + vcirc_bdm(nmax),npar_bdm(nmax), + npar_bdm1(nmax),npar_bdm2(nmax),npar_bdm3(nmax), + indx(nmax), rmax(nmax),conc(nmax) c box = 50. write(*,*) 'filename of the catshort ' READ (*,'(A)') f_tmp c f_tmp = 'Catshort.DAT' CALL get_name(f_tmp,i,j) name2 = f_tmp(i:j) name1 = f_tmp(i:j)//'_ord' CALL get_name(name1,in,jn) OPEN(11,FILE= name1(in:jn),FORM='FORMATTED', + STATUS='UNKNOWN') CALL get_name(name2,in,jn) OPEN(10,FILE= name2(in:jn),FORM='FORMATTED', + STATUS='OLD') i = 1 Do k = 1,21 READ(10,'(A)') trash write(11,'(A)') trash ENDDO READ(10,'(A)') line WRITE(11,123) line write(*,123) line 123 Format(' # ',A120) DO WHILE (.TRUE.) READ(10,*,END=260) + r_bdm(1,i),r_bdm(2,i),r_bdm(3,i), + v_bdm3(1,i),v_bdm3(2,i),v_bdm3(3,i), + total_mass_bdm(i),rh_bdm(i), Vrms_bdm(i), + vcirc_bdm(i),npar_bdm(i), + rmax(i),conc(i) i = i + 1 ENDDO 260 n_bdm = i - 1 write(*,*) n_bdm,' objects read from the bdm_catalog' CALL indexx(n_bdm,vcirc_bdm,indx) n_40 = 0 n_50 = 0 n_80 = 0 n_100 = 0 kk = 0 DO k = n_bdm,1,-1 kk = kk + 1 i = indx(k) write(11,200) kk, + r_bdm(1,i),r_bdm(2,i),r_bdm(3,i), + v_bdm3(1,i),v_bdm3(2,i),v_bdm3(3,i), + total_mass_bdm(i),rh_bdm(i), Vrms_bdm(i), + vcirc_bdm(i),npar_bdm(i), + rmax(i),conc(i) IF(vcirc_bdm(i).GT.100.) n_100 = n_100 + npar_bdm(i) IF(vcirc_bdm(i).GT.80.) n_80 = n_80 + npar_bdm(i) IF(vcirc_bdm(i).GT.50.) n_50 = n_50 + npar_bdm(i) IF(vcirc_bdm(i).GT.40.) n_40 = n_40 + npar_bdm(i) ENDDO write(*,*) n_100,' particles in halos with v_circ > 100 km/s' write(*,*) n_80,' particles in halos with v_circ > 80 km/s' write(*,*) n_50,' particles in halos with v_circ > 50 km/s' write(*,*) n_40,' particles in halos with v_circ > 40 km/s' 200 FORMAT(I8,3F8.3,3X,3F8.1,G12.3,3F8.2,I8,2F8.1) STOP END ******************************************************************* c distance with periodic boundary conditions FUNCTION DIST(boxsize,XI,XJ,YI,YJ,ZI,ZJ) DX = ABS(XI - XJ) DDX =ABS(DX - boxsize) DX = MIN(DX,DDX) DY = ABS(YI - YJ) DDY =ABS(DY - boxsize) DY = MIN(DY,DDY) DZ = ABS(ZI - ZJ) DDZ =ABS(DZ - boxsize) DZ = MIN(DZ,DDZ) DIST = SQRT(DX**2 + DY**2 + DZ**2) RETURN END ******************************************************************* C Ordering of a real array C C Indexing Sect. 8.4. of Numerical Recipis C Changed to integer field SUBROUTINE indexx(n,arr,indx) INTEGER n,indx(n),M,NSTACK REAL arr(n) PARAMETER (M=7,NSTACK=50) INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK) REAL a do 11 j=1,n indx(j)=j 11 continue jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 13 j=l+1,ir indxt=indx(j) a=arr(indxt) do 12 i=j-1,1,-1 if(arr(indx(i)).le.a)goto 2 indx(i+1)=indx(i) 12 continue i=0 2 indx(i+1)=indxt 13 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 itemp=indx(k) indx(k)=indx(l+1) indx(l+1)=itemp if(arr(indx(l+1)).gt.arr(indx(ir)))then itemp=indx(l+1) indx(l+1)=indx(ir) indx(ir)=itemp endif if(arr(indx(l)).gt.arr(indx(ir)))then itemp=indx(l) indx(l)=indx(ir) indx(ir)=itemp endif if(arr(indx(l+1)).gt.arr(indx(l)))then itemp=indx(l+1) indx(l+1)=indx(l) indx(l)=itemp endif i=l+1 j=ir indxt=indx(l) a=arr(indxt) 3 continue i=i+1 if(arr(indx(i)).lt.a)goto 3 4 continue j=j-1 if(arr(indx(j)).gt.a)goto 4 if(j.lt.i)goto 5 itemp=indx(i) indx(i)=indx(j) indx(j)=itemp goto 3 5 indx(l)=indx(j) indx(j)=indxt jstack=jstack+2 if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx' if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 END C (C) Copr. 1986-92 Numerical Recipes Software .-35?421.1-9. *********************************************************** SUBROUTINE get_name(name,i,j) CHARACTER*120 name DO i=1,120 IF(name(i:i).NE.' ') GOTO 1 ENDDO i=1 j=1 RETURN 1 CONTINUE DO j=i,120 IF(name(j:j).EQ.' ') GOTO 2 ENDDO j=121 2 j=j-1 IF(name(i:i).EQ.'''') i=i+1 IF(name(j:j).EQ.'''') j=j-1 END