SUBROUTINE reconstruct(params,x,y,z,Kc)

  USE alias
  USE input_exterior
  use input_params

  IMPLICIT NONE

  INTERFACE

     SUBROUTINE harmonics(l,m,z,south,ispole,v)
       INTEGER, INTENT(in) :: l,m
       COMPLEX*16, INTENT(in) ::  z
       COMPLEX*16, DIMENSION(5), INTENT(out) :: v(5)
       LOGICAL, INTENT(in) :: south,ispole
     END SUBROUTINE harmonics

  END INTERFACE

  real*8, dimension(100), intent(in)    :: params
  REAL*8, INTENT(in)                    :: x,y,z
  REAL*8, DIMENSION(3,3),INTENT(out)    :: Kc

  LOGICAL                  :: south, ispole
  REAL*8                   :: r,theta,phi,cost,sint,eps,del,zero,one,two,&
       cosp,sinp,cos2t,sin2t,cos2p,sin2p,Nt2
                              
  real*8, dimension(2)     :: ax,ax_r,ap,ap_r,ap_2r,hp,hp_r,hp_2r,       &
                              bx,bp,cp,dp

  REAL*8, DIMENSION(2,3,3) :: Xa,Xb,Pa,Pb,Pc,Pd
  INTEGER                  :: i,j,k,l,ip,jp,lp,m,n,ml
  COMPLEX*16               :: xi
  COMPLEX*16, DIMENSION(5) :: tmp1
  REAL*8, DIMENSION(2,5)   :: v


!!$ Spherical coordinates


  zero = 0.0d0
  one  = 1.0d0
  two  = 2.0d0

  ispole = .FALSE.
  IF(z.GE.zero) THEN
     south = .FALSE.
  ELSE 
     south = .TRUE.
  ENDIF

  r = SQRT(x**2 + y**2 + z**2)
  theta = ACOS(z/r)
  phi = ATAN2(y,x)

  Nt2 = (1 - params(P_HOLE1_RAD)/r)

  cost = z/r
  sint = SIN(theta)
  cosp = COS(phi)
  sinp = SIN(phi)

  cos2t = cost**2 - sint**2
  sin2t = 2*sint*cost
  cos2p = cosp**2 - sinp**2
  sin2p = 2*sinp*cosp

! Index point (grid point where r_i < r < r_i+1

  i = INT((r-rr(1))/dr)+1

  eps = (r-rr(i))/dr
  del = one - eps

  Kc = zero

!!$ Stereographic coordinates, and forward/backward values

  xi = stereo(theta,phi,south)

!!$ Loop over l,m modes to reconstruct K,L,M in a spherical basis


  DO l= lmin,lmax
     ml = min(l,mmax)
     DO m=0,ml

! Interpolate moments for construction of K

        ax(1) = del*ao_p1(l,m,i)  + eps*ao_p1(l,m,i+1)
        ax(2) = del*ao_p1(l,-m,i) + eps*ao_p1(l,-m,i+1)
        ap(1) = del*ae_p1(l,m,i)  + eps*ae_p1(l,m,i+1)
        ap(2) = del*ae_p1(l,-m,i) + eps*ae_p1(l,-m,i+1)
        hp(1) = del*he_p1(l,m,i)   + eps*he_p1(l,m,i+1)
        hp(2) = del*he_p1(l,-m,i)  + eps*he_p1(l,-m,i+1)

! Interpolate numerical derivatives of these moments:
!             ax', ax'', a+', a+'', a+''', h', h'', h'''

        ax_r(1) = ( eps*ao_p1(l,m,i+2)  + del*ao_p1(l,m,i+1)               &
                  - eps*ao_p1(l,m,i)    - del*ao_p1(l,m,i-1)  )/(2*dr)
        ax_r(2) = ( eps*ao_p1(l,-m,i+2) + del*ao_p1(l,-m,i+1)              &
                  - eps*ao_p1(l,-m,i)   - del*ao_p1(l,-m,i-1) )/(2*dr)

        ap_r(1) = ( eps*ae_p1(l,m,i+2)  + del*ae_p1(l,m,i+1)               &
                  - eps*ae_p1(l,m,i)    - del*ae_p1(l,m,i-1)  )/(2*dr)
        ap_r(2) = ( eps*ae_p1(l,-m,i+2) + del*ae_p1(l,-m,i+1)              &
                  - eps*ae_p1(l,-m,i)   - del*ae_p1(l,-m,i-1) )/(2*dr)

        ap_2r(1) = (   eps*ae_p1(l,m,i+2)  + (del-2*eps)*ae_p1(l,m,i+1)     &
             + (eps-2*del)*ae_p1(l,m,i)    +         del*ae_p1(l,m,i-1) ) / &
             (dr**2)
        ap_2r(2) = (   eps*ae_p1(l,-m,i+2) + (del-2*eps)*ae_p1(l,-m,i+1)    &
             + (eps-2*del)*ae_p1(l,-m,i)   +         del*ae_p1(l,-m,i-1) )/ &
             (dr**2)

        hp_r(1) = ( eps*he_p1(l,m,i+2)  + del*he_p1(l,m,i+1)               &
                  - eps*he_p1(l,m,i)    - del*he_p1(l,m,i-1)  )/(2*dr)
        hp_r(2) = ( eps*he_p1(l,-m,i+2) + del*he_p1(l,-m,i+1)              &
                  - eps*he_p1(l,-m,i)   - del*he_p1(l,-m,i-1) )/(2*dr)

        hp_2r(1) = (   eps*he_p1(l,m,i+2)  + (del-2*eps)*he_p1(l,m,i+1)     &
             + (eps-2*del)*he_p1(l,m,i)    +         del*he_p1(l,m,i-1) ) / &
             (dr**2)
        hp_2r(2) = (   eps*he_p1(l,-m,i+2) + (del-2*eps)*he_p1(l,-m,i+1)    &
             + (eps-2*del)*he_p1(l,-m,i)   +         del*he_p1(l,-m,i-1) )/ &
             (dr**2)


! Compute other variables from constraints:

        if(l.eq.0) then
           ax = zero
           bx = zero
           bp = zero
           dp = zero
           cp = (hp - ap)/2
        else if(l.eq.1) then
           bx = zero
           bp = (3*ap + r*ap_r - hp - r*hp_r)/(l*l+l)
           cp = zero
           dp = (ap - hp)/2
        else 
           bx = -( (1+3*Nt2)*ax + 2*Nt2*r*ax_r )/(l*l+l-2)
           bp =  ( 3*ap + r*ap_r - hp - r*hp_r )/(l*l+l)
           cp =  ( ( 5*(1+3*Nt2) - 2*(l*l+l) )*ap - (3+5*Nt2)*hp       &
                 + (1+13*Nt2)*r*ap_r - (1+9*Nt2)*r*hp_r                &
                 + 2*Nt2*r*r*(ap_2r - hp_2r) )/(2*(l*l+l-2))
           dp =  ( ap + 2*cp - hp )/(l*l+l)
        end if

! Angular dependence

        CALL harmonics(l,m,xi,south,ispole,tmp1)

        v(1,:) = REAL(tmp1)
        v(2,:) = AIMAG(tmp1)


! Spherical basis tensors for K
! (first index denotes real or imaginary part)
! Note that while there is an r-dependence in these tensors, there
!      is no r-dependence in the cartesian basis

        Xa(:,1,1) = -2*sint*cosp*( v(:,2)*sinp + v(:,4)*cost*cosp )
        Xa(:,1,2) =  sint*( v(:,2)*cos2p - v(:,4)*cost*sin2p )
        Xa(:,1,3) = - ( v(:,2)*cost*sinp + v(:,4)*cosp*cos2t )
        Xa(:,2,2) =  2*sint*sinp*( v(:,2)*cosp - v(:,4)*cost*sinp )
        Xa(:,2,3) =  v(:,2)*cost*cosp - v(:,4)*cos2t*sinp
        Xa(:,3,3) =  v(:,4)*sin2t
        Xa(:,2,1) = Xa(:,1,2)
        Xa(:,3,1) = Xa(:,1,3)
        Xa(:,3,2) = Xa(:,2,3)


        Xb(:,1,1) = v(:,5)*( (cost**2)*(cosp**2) - sinp**2 )           &
                  + (l*(l+1)/2 * v(:,1) + v(:,3)) * cost*sin2p
        Xb(:,1,2) = v(:,5)*(1+cost**2)*cosp*sinp                       &
                  - (l*(l+1)/2 * v(:,1) + v(:,3)) * cost*sin2p
        Xb(:,1,3) = -v(:,5) * sint*cost*cosp                           &
                  - (l*(l+1)/2 * v(:,1) + v(:,3)) * sint*sinp
        Xb(:,2,2) = v(:,5)*( (cost**2)*(sinp**2) - cosp**2 )           &
                  - (l*(l+1)/2 * v(:,1) + v(:,3)) * cost*sin2p
        Xb(:,2,3) = -v(:,5) * sint*cost*sinp                           &
                  + (l*(l+1)/2 * v(:,1) + v(:,3)) * sint*cosp
        Xb(:,3,3) = v(:,5) * sint**2
        Xb(:,2,1) = Xb(:,1,2)
        Xb(:,3,1) = Xb(:,1,3)
        Xb(:,3,2) = Xb(:,2,3)

        Pa(:,1,1) = v(:,1) * (sint**2) * (cosp**2)
        Pa(:,1,2) = v(:,1) * (sint**2) * (cosp*sinp)
        Pa(:,1,3) = v(:,1) * (sint*cost) * cosp
        Pa(:,2,2) = v(:,1) * (sint**2) * (sinp**2)
        Pa(:,2,3) = v(:,1) * (sint*cost) * sinp
        Pa(:,3,3) = v(:,1) * (cost**2)
        Pa(:,2,1) = Pa(:,1,2)
        Pa(:,3,1) = Pa(:,1,3)
        Pa(:,3,2) = Pa(:,2,3)
        
        Pb(:,1,1) = 2*sint*cosp * ( v(:,2)*cost*cosp - v(:,4)*sinp )
        Pb(:,1,2) = sint* ( v(:,2)*cost*sin2p + v(:,4)*cos2p )
        Pb(:,1,3) = v(:,2)*cos2t*cosp - v(:,4)*cost*sinp
        Pb(:,2,2) = 2*sint*sinp * ( v(:,2)*cost*sinp + v(:,4)*cosp )
        Pb(:,2,3) = v(:,2)*cos2t*sinp + v(:,4)*cost*cosp
        Pb(:,3,3) = -v(:,2) * sin2t
        Pb(:,2,1) = Pb(:,1,2)
        Pb(:,3,1) = Pb(:,1,3)
        Pb(:,3,2) = Pb(:,2,3)

        Pc(:,1,1) =  v(:,1) * ( (cost**2)*(cosp**2) + sinp**2 )
        Pc(:,1,2) = -v(:,1) * (sint**2)*cosp*sinp
        Pc(:,1,3) = -v(:,1) * cost*sint*cosp
        Pc(:,2,2) =  v(:,1) * ( (cost**2)*(sinp**2) + cosp**2 )
        Pc(:,2,3) = -v(:,1) * cost*sint*sinp
        Pc(:,3,3) =  v(:,1) * sint**2
        Pc(:,2,1) = Pc(:,1,2)
        Pc(:,3,1) = Pc(:,1,3)
        Pc(:,3,2) = Pc(:,2,3)
        
        Pd(:,1,1) = v(:,3) * (cost**2)*(cosp**2) - v(:,5)*cost*sin2p    &
                  - (l*(l+1)*v(:,1) + v(:,3)) * sinp**2
        Pd(:,1,2) = (l*(l+1)*v(:,1) + v(:,3)*(1 + cost**2)) * cosp*sinp &
                  + v(:,5) * cost*cos2p
        Pd(:,1,3) = sint * ( -v(:,3)*cost*cosp + v(:,5)*sinp )
        Pd(:,2,2) = v(:,3) * (cost**2)*(sinp**2) + v(:,5)*cost*sin2p    &
                  - (l*(l+1)*v(:,1) + v(:,3)) * cosp**2
        Pd(:,2,3) = sint * (  v(:,3)*cost*sinp - v(:,5)*cosp )
        Pd(:,3,3) = v(:,3) * sint**2
        Pd(:,2,1) = Pd(:,1,2)
        Pd(:,3,1) = Pd(:,1,3)
        Pd(:,3,2) = Pd(:,2,3)

! Construct K 

        DO n = 1,2
           Kc = Kc +   ax(n)/r * Xa(n,:,:) + bx(n)/r * Xb(n,:,:)  &
                   + ap(n)/Nt2 * Pa(n,:,:) +   bp(n) * Pb(n,:,:)  &
                   +     cp(n) * Pc(n,:,:) +   dp(n) * Pd(n,:,:)
        ENDDO

     ENDDO
  ENDDO

  RETURN

CONTAINS

  COMPLEX*16 FUNCTION stereo(theta,phi,south)

    REAL*8, INTENT(in)      :: theta, phi
    LOGICAL, INTENT(in)     :: south
    REAL*8                  :: cost

    cost = COS(theta)

    IF(south) THEN
       stereo = SQRT((1+cost)/(1-cost))*EXP(CMPLX(zero,one)*phi)
    ELSE
       stereo = SQRT((1-cost)/(1+cost))*EXP(CMPLX(zero,one)*phi)
    ENDIF

    RETURN

  END FUNCTION stereo

END SUBROUTINE reconstruct
