module XctCrossCalc_M
  use, intrinsic :: ISO_C_BINDING
  use CrossSectionCalculator_M
  use DerivativeHandler_M
  use SammyRMatrixParameters_M
  use CovarianceData_M
  use AdjustedRadiusData_M
  use AllocateFunctions_m
  implicit none

  type, extends(CrossSectionCalculator) :: XctCrossCalc
     ! do we need the cross sections for a given channel
     ! this is assuming first width is gamma and is not counted, and all other width are as in normal sammy definition
     logical,allocatable,dimension(:)::Ifcros(:)
     logical::separateIso

      type(DerivativeHandler)::crossDataSelf     ! the self indicated cross section and derivative calculated for a given energy
      logical,allocatable,dimension(:)::crossSelfWhy  ! reproduce a SAMMY bug for self-indication experiments. To Do fix the bug instead

      integer::lllmax = 0                  ! maximum number of Clebsch-Gordon coefficients
      logical::wantSelfIndicate = .false.  ! do we need to calculate self-indicated cross section data

      integer::Kpolar = 0                  ! are fission width given in polar form

      integer::ntotc    ! the maximum number of channels (excluding gamma width) in any spin group
      integer::ntriag   ! size of lower triangular matrix over width, i.e. (ntotc*(ntotc+1))/2

      logical::debug   ! do we want debug print out

      integer::Kslow = 0     ! should we use the slow or fast version of Clebsch-Gordon calculation
      integer::C_G_Kxlmn = 1

      integer::Ifcoul = 0   ! do we need to calculate coulomb data
      logical::addElimKapt  !ADD ELIMINATED CAPTURE CHANNEL TO FINA

      real(kind=8),allocatable,dimension(:,:)::Alj   ! Eqs. (II B1 b.13) and (II B1 b.4) for 1- and 2-channel case, respectively [SAMMY 8 Users' Guide]

! *** Alj = (2*J+1) * (2*el+1) *del2 (el, s, J) in one-channel case
! *** Alj = (2*J+1) * sqrt { (2*el +1) *del2 (el , s , J) } * 2
! ***               * sqrt { (2*el'+1) *del2 (el', s', J) } in 2-channel

      real(kind=8),allocatable,dimension(:)::Xx      ! 0.0 or if SHIFT RESONANCE ENERGIES VIA SHIFT FACTOR, the factor
      real(kind=8),allocatable,dimension(:)::XxHelper    ! helper array to calculate Xx
      real(kind=8),allocatable,dimension(:)::Xlmn
      integer,allocatable,dimension(:,:)::Ixlmn

     real(kind=8),allocatable,dimension(:,:)::Br, Bi    ! energy independent part of derivatives for resonance parameters
     real(kind=8),allocatable,dimension(:,:,:)::Bga !E ENERGY-INDEPENDENT PORTION OF PARTIAL of R wrt reduced-width-amplitudes (for those which are not varied)
     real(kind=8),allocatable,dimension(:,:,:)::Pgar, Pgai  ! ENERGY-INDEPENDENT PORTION OF PARTIAL of R wrt reduced-width-amplitudes (for those which are not varied)

     real(kind=8),allocatable,dimension(:)::Alphar  !  (DEL E) / ( (DEL E)**2 + (Gamgam/2)**2 )
     real(kind=8),allocatable,dimension(:)::Alphai  ! Gamgam/2 / ( (DEL E)**2 + (Gamgam/2)**2 )
     logical,allocatable,dimension(:)::needAlphai   ! is Alphai big enough that it needs to be calculated and added
     real(kind=8),allocatable,dimension(:)::Difen   ! resonance energy minus currrent energy or 0 if resonance is not included in calculation


     real(kind=8),allocatable,dimension(:,:)::Pi, Pr ! Penetrability, P, real and imaginary part thereof.

     real(kind=8),allocatable,dimension(:,:)::Cscs  ! Looks like Cos(a + b), Sin(a + b), where b may be phi
     real(kind=8),allocatable,dimension(:)::Sinsqr  ! sin^2( phase shift )
     real(kind=8),allocatable,dimension(:)::Sin2ph  ! sin( 2 * phase shift )
     real(kind=8),allocatable,dimension(:)::Dphi ! partial derivative of Phi wrt Rho
     real(kind=8),allocatable,dimension(:)::Sinphi, Cosphi ! ( sin, cos ) ( phase shift )
     real(kind=8),allocatable,dimension(:)::Dpdr ! partial P wrt Rho
     real(kind=8),allocatable,dimension(:)::Dsdr ! partial S wrt Rho

     real(kind=8),allocatable,dimension(:,:)::Rmat ! Rmat = SUM Beta*Beta/((DEL E)-i(Gamgam/2))
     real(kind=8),allocatable,dimension(:,:)::Ymat ! Ymat = Linv - Rmat , Linv = 1/(S-B+IP)
     real(kind=8),allocatable,dimension(:,:)::Yinv ! Yinv = 1 / Ymat
     real(kind=8),allocatable,dimension(:)::Rootp  ! Rmat = SUM Beta*Beta/((DEL E)-i(Gamgam/2))
     real(kind=8),allocatable,dimension(:)::Elinvr,Elinvi ! Real and imaginary part of Linv
     real(kind=8),allocatable,dimension(:)::Psmall ! Is Psmall used?

! The following comments were copied from the preamble of the Subroutine Setxqx in mxct09.f90: 
! ***            XQ   = Yinv * Rmat       and
! ***            XXXX = P/L + sqrt(P)/L         (1/L-R)**-1          sqrt(P)/L
! ***                 =       sqrt(P)/(S-B+IP) * Yinv       * Rmat * sqrt(P)
! ***                 =       sqrt(P)/L        * XQ                * sqrt(P)
!
! ***            Note that the matrix W defined in SAMMY manual is given
! ***                 by W(c,c') = delta(c,c') + 2i XXXX(c,c')
! ***                 as in Eq. (II.B1.3) in SAMMY manual R7
!
! ***         ie W    = I + 2i XXXX

     real(kind=8),allocatable,dimension(:)::Xxxxr, Xxxxi ! *** Xxxx = sqrt(P)/L  * xq * sqrt(P) ... symmetric
     real(kind=8),allocatable,dimension(:,:)::Xqr,Xqi    ! *** Xq   = (L**-1-R)**-1 * R ... note asymmetry
     real(kind=8),allocatable,dimension(:,:)::Pxrhor, Pxrhoi ! Pxrho_ = partial(Xxxx_) wrt (Rho)

! *** Purpose -- Generate QR,QI =
! ***                   SQRT(P)/(S-B+IP) * Yinv*Yinv * SQRT(P)/(S-B+IP)
! ***
! *** That is, Qr(KL,Ij) is (real part of) partial of XXXX(Kl) wrt R(Ij)

     real(kind=8),allocatable,dimension(:,:)::Qr, Qi ! see above

! ***    For angular distributions
! ***      prtl (1-U)(kl) wrt R(ij) = prtl(1-U)(kl) wrt X(kl) * prtl X wrt R
! ***      Tx(1,ij,kl) = prtl Re(1-U) wrt ReR = -prtl Im (1-U) wrt ImR
! ***      Tx(2,ij,kl) = prtl Re(1-U) wrt ImR =  prtl Re (1-U) wrt ReR

     real(kind=8),allocatable,dimension(:,:,:)::Tx ! see above

! *** generate Deriv(k,Itzero) & Deriv(k,Ilzero); ditto Derivx
! *** ie Deriv(k,i) = partial Crss(k) wrt either Tzero or LZero
! ***    Calculating Deriv(K,I?z) = (2T)(Pre)(Z)(E?z), where
! ***               Tr   + i Ti   = (partial sigmaX  wrt R) * Half
! ***               Prer + i Prei =  partial R       wrt E
! ***               Z             =  partial E       wrt sqrt(E) = 2*sqrt(E)
! ***               Etz           =  partial sqrt(E) wrt  tzero
! ***               ELz           =  partial sqrt(E) wrt eLzero
! ***    but do not include the E?z part yet, or missing {dgoj * 1/E^2},
! ***    and store in Ddddtl

     real(kind=8),allocatable,dimension(:,:)::Tr, Ti ! see above
     real(kind=8),allocatable,dimension(:,:)::Prei, Prer ! see above
     real(kind=8),allocatable,dimension(:)::Ddddd, Ddddtl ! see above

     logical,allocatable, dimension(:,:)::useChannel 

! *** generate derivatives of Crss & Crssx wrt rho
! *** Dsf (   I) = Deriv of Crss(1   ) wrt rho via phi(I)
! *** Dst (J ,I) = Deriv of Crss(J   ) wrt rho via P(I) & S(I)
! *** Dstt(Jj,I) = Deriv of Crss(Jj+2) wrt rho via P(I) & S(I)

     real(kind=8),allocatable,dimension(:)::Dsf  ! see above
     real(kind=8),allocatable,dimension(:,:)::Dstt, Dst ! see above
     real(kind=8),allocatable,dimension(:,:,:)::Dsfx, Dstx ! ***          Note that the "f" in Dsf is for "effective" rho, i.e., radius.
     real(kind=8),allocatable,dimension(:)::Xden ! The inverse of Eq. (II D1 b.7) = 1/[(E-E_resonance)^2 + (gamma/2)^2]

     real(kind=8),allocatable,dimension(:)::termf, termfx ! termf(1) = elastic, (2) = absorption, (3) = reaction channels; "x" number of parameters (Numdet) of the ETA detector efficiency is greater than 0. 

     real(kind=8),allocatable,dimension(:,:,:)::crossInternal           ! used to keep track of cross section+deriv for one row and all channels
     real(kind=8),allocatable,dimension(:,:,:,:,:)::angInternal         ! used to keep track of angula section+deriv for one row and all channels
     real(kind=8),allocatable,dimension(:)::uniqueEchan  ! the unique value of echan for a given isotope, used as a scratch array during calculation
     integer,allocatable,dimension(:)::iradIndex  ! map radius id

     real(kind=8),allocatable,dimension(:,:)::Ccoulx

     real(kind=8):: Etz   ! partial sqrt(E) wrt  tzero, energy dependent
     real(kind=8):: Elz   ! partial sqrt(E) wrt eLzero, energy dependent

     ! direct capture arrays. Todo update the direct capture to read and store data in C++ class
     logical::hasDirectCapture
     real(kind=8),allocatable,dimension(:,:)::Edrcpt 
     real(kind=8),allocatable,dimension(:,:)::Cdrcpt
     real(kind=8),allocatable,dimension(:)::Xdrcpt
     integer,allocatable,dimension(:)::Ndrcpt
     real(kind=8),allocatable,dimension(:)::DrcpCoef ! the coefficent for a given nuclde (Coff from table 11.11 in input)
     integer,allocatable,dimension(:)::IflDrcpCoef  ! should the coefficient be varied
     integer::numdrc  !  number of nuclides in the direct capture file
     integer::Nucdrc  ! max number of energy/value pairs
     ! end direct capture data
     contains
     procedure, pass(this) :: setUpDerivativeList => XctCrossCalc_setUpDerivativeList    ! set up  crossData, depending on number of isotopes
     procedure, pass(this) :: setAddtionalParams => XctCrossCalc_setAddtionalParams
     procedure, pass(this) :: setEnergyIndependent => XctCrossCalc_setEnergyIndependent  ! set energy independent values using current parameter values
     procedure, pass(this) :: getParamPerSpinGroup => XctCrossCalc_getParamPerSpinGroup  ! also fill useChannel
     procedure, pass(this) :: calcCross => XctCrossCalc_calcCross
     procedure, pass(this) :: initialize => XctCrossCalc_initialize
     procedure, pass(this) :: destroy => XctCrossCalc_destroy
  end type XctCrossCalc

contains
subroutine XctCrossCalc_setAddtionalParams(this,  lllmax, Kfinit, wantSelfIndicate, Kssmsc, debug)
     use mxct25_m
     use mxct30_m

     class(XctCrossCalc)::this
     integer::lllmax            ! maximum number of Clebsch-Gordon coefficients
     integer::Kfinit            ! finite-size corrections flag
     integer::Kssmsc            ! flag for scattering correction
     logical::wantSelfIndicate  ! do we need to calculate self-indicated cross section data
     integer::numIso
     logical::debug

     integer::ntotTrig, ng, dd

    this%debug = debug


     this%wantSelfIndicate = wantSelfIndicate
     if ( wantSelfIndicate) then
         numIso = this%crossData%getNumberIsotopes()
         if (this%crossDataSelf%getNumberIsotopes().ne.numIso) then
            call this%crossDataSelf%setUpList(this%resData, this%radiusData, numIso)
         end if
     end if


     if (Kfinit.ne.0) this%Ifcros = .true. ! need all cross sections
     this%lllmax = lllmax

     this%Kslow = 1
     this%C_G_Kxlmn = 1
     if (this%lllmax.gt.0) then
        IF (this%reactType.EQ.7 .OR. this%reactType.EQ.11 .OR. Kssmsc.GT.0) THEN
           if (.not.allocated(this%Alj)) then
               ng = this%resData%getNumSpinGroups()
               ntotTrig = (this%ntotc*(this%ntotc+1))/2
               allocate(this%Alj(ntotTrig, ng))
           end if
           CALL Kount_Clbsch (this%Zke, this%Alj, this%Kslow, this%resData, this%C_G_Kxlmn, this%Lllmax, this%ntotc)

           call allocate_real_data(this%Xlmn, this%C_G_Kxlmn)
           if ( allocated(this%Ixlmn)) then
                if (size(this%Ixlmn, dim=1).lt.this%C_G_Kxlmn.or. &
                    size(this%Ixlmn, dim=2).le.dd)  deallocate(this%Ixlmn)
           end if

           dd = 3
           IF (this%Kslow.EQ.0) dd = 1
           if (.not.allocated(this%Ixlmn)) then
               allocate(this%Ixlmn(this%C_G_Kxlmn,dd))
           end if
           IF (this%Kslow.EQ.0) THEN
              CALL Clbsch  ( this%resData, this%Xlmn, this%Alj, this%Ixlmn(:,1), this%C_G_Kxlmn, this%Lllmax, this%Ntotc,  this%Debug)
           ELSE
              CALL Clbsch_Slow (this%resData, this%Xlmn, this%Alj, this%Ixlmn, this%C_G_Kxlmn, this%Lllmax, this%Ntotc,  this%Debug)
           END IF
        END IF
    end if
end subroutine
subroutine XctCrossCalc_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)
    use mcro8_m
    use mxct26_m
    class(XctCrossCalc) :: this
    integer::kwcoul, reactType
    real(kind=8)::Twomhb, Etac

    call CrossSectionCalculator_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)

    ! set bounds
    call this%setBound(Twomhb, kwcoul, Etac)

    ! calculate energy shift if needed
    if (this%doResShift) then
         call Fixx (this, this%XX,  this%XxHelper)
    end if

    ! determine whether we need coulomb interaction
    call Find_If_Coulomb (this, this%IfCoul)

    if (this%Ifcoul.gt.0) then
        if (.not.allocated(this%Ccoulx)) then
           allocate(this%Ccoulx(this%ntotc, this%resData%getNumSpinGroups()))
        end if
        call Start_Coul (this%Zke, this%Ccoulx, this%resData)
    end if
end subroutine
subroutine XctCrossCalc_calcCross(this, ener, Ipoten)
    class(XctCrossCalc):: this
    real(kind=8)::ener
    integer::ipoten

    integer::itot

    call CrossSectionCalculator_calcCross(this, ener, Ipoten)
    if (this%wantSelfIndicate) then
       itot = this%covariance%getNumTotalParam() + 1
       if ( .not.this%wantDerivs) itot = 1
       call this%crossDataSelf%reserveColumns(this%row, itot)
    end if
end subroutine
subroutine XctCrossCalc_setUpDerivativeList(this)
   class(XctCrossCalc) :: this
   integer::itot, i, ii, Nnnsig,numIso
   type(SammySpinGroupInfo)::spinInfo

   call CrossSectionCalculator_setUpDerivativeList(this)

    if (this%wantSelfIndicate) then
        call this%crossDataSelf%nullify()
        call this%crossDataSelf%setNnsig(1)
    end if
end subroutine
subroutine XctCrossCalc_getParamPerSpinGroup(this, ires, igr,  fillIt)
    class(XctCrossCalc) :: this
    integer,intent(inout)::ires  ! starting resonance for this group
    integer,intent(in)::igr     ! spin group for which to get the number of varied parameters
    logical,intent(in)::fillIt

    type(SammySpinGroupInfo):: spinInfo
    type(SammyChannelInfo)::channelInfo
    integer::Ichan
    integer::k, nent
    logical::inc, fillLocal
    optional fillIt

    fillLocal = .true.
    if (present(fillIt)) fillLocal = fillIt
    call CrossSectionCalculator_getParamPerSpinGroup(this, ires, igr,  fillLocal)
    if (.not.fillLocal) return ! for getting the number of parameters, not for filling data

    this%useChannel(:,1) = .false.
    call this%resData%getSpinGroupInfo(spinInfo, igr)
    nent = spinInfo%getNumEntryChannels()

    ! 1-d data
    do ichan = 1, this%ntotc + 1
       if (.not.this%Ifcros(ichan)) cycle

       if (ichan.le.2) then
          this%useChannel(ichan, 1) = .true.
          cycle
        end if

        k = Ichan-2+Nent
        inc = .true.
        if (k.le.spinInfo%getNumChannels()) then
           call spinInfo%getChannelInfo(channelInfo, k)
           inc =  channelInfo%getIncludeInCalc()
        end if
        if (inc) then
           if (.not.this%addElimKapt) then
               this%useChannel(ichan, 1) = .true.
           end if
        else
           if (this%addElimKapt) then
              this%useChannel(ichan, 1) = .true.
           end if
        end if
    end do

    ! 2-d data
    if (.not.this%needAngular) return
    this%useChannel(:,2) = .true.
    do ichan = 1, this%ntotc + 1
        IF (this%reactType.ne.11) then  ! elastic
            if(Ichan.GT.Nent) this%useChannel(ichan, 2)  = .false. ! but not an elastic channel
       ELSE   !  reaction of some kind
          inc = .true.
          if (ichan.le.spinInfo%getNumChannels()) then
             call spinInfo%getChannelInfo(channelInfo, ichan)
             inc =  channelInfo%getIncludeInCalc()
             if (channelInfo%getExcludeCompletely()) then
                this%useChannel(ichan, 2) = .false. ! Do not want it anywhere
                cycle
             end if
          end if
          IF (Ichan.LE.Nent) THEN
              this%useChannel(ichan, 2) = .false.  ! Do not want elastic
          else if (channelInfo%getIncludeInCalc()) then
              if ( this%addElimKapt) this%useChannel(ichan, 2) = .false. !  Will subtract only excluded channels from absorption
          else  ! normally excclude but calc%addElimKapt.eq.true overrides exclusion
               IF (.not.this%addElimKapt)  this%useChannel(ichan, 2) = .false. !  Do not want excluded channel in final state
          END IF
       END IF
    end do
end subroutine
subroutine XctCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
    class(XctCrossCalc) :: this
    type(SammyRMatrixParameters)::pars
    type(CovarianceData)::cov
    type(AdjustedRadiusData)::rad
    integer,intent(in)::niso, Itzero, Ilzero
    type(VariedParameterInfo)::info
    integer::ntot
    logical,intent(in)::needAngular, doShiftRes

    integer::itot,  ipup, irad, nres, ngroup, igr, ichan, ifl
    logical(C_BOOL)::countCombined
    type(SammySpinGroupInfo)::spinInfo

    call CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
    this%separateIso = .true.
    this%addElimKapt = .false.
    if (niso.eq.1) this%separateIso = .false.
    call info%initialize(pars, cov, rad)

     ntot = info%getMaxChannel()
     if (ntot.eq.0) ntot = 1
     allocate(this%Ifcros(ntot))
     this%Ifcros = .true.  ! assume cross sections are to be calculated for all channels

    call allocate_real_data(this%termf, max(2,ntot))
    call allocate_real_data(this%termfx, max(2,ntot))

    this%ntotc = info%getMaxChannel() -1  ! for most of the internal arrays the gamma channel is not used in the internal arrays
    this%ntriag = (this%ntotc*(this%ntotc+1))/2
    countCombined = .true.
    call info%getNumAll(itot, ipup, countCombined)
    if(this%Itzero.gt.0 .or. this%Ilzero.gt.0) itot = itot + 2
    irad = this%radiusData%getNumTotalVaried()
    nres = this%resData%getNumResonances()
    ngroup = this%resData%getNumSpinGroups()
    call info%destroy()

    call allocate_real_data(this%Alphai, nres)
    call allocate_real_data(this%Alphar, nres)
    allocate(this%needAlphai(nres))
    call allocate_real_data(this%Difen, nres)
    call allocate_real_data(this%Xden, nres)

    if (this%Itzero.gt.0 .or. this%Ilzero.gt.0) then
       call reallocate_real_data_2d(this%Prer, this%ntriag, 0, ngroup, 0)
       call reallocate_real_data_2d(this%Prei, this%ntriag, 0, ngroup, 0)
    end if
    if (this%Itzero.gt.0 .or. this%Ilzero.gt.0.or.irad.gt.0) then
        allocate(this%Dstt(this%ntotc, this%ntotc))
        allocate(this%Dst(2, this%ntotc))
        call allocate_real_data(this%Dsf, this%ntotc+1) ! for all reacttions
        if (this%needAngular) then
           allocate(this%Dsfx(2, this%ntotc, this%ntotc))
           allocate(this%Dstx(2, this%ntotc, this%ntotc))
        end if
    end if

    if (doShiftRes) then
       call allocate_real_data(this%Xx, pars%getNumResonances())
       call allocate_real_data(this%XxHelper, pars%getNumResonances())
    end if

    if (itot.gt.0) then
       call reallocate_real_data_2d(this%Br, this%ntriag, 0, itot, 0)
       call reallocate_real_data_2d(this%Bi, this%ntriag, 0, itot, 0)
       call reallocate_real_data_2d(this%Pr, this%ntriag, 0, itot, 0)
       call reallocate_real_data_2d(this%Pi, this%ntriag, 0, itot, 0)
       call reallocate_real_data_2d(this%Qr, this%ntriag, 0, this%ntriag, 0)
       call reallocate_real_data_2d(this%Qi, this%ntriag, 0, this%ntriag, 0)
       call reallocate_real_data_2d(this%Tr, this%ntotc+1, 0, this%ntriag, 0)
       call reallocate_real_data_2d(this%Ti, this%ntotc+1, 0, this%ntriag, 0)
       call allocate_real_data(this%Ddddd, this%ntotc+1)  ! derivative calculated for all entry channels
       call allocate_real_data(this%Ddddtl, this%ntotc+1)
       if (this%needAngular) then
          allocate(this%Tx(2, this%ntriag, this%ntriag))
       end if
    end if
    allocate(this%useChannel(this%ntotc+1,2))
    allocate(this%crossSelfWhy(this%covariance%getNumTotalParam()+1))


    call reallocate_real_data_2d(this%Rmat, 2, 0, this%ntriag, 0)
    call reallocate_real_data_2d(this%Ymat, 2, 0, this%ntriag, 0)
    call reallocate_real_data_2d(this%Yinv, 2, 0, this%ntriag, 0)

     call allocate_real_data(this%Xxxxr, this%ntriag)
     call allocate_real_data(this%Xxxxi, this%ntriag)
     call reallocate_real_data_2d(this%Xqr, this%ntotc, 0, this%ntotc, 0)
     call reallocate_real_data_2d(this%Xqi, this%ntotc, 0, this%ntotc, 0)


     call reallocate_real_data_2d(this%Bound, this%ntotc, 0, ngroup, 0)

     if (this%needAngular) then
        call reallocate_real_data_2d(this%Cscs, 2, 0, this%ntriag, 0)
      end if

     call allocate_real_data(this%Sinsqr, this%ntotc)
     call allocate_real_data(this%Sin2ph, this%ntotc)
     call allocate_real_data(this%Dphi, this%ntotc)
     call allocate_real_data(this%Sinphi, this%ntotc)
     call allocate_real_data(this%Cosphi, this%ntotc)
     call allocate_real_data(this%Rootp, this%ntotc)
     call allocate_real_data(this%Elinvr, this%ntotc)
     call allocate_real_data(this%Elinvi, this%ntotc)
     call allocate_real_data(this%Psmall, this%ntotc)


    if (this%Itzero.ne.0 .or. this%Ilzero.ne.0 .OR. irad.gt.0) then
       call allocate_real_data(this%Dpdr, this%ntotc)
       call allocate_real_data(this%Dsdr, this%ntotc)
       call reallocate_real_data_2d(this%Pxrhor, this%ntriag, 0, this%ntotc, 0)
       call reallocate_real_data_2d(this%Pxrhoi, this%ntriag, 0, this%ntotc, 0)
    end if

    if (irad.gt.0) then
       if( nres.gt.0) allocate(this%Bga(this%ntriag, this%ntotc, nres))

       ! collect the unique fit flags for the true radius
       ! in iradIndex
       allocate(this%iradIndex(irad))
       this%iradIndex = 0
       irad = 0
       do igr = 1, ngroup
          call this%resData%getSpinGroupInfo(spinInfo,igr)
          DO ichan=1, spinInfo%getNumChannels()
             ifl = this%radiusData%getTrueFitFlag(Igr, Ichan)
             if (ifl.eq.0) cycle
             ! check whether we need to add it
             if (.not.any(this%iradIndex.eq.ifl)) then
                irad = irad + 1
                this%iradIndex(irad) = ifl
             end if
          end do
       end do
       ! the new irad might be smaller than the first irad in case some of the true radii are linked
       if (irad.ne.0) then ! none of the true radii is fitted
          allocate(this%Pgar(this%ntriag, irad, ngroup))
          allocate(this%Pgai(this%ntriag, irad, ngroup))
       end if
    end if

    allocate(this%crossInternal(ntot, this%resData%getNumSpinGroups(), 0:this%covariance%getNumTotalParam()))

    if (this%needAngular) then
       allocate(this%angInternal(2, ntot, ntot, this%resData%getNumSpinGroups(), 0:this%covariance%getNumTotalParam()))
       allocate(this%uniqueEchan(this%numIso))
    end if


    ! direct capture arrays. Todo update the direct capture to read and store data in C++ class
    this%numdrc = 0
    this%Nucdrc = 0
    this%hasDirectCapture = .false.
    ! end direct capture data

end subroutine
subroutine XctCrossCalc_destroy(this)
    class(XctCrossCalc) :: this

    call CrossSectionCalculator_destroy(this)
    deallocate(this%Ifcros)
    deallocate(this%termf)
    deallocate(this%termfx)
    deallocate(this%crossInternal)
    if (allocated(this%angInternal)) deallocate(this%angInternal)
    if (allocated(this%uniqueEchan)) deallocate(this%uniqueEchan)
    deallocate(this%Bound)

    if (allocated(this%Alj))      deallocate(this%Alj)
    if (allocated(this%Xx))       deallocate(this%Xx)
    if (allocated(this%XxHelper)) deallocate(this%XxHelper)
    if (allocated(this%Xlmn))     deallocate(this%Xlmn)
    if (allocated(this%Ixlmn))    deallocate(this%Ixlmn)
    if (allocated(this%Br))       deallocate(this%Br)
    if (allocated(this%Bi))       deallocate(this%Bi)
    if (allocated(this%Pr))       deallocate(this%Pr)
    if (allocated(this%Pi))       deallocate(this%Pi)
    if (allocated(this%Bga))      deallocate(this%Bga)
    if(allocated(this%Pgar))      deallocate(this%Pgar)
    if(allocated(this%Pgai))      deallocate(this%Pgai)
    if (allocated(this%Alphai))   deallocate(this%Alphai)
    if (allocated(this%Alphar))   deallocate(this%Alphar)
    if (allocated(this%needAlphai))     deallocate(this%needAlphai)
    if (allocated(this%Difen))    deallocate(this%Difen)
    if (allocated(this%iradIndex)) deallocate(this%iradIndex)
    if (allocated(this%Cscs)) deallocate(this%Cscs)
    deallocate(this%Sinsqr)
    deallocate(this%Sin2ph)
    deallocate(this%Dphi)
    deallocate(this%Sinphi)
    deallocate(this%Cosphi)
    if (allocated(this%Dpdr)) deallocate(this%Dpdr)
    if (allocated(this%Dsdr)) deallocate(this%Dsdr)
    if( allocated(this%Ddddd)) deallocate(this%Ddddd)
    if( allocated(this%Ddddtl)) deallocate(this%Ddddtl)
    if( allocated(this%useChannel)) deallocate(this%useChannel)
    if(allocated(this%Prer)) deallocate(this%Prer)
    if(allocated(this%Prei)) deallocate(this%Prei)
    if (allocated(this%Dsf)) deallocate(this%Dsf)
    if (allocated(this%Dsfx)) deallocate(this%Dsfx)
    if (allocated(this%Dsf)) deallocate(this%Dsf)
    if (allocated(this%Dstx)) deallocate(this%Dstx)
    if (allocated(this%Dstt)) deallocate(this%Dstt)
     if (allocated(this%Dst)) deallocate(this%Dst)
    deallocate(this%Rmat)
    deallocate(this%Ymat)
    deallocate(this%Yinv)
    deallocate(this%Rootp)
    deallocate(this%Elinvr)
    deallocate(this%Elinvi)
    deallocate(this%Psmall)
    deallocate(this%Xxxxr)
    deallocate(this%Xxxxi)
    deallocate(this%Xqr)
    deallocate(this%Xqi)
    deallocate(this%Xden)
    if(allocated(this%Qr))deallocate(this%Qr)
    if(allocated(this%Qi)) deallocate(this%Qi)
    if(allocated(this%Tr)) deallocate(this%Tr)
    if(allocated(this%Ti)) deallocate(this%Ti)
    if(allocated(this%Tx)) deallocate(this%Tx)
    if(allocated(this%Pxrhor)) deallocate(this%Pxrhor)
    if(allocated(this%Pxrhoi)) deallocate(this%Pxrhoi)
    if(allocated(this%Ccoulx)) deallocate(this%Ccoulx)

    if (allocated(this%crossSelfWhy)) deallocate(this%crossSelfWhy)

    ! direct capture arrays. Todo update the direct capture to read and store data in C++ class
    if(allocated(this%Edrcpt)) deallocate(this%Edrcpt)
    if(allocated(this%Cdrcpt)) deallocate(this%Cdrcpt)
    if(allocated(this%Xdrcpt)) deallocate(this%Xdrcpt)
    if(allocated(this%Ndrcpt)) deallocate(this%Ndrcpt)
    if(allocated(this%DrcpCoef)) deallocate(this%DrcpCoef)
    if(allocated(this%IflDrcpCoef)) deallocate(this%IflDrcpCoef)
    ! end direct capture data
end subroutine

end module XctCrossCalc_M
