Reputation: 17
I have questions : Is it possible to define final-value of do loop(ngal) as a free parameter with out exact value or not when I don't know the exact size of my catalog either I don't know the size of nm and nz? ( the whole code as you wished )
integer i,iz,im,ngal,n_thresh,nmax,nz,nm,ok,total,ii,j
integer, dimension(:,:),allocatable :: ngalb
double precision ramin,ramax,decmin,decmax,zmin,zmax,mag_min,mag_max,step,bin,S,dix,fi,Vmax,zbin,magbin
double precision, dimension (:),allocatable :: ra,dec,mag_g,mag_r,mag_i,redshift,zup
double precision, dimension (:,:,:),allocatable :: magb,zb
character plotname*100,zbin_str*4
open(unit=1,file='Lagos12.i.asc')
open(2,file='bin.asc')
open(3,file='mag_vol.asc')
open(4,file='luminosity_func.asc')
ramin=0.0
ramax=2.0
decmin=-2.0
decmax=2.0
zmin=0.0
zmax=3.2
mag_min=13
mag_max=26
step=0.1
bin=0.5
n_thresh=20
allocate (ra(ngal),stat=ok)
allocate (dec(ngal),stat=ok)
allocate (mag_g(ngal),stat=ok)
allocate (mag_r(ngal),stat=ok)
allocate (mag_i(ngal),stat=ok)
allocate (redshift(ngal),stat=ok)
allocate(ngalb(nz,nm),stat=ok)
ngalb=0
do i=1,ngal
read(1,*) ra(i),dec(i),mag_g(i),mag_r(i),mag_i(i),redshift(i)
iz=int((redshift(i)-zmin)/step)+1
im=int((mag_i(i)-mag_min)/bin)+1
ngalb(iz,im)=ngalb(iz,im)+1
enddo
nmax=0
do iz=1,nz
do im=1,nm
nmax=max(nmax,ngalb(iz,im))
enddo
enddo
allocate(magb(nz,nm,nmax),stat=ok)
allocate(zb(nz,nm,nmax),stat=ok)
allocate(zup(nz),stat=ok)
ngalb=0
do i=1,ngal
read(1,*)ra(i),dec(i),mag_g(i),mag_r(i),mag_i(i),redshift(i)
iz=int((redshift(i)-zmin)/step)+1
im=int((mag_i(i)-mag_min)/bin)+1
ngalb(iz,im)=ngalb(iz,im)+1
j=ngalb(iz,im)
magb(iz,im,j)=mag_i(i)
zb(iz,im,j)=redshift(i)
write(2,'(2x,3i5,2x,i5,2x,f10.7,2x,f10.5)')iz,im,j,ngalb(im,iz),magb(im,iz,j),zb(im,iz,j)
enddo
deallocate(magb,zb)
deallocate(ngalb)
deallocate (ra,dec,mag_g,mag_r,mag_i,redshift,zup)
Upvotes: 0
Views: 233
Reputation: 18098
To sum up the comments by High Performance Mark, steabert, Vladimir F, and myself:
! First pass: Determine sizes
nm = 0 ; nz = 0 ; ngal = 0
do
read(1,*,iostat=ierr) dummy_ra, dummy_dec, dummy_mag_g, dummy_mag_r, &
dummy_mag_i, dummy_redshift
if ( ierr < 0 ) exit ! End of file
if ( ierr > 0 ) stop 'An error occured reading from unit 1'
nz = max(nz,int((dummy_redshift-zmin)/step)+1)
nm = max(nm,int((dummy_mag_i-mag_min)/bin)+1)
ngal = ngal + 1
enddo
! Allocate arrays
allocate (ra(ngal), dec(ngal), mag_g(ngal), mag_r(ngal), mag_i(ngal), &
redshift(ngal), ngalb(nm,nz), stat=ierr)
if ( ierr /= 0 ) stop 'Cannot allocate memory'
ngalb=0
! Second pass: Read in
rewind(1)
do i=1,ngal
read(1,*) ra(i), dec(i), mag_g(i), mag_r(i), mag_i(i), redshift(i)
iz = int((redshift(i)-zmin)/step)+1
im = int((mag_i(i)-mag_min)/bin)+1
ngalb(im,iz) = ngalb(im,iz)+1
enddo
Upvotes: 4