*==============================================================================
* Program to find the completeness at a given z and b:
*
* compile: f77 -o complete_interp complete_interp.f
*
* written by S.M. Croom 6/7/99
*
* version 1.1 (06/06/03) updated for larger mag range and different
*                        format of input file.
*
*###############################################################################
*                                                                              #
*          LICENCE                                                             #
*                                                                              #
*          This program is free software; you can redistribute it and/or       #
*          modify it under the terms of the GNU General Public License         #
*          as published by the Free Software Foundation; either version 2      #
*          of the License, or (at your option) any later version.              #
*                                                                              #
*          This program is distributed in the hope that it will be useful,     #
*          but WITHOUT ANY WARRANTY; without even the implied warranty of      #
*          MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       #
*          GNU General Public License for more details.                        #
*                                                                              #
*          You should have received a copy of the GNU General Public License   #
*          along with this program; if not, write to the Free Software         #
*          Foundation, Inc., 59 Temple Place - Suite 330,                      #
*          Boston, MA  02111-1307, USA.                                        #
*                                                                              #
*###############################################################################
*
      program comp_intr
*
* Definition of variables:     
      implicit none
      integer i,j,mbin,zbin,mb,zb
      parameter (mbin=50,zbin=30)
      real m,z,c,compar(mbin,zbin)
      real x1(mbin),x2(zbin),comp,compval

* common and data blocks:
      common /complete/ compar,x1,x2
*------------------------------------------------------------------------------
*
* read in completeness array:
      open(20,file='col_comp_data.txt',status='old')
      do i=1,100000
         read(20,*,end=100,err=999) z,m,c
         if (x2(zb).lt.z) zb=zb+1
         if (x1(mb).gt.m) mb=0
         mb=mb+1
         x1(mb)=m
         x2(zb)=z
         compar(mb,zb)=c*100.0
      end do
 100  close(20)      
*
* write out array to screen:
      write(*,*) 'completeness array:'
      write(*,10) (x2(j),j=1,zbin)
 10   format(' b\\z',30f5.1)
      do i=1,mbin
         m=x1(i)
         write(*,'(f4.1,30f5.1)') m,(compar(i,j),j=1,zbin)
      end do
*
* read in m,z values from user and give completeness value at that point:
      write(*,*) 'input b_J and z:'
      read(*,*) m,z
      comp=compval(m,z)
      write(*,*) 'completeness=',comp,' percent'

      stop
 999  stop 'ERROR: cant read file'
      end
*==============================================================================
* function to do simple interpolation of completeness grid
* will work over the range Bj=16.0-20.9 and z=0.1-3.0
*
      function compval(b,z)
      implicit none
      integer mbin,zbin,mb,zb
      real b,z,boff,zoff,bstep,zstep
      parameter (mbin=50,zbin=30)
      real c(mbin,zbin),compval,t,u
      real x1(mbin),x2(zbin)

* common and data blocks:
      common /complete/ c,x1,x2
*------------------------------------------------------------------------------
* find which array element to use.  need, mb, mb+1 and zb, zb+1
* assumes b and z increase, and have constant binning:
      boff=x1(1)
      zoff=x2(1)
      bstep=x1(2)-x1(1)
      zstep=x2(2)-x2(1)
      mb=int((b-boff)/bstep)+1
      zb=int((z-zoff)/zstep)+1
* do a simple bilinear interpolation:
      t=(b-x1(mb))/(x1(mb+1)-x1(mb))
      u=(z-x2(zb))/(x2(zb+1)-x2(zb))
      compval=(1.0-t)*(1.0-u)*c(mb,zb)+t*(1.0-u)*c(mb+1,zb)+t*u*c(mb+1
     &     ,zb+1)+(1-t)*u*c(mb,zb+1)
      return
      end
*==============================================================================