    !First version AL October 2010

    !Calculates local fnl and CMB lensing bispectra
    !CMB lensing bispectra are calculated in the approximation in which the first-order
    !result is used replacing the unlensed with the lensing power spectra
    !This is non-perturbatively correct to about 1%
    !Note the lensing bispectrum only includes the linear potentials, no Rees-Sciama or SZ
    !See Lewis, Challinor & Hanson 2010 for details

    !Note that the primordial local bispectra are unlensed (see arXiv: 0905.4732)

    !Compile with LAPACK and -DFISHER if you want to get the fisher matrix outputs
    !This is disabled by default in order not to require LAPACK

    !Fisher results are with and without the cosmic variance from low L_1

    module Bispectrum
    use ModelParams
    use ModelData
    use InitialPower
    use SpherBessels
    use IniFile
    implicit none

    integer, parameter :: max_bispectrum_deltas = 5, max_bispectrum_fields=3

    Type TBispectrumParams
        logical do_lensing_bispectrum
        logical do_primordial_bispectrum
        integer nfields
        integer Slice_Base_L, ndelta, deltas(max_bispectrum_deltas)
        logical do_parity_odd
        logical DoFisher
        logical export_alpha_beta
        real(dl) FisherNoise, FisherNoisePol, FisherNoiseFwhmArcmin
        character(LEN=Ini_max_string_len)  FullOutputFile
        logical SparseFullOutput
    end Type

    !global parameter for now, only intend for this module to be used interactively for the moment
    Type(TBispectrumParams)  :: BispectrumParams

    Type TBispectrum
        real(dl), pointer :: b(:,:)
    end Type
    Type TCov
        real(dl), pointer :: C(:,:)
    end type TCov

    Type TCov2
        real(dl) :: C(2,2)
    end type TCov2

    Type TCov3
        real(dl) :: C(3,3)
    end type TCov3

    real(dl), allocatable :: dJl(:,:), dddJl(:,:)
    real(dl), parameter :: InternalScale = 1d10
    character(LEN=1024) :: output_root =''
    integer, parameter :: shape_local = 1, shape_warm = 2, shape_warm2 = 3
    !integer, parameter :: shape = shape_local
    !PDM 2016 adding shapes
    integer, parameter :: shape_equil = 4, shape_ortho = 5, shape_flat = 6
    integer ::  shape = shape_local

    real(dl), allocatable :: TransferPolFac(:)      !sqrt((l+2)!/(l-2)!)
    contains

    subroutine InitBesselDerivs(CTrans)
    ! j_l' array for interpolation if needed; not needed for local fnl
    Type(ClTransferData) :: CTrans
    integer i,l1,j
    real(dl) Jm, Jp

    if (allocated(dJl)) then
        deallocate(dJL, dddJl)
    end if
    allocate(dJl(BessRanges%npoints,CTrans%ls%l0),dddJl(BessRanges%npoints,CTrans%ls%l0))

    do i=1, CTrans%ls%l0
        !Spline agrees well
        !  call spline_deriv(BessRanges%points,ajl(1,i),ajlpr(1,i),dJl(1,i),BessRanges%npoints)
        !  call spline(BessRanges%points,dJl(1,i),BessRanges%npoints,spl_large,spl_large,dddJl(1,i))

        l1 = CTrans%ls%l(i)
        do j=1, BessRanges%npoints
            call BJL(l1-1,BessRanges%points(j),Jm)
            call BJL(l1+1,BessRanges%points(j),Jp)
            dJl(j,i)=  ( l1*Jm - (l1+1)*Jp)/(2*l1+1)
        end do
        call spline(BessRanges%points,dJl(1,i),BessRanges%npoints,spl_large,spl_large,dddJl(1,i))

    end do

    end subroutine InitBesselDerivs

    subroutine NonGauss_l_r_localOpt(CTrans, ind, indP, res, resP, nfields, r)
    !functions of the form int dk k^2 k^i j_l(kr) Delta_l(k) [P]
    !ind and indP are arrays of required k^i powers
    !res and resP are the results without and with the power spectrum P in the integrand
    Type(ClTransferData) :: CTrans
    integer, intent(in) :: ind(:), indP(:)
    integer :: nfields
    real(dl) res(CTrans%ls%l0,size(ind),nfields), resP(CTrans%ls%l0,size(indP),nfields)
    real(dl), intent(in) :: r
    integer q_ix, j, bes_ix
    integer n, nP, ellmax
    real(dl) xf , J_l, fac, a2, k, dlnk, term, P, kpow, kpowP

    n = size(ind)
    nP =size(indP)
    res=0
    resP = 0
    do q_ix = 1, CTrans%q%npoints
        k = CTrans%q%points(q_ix)
        xf = k*r
        bes_ix=Ranges_indexOf(BessRanges,xf)
        fac=BessRanges%points(bes_ix+1)-BessRanges%points(bes_ix)
        a2=(BessRanges%points(bes_ix+1)-xf)/fac
        fac=fac**2*a2/6
        dlnk = CTrans%q%dpoints(q_ix) /k
        P = ScalarPower(k, 1)*InternalScale  !!only first index for now

        ellmax = max(xf/(1-xlimfrac), xf + xlimmin) * AccuracyBoost
        kpow =  k**(ind(1)+3)
        kpowP = k**indP(1) * P
        do j=1,CTrans%ls%l0
            if (CTrans%ls%l(j) <= ellmax) then
                J_l=a2*ajl(bes_ix,j)+(1-a2)*(ajl(bes_ix+1,j) - ((a2+1) &
                    *ajlpr(bes_ix,j)+(2-a2)*ajlpr(bes_ix+1,j))* fac) !cubic spline
                term = CTrans%Delta_p_l_k(1,j,q_ix)*J_l*dlnk
                res(j,1,1) = res(j,1,1) + term * kpow
                resP(j,1,1) = resP(j,1,1) + term * kpowP
                if (nfields>1) then
                    !E pol
                    term = CTrans%Delta_p_l_k(2,j,q_ix)*J_l*dlnk * TransferPolFac(CTrans%ls%l(j))
                    res(j,1,2) = res(j,1,2) + term * kpow
                    resP(j,1,2) = resP(j,1,2) + term * kpowP
                    if (nfields>2) then
                        !lensing potential
                        term = CTrans%Delta_p_l_k(3,j,q_ix)*J_l*dlnk
                        res(j,1,3) = res(j,1,3) + term * kpow
                        resP(j,1,3) = resP(j,1,3) + term * kpowP
                    end if
                end if

            end if
        end do
    end do
    resP = resP * fourpi
    res = res * 2/pi

    end subroutine NonGauss_l_r_localOpt

    subroutine NonGauss_l_r(CTrans, ind, indP,res, resP,nfields, r)
    !functions of the form int dk k^2 k^i j_l(kr) Delta_l(k) [P]
    !ind and indP are arrays of required k^i powers
    !res and resP are the results without and with the power spectrum P in the integrand
    !Output of P scaled by 1d10 (so bispectrum by 1d20)
    Type(ClTransferData) :: CTrans
    integer:: nfields
    integer, intent(in) :: ind(:), indP(:)
    real(dl) res(CTrans%ls%l0,size(ind),nfields), resP(CTrans%ls%l0,size(indP),nfields)
    real(dl), intent(in) :: r
    integer q_ix, j, bes_ix, i
    integer n, nP, ellmax
    real(dl) xf , J_l, fac, a2, k, dlnk, term, P, kpow(size(ind)), kpow2(size(indP))

    if (shape == shape_local) then
        call NonGauss_l_r_localOpt(CTrans, ind, indP,res, resP, nfields, r)
        return
    end if

    n = size(ind)
    nP =size(indP)
    res=0
    resP = 0
    do q_ix = 1, CTrans%q%npoints
        k = CTrans%q%points(q_ix)
        xf = k*r
        bes_ix=Ranges_indexOf(BessRanges,xf)
        fac=BessRanges%points(bes_ix+1)-BessRanges%points(bes_ix)
        a2=(BessRanges%points(bes_ix+1)-xf)/fac
        fac=fac**2*a2/6
        dlnk = CTrans%q%dpoints(q_ix) /k
        P = ScalarPower(k, 1)*InternalScale  !!only first index for now

        ellmax = max(xf/(1-xlimfrac), xf + xlimmin) * AccuracyBoost
        do i=1,n
            kpow(i)=k**(ind(i)+3)
        end do
        do i=1,nP
            kpow2(i)=k**indP(i) * P
        end do

        do j=1,CTrans%ls%l0
            if (CTrans%ls%l(j) <= ellmax) then

                J_l=a2*ajl(bes_ix,j)+(1-a2)*(ajl(bes_ix+1,j) - ((a2+1) &
                    *ajlpr(bes_ix,j)+(2-a2)*ajlpr(bes_ix+1,j))* fac) !cubic spline
                !call BJL(CTrans%ls%l(j), xf, J_l)
                term = CTrans%Delta_p_l_k(1,j,q_ix)*J_l*dlnk
                do i=1,n
                    res(j,i,1) = res(j,i,1) + term *kpow(i)
                end do
                do i=1,nP
                    resP(j,i,1) = resP(j,i,1) + term * kpow2(i)
                end do
                !     if (CTrans%ls%l(j)==8) write (1,'(9D20.7)') &
                !      k, xf, real(term * k**3/dlnk), real(term * k**indP(1) * P), &
                !      real(res(j,1)),real(resP(j,1)), J_l, real(term), real(CTrans%Delta_p_l_k(1,j,q_ix))
                if (nfields>1) then
                    !E pol
                    term = CTrans%Delta_p_l_k(2,j,q_ix)*J_l*dlnk* TransferPolFac(CTrans%ls%l(j))
                    do i=1,n
                        res(j,i,2) = res(j,i,2) + term *kpow(i)
                    end do
                    do i=1,nP
                        resP(j,i,2) = resP(j,i,2) + term * kpow2(i)
                    end do
                    if (nfields>2) then
                        !lensing potential
                        term = CTrans%Delta_p_l_k(3,j,q_ix)*J_l*dlnk
                        do i=1,n
                            res(j,i,3) = res(j,i,3) + term *kpow(i)
                        end do
                        do i=1,nP
                            resP(j,i,3) = resP(j,i,3) + term * kpow2(i)
                        end do
                    end if
                end if

            end if
        end do
    end do
    !resP = resP * fourpi
    !res = res * 2/pi

    end subroutine NonGauss_l_r


    subroutine GetBispectrum(CTrans)
      !Note: may need high maxetak to make sure oscillatory k integrals cancel correctly
      !for accurate alpha(r), beta(r), e.g. 8000; not so important for bispectrum
      !increase accuracy_boost
      use lensing
      use lvalues
      use constants
      use Ranges
#ifdef FISHER
      use MatrixUtils
#endif          
      integer, parameter :: max_bispectra = 2  !fnl, lensing
      Type(ClTransferData) :: CTrans
      Type(Regions) :: TimeStepsNongauss
      integer, allocatable ::  ind(:), indP(:), indPd(:)
      real(dl), allocatable :: res(:,:,:), resP(:,:,:), resPd(:,:)
      real(dl), allocatable :: res_l(:,:,:), resP_l(:,:,:), resPd_l(:,:)
      real(dl) r, term
      Type(TBispectrum), target,allocatable :: Bispectra(:,:,:,:), OddBispectra (:)
      !TTT, TTE, etc; last index is bispectrum kind, default 1=fnl, 2=lensing
      !OddBispectra are parity odd terms like TEB (if do_parity_odd requested for lensing)
      Type(TBispectrum), pointer :: Bispectrum
      !For use in Fisher approximations
      real(dl) test(lmin:CTrans%ls%l(CTrans%ls%l0))
      integer i, j, l1,l2,l3, il1, n,np, npd
      integer min_l, max_l, lmax
      real(dl) tmp, tmp1, tmp2, tmp3
      !PDM 2016: more temps!
      real(dl) tmp4, tmp5, tmp6, tmp7
      integer :: ellmaxl, ellminl
      !PDM 2016
      real(dl) a3j(0:CTrans%ls%l(CTrans%ls%l0)*2+1)
      real(dl) a3j2(0:CTrans%ls%l(CTrans%ls%l0)*2+1,4,2)
      real(dl) CLForLensingIn(4,lmin:CTrans%ls%l(CTrans%ls%l0)),CPhi(3,lmin:CTrans%ls%l(CTrans%ls%l0))
      Type(lSamples) :: SampleL
      real starttime
      real(dl) Bscale
      integer field, field1,field2,field3, bi_ix,bix
      Type(TCov2), allocatable :: CForLensing(:)
      integer nfields, nbispectra, bispectrum_type
      integer :: fnl_bispectrum_ix = 1
      integer :: lens_bispectrum_ix = 2
      character(LEN=256) ::  file_tag = ''
      integer idelta, fileid
      character(LEN=26) :: BispectrumNames(max_bispectra)
      integer :: parities(3), oddix
      integer, parameter :: lmax_lensing_corrT = 300
      !assume C^{T\psi} zero above this for CMB lensing; also neglect lensing contributions to variance
      integer, parameter :: lmax_lensing_corrE = 40 !assume C^{E\psi} zero above this for CMB lensing

      integer, parameter :: first_order_unlensed = 1, first_order_lensed = 2
      integer, parameter :: lens_bispectrum_approx = first_order_lensed
#ifdef FISHER
      Type(TBispectrum), pointer :: Bispectrum2
      real(dl) Cl(4,lmin:CTrans%ls%l(CTrans%ls%l0))
      real(dl)  a3j_00(0:CTrans%ls%l(CTrans%ls%l0)*2+1)
      integer lstart
      real(dl) Noise, NoiseP, bias
      real(dl), allocatable:: fish_contribs(:,:,:)
      real(dl) fish_contribs_sig(lmin:CTrans%ls%l(CTrans%ls%l0))
      real(dl), allocatable :: ifish_contribs(:,:,:,:,:), Fisher(:,:), tmpFisher(:,:),OptimalFisher(:,:)
      real(dl), allocatable :: tmpBigFisher(:,:), Fisher_L1(:,:,:),tmpProjFisher(:,:)
      real(dl), allocatable :: fish_l1(:,:,:,:), fish_L_ij(:,:), fish_L_noise(:,:)
      Type(TBispectrum), target,allocatable :: SqueezedLensingKernel(:,:)
      real(dl) sigma2, xlc, tmpf(3)
      integer  f1,f2,f3, minl2,bigi,bigj, bispectrum_type2, lmaxcuti
      integer sz,corrsize
      Type(TCov), allocatable :: InvC(:)
      integer ix1,ix2
      real(dl) tmpArr(lmin:CTrans%ls%l(CTrans%ls%l0))
#endif

      parities(1)=1  !T
      parities(2)=1  !E
      parities(3)=-1 !B


      if (BispectrumParams%do_primordial_bispectrum) then
         fnl_bispectrum_ix = 1
         nbispectra=1
         BispectrumNames(fnl_bispectrum_ix)='fnl'
      else
         fnl_bispectrum_ix = 0
         nbispectra=0
      end if
      if (BispectrumParams%do_lensing_bispectrum) then
         lens_bispectrum_ix = fnl_bispectrum_ix+1
         nbispectra=nbispectra+1
         BispectrumNames(lens_bispectrum_ix)='lensing'
      end if
      if (nbispectra>max_bispectra) call MpiStop('check max_bispectra')

      if (CP%InitPower%nn>1) call MpiStop('Bispectrum: multiple initial power spectra not supported')

      nfields=BispectrumParams%nfields

      if (lSampleBoost <50) call MpiStop('Bispectrum assumes lSampleBoost=50 (all L sampled)')

      if (lens_bispectrum_approx == first_order_unlensed) file_tag='_unlens'

      lmax = CTrans%ls%l(CTrans%ls%l0)
      if (CP%DoLensing) lmax = lmax_lensed
      SampleL%l0=0
      l1=1
      do
         if (l1<=lmax_lensing_corrE) then
            l1 = l1+1
         else if (l1<120) then
            l1 =l1+nint(7/AccuracyBoost)
         else
            l1 =l1+nint(50/AccuracyBoost)
         end if
         if (l1>lmax) then
            l1 =lmax
         end if
         if (BispectrumParams%Slice_Base_L>0 .and. SampleL%l0>0) then
            !Make sure requested slice base is actually calculated
            if ( BispectrumParams%Slice_Base_L <l1 .and. BispectrumParams%Slice_Base_L>SampleL%l(SampleL%l0)) then
               SampleL%l0= SampleL%l0 + 1
               SampleL%l(SampleL%l0) = BispectrumParams%Slice_Base_L
            end if
         end if
         SampleL%l0= SampleL%l0 + 1
         ! print *,l1
         SampleL%l(SampleL%l0) = l1
         if (l1 == lmax) exit
      end do

      allocate(Bispectra(nfields,nfields,nfields,nbispectra))
      do field1=1,nfields
         do field2=1,nfields
            do field3=1,nfields
               !Only store l2,l3 that are non-zero, array size is approx
               do bispectrum_type=1,nbispectra
                  allocate(Bispectra(field1,field2,field3,bispectrum_type)%b((lmax*(lmax+1))/4,SampleL%l0))
                  Bispectra(field1,field2,field3,bispectrum_type)%b=0
               end do
            end do
         end do
      end do

      if (BispectrumParams%do_lensing_bispectrum) then

         if (.not. CP%DoLensing) call MpiStop('Must turn on lensing to get lensing bispectra')
         print *,'Getting lensing reduced bispectra'

         allocate(CForLensing(lmax))

         CPhi=0
         do i=lmin,lmax
            CPhi(1,i) = Cl_scalar(i,1,C_Phi)/real(i,dl)**4 * InternalScale
            !set correlations to zero where very small to avoid numerical issues
            if (i<=lmax_lensing_corrT) then
               CPhi(2,i) = Cl_scalar(i,1,C_PhiTemp) /real(i,dl)**3 * InternalScale
            end if
            if (i<=lmax_lensing_corrE) then
               CPhi(3,i) = Cl_scalar(i,1,C_PhiE) /real(i,dl)**3 * InternalScale
            end if
            tmp = i*(i+1)/(2*pi)
            CLForLensingIn(:,i) = CL_lensed(i,1,CT_Temp:CT_Cross) * InternalScale/tmp
            ! CForLensing(i)%C=0
            CForLensing(i)%C(1,1)=CLForLensingIn(1,i)
            CForLensing(i)%C(1,2)=CLForLensingIn(4,i)
            CForLensing(i)%C(2,1)=CLForLensingIn(4,i)
            CForLensing(i)%C(2,2)=CLForLensingIn(2,i)
            ! CForLensing(i)%C(3,3)=CL_lensed(i,1,CT_B) * InternalScale/tmp
         end do

#ifdef FISHER             
         allocate(SqueezedLensingKernel(nfields,nfields))
         do field2=1,nfields
            do field3=1,nfields
               allocate(SqueezedLensingKernel(field2,field3)%b((lmax*(lmax+1))/4,SampleL%l0))
               SqueezedLensingKernel(field2,field3)%b=0
            end do
         end do
#endif
         if (DebugMsgs) starttime=GetTestTime()


         !$OMP PARAllEl DO DEFAUlT(SHARED),SCHEDULE(STATIC,3) &
         !$OMP PRIVATE(il1,l1,l2,l3,max_l,min_l,bix,bi_ix, tmp1,tmp2,tmp3), &
         !$OMP PRIVATE(field1,field2,field3, Bispectrum, a3j,a3j2)

         do il1= 1, SampleL%l0
            l1 = SampleL%l(il1)
            if (l1 > lmax_lensing_corrT) cycle !no exit in parallel loop
            tmp1=l1*(l1+1)
            bi_ix=0
            do l2= max(lmin,l1), lmax
               tmp2=l2*(l2+1)
               min_l = max(abs(l1-l2),l2)
               if (mod(l1+l2+min_l,2)/=0) then
                  min_l = min_l+1
               end if
               max_l = min(lmax,l1+l2)
               bix=bi_ix
               a3j2(:,:,1)=0.5d0
               if (nfields>1) then
                  call GetThreeJs(a3j(abs(l2-l1)),l1,l2,0,0)
                  call GetThreeJs(a3j2(max(2,abs(l2-l1)),1,2),l1,l2,2,0)
                  call GetThreeJs(a3j2(max(2,abs(l2-l1)),2,2),l1,l2,0,2)
                  call GetThreeJs(a3j2(max(0,abs(l2-l1)),3,2),l1,l2,2,-2)
                  do l3=min_l,max_l ,2
                     a3j2(l3,:,2) = a3j2(l3,:,2)/a3j(l3)*0.5d0
                  end do
               end if
               do field1=1,nfields
                  do field2=1,nfields
                     do field3=1,nfields
                        Bispectrum=> Bispectra(field1,field2,field3, lens_bispectrum_ix)
                        bi_ix=bix
                        do l3=min_l,max_l ,2
                           bi_ix=bi_ix+1
                           tmp3=l3*(l3+1)
                           !bispectrum is the reduced bispectrum
#ifdef FISHER             
                           if (field1==1) then
                              SqueezedLensingKernel(field2,field3)%b(bi_ix,il1)=  &
                                   (-tmp2+tmp3+tmp1)*(a3j2(l3,2,field2)*CForLensing(l3)%C(field2,field3)) + &
                                   (-tmp3+tmp1+tmp2)*(a3j2(l3,2,field3)*CForLensing(l2)%C(field3,field2))
                           end if
#endif            

                           Bispectrum%b(bi_ix,il1)=  &
                                (-tmp1+tmp2+tmp3) *  &
                                (a3j2(l3,1,field1)*CPhi(1+field2,l2)*CForLensing(l3)%C(field1,field3) + &
                                a3j2(l3,3,field1)*CPhi(1+field3,l3)*CForLensing(l2)%C(field1,field2) ) + &
                                (-tmp2+tmp3+tmp1)* &
                                (a3j2(l3,3,field2)*CPhi(1+field3,l3)*CForLensing(l1)%C(field2,field1) + &
                                a3j2(l3,2,field2)*CPhi(1+field1,l1)*CForLensing(l3)%C(field2,field3) ) + &
                                (-tmp3+tmp1+tmp2)* &
                                (a3j2(l3,2,field3)*CPhi(1+field1,l1)*CForLensing(l2)%C(field3,field2) + &
                                a3j2(l3,1,field3)*CPhi(1+field2,l2)*CForLensing(l1)%C(field3,field1) )

                        end do

                     end do
                  end do
               end do

            end do
         end do
         !$OMP END PARAllEl DO
         if (DebugMsgs) print *,'Time for lensing:', GetTestTime()-starttime

         if (nfields==1) BispectrumParams%do_parity_odd=.false.

         if (BispectrumParams%do_parity_odd) then

            allocate(OddBispectra(12))
            oddix=0
            do field1=1,3
               do field2=1,3
                  do field3=1,3
                     if (parities(field1)+parities(field2)+parities(field3)/=1) cycle
                     oddix=oddix+1
                     !Only store l2,l3 that are non-zero, array size is approx
                     allocate(OddBispectra(oddix)%b((lmax*(lmax+1))/4,SampleL%l0))
                     OddBispectra(oddix)%b=0
                  end do
               end do
            end do

            !$OMP PARAllEl DO DEFAUlT(SHARED),SCHEDULE(STATIC,3) &
            !$OMP PRIVATE(il1,l1,l2,l3,max_l,min_l,bix,bi_ix, tmp1,tmp2,tmp3), &
            !$OMP PRIVATE(field1,field2,field3, Bispectrum, a3j2, oddix)

            do il1= 1, SampleL%l0
               l1 = SampleL%l(il1)
               if (l1 > lmax_lensing_corrT) cycle !no exit in parallel loop
               tmp1=l1*(l1+1)
               bi_ix=0
               do l2= max(lmin,l1), lmax
                  tmp2=l2*(l2+1)
                  min_l = max(abs(l1-l2),l2)
                  if (mod(l1+l2+min_l,2)/=1) then
                     min_l = min_l+1
                  end if
                  max_l = min(lmax,l1+l2)
                  bix=bi_ix

                  a3j2(:,:,1)=0.5d0
                  call GetThreeJs(a3j2(max(2,abs(l2-l1)),1,2),l1,l2,2,0)
                  call GetThreeJs(a3j2(max(2,abs(l2-l1)),2,2),l1,l2,0,2)
                  call GetThreeJs(a3j2(max(0,abs(l2-l1)),3,2),l1,l2,2,-2)
                  do l3=min_l,max_l ,2
                     a3j2(l3,:,2) = a3j2(l3,:,2)*0.5d0*sqrt(real((2*L1+1)*(2*L2+1),dl)*(2*L3+1)/(3.1415926535*4))
                  end do


                  oddix=0
                  do field1=1,3
                     do field2=1,3
                        do field3=1,3
                           !Only calculate terms with one B
                           if (parities(field1)+parities(field2)+parities(field3)/=1) cycle
                           oddix=oddix+1
                           Bispectrum=> OddBispectra(oddix)
                           bi_ix=bix
                           do l3=min_l,max_l ,2
                              bi_ix=bi_ix+1

                              tmp3=l3*(l3+1)
                              !bispectrum the non-reduced bispectrum without the i
                              if (parities(field1)==-1) then
                                 Bispectrum%b(bi_ix,il1)=  &
                                      (-tmp1+tmp2+tmp3) *  &
                                      (a3j2(l3,1,min(2,field1))*CPhi(1+field2,l2)*CForLensing(l3)%C(2,field3)  &
                                      -a3j2(l3,3,min(2,field1))*CPhi(1+field3,l3)*CForLensing(l2)%C(2,field2) )
                              elseif (parities(field2)==-1) then
                                 Bispectrum%b(bi_ix,il1)=  &
                                      (-tmp2+tmp3+tmp1)* &
                                      (-a3j2(l3,3,min(2,field2))*CPhi(1+field3,l3)*CForLensing(l1)%C(2,field1)  &
                                      -a3j2(l3,2,min(2,field2))*CPhi(1+field1,l1)*CForLensing(l3)%C(2,field3) )
                              else if (parities(field3)==-1) then
                                 Bispectrum%b(bi_ix,il1)=  &
                                      (-tmp3+tmp1+tmp2)* &
                                      (-a3j2(l3,2,min(2,field3))*CPhi(1+field1,l1)*CForLensing(l2)%C(2,field2) + &
                                      a3j2(l3,1,min(2,field3))*CPhi(1+field2,l2)*CForLensing(l1)%C(2,field1) )
                              end if

                           end do

                        end do
                     end do
                  end do

               end do
            end do
            !$OMP END PARAllEl DO
         end if

      end if

      if (BispectrumParams%do_primordial_bispectrum) then

         !print *,'getting reduced local fnl bispectra'

         allocate(TransferPolFac(lmax))
         do i=2,lmax
            TransferPolFac(i) =sqrt( real((i+1)*i,dl)*(i+2)*(i-1))
         end do

         !if (shape /= shape_local) call MpiStop('Non-local shapes not working')

         if (shape == shape_local) then
            print *,'getting reduced local bispectra'
            n=1
            np=1
            npd=0 !derivatives of function
         else if (shape == shape_warm) then
            n=2
            np=3
            npd=0
         else if (shape == shape_warm2) then
            n=1
            np=2
            npd=2
            !PDM 2016: adding 3 other shapes.
            !for now I am doing ns = 1 forecasting
            !Can not get code to work differently for now.
            !should be fine for rough numbers. 
         elseif (shape == shape_equil) then
            print *,'getting reduced equilateral bispectra'
            n = 4
            np = 0
            npd = 0
         elseif (shape == shape_ortho) then
            print *,'getting reduced orthogonal bispectra'
            n = 4
            np = 0
            npd = 0
         elseif (shape == shape_flat) then
            print *,'getting reduced flat bispectra'
            n = 4
            np = 0
            npd = 0
         else
            call MpiStop('unknown shape')
         end if

!!$         if (shape == shape_local) then
!!$            n=1
!!$            np=1
!!$            npd=0 !derivatives of function
!!$         else if (shape == shape_warm) then
!!$            n=2
!!$            np=3
!!$            npd=0
!!$         else if (shape == shape_warm2) then
!!$            n=1
!!$            np=2
!!$            npd=2
!!$         else
!!$            call MpiStop('unknown shape')
!!$         end if

         allocate(ind(n))
         allocate(indP(np))

         if (npd>0) then
            call InitBesselDerivs(CTrans)
            allocate(indPd(npd))
         end if

         if (shape==shape_warm) then
            !Separable form is very unstable and unworkable probably
            ind(1) = 0
            ind(2) = 2
            indP(1) = 0
            indP(2) = 2           
            indP(3) = -2
         else if (shape==shape_warm2) then
            ind(1) = 0
            indP(1) = 0
            indP(2) = -2
            indPd(1) = 0
            indPd(2) = -2
            !PDM 2016 Including more shapes 
         elseif (shape==shape_local) then
            ind(1) = 0
            indP(1) = 0 !k^0 P = P^1 = k^(-4+ns)
         elseif (shape==shape_equil) then
            !Should eventually code this up, now use below (ns = 1):
            !ind(1) = 0  !k^3 dlnk =  k^2 dk
            !indP(1) = 3 !k^3 dlnk P^(3/3) = k^2 dk P^1 = k^(-4+ns) k^2 dk =  k^(-2+ns) dk
            !indP(2) = 1 !k^3 dlnk P^(1/3) = k^2 dk P^1/3 = k^((-4+ns)/3) k^2 dk
            !indP(3) = 2 !k^3 dlnk P^(2/3) = k^2 dk P^1/3 = k^(2(-4+ns)/3) k^2 dk
            ind(1) =  0
            ind(2) = -1
            ind(3) = -2
            ind(4) = -3
         elseif (shape==shape_ortho) then
            ind(1) =  0
            ind(2) = -1
            ind(3) = -2
            ind(4) = -3
         elseif (shape==shape_flat) then
            ind(1) =  0
            ind(2) = -1
            ind(3) = -2
            ind(4) = -3
         end if

!!$         if (shape==shape_warm) then
!!$            !Separable form is very unstable and unworkable probably
!!$            ind(1) = 0
!!$            ind(2) = 2
!!$            indP(1) = 0
!!$            indP(2) = 2
!!$            indP(3) = -2
!!$         else if (shape==shape_warm2) then
!!$            ind(1) = 0
!!$            indP(1) = 0
!!$            indP(2) = -2
!!$            indPd(1) = 0
!!$            indPd(2) = -2
!!$         else
!!$            ind(1) = 0
!!$            indP(1) = 0
!!$         end if

         test=0
         call Ranges_Nullify(TimeStepsNongauss)
         call Ranges_Assign(TimeStepsNongauss,TimeSteps)
         call Ranges_Add_delta(TimeStepsNongauss, -taurst*10*AccuracyBoost, taurst, dtaurec)
         call Ranges_getArray(TimeStepsNongauss, .true.)

         !$        if (BispectrumParams%export_alpha_beta) call OMP_SET_NUM_THREADS(1)
         if (BispectrumParams%export_alpha_beta) then
            !Note that all the points outside recombination are not really needed
            !And these are for curvature perturbation, so do not include 3/5 factor
            call CreateTxtFile(trim(output_root)//'_alpha.txt',100)
            call CreateTxtFile(trim(output_root)//'_beta.txt',101)
            call CreateTxtFile(trim(output_root)//'_alpha_beta_r.txt',102)
         end if

         if (DebugMsgs) starttime=GetTestTime()

!!$         !$OMP PARALLEL DO DEFAUlT(SHARED),SCHEDULE(STATIC,3) &
!!$         !$OMP PRIVATE(i,r,res,resP,resPd,res_l,resP_l,resPd_l,term,j), &
!!$         !$OMP PRIVATE(il1,l1,l2,l3,min_l,max_l,tmp,tmp1,tmp2,Bispectrum), &
!!$         !$OMP PRIVATE(bi_ix,bix,field1,field2,field3,field)
         
         !$OMP PARALLEL DO DEFAUlT(SHARED),SCHEDULE(STATIC,3) &
         !$OMP PRIVATE(i,r,res,resP,resPd,res_l,resP_l,resPd_l,term,j), &
         !$OMP PRIVATE(il1,l1,l2,l3,min_l,max_l,tmp,tmp1,tmp2,tmp4,tmp5,tmp6,tmp7,Bispectrum), &
         !$OMP PRIVATE(bi_ix,bix,field1,field2,field3,field)

         do i= TimeStepsNongauss%npoints-1, 2,-1
            r=(CP%tau0-TimeStepsNongauss%points(i))

            allocate(res(CTrans%ls%l0,n,nfields))
            allocate(resP(CTrans%ls%l0,np,nfields))

            allocate(res_l(1:CTrans%ls%l(CTrans%ls%l0),n,nfields))
            allocate(resP_l(1:CTrans%ls%l(CTrans%ls%l0),np,nfields))
            if (npd>0) then
               allocate(resPd(CTrans%ls%l0,npd))
               allocate(resPd_l(1:CTrans%ls%l(CTrans%ls%l0),npd))
            end if

            call NonGauss_l_r(CTrans, ind, indP,res, resP, nfields, r)
            if (npd>0) call NonGauss_deriv_l_r(CTrans, indPd,resPd, r, dJl,dddJl)

            do field=1,nfields
               do j=1,n
                  call InterpolateClArr(CTransScal%ls,res(1,j,field),res_l(lmin,j,field),CTransScal%ls%l0)
               end do
               do j=1,np
                  call InterpolateClArr(CTransScal%ls,resP(1,j,field),resP_l(lmin,j,field),CTransScal%ls%l0)
               end do
            end do
            deallocate(res,resP)

            if (BispectrumParams%export_alpha_beta) then
               write(100,concat('(',lmax-lmin+1 ,'E15.5)')) res_l(lmin:lmax,1,1)
               write(101,concat('(',lmax-lmin+1 ,'E15.5)')) resP_l(lmin:lmax,1,1)
               write(102,'(1E15.5)') r
            end if

            if (npd>0) then
               do j=1,npd
                  call InterpolateClArr(CTransScal%ls,resPd(1,j),resPd_l(lmin,j),CTransScal%ls%l0)
               end do
               deallocate(resPd)
            end if

            term = r**2 * TimeStepsNongauss%dpoints(i) * (3./5)


            !Restrict to l1<=l2<=l3
            do il1= 1, SampleL%l0
               l1 = SampleL%l(il1)
               bi_ix=0
               do l2= max(lmin,l1), lmax
                  min_l = max(abs(l1-l2),l2)
                  if (mod(l1+l2+min_l,2)/=0) then
                     min_l = min_l+1
                  end if
                  max_l = min(lmax,l1+l2)
                  do field1=1,nfields
                     do field2=1,nfields
!!$                        tmp1 = 2*term*(res_l(l1,1,field1)*resP_l(l2,1,field2) + &
!!$                             res_l(l2,1,field2)*resP_l(l1,1,field1))
!!$                        tmp2 = 2*term*(resP_l(l1,1,field1)*resP_l(l2,1,field2))
!!$                        do field3=1,nfields
!!$                           Bispectrum => Bispectra(field1,field2,field3,fnl_bispectrum_ix)
!!$                           bix=bi_ix
!!$                           do l3=min_l,max_l ,2
!!$                              bix=bix+1
!!$                              Bispectrum%b(bix,il1) = Bispectrum%b(bix,il1) + &
!!$                                   (tmp1*resP_l(l3,1,field3) +   tmp2*res_l(l3,1,field3))
!!$                           end do
!!$                        end do
                        !PDM 2016: adding other shapes
                        if(shape == shape_local) then 
                           tmp1 = 2*term*(res_l(l1,1,field1)*resP_l(l2,1,field2) + &
                                res_l(l2,1,field2)*resP_l(l1,1,field1))
                           tmp2 = 2*term*(resP_l(l1,1,field1)*resP_l(l2,1,field2))
                        elseif(shape == shape_equil) then
                           !extra 2 because of equileteral normalization
                           tmp1 = -6*term*(res_l(l1,1,field1)*res_l(l2,4,field2) + &
                                res_l(l2,1,field2)*res_l(l1,4,field1))
                           tmp2 = -6*term*(res_l(l1,4,field1)*res_l(l2,4,field2))
                           !multiplicity 2 (since it can be T or E). All terms are the same
                           !so no two mirrored terms as in local as above
                           !tmp3 already used
                           tmp4 = -12*term*(res_l(l1,3,field1)*res_l(l2,3,field2))
                           tmp5 = 6*term*(res_l(l1,2,field1)*res_l(l2,3,field2) + &
                                res_l(l2,2,field2)*res_l(l1,3,field1)) 
                           tmp6 = 6*term*(res_l(l1,3,field1)*res_l(l2,4,field2) + &
                                res_l(l2,3,field2)*res_l(l1,4,field1)) 
                           tmp7 = 6*term*(res_l(l1,2,field1)*res_l(l2,4,field2) + &
                                res_l(l2,2,field2)*res_l(l1,4,field1))
                        elseif(shape == shape_flat) then
                           tmp1 = 6*term*(res_l(l1,1,field1)*res_l(l2,4,field2) + &
                                res_l(l2,1,field2)*res_l(l1,4,field1))
                           tmp2 = 6*term*(res_l(l1,4,field1)*res_l(l2,4,field2))

                           tmp4 = 18*term*(res_l(l1,3,field1)*res_l(l2,3,field2))
                           tmp5 = -6*term*(res_l(l1,2,field1)*res_l(l2,3,field2) + &
                                res_l(l2,2,field2)*res_l(l1,3,field1)) 
                           tmp6 = -6*term*(res_l(l1,3,field1)*res_l(l2,4,field2) + &
                                res_l(l2,3,field2)*res_l(l1,4,field1)) 
                           tmp7 = -6*term*(res_l(l1,2,field1)*res_l(l2,4,field2) + &
                                res_l(l2,2,field2)*res_l(l1,4,field1))                     
                        elseif(shape == shape_ortho) then
                           tmp1 = -18*term*(res_l(l1,1,field1)*res_l(l2,4,field2) + &
                                res_l(l2,1,field2)*res_l(l1,4,field1))
                           tmp2 = -18*term*(res_l(l1,4,field1)*res_l(l2,4,field2))

                           tmp4 = -48*term*(res_l(l1,3,field1)*res_l(l2,3,field2))
                           tmp5 = 18*term*(res_l(l1,2,field1)*res_l(l2,3,field2) + &
                                res_l(l2,2,field2)*res_l(l1,3,field1)) 
                           tmp6 = 18*term*(res_l(l1,3,field1)*res_l(l2,4,field2) + &
                                res_l(l2,3,field2)*res_l(l1,4,field1)) 
                           tmp7 = 18*term*(res_l(l1,2,field1)*res_l(l2,4,field2) + &
                                res_l(l2,2,field2)*res_l(l1,4,field1))
                        endif
                        do field3=1,nfields
                           Bispectrum => Bispectra(field1,field2,field3,fnl_bispectrum_ix)
                           bix=bi_ix    
                           do l3=min_l,max_l ,2
                              bix=bix+1
                              !PDM 2016: adding other shapes
                              if(shape == shape_local) then
                                 Bispectrum%b(bix,il1) = Bispectrum%b(bix,il1) + &
                                      (tmp1*resP_l(l3,1,field3) +   tmp2*res_l(l3,1,field3))
                              else!equilateral, orthogonal and flat
                                 Bispectrum%b(bix,il1) = Bispectrum%b(bix,il1) + &
                                      (21)**2*32*pi*(tmp1*res_l(l3,4,field3) +   tmp2*res_l(l3,1,field3)  + &
                                      tmp4*res_l(l3,3,field3)  +   tmp5*res_l(l3,4,field3) + &
                                      tmp6*res_l(l3,2,field3)  +   tmp7*res_l(l3,3,field3))
                              endif
                           end do !l3
                        end do !field3
                     end do !field2
                  end do !field1
                  bi_ix=bix

               end do !l2
            end do !il1

            deallocate(res_l,resP_l)
            if (npd>0) deallocate(resPd_l)
         end do !TimeStepsNongauss
         !$OMP END PARAllEl DO
         if (BispectrumParams%export_alpha_beta) then
            close(100)
            close(101)
            close(102)
         end if
         deallocate(TransferPolFac)
         call Ranges_Free(TimeStepsNongauss)

         if (DebugMsgs) print *,'Time for fnl bispectrum:', GetTestTime()-starttime

      end if !DoPrimordial

      if (BispectrumParams%Slice_Base_L>0 .or. BispectrumParams%FullOutputFile/='') then
         !write out slice in (muK)^3 units
         Bscale=(COBE_CMBTemp*1d6)**3/InternalScale**2;
         do bispectrum_type=1,nbispectra
            if (BispectrumParams%Slice_Base_L>0) then
               do idelta=1,BispectrumParams%ndelta
                  if (mod(BispectrumParams%Slice_Base_L + BispectrumParams%deltas(idelta),2)==1 &
                       .and. bispectrum_type/=lens_bispectrum_ix) cycle
                  call CreateTxtFile(concat(trim(output_root)//'bispectrum_'//trim(BispectrumNames(bispectrum_type))//'_base_', &
                       BispectrumParams%Slice_Base_L,'_delta_',BispectrumParams%deltas(idelta),trim(file_tag)//'.dat'),&
                       nbispectra +BispectrumParams%ndelta*(bispectrum_type-1)+idelta)
               end do
            end if
            if (BispectrumParams%FullOutputFile/='') then
               call CreateTxtFile(concat(output_root,BispectrumParams%FullOutputFile, &
                    '_', BispectrumNames(bispectrum_type), file_tag, '.dat'),bispectrum_type)
            end if
         end do
         do il1= 1, SampleL%l0
            l1 = SampleL%l(il1)
            bi_ix=0
            do l2= max(lmin,l1), lmax
               min_l = max(abs(l1-l2),l2)
               if (mod(l1+l2+min_l,2)/=0) then
                  min_l = min_l+1
               end if
               max_l = min(lmax,l1+l2)
               do l3=min_l, max_l ,2
                  bi_ix=bi_ix+1
                  if (l1==BispectrumParams%Slice_Base_L &
                       .and. any(l3-l2==BispectrumParams%deltas(1:BispectrumParams%ndelta))) then
                     !Particular slice
                     idelta=IndexOf(l3-l2,BispectrumParams%deltas,BispectrumParams%ndelta)
                     do bispectrum_type=1,nbispectra
                        fileid=nbispectra +BispectrumParams%ndelta*(bispectrum_type-1)+idelta
                        write (fileid,'(1I5)', advance='NO') L2
                        do field1=1,nfields
                           do field2=1,nfields
                              do field3=1,nfields
                                 write(fileid,'(1E15.5)', advance='NO') &
                                      Bispectra(field1,field2,field3,bispectrum_type)%b(bi_ix,il1)*Bscale
                              end do
                           end do
                        end do
                        write (fileid,'(a)') ''
                     end do
                  end if !slice
                  if (BispectrumParams%FullOutputFile/='') then
                     if (BispectrumParams%SparseFullOutput .and. .not. any( SampleL%l(1:SampleL%l0)==L2) .or. &
                          l1 > 30 .and. mod(l3-min_l,10)/=0 .and. l3 /= max_l) cycle

                     do bispectrum_type=1,nbispectra
                        if (bispectrum_type==lens_bispectrum_ix .and. L1 > lmax_lensing_corrT) cycle
                        write(bispectrum_type,'(3I5)', advance='NO') L1, L2, L3
                        do field1=1,nfields
                           do field2=1,nfields
                              do field3=1,nfields
                                 write(bispectrum_type,'(1E14.5)', advance='NO') &
                                      Bispectra(field1,field2,field3,bispectrum_type)%b(bi_ix,il1)*Bscale
                              end do
                           end do
                        end do
                        write (bispectrum_type,'(a)') ''
                     end do
                  end if

               end do
            end do
         end do
         if (BispectrumParams%do_parity_odd) then
            do il1= 1, SampleL%l0
               l1 = SampleL%l(il1)
               bi_ix=0
               do l2= max(lmin,l1), lmax
                  min_l = max(abs(l1-l2),l2)
                  if (mod(l1+l2+min_l,2)/=1) then
                     min_l = min_l+1
                  end if
                  max_l = min(lmax,l1+l2)
                  do l3=min_l, max_l ,2
                     bi_ix=bi_ix+1
                     if (l1==BispectrumParams%Slice_Base_L &
                          .and. any(l3-l2==BispectrumParams%deltas(1:BispectrumParams%ndelta))) then
                        !Particular slice
                        idelta=IndexOf(l3-l2,BispectrumParams%deltas,BispectrumParams%ndelta)
                        do bispectrum_type=1,nbispectra
                           if (bispectrum_type/=lens_bispectrum_ix) cycle
                           fileid=nbispectra +BispectrumParams%ndelta*(bispectrum_type-1)+idelta
                           write (fileid,'(1I5)', advance='NO') L2
                           oddix=0
                           do field1=1,3
                              do field2=1,3
                                 do field3=1,3
                                    if (parities(field1)+parities(field2)+parities(field3)/=1) cycle
                                    oddix=oddix+1
                                    write(fileid,'(1E15.5)', advance='NO') OddBispectra(oddix)%b(bi_ix,il1)*Bscale
                                 end do
                              end do
                           end do
                           write (fileid,'(a)') ''
                        end do
                     end if !slice

                  end do
               end do
            end do
         end if

         do bispectrum_type=1,nbispectra
            if (BispectrumParams%Slice_Base_L>0) then
               do idelta=1,BispectrumParams%ndelta
                  if (mod(BispectrumParams%Slice_Base_L + BispectrumParams%deltas(idelta),2)==1 &
                       .and. bispectrum_type/=lens_bispectrum_ix) cycle
                  close(nbispectra +BispectrumParams%ndelta*(bispectrum_type-1)+idelta)
               end do
            end if
            if (BispectrumParams%FullOutputFile/='') close(bispectrum_type)
         end do

      end if

#ifdef FISHER
      if (BispectrumParams%DoFisher) then
         !Get stuff for Fisher etc.

         print *,'Getting Fisher for lmax = ', lmax

         Noise = BispectrumParams%FisherNoise/ (COBE_CMBTemp*1e6)**2  !Planckish, dimensionless units
         NoiseP = BispectrumParams%FisherNoisePol/ (COBE_CMBTemp*1e6)**2

         do i=lmin,lmax
            if (CP%DoLensing) then
               cl(:,i) = CL_lensed(i,1,CT_Temp:CT_Cross)
            else
               cl(1,i) = CL_Scalar(i,1,C_Temp)
               cl(2,i) = CL_Scalar(i,1,C_E)
               cl(4,i) = CL_Scalar(i,1,C_Cross)
               cl(3,i) = 0
            end if
            if (CP%WantTensors .and. i<= CP%Max_l_tensor .and. i>=2) then
               cl(:,i) = cl(:,i) + Cl_tensor(i,1,CT_Temp:CT_Cross)
            end if
         end do
         if (.false.) then
            call OpenTxtFile('CAMBdefault_lensedCls.dat',3)
            do i=lmin,lmax
               !Assume T,E,B,X ordering
               read(3,*) j, cl(1:4,i)
               if (j<lmin) read(3,*) j, cl(1:4,i)
               cl(:,i)=cl(:,i)/(COBE_CMBTemp*1e6)**2
            end do
            close(3)
         end if

         if (Noise >0) then
            file_tag = concat(file_tag,'_noise')
         end if
         xlc= 180*sqrt(8.*log(2.))/3.14159
         sigma2 = (BispectrumParams%FisherNoiseFwhmArcmin/60/xlc)**2
         allocate(InvC(lmax))
         do l1= lmin, lmax
            tmp = l1*(l1+1)/(2*pi)
            !PDM adding noise to fake foregrounds to TT:
            if (Noise .eq. 0 .and. l1 .ge. 3000) then
               Noise = 1e-3/(COBE_CMBTemp*1e6)**2
               sigma2 = (1/60/xlc)**2
            endif
            Cl(1,l1) = Cl(1,l1)/tmp + Noise*exp(l1*(l1+1)*sigma2)

            Cl(2:3,l1) = Cl(2:3,l1)/tmp + NoiseP*exp(l1*(l1+1)*sigma2)
            Cl(4,l1) = Cl(4,l1)/tmp
            allocate(InvC(l1)%C(nfields,nfields))
            if (nfields > 2) call MpiStop('Not implemented nfields>2 in detail')
            if (nfields==1) then
               InvC(l1)%C(1,1)=(2*l1+1)/cl(1,l1)/InternalScale
            else
               InvC(l1)%C(1,1)=cl(2,l1)
               InvC(l1)%C(1,2)=-cl(4,l1)
               InvC(l1)%C(2,1)=-cl(4,l1)
               InvC(l1)%C(2,2)=cl(1,l1)
               InvC(l1)%C= InvC(l1)%C * (2*l1+1)/(cl(1,l1)*cl(2,l1)-cl(4,l1)**2)/InternalScale
            end if
         end do

         if (debugMsgs) starttime=GetTestTime()
         allocate(ifish_contribs(SampleL%l0,nbispectra,nbispectra,nfields,nfields) )
         !This loop is just in case want to plot out lmax dependence
         do lmaxcuti=SampleL%l0, SampleL%l0
            !    call CreateTxtFile('TE-Planck-LensFish.txt',20)
            !   do lmaxcuti=1, SampleL%l0
            !   if (SampleL%l(lmaxcuti) < 425) cycle

            lmax= SampleL%l(lmaxcuti)

            ifish_contribs=0
            lstart = 2 !lmin
            !$OMP PARAllEl DO DEFAUlT(SHARED),SCHEDULE(STATIC,3) &
            !$OMP PRIVATE(il1,l1,l2,l3,fish_l1,bi_ix,min_l,max_l,a3j_00,a3j), &
            !$OMP PRIVATE(Bispectrum,Bispectrum2,minl2,bix,tmp,tmp1,tmp2,tmpf), &
            !$OMP PRIVATE(field1,field2,field3,f1,f2,f3,bispectrum_type,bispectrum_type2)

            do il1= 1,  lmaxcuti !!!SampleL%l0
               allocate(fish_l1(nbispectra,nbispectra,nfields,nfields)) !last indices are field1,f1
               l1 = SampleL%l(il1)
               if (l1< lstart)  then
                  deallocate(fish_L1)
                  cycle
               endif
               fish_l1=0
               bi_ix=0
               do l2 = l1,lmax
                  if (l2< lstart) cycle
                  min_l = max(lstart,max(abs(l1-l2),l2))
                  if (mod(l1+l2+min_l,2)/=0) then
                     min_l = min_l+1
                  end if
                  max_l = min(lmax,l1+l2)
                  call GetThreeJs(a3j(abs(l2-l1)),l1,l2,0,0)
                  do l3=min_l,max_l ,2
                     a3j_00(l3)=a3j(l3)**2
                  end do

                  tmp1= 1.d0/(4*pi)  !(2l+1) factors included in InvC
                  minl2=min_l
                  bix=bi_ix
                  do field1=1,nfields
                     do f1=1,nfields
                        tmpf(1)= InvC(l1)%C(field1,f1)*tmp1
                        do field2=1,nfields
                           do f2=1,nfields
                              tmpf(2)= InvC(l2)%C(field2,f2)*tmpf(1)
                              do field3=1,nfields
                                 do bispectrum_type=1,nbispectra
                                    if (bispectrum_type==lens_bispectrum_ix) then
                                       Bispectrum=>SqueezedLensingKernel(field2,field3)
                                    else
                                       Bispectrum=>bispectra(field1,field2,field3,bispectrum_type)
                                    end if
                                    do f3=1,nfields

                                       do bispectrum_type2=bispectrum_type,nbispectra
                                          if (bispectrum_type2==lens_bispectrum_ix) then
                                             Bispectrum2=>SqueezedLensingKernel(f2,f3)
                                          else
                                             Bispectrum2=>Bispectra(f1,f2,f3,bispectrum_type2)
                                          end if


                                          min_l=minl2
                                          bi_ix=bix
                                          if (min_l==l2) then
                                             !Symmetry factors
                                             bi_ix=bi_ix+1
                                             l3=l2
                                             if (l2==l1) then
                                                !l1=l2=l3
                                                tmp = Bispectrum%b(bi_ix,il1)*tmpf(2)*Bispectrum2%b(bi_ix,il1) &
                                                     *InvC(l3)%C(field3,f3)*a3j_00(l3)/6
                                             else
                                                !l3=l2 (l3=l1<>l2 can't happen because l1<=l2<=l3)
                                                tmp = Bispectrum%b(bi_ix,il1)*tmpf(2)*Bispectrum2%b(bi_ix,il1) &
                                                     * InvC(l3)%C(field3,f3)*a3j_00(l3)/2
                                             end if
                                             min_l = min_l+2
                                          else
                                             tmp=0
                                          end if
                                          tmp2=0
                                          do l3=min_l,max_l ,2
                                             bi_ix=bi_ix+1
                                             tmp2 = tmp2 + Bispectrum%b(bi_ix,il1)*Bispectrum2%b(bi_ix,il1) &
                                                  * InvC(l3)%C(field3,f3)*a3j_00(l3)
                                          end do
                                          if (l2==l1) then
                                             tmp2=tmp2*tmpf(2)/2
                                          else
                                             tmp2=tmp2*tmpf(2)
                                          end if
                                          fish_l1(bispectrum_type,bispectrum_type2,field1,f1)= &
                                               fish_l1(bispectrum_type,bispectrum_type2,field1,f1)+(tmp+tmp2)

                                       end do !bispectrum_type2

                                    end do
                                 end do !bispectrum_type
                              end do !field3
                           end do !f2
                        end do !field2
                     end do !f1
                  end do  !field1

               end do !l2

               ifish_contribs(il1,:,:,:,:)=fish_l1(1:nbispectra,1:nbispectra,:,:) /InternalScale
               deallocate(fish_L1)

            end do
            !$OMP END PARAllEl DO
            if (DebugMsgs) print *,'Time for Fisher:', GetTestTime()-starttime


            allocate(fish_contribs(lmin:CTrans%ls%l(CTrans%ls%l0),nfields,nfields))
            allocate(Fisher(nbispectra,nbispectra))
            allocate(tmpFisher(nbispectra,nbispectra))
            allocate(Fisher_L1(lmin:CTrans%ls%l(CTrans%ls%l0),nbispectra*nfields,nbispectra*nfields))

            do bispectrum_type=1,nbispectra
               do bispectrum_type2=bispectrum_type,nbispectra

                  fish_contribs=0
                  do field1=1,nfields
                     do f1=1,nfields
                        call InterpolateClArr(SampleL,ifish_contribs(1,bispectrum_type,bispectrum_type2,field1,f1), &
                             fish_contribs(lmin,field1,f1),lmaxcuti)  !SampleL%l0)
                     end do
                  end do
                  Fisher(bispectrum_type,bispectrum_type2)=0
                  do i=lmin, CTrans%ls%l(CTrans%ls%l0)
                     do field1=1,nfields
                        do f1=1,nfields
                           Fisher_L1(i,(bispectrum_type-1)*nfields+field1,(bispectrum_type2-1)*nfields+f1) = &
                                fish_contribs(i,field1,f1)
                           Fisher_L1(i,(bispectrum_type2-1)*nfields+f1,(bispectrum_type-1)*nfields+field1) = &
                                fish_contribs(i,field1,f1)
                           tmp=fish_contribs(i,field1,f1)
                           if (bispectrum_type==lens_bispectrum_ix) then
                              tmp = tmp * CPhi(1+field1,i)
                           end if
                           if (bispectrum_type2==lens_bispectrum_ix) then
                              tmp = tmp * CPhi(1+f1,i)
                           end if
                           Fisher(bispectrum_type,bispectrum_type2)=Fisher(bispectrum_type,bispectrum_type2)+ tmp
                        end do
                     end do
                  end do
                  Fisher(bispectrum_type2,bispectrum_type) = Fisher(bispectrum_type,bispectrum_type2)
                  do il1= 1, lmaxcuti
                     do field1=1,nfields
                        do f1=1,nfields
                           write(*,*)  il1,field1,f1, fish_contribs(SampleL%l(il1),field1,f1), &
                                ifish_contribs(il1,bispectrum_type,bispectrum_type2,field1,f1)
                        enddo
                     enddo
                  enddo
                  print *,'Zero-signal Fisher ',trim(BispectrumNames(bispectrum_type)) &
                       //'-'//trim(BispectrumNames(bispectrum_type2)), &
                       ':', Fisher(bispectrum_type2,bispectrum_type)

!!!! contribution of lensing to the fnl variance for temperature:
                  !           if (bispectrum_type == fnl_bispectrum_ix .and. bispectrum_type2 == fnl_bispectrum_ix ) then
                  !              fish_contribs_sig=0
                  !              tmpArr=0
                  !              call InterpolateClArr(SampleL,ifish_contribs(1,1,2,1,1),tmpArr(lmin),lmaxcuti)
                  !               do i=lmin, lmax_lensing_corrT
                  !               if (CPhi(1+1,i)/=0) then
                  !                 fish_contribs_sig(i) = &
                  !                  tmpArr(i)**2*(1+ CPhi(1,i)*CForLensing(i)%C(1,1)/(CPhi(1+1,i)*CPhi(1+1,i)))/(2*i+1)
                  !               end if
                  !              end do
                  !              print *,'signal contribution to fnl variance', sum(fish_contribs_sig)/Fisher(bispectrum_type,bispectrum_type2)**2
                  !           end if
!!! same with polarization
                  !  if (bispectrum_type == fnl_bispectrum_ix .and. bispectrum_type2 == lens_bispectrum_ix ) then
                  !            fish_contribs_sig=0
                  !            do i=lmin, lmax_lensing_corrT
                  !             corrsize = count(CPhi(2:1+nfields,i)/=0)
                  !             allocate(fish_L_ij(corrsize,corrsize))
                  !             allocate(fish_L_noise(corrsize,corrsize))
                  !             fish_L_ij=0
                  !             ix1=0
                  !             do field1=1,nfields
                  !              ix1=ix1+1
                  !              ix2=0
                  !              do f1=1,nfields
                  !                 ix2=ix2+1
                  !                 fish_L_noise(ix1,ix2)= fish_contribs(i,field1, f1)
                  !                 fish_L_ij(ix1,ix2) = (CPhi(1+f1,i)*CPhi(1+field1,i) + &
                  !                     CPhi(1,i)*CForLensing(i)%C(field1,f1) )/(2*i+1)
                  !              end do
                  !             end do
                  !             fish_L_ij=matmul(matmul((fish_L_noise(:,:)),fish_L_ij),transpose(fish_L_noise(:,:)))
                  !             fish_contribs_sig(i) = sum(fish_L_ij)!Matrix_trace(fish_L_ij)
                  !             deallocate(fish_L_ij)
                  !             deallocate(fish_L_noise)
                  !            end do
                  !            print *,'fnl signal factor', sum(fish_contribs_sig)
                  !       end if


                  if (bispectrum_type == lens_bispectrum_ix .and. bispectrum_type2 == lens_bispectrum_ix ) then
                     print *,'doing signal part of the lensing variance'
                     fish_contribs_sig=0
                     do i=lmin, lmax_lensing_corrT
                        corrsize = count(CPhi(2:1+nfields,i)/=0)
                        allocate(fish_L_ij(corrsize,corrsize))
                        allocate(fish_L_noise(corrsize,corrsize))
                        fish_L_ij=0
                        ix1=0
                        do field1=1,nfields
                           if (CPhi(1+field1,i)/=0) then
                              ix1=ix1+1
                              ix2=0
                              do f1=1,nfields
                                 if (CPhi(1+f1,i)/=0) then
                                    ix2=ix2+1
                                    fish_L_noise(ix1,ix2)= fish_contribs(i,field1, f1)*CPhi(1+f1,i)*CPhi(1+field1,i)
                                    fish_L_ij(ix1,ix2) =  (1 + &
                                         CPhi(1,i)*CForLensing(i)%C(field1,f1)/(CPhi(1+f1,i)*CPhi(1+field1,i)))/(2*i+1)

                                 end if
                              end do
                           end if
                        end do
                        call Matrix_Inverse(fish_L_noise)
                        fish_L_ij=fish_L_ij+fish_L_noise
                        call Matrix_Inverse(fish_L_ij)
                        fish_contribs_sig(i) = sum(fish_L_ij)
                        deallocate(fish_L_ij)
                        deallocate(fish_L_noise)
                     end do
                     print *,'Lensing Fisher including lensing variance ', sum(fish_contribs_sig)

                  end if

               end do
            end do
            deallocate(fish_contribs)

            print *, 'Results assuming zero fiducial bispectra'
            do bispectrum_type=1,nbispectra
               print *,trim(IntToStr(bispectrum_type))//'-'//trim(BispectrumNames(bispectrum_type)), &
                    ': 1/sqrt(F_ii) = ',1/sqrt(Fisher(bispectrum_type,bispectrum_type))
            end do

            do bispectrum_type=1,nbispectra
               do bispectrum_type2=bispectrum_type+1,nbispectra
                  bias = Fisher(bispectrum_type2,bispectrum_type)/Fisher(bispectrum_type,bispectrum_type)
                  print *,'Bias of ',trim(BispectrumNames(bispectrum_type2)),' on ', &
                       trim(BispectrumNames(bispectrum_type)),':', bias
               end do
               !           if (bispectrum_type==1) write(1,concat('(1I10,',nbispectra,'E15.5)')) lmax, Fisher(1,1), &
               !                                  Fisher(2:nbispectra,1)/Fisher(1,1)
            end do

            tmpFisher=Fisher
            do bispectrum_type=1,nbispectra
               tmp = sqrt(tmpFisher(bispectrum_type,bispectrum_type))
               tmpFisher(bispectrum_type,:)=tmpFisher(bispectrum_type,:)/tmp
               tmpFisher(:,bispectrum_type)=tmpFisher(:,bispectrum_type)/tmp
            end do
            if (nbispectra>1) then
               print *,'Zero-signal Bispectrum correlation matrix:'
               do bispectrum_type=1,nbispectra
                  print *,tmpFisher(:,bispectrum_type)
               end do

               if (nbispectra > 1) then
                  allocate(OptimalFisher(nbispectra,nbispectra))
                  allocate(tmpBigFisher(nbispectra*nfields,nbispectra*nfields))
                  !              allocate(diag(nbispectra*nfields))
                  OptimalFisher=0
                  do i=lmin, lmax
                     tmpBigFisher=Fisher_L1(i,:,:)
                     if (lens_bispectrum_ix/=0 .and. i<=lmax_lensing_corrT) then

                        call Matrix_Inverse(tmpBigFisher)
                        do field1 = 1,nfields
                           do field2 =1, nfields
                              tmpBigFisher((lens_bispectrum_ix-1)*nfields+field1, &
                                   (lens_bispectrum_ix-1)*nfields+field2) = &
                                   tmpBigFisher((lens_bispectrum_ix-1)*nfields+field1, &
                                   (lens_bispectrum_ix-1)*nfields+field2) + &
                                   (CPhi(1+field1,i)*CPhi(1+field2,i) + CPhi(1,i)*CForLensing(i)%C(field1,field2))/(2*i+1)
                           end do
                        end do
                        call Matrix_Inverse(tmpBigFisher)

                     end if

                     do field1=1,nfields
                        do field2=1,nfields
                           do bispectrum_type=1,nbispectra
                              do bispectrum_type2=1,nbispectra
                                 if (bispectrum_type==lens_bispectrum_ix) then
                                    tmp=CPhi(1+field1,i)
                                 else
                                    tmp=1
                                 end if
                                 if (bispectrum_type2==lens_bispectrum_ix) then
                                    tmp=tmp*CPhi(1+field2,i)
                                 end if
                                 OptimalFisher(bispectrum_type,bispectrum_type2)=OptimalFisher(bispectrum_type,bispectrum_type2)+ &
                                      tmp*tmpBigFisher((bispectrum_type-1)*nfields+field1,(bispectrum_type2-1)*nfields+field2)
                              end do
                           end do
                        end do
                     end do
                  end do
                  deallocate(tmpBigFisher)
               end if

               do bispectrum_type=1,nbispectra
                  print *,'Optimal Inc. lensing:', trim(IntToStr(bispectrum_type))//'-'//trim(BispectrumNames(bispectrum_type)), &
                       ': 1/sqrt(F_ii) = ',1/sqrt(OptimalFisher(bispectrum_type,bispectrum_type))
               end do
               tmpFisher=OptimalFisher
               call Matrix_Inverse(tmpFIsher)
               do bispectrum_type=1,nbispectra
                  print *,'Optimal Inc. lensing:', trim(IntToStr(bispectrum_type))//'-'//trim(BispectrumNames(bispectrum_type)), &
                       ': Cov_ii = ', sqrt(tmpFIsher(bispectrum_type,bispectrum_type))
               end do
               tmpFisher=OptimalFisher
               do bispectrum_type=1,nbispectra
                  tmp = sqrt(tmpFisher(bispectrum_type,bispectrum_type))
                  tmpFisher(bispectrum_type,:)=tmpFisher(bispectrum_type,:)/tmp
                  tmpFisher(:,bispectrum_type)=tmpFisher(:,bispectrum_type)/tmp
               end do
               print *,'Optimal Bispectrum correlation matrix with lensing variance:'
               do bispectrum_type=1,nbispectra
                  print *,tmpFisher(:,bispectrum_type)
               end do
            end if

            deallocate(Fisher, tmpFisher)

         end do

         deallocate(ifish_contribs)
         deallocate(InvC)

      end if !DoFIsher
#else
      if (BispectrumParams%DoFisher) call MpiStop('compile with FISHER defined')
#endif                    

      !Tidy up a bit
      do field1=1,nfields
         do field2=1,nfields
            do field3=1,nfields
               do bispectrum_type=1,nbispectra
                  deallocate(Bispectra(field1,field2,field3,bispectrum_type)%b)
               end do
            end do
         end do
      end do
      deallocate(Bispectra)
      if (allocated(CForLensing)) deallocate(CForLensing)

      if (allocated(OddBispectra)) then
         oddix=0
         do field1=1,3
            do field2=1,3
               do field3=1,3
                  if (parities(field1)+parities(field2)+parities(field3)/=1) cycle
                  oddix=oddix+1
                  deallocate(OddBispectra(oddix)%b)
               end do
            end do
         end do
         deallocate(OddBispectra)
      end if

    end subroutine GetBispectrum


    !not needed for local NG
    subroutine NonGauss_deriv_l_r(CTrans, indP,resP, r, dJl, dddJl)
    !As above, but integral against derivative of bessel function to get derivative of function
    Type(ClTransferData) :: CTrans
    real(dl), intent(in) :: dJl(BessRanges%npoints,CTrans%ls%l0), dddJl(BessRanges%npoints,CTrans%ls%l0)
    integer, intent(in) :: indP(:)
    real(dl) resP(CTrans%ls%l0,size(indP))
    real(dl), intent(in) :: r
    integer q_ix, j, bes_ix, i
    integer nP, ellmax
    real(dl) xf , dJ_l, fac, a2,  k, dlnk, term, P

    nP =size(indP)
    resP = 0
    do q_ix = 1, CTrans%q%npoints

        k = CTrans%q%points(q_ix)
        xf = k*r  !kr
        bes_ix=Ranges_indexOf(BessRanges,xf)
        fac=BessRanges%points(bes_ix+1)-BessRanges%points(bes_ix)
        a2=(BessRanges%points(bes_ix+1)-xf)/fac
        fac=fac**2*a2/6
        dlnk = CTrans%q%dpoints(q_ix) /k
        P = ScalarPower(k, 1)  !!only first index for now
        ellmax = max(xf/(1-xlimfrac), xf + xlimmin) * AccuracyBoost

        do j=1,CTrans%ls%l0
            if (CTrans%ls%l(j) <= ellmax) then
                dJ_l=a2*djl(bes_ix,j)+(1-a2)*(djl(bes_ix+1,j) - ((a2+1) &
                    *dddjl(bes_ix,j)+(2-a2)*dddjl(bes_ix+1,j))* fac) !cubic spline

                term = CTrans%Delta_p_l_k(1,j,q_ix)*dJ_l*dlnk *k
                do i=1,nP
                    resP(j,i) = resP(j,i) + term * k**indP(i) * P
                end do

            end if
        end do
    end do
    resP = resP * fourpi

    end subroutine NonGauss_deriv_l_r

    subroutine Bispectrum_SetDefParams(B)
    Type(TBispectrumParams) :: B

    B%nfields=2
    B%Slice_Base_L=0
    B%ndelta=0
    B%DoFisher = .true.
    B%FisherNoise = 0 !2d-4 !Planckish
    B%FisherNoisePol = 4*B%FisherNoise
    B%FisherNoiseFwhmArcmin = 7
    B%FullOutputFile='' !quite large
    B%SparseFullOutput = .false.
    B%do_lensing_bispectrum = .true.
    B%do_primordial_bispectrum = .true.
    B%do_parity_odd = .false.
    B%export_alpha_beta = .false.

    end subroutine Bispectrum_SetDefParams


    subroutine Bispectrum_ReadParams(B, Ini, outroot)
    use IniFile
    Type(TBispectrumParams) :: B
    character(LEN=*), intent(in) :: outroot
    Type(TIniFile) :: Ini
    integer i

    call Bispectrum_SetDefParams(B)

    B%do_lensing_bispectrum = Ini_Read_Logical_File(Ini,'do_lensing_bispectrum',.false.)
    B%do_primordial_bispectrum = Ini_Read_Logical_File(Ini,'do_primordial_bispectrum',.false.)

    do_bispectrum= B%do_lensing_bispectrum .or. B%do_primordial_bispectrum

    if (do_bispectrum) then

        output_root = outroot

        B%nfields = Ini_Read_Int_File(Ini,'bispectrum_nfields',B%nfields)
        if (B%nfields /= 2 .and. B%nfields/=1) call MpiStop('Bispectrum: nfields=1 for T only or 2 for polarization')
        B%do_parity_odd = Ini_Read_Logical_File(Ini,'do_parity_odd',.false.)
        if (B%do_parity_odd .and. (.not.  B%do_lensing_bispectrum .or. B%nfields==1)) then
            B%do_parity_odd = .false.
            write(*,*) 'Ignoring do_parity_odd since do_lensing_bispectrum=F or no polarization'
        end if
        B%Slice_Base_L = Ini_Read_Int_File(Ini,'bispectrum_slice_base_L',B%Slice_Base_L)
        if (B%Slice_Base_L>0) then
            B%ndelta = Ini_Read_Int_File(Ini,'bispectrum_ndelta',B%ndelta)
            if (B%ndelta > max_bispectrum_deltas) call MpiStop('Bispectrum : increase max_bispectrum_deltas')
            do i=1, B%ndelta
                B%deltas(i) = Ini_Read_Int_Array_File(Ini,'bispectrum_delta', i)
            end do
            if (.not. B%do_parity_odd .and. B%Slice_Base_L>0 .and. &
                any(mod(B%Slice_Base_L + B%deltas(1:B%ndelta),2) /= 0)) &
                call MpiStop('Slice is zero for even parity with L1+L2+L3 odd, i.e. Base+DeltaL3 odd')

        end if
        B%DoFisher = Ini_Read_Logical_File(Ini,'bispectrum_do_fisher',B%DoFisher)
        if (B%DoFisher) then
            B%FisherNoise = Ini_Read_Double_File(Ini,'bispectrum_fisher_noise',B%FisherNoise)
            B%FisherNoisePol = Ini_Read_Double_File(Ini,'bispectrum_fisher_noise_pol',B%FisherNoisePol )
            B%FisherNoiseFwhmArcmin = Ini_Read_Double_File(Ini,'bispectrum_fisher_fwhm_arcmin',B%FisherNoiseFwhmArcmin)
        end if
        B%FullOutputFile = Ini_Read_String_File(Ini,'bispectrum_full_output_file')
        if (B%FullOutputFile /='') then
            B%SparseFullOutput = Ini_Read_Logical_file(Ini,'bispectrum_full_output_sparse',B%SparseFullOutput)
        end if

        if (B%do_primordial_bispectrum) then
            B%export_alpha_beta = Ini_Read_Logical_file(Ini,'bispectrum_export_alpha_beta', B%export_alpha_beta)
        end if

    end if

    end subroutine Bispectrum_ReadParams

    end module Bispectrum
