c  *********************************************************
      subroutine TAPER(fpar,xyz,b)
c
c  Inverse cubic fit for tapering magnetic field
C  From ICOOL subtroutine SOLENOID for fpar = 11
c   http://www.hep.princeton.edu/~mcdonald/mumu/target/taper.pdf
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      dimension zsol(1000),bsol(1000)
      dimension ds0(8),c(8)
      save bsc,bsd,per,fns,cs,sn,iord
      save zsol,bsol,npts
      save sold
!   
!     real*8 cs(0:199),sn(0:199),bsc(0:199),bsd(0:199)
      real*8 xyz(3),b(3)
      real*8 fpar(15)
!      character*10 fname
!      data bnorm/2.e-7/
      data i11/0/    ! to test for first entry to model 11 (KTM)
!
      do i=1,3
      b(i) = 0.
      end do
!
!      model = fpar(1)

      rp = SQRT( xyz(1)**2 + xyz(2)**2 )
!
      if( rp .gt. 0. ) then
        phi = ATAN2( xyz(2),xyz(1) )
        sphi = SIN(phi)
        cphi = COS(phi)
      else
        phi = 0.
        sphi = 0.
        cphi = 1.
      end if
!
!
! ---------------------------------------------------------
!      else if( model .eq. 11) then    ! field from inverse cubic
!
        if (i11 .eq. 0) then      ! initialize
	  i11 = 1
	  pow = fpar(2)    ! pow = 1 is basic inverse cubit fit
	  if (pow .eq. 0.) then
	    write(2,*) 'SOL model 11 power cannot be 0'
	    stop
	  end if
	  iord = fpar(3)   ! iord = 0 = 1st order, = 1 = 3rd order
	  z1 = fpar(4)     !  z of beginning of taper region
	  b1 = fpar(5)     !  b_z(0,z1) = axial field at z1
	  db1 = fpar(6)    ! = d b_z(0,z1) / dz  at z1
	  z2 = fpar(7)     !  z at end of taper
	  b2 = fpar(8)     !  b_z(0,z2) = axial field at z2; d b_z(0,z2) = 0.
	  if (b2 .le. 0.) then
	    write(2,*) 'SOL model 11 field at z2 must be positive'
	    stop
	  end if
	  db2 = fpar(9)
	  a1 = - db1/b1/pow
	  z21 = z2 - z1
	  temp1 = ((b1/b2)**(1./pow) - 1.) / z21**2
	  temp2 = a1 / z21
	  a2 = 3.*temp1 - 2.*temp2
	  a3 = (-2.*temp1 + temp2) / z21
	end if
!
! The inverse cubic fit
!
        dz = xyz(3) - z1
	cubic = 1. + a1*dz + a2*dz**2 + a3*dz**3
	if (cubic .le. 0.) then
	  write(2,*) 'SOL model 11: cubic went negative'
	  stop
	end if
	bz = b1 / cubic**pow
	temp3 = - pow*bz / cubic
	temp4 = a1 + 2.*a2*dz + 3.*a3*dz**2
	dbdz1 = temp3*temp4
        br = -rp / 2. * dbz1       
!
        if( iord .gt. 1 ) then           ! 3rd order
          temp5 = 2.*a2 + 6.*a3*dz
	  temp6 = -(pow + 1.)*temp3 / cubic
	  dbz2 = temp3*temp5 + temp6*temp4**2
          bz = bz - rp**2 / 4. * dbz2
	  temp7 = -(pow + 2.)*temp6 / cubic
	  dbz3 = 6.*a3*temp3 + 2.*temp6*temp4*temp5 + temp7*temp4**3
          br = br + rp**3 / 16. * dbz3
        end if
		
        b(1) = br * cphi
        b(2) = br * sphi
        b(3) = bz	
!
!     end if      ! end of test on models
!  --------------------------------------------------------
900   continue
      end