The IDL procedure calcdata

pro calcdata, file_in, ndata, a, maxv, $
	      z_low, z_high, a_low, a_high, d_low, d_high, new_n

;   Bo Milvang-Jensen, May 21, 1994.

    if n_params() eq (0) then begin
      print
      print, 'Procedure calcdata'
      print
      print, 'Usage:'
      print, 'calcdata, file_in, ndata, a, maxv, $'
      print, '          z_low, z_high, a_low, a_high, d_low, d_high, new_n'
      print
      print, 'Description:'
      print, 'This procedure converts an ascii data file, with alpha, delta, and'
      print, 'redshift, to cartesian coordinates, stored in the array a'
      print
      print, 'file_in : name of data ascii input file with alpha, delta, redshift'
      print, 'file_in: name of data ascii file with alpha, delta, redshift'
      print, 'ndata  : number of datapoint in file_in'
      print, 'a      : array which returns catesian coordinates'
      print, 'mav    : array which return maximum absolute values in x, y and z'
      print, 'z_low  : minimum redshift value'
      print, 'z_high : maximum redshift value'
      print, 'a_low  : minimum alpha    value'
      print, 'a_high : maximum alpha    value'
      print, 'd_low  : mimimum delta    value'
      print, 'd_high : maximum delta    value'
      print, 'new_n  : number of data points left'
      print
      return
    end

    diam = 1.     ; ball diameter 
    
    a=fltarr(3,ndata)

    openr,    lun, file_in, /GET_LUN
    readf,    lun, a
    close,    lun
    free_lun, lun
    
    a=shift(a,1,0) ; before: alpha,    delta, redshift
		   ; after : redshift, alpha, delta  

    b=reform(a,ndata,3) ; transpose 
    b(*,0) = a(0,*)
    b(*,1) = a(1,*)
    b(*,2) = a(2,*)
    a=b
    
    c = where(a(*,0) gt z_low, ndata2)
    if ndata2 ne 0 then begin
      d = a(c,*)
      a = d
      print, 'z     low  cut: ','before = ', ndata, '  after = ', ndata2
    endif
    ndata = ndata2

    c = where(a(*,0) lt z_high, ndata2)
    if ndata2 ne 0 then begin
      d = a(c,*)
      a = d
      print, 'z     high cut: ','before = ', ndata, '  after = ', ndata2
    endif
    ndata = ndata2

    c = where(a(*,1) gt a_low, ndata2)
    if ndata2 ne 0 then begin
      d = a(c,*)
      a = d
      print, 'alpha low  cut: ','before = ', ndata, '  after = ', ndata2
    endif
    ndata = ndata2

    c = where(a(*,1) lt a_high, ndata2)
    if ndata2 ne 0 then begin
      d = a(c,*)
      a = d
      print, 'alpha high cut: ','before = ', ndata, '  after = ', ndata2
    endif
    ndata = ndata2

    c = where(a(*,2) gt d_low, ndata2)
    if ndata2 ne 0 then begin
      d = a(c,*)
      a = d
      print, 'delta low  cut: ','before = ', ndata, '  after = ', ndata2
    endif
    ndata = ndata2

    c = where(a(*,2) lt d_high, ndata2)
    if ndata2 ne 0 then begin
      d = a(c,*)
      a = d
      print, 'delta high cut: ','before = ', ndata, '  after = ', ndata2
    endif
    ndata = ndata2

			      ; NO scaling of redshift!!!
    a(*,1) = (a(*,1)/180)*!pi ; conversion to radians
    a(*,2) = (a(*,2)/180)*!pi ; conversion to radians
     
    a=new_pol_to_cart(a) ; conversion to cartesian coordinates
    
    maxv=fltarr(3)
    maxv(0)=max(abs(a(*,0)),min=minx)
    maxv(1)=max(abs(a(*,1)),min=miny)
    maxv(2)=max(abs(a(*,2)),min=minz)
    print
    print, 'Maximum absolute values (maxv) are: ' ,maxv

    new_n = ndata2
    print, 'Number of data points left (new_n) = ', new_n

end