MODULE nonstatistical_mod
USE gem_kind
USE angle
USE functs
USE trig
IMPLICIT NONE
PRIVATE


PUBLIC :: non_statistical 
PRIVATE :: prepare_nonstatistical

REAL (KIND=r4), DIMENSION(68), SAVE :: BE6_1 = &
(/2.0270086E-09, 5.8145673E-07, 7.1854870E-06, 3.2172113E-05,  &
  8.9211062E-05, 1.8850228E-04, 3.3525753E-04, 5.3003017E-04,  &
  7.6970848E-04, 1.0486309E-03, 1.3595871E-03, 1.6944309E-03,  &
  2.0449837E-03, 2.4032623E-03, 2.7617388E-03, 3.1135539E-03,  &
  3.4528391E-03, 3.7744523E-03, 4.0741344E-03, 4.3484448E-03,  &
  4.5947176E-03, 4.8111109E-03, 4.9964017E-03, 5.1498408E-03,  &
  5.2712206E-03, 5.3608082E-03, 5.4191947E-03, 5.4473905E-03,  &
  5.4465602E-03, 5.4181116E-03, 5.3636283E-03, 5.2847988E-03,  &
  5.1833903E-03, 5.0613396E-03, 4.9207220E-03, 4.7632651E-03,  &
  4.5908946E-03, 4.4055441E-03, 4.2090644E-03, 4.0033292E-03,  &
  3.7900696E-03, 3.5710568E-03, 3.3479610E-03, 3.1227483E-03,  &
  2.8964968E-03, 2.6708802E-03, 2.4473441E-03, 2.2272691E-03,  &
  2.0119594E-03, 1.8026819E-03, 1.6005816E-03, 1.4068034E-03,  &
  1.2223840E-03, 1.0483334E-03, 8.8556873E-04, 7.3493982E-04,  &
  5.9724285E-04, 4.7317994E-04, 3.6334124E-04, 2.6820163E-04,  &
  1.8810487E-04, 1.2316174E-04, 7.3201008E-05, 3.7638565E-05,  &
  1.5269745E-05, 3.9710330E-06, 3.5455668E-07, 2.4267577E-10/)
REAL (KIND=r4), DIMENSION(41), SAVE :: he6_1 = &
(/5.9025350E-04, 1.4816680E-03, 2.4142121E-03, 3.2942002E-03,  &
  4.0769968E-03, 4.7419108E-03, 5.2821259E-03, 5.6987219E-03,  &
  5.9972950E-03, 6.1864676E-03, 6.2764324E-03, 6.2776804E-03,  &
  6.2013729E-03, 6.0580783E-03, 5.8579575E-03, 5.6113214E-03,  &
  5.3261835E-03, 5.0112591E-03, 4.6740063E-03, 4.3213200E-03,  &
  3.9594155E-03, 3.5939862E-03, 3.2301033E-03, 2.8722088E-03,  &
  2.5243657E-03, 2.1901692E-03, 1.8727884E-03, 1.5749880E-03,  &
  1.2991506E-03, 1.0474186E-03, 8.2146755E-04, 6.2264630E-04,  &
  4.5202073E-04, 3.1013964E-04, 1.9707170E-04, 1.1218890E-04,  &
  5.3955449E-05, 1.9486508E-05, 3.9684455E-06, 1.6708732E-07,  &
  8.2843659E-15/)


REAL (KIND=r4), DIMENSION(85), SAVE :: be9_2 = &
(/8.2305390E-26, 2.8723518E-18, 6.4047232E-15, 6.3987297E-13, &
  1.4879972E-11, 1.5212331E-10, 9.2700458E-10, 3.9760013E-09, &
  1.3261308E-08, 3.6661554E-08, 8.7709815E-08, 1.8716629E-07, &
  3.6419763E-07, 6.5686294E-07, 1.1119278E-06, 1.7841808E-06, &
  2.7345291E-06, 4.0288242E-06, 5.7347315E-06, 7.9200827E-06, &
  1.0650383E-05, 1.3986260E-05, 1.7982296E-05, 2.2684633E-05, &
  2.8128959E-05, 3.4341094E-05, 4.1336891E-05, 4.9120306E-05, &
  5.7680518E-05, 6.7000161E-05, 7.7047647E-05, 8.7784290E-05, &
  9.9158489E-05, 1.1111212E-04, 1.2357711E-04, 1.3648247E-04, &
  1.4974950E-04, 1.6329350E-04, 1.7702780E-04, 1.9086317E-04, &
  2.0470416E-04, 2.1845852E-04, 2.3203198E-04, 2.4532984E-04, &
  2.5825755E-04, 2.7072892E-04, 2.8265314E-04, 2.9394246E-04, &
  3.0451489E-04, 3.1429125E-04, 3.2319827E-04, 3.3116009E-04, &
  3.3811448E-04, 3.4400052E-04, 3.4876211E-04, 3.5234989E-04, &
  3.5471984E-04, 3.5583484E-04, 3.5566138E-04, 3.5417802E-04, &
  3.5136752E-04, 3.4721725E-04, 3.4172178E-04, 3.3488654E-04, &
  3.2672478E-04, 3.1725483E-04, 3.0650699E-04, 2.9451671E-04, &
  2.8132813E-04, 2.6699918E-04, 2.5159851E-04, 2.3520501E-04, &
  2.1790304E-04, 1.9980530E-04, 1.8102973E-04, 1.6171398E-04, &
  1.4202381E-04, 1.2212848E-04, 1.0225191E-04, 8.2647246E-05, &
  6.3622749E-05, 4.5560981E-05, 2.8965014E-05, 1.4560422E-05, &
  3.6027143E-06/)


  TYPE (vell), PARAMETER :: blank_vell = vell(0.0,0.0,0.0)
  TYPE (spinn), PARAMETER :: blank_spinn = spinn(0.0,0.0)
  TYPE (identity), PARAMETER :: neutron = &
  identity(0,1,0,0,0.0,0.0,0.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: proton = &
  identity(1,1,0,0,0.0,0.0,0.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: deuteron = &
  identity(1,2,0,0,0.0,0.0,1.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: triton = &
  identity(1,3,0,0,0.0,0.0,0.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: he3 = &
  identity(2,3,0,0,0.0,0.0,0.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: alpha = &
  identity(2,4,0,0,0.0,0.0,0.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: M1 = &
  identity(-1,11,0,0,0.0,0.0,1.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: E2 = &
  identity(-1,2,0,0,0.0,0.0,2.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: E1 = &
  identity(-1,1,0,0,0.0,0.0,1.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: he6 = &
  identity(2,6,0,0,0.0,0.0,0.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: he7 = &
  identity(2,7,0,0,0.0,0.0,1.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: li6 = &
  identity(3,6,0,0,0.0,0.0,1.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: li7 = &
  identity(3,7,0,0,0.0,0.0,1.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: li8 = &
  identity(3,8,0,0,0.0,0.0,2.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: be7 = &
  identity(4,7,0,0,0.0,0.0,1.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: be8 = &
  identity(4,8,0,1,0.0,0.0,0.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: be9 = &
  identity(4,9,0,0,0.0,0.0,1.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: be10 = &
  identity(4,10,0,0,0.0,0.0,0.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: li5 = &
  identity(3,5,0,1,0.0,0.0,1.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: he5 = &
  identity(2,5,0,1,0.0,0.0,1.5,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: n2 = &
  identity(0,2,0,1,0.0,0.0,0.0,0.0,0.0,0.0,blank_spinn,blank_vell)

  TYPE (identity), PARAMETER :: n3 = &
  identity(0,3,0,1,0.0,0.0,1.5,0.0,0.0,0.0,blank_spinn,blank_vell)

CONTAINS

SUBROUTINE prepare_nonstatistical ()
    REAL (KIND=r4) :: tot
    INTEGER :: i
    
    tot = SUM(be6_1)
    DO I=1,68
        be6_1(i) = be6_1(i)/tot
        IF (i > 2) be6_1(i) = be6_1(i) + be6_1(i-1)
    END DO

    tot = SUM(be9_2)
    DO I=1,85
        be9_2(i) = be9_2(i)/tot
        IF (i > 2) be9_2(i) = be9_2(i) + be9_2(i-1)
    END DO

    tot = SUM(he6_1)
    DO I=1,41
        He6_1(i) = He6_1(i)/tot
        IF (i > 2) He6_1(i) = He6_1(i) + He6_1(i-1)
    END DO
END SUBROUTINE prepare_nonstatistical 


SUBROUTINE non_statistical(parent,daughter,n_extra)
        
!+
! 
! Does the decay of specific know levels of light nuclei.
!
! 3-body decays, He6(1.8 Mev), Be6, Li6(5.36) treats as
! two 2-body decays following Lane + Thomas,
! Rev Modern Physics, vol 30 1958 page=257
! R matrix parameters for He5 and Li5 intermediates
! from Stammbach and Walters, nucl phys A180 (1960) page=225.
! FUNCTIONAL DESCRIPTION: 
! 
!-
TYPE (identity), INTENT(IN) :: parent
TYPE (identity), DIMENSION(3), INTENT(out) :: daughter
INTEGER, INTENT(out) :: n_extra
INTEGER :: i,j,ll
LOGICAL, save :: first = .TRUE.
REAL (KIND=r4) :: gamma_width,decay_time,vrel,v1,v2,ared,theta,phi
REAL (KIND=r4) :: e,prob,ran
REAL (KIND=r4), DIMENSION(3,3) :: v_f
REAL (KIND=r4), DIMENSION(3) :: vv_f,p_f
REAL (KIND=r4) :: e_kinetic,l_p_s

IF (first) THEN
    CALL prepare_nonstatistical ()
    first = .FALSE.
END IF

IF (parent%level <= 0) RETURN


frag_Z: SELECT CASE (parent%Z)  ! select z
CASE (0) ! dineutrons etc
    frag_A_n: SELECT CASE (parent%A)
    CASE (2) !dineutron dineutron dineutron
        daughter(1) = neutron 
        daughter(2) = neutron
        gamma_width = 1
        e_kinetic = 0.1
        n_extra = 2
        ll = 0
        l_p_s = 0.5
    CASE (3) !trineutron trineutron trineutron
        daughter(1) = neutron
        daughter(2) = n2
        gamma_width = 1
        e_kinetic = 0.1
        n_extra = 2
        ll = 1
        l_p_s = 1.5
    CASE (4) ! 4n 4n 4n 4n 4n 4n 4n 4n 4n 4n 4n
        daughter(1) = neutron
        daughter(2) = n3
        gamma_width = 1
        e_kinetic = 0.1
        n_extra = 2
        ll = 1
        l_p_s = 1.5
    END SELECT frag_a_n
!HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
CASE (1)
! 4H  4H  4H  4H  4H  4H  4H  4H  4H  4H  4H  4H  4H
    daughter(1) = neutron
    daughter(2) = triton
    gamma_width = 1
    e_kinetic = 2.9 + parent%ex
    n_extra = 2
    ll = 1
    l_p_s = 1.5

!HEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEHEH
CASE (2)  !He fragments
    frag_a_He: SELECT CASE (parent%A)

    CASE (4) ! 4He 4He 4He 4He 4He 4He 4He 4He 4He 4He 4He 
        SELECT CASE (parent%level)
        CASE (1)
            daughter(1) = proton
            daughter(2) = triton
            gamma_width = 0.5
            e_kinetic = 0.3950 + parent%ex
            n_extra = 2
            ll = 0
            l_p_s = 0.5
        CASE (2)
            n_extra = 2
            gamma_width = 0.8
            ll = 0
            l_p_s = 0.5
            e_kinetic = 0.0
            CALL RANDOM_NUMBER(ran)
            IF (RAN <= 0.2381) THEN
                e_kinetic = 0.4320 + parent%ex
            END IF
            IF (e_kinetic > 0) THEN
                daughter(1) = neutron
                daughter(2) = he3
            ELSE
                daughter(1) = proton
                daughter(2) = triton
                e_kinetic = 1.195 + parent%ex
            END IF
        CASE (3)
            n_extra = 2
            gamma_width = 2.01
            ll = 1
            l_p_s = 1.5
            e_kinetic = 0.0
            CALL RANDOM_NUMBER(ran)
            IF (RAN <= 0.3731) THEN
                e_kinetic = 1.262 + parent%ex
            END IF
            IF (e_kinetic > 0) THEN
                daughter(1) = neutron
                daughter(2) = he3
            ELSE
                daughter(1) = proton
                daughter(2) = triton
                e_kinetic = 2.025 + parent%ex
            END IF
        CASE (4)
            n_extra = 2
            gamma_width = 5.01
            ll = 1
            l_p_s = 1.5
            e_kinetic = 0.0
            CALL RANDOM_NUMBER(ran)
            IF (RAN <= 0.4741) THEN
                e_kinetic = 2.75 + parent%ex
            END IF
            IF (e_kinetic > 0) THEN
                daughter(1) = neutron
                daughter(2) = he3
            ELSE
                daughter(1) = proton
                daughter(2) = triton
                e_kinetic = 3.515 + parent%ex
            END IF

        END SELECT
    CASE (5) ! 5He 5He 5He 5He 5He 5He 5He 5He 5He 5He 5He 5He 5He 5He 5He 
        SELECT CASE (parent%level)
        CASE (1)
            daughter(1) = neutron
            daughter(2) = alpha
            gamma_width = 0.6
            e_kinetic = parent%ex + 0.89
            if (e_kinetic <= 0) WRITE (UNIT=*, FMT=*) "he5",parent%ex
            n_extra = 2
            ll = 1
            l_p_s = 1.5
        END SELECT
    CASE (6) ! 6He 6He 6He 6He 6He 6He 6He 6He 6He 6He 6He 6He 6He 6He 6He 6He 
        SELECT CASE (parent%level)
        CASE (1)
            CALL RANDOM_NUMBER(prob)
            i = 1
            DO
                IF (prob <= he6_1(i)) EXIT
                i = i + 1
                IF (i >= 41) EXIT
            END DO
            e_kinetic = REAL(i)*0.02
            daughter(1) = neutron
            daughter(2) = he5
            daughter(2)%ex = 0.8245 - e_kinetic -0.89
            gamma_width = 0.113
            ll = 1
            l_p_s = 0.5 
            n_extra = 2
        END SELECT
    CASE (7) ! 7He 7He 7He 7He 7He 7He 7He 7He 7He 7He 7He 7He 7He 7He 7He 7He 
        daughter(1) = neutron
        daughter(2) = he6
        gamma_width = 0.16
        e_kinetic = 0.445 + parent%ex
        n_extra = 2
        ll = 1
        l_p_s = 1.5
    CASE (8) ! 8He 8He 8He 8He 8He 8He 8He 8He 8He 8He 8He 8He 8He 8He 8He 8He 
        daughter(1) = neutron
        daughter(2) = he7
        gamma_width = 0.8
        e_kinetic = 1.01 + parent%ex
        n_extra = 2
        ll = 1
        l_p_s = 1.5
    END SELECT frag_A_He
!LILILILILILILILILILILILILILILILILILILILILILILILILILILI
CASE (3)  !li fragments                     
    frag_A_Li: SELECT CASE (parent%A)
    CASE (5) ! 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 5Li 
        SELECT CASE (parent%level)
        CASE (1)
            daughter(1) = proton
            daughter(2) = alpha
            gamma_width = 1.5
            e_kinetic = 1.96 + parent%ex
            n_extra = 2
            ll = 1
            l_p_s = 1.5
        END SELECT 
    CASE (6) ! 6LI 6Li 6LI 6Li 6LI 6Li 6LI 6Li 6LI 6Li 6LI 6Li 6LI 6Li 6LI 6Li 
        SELECT CASE(parent%level)
        CASE(2)
            daughter(1) = E2
            daughter(2) = li6
            gamma_width = 8.2e-6
            e_kinetic = 0.0
            n_extra = 2
            ll = 2
            l_p_s = 2
        CASE (4)
    !       the distribution of kinetic energies between the
    !       protons and the 5He is well approximated by a gaussian  
            DO
                e_kinetic = 0.8875 + 0.23872*Gauss()
                IF (e_kinetic > 0.0 .AND.&
                e_kinetic <= 1.66) EXIT
            END DO
            
            daughter(1) = proton
            daughter(2) = he5
            daughter(2)%ex = 1.667 - e_kinetic - 0.89
            ll = 1
            l_p_s = 0.5
            gamma_width = 0.54
            n_extra = 2
        CASE DEFAULT
            daughter(1) = deuteron
            daughter(2) = alpha
            IF (parent%level == 1) THEN
                gamma_width = 0.024
                e_kinetic = 2.186 +parent%ex - 1.4753
                ll = 2
                l_p_s = 3 
            ELSE IF (parent%level == 3) THEN
                gamma_width = 1.7
                e_kinetic = 4.31 + parent%ex - 1.4753
                ll = 2
                l_p_s = 2
            ELSE IF (parent%level == 5) THEN
                gamma_width = 1.5
                e_kinetic = 5.65 + parent%ex - 1.4753
                ll = 0
                l_p_s = 1
            END IF
            n_extra = 2
        END SELECT
    CASE (7) ! 7Li 7Li 7Li 7Li 7Li 7Li 7Li 7Li 7Li 7Li 7Li 7Li 7Li 7Li 7Li 
        SELECT CASE (parent%level)
        CASE (1)
            daughter(1) = M1
            daughter(2) = li7
            gamma_width = 0.00626e-6
            e_kinetic = 0.0
            n_extra = 2
            ll = 1
            l_p_s = 1
        CASE DEFAULT
            daughter(1) = triton
            daughter(2) = alpha
            gamma_width = 0.093
            e_kinetic = 4.63 - 2.4681 + parent%ex
            n_extra = 2
            ll = 3
            l_p_s = 3.5
        END SELECT
    CASE (8) ! 8Li 8Li 8Li 8Li 8Li 8Li 8Li 8Li 8Li 8Li 8Li 8Li 8Li 8Li 
        SELECT CASE (parent%level)
        CASE (1)
            daughter(1) = M1
            daughter(2) = li8
            gamma_width = 3.8e-8
            e_kinetic = 0.0
            n_extra = 2
            ll = 1
            l_p_s = 1
        CASE (2)
            daughter(1) = neutron
            daughter(2) = li7
            gamma_width = 0.033
            e_kinetic = 0.22
            n_extra = 2
            ll = 1
            l_p_s = 1.5
        END SELECT
    END SELECT frag_A_Li
!BEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBEBE
CASE (4) ! Be fragments
    frag_a_Be: SELECT CASE (parent%A)
    CASE (6) ! 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 6Be 
        SELECT CASE (parent%level)
        CASE (1) 
            CALL RANDOM_NUMBER(prob)
            i = 1
            DO
                IF (prob <= be6_1(i)) EXIT
                i = i + 1
            END DO
            e_kinetic = REAL(i)*0.02
            daughter(1) = proton
            daughter(2) = li5
            daughter(2)%ex = 1.3706 - e_kinetic - 1.96
            gamma_width = 0.096
            ll = 1
            l_p_s = 1.5 
            n_extra = 2
        CASE (2)
            daughter(1) = proton
            daughter(2) = proton
            daughter(3) = alpha
            gamma_width = 1.16
            e_kinetic = 3.041
            n_extra = 3
        END SELECT
    CASE (7) ! 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 7Be 
        SELECT CASE (parent%level)
        CASE (1)
            daughter(1) = M1
            daughter(2) = be7
            gamma_width = 2.4e-9
            e_kinetic = 0.0
            n_extra = 2
            ll = 1
            l_p_s = 1
        CASE (2)
            daughter(1) = he3
            daughter(2) = alpha
            gamma_width = 0.175
            e_kinetic = 2.9824 + parent%ex
            n_extra = 2
            ll = 3
            l_p_s = 3.5
         END SELECT
    CASE (8) ! 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 8Be 
        SELECT CASE (parent%level)
        CASE (1)
            daughter(1) = alpha
            daughter(2) = alpha
            gamma_width = 6.8e-6
            e_kinetic = 0.09178 + parent%ex
            n_extra = 2
            ll = 0
            l_p_s = 0
        CASE (2) 
            daughter(1) = alpha
            daughter(2) = alpha
            gamma_width = 1.5
            e_kinetic = 3.1318 + parent%ex
            n_extra = 2
            ll = 2
            l_p_s = 2
        END SELECT
    CASE (9) ! 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 9Be 
        SELECT CASE (parent%level)
        CASE (1) 
            daughter(1) = E1
            daughter(2) = Be9
            gamma_width = 0.15
            e_kinetic = 0.0
            n_extra = 2
            ll = 1
            l_p_s = 1
        CASE (3:6) 
            daughter(1) = neutron
            daughter(2) = be8
            SELECT CASE (parent%level)
            CASE (3)
                gamma_width = 1.08
                e_kinetic = 1.1345 + parent%ex
                ll = 1
                l_p_s = 0.5
            CASE (4)
                gamma_width = 0.282
                e_kinetic = 1.3845 + parent%ex
                ll = 2
                l_p_s = 2.5
            CASE (5)
                gamma_width = 0.743
                e_kinetic = 4.7 - 1.6655 + parent%ex
                ll = 1
                l_p_s = 1.5
            CASE (6)
                gamma_width = 1.54
                e_kinetic = 6.76 - 1.6655 + parent%ex
                ll = 3
                l_p_s = 3.5
            END SELECT
            n_extra = 2
        CASE (2)
            daughter(1) = alpha
            daughter(2) = he5
            gamma_width = 7.7e-4
            CALL RANDOM_NUMBER(prob)
            i = 1
            DO
                IF (prob <= be9_2(i)) EXIT
                i = i + 1
            END DO
            e_kinetic = REAL(i)*0.01
            daughter(2)%ex = 0.8563 - e_kinetic - 0.89
            ll = 1
            l_p_s = 1
            n_extra = 2
        END SELECT
    CASE (10) ! 10Be 10Be 10Be 10Be 10Be 10Be 10Be 10Be 10Be 10Be 10Be 10Be 
        SELECT CASE (parent%level)
        CASE (1)
            daughter(1) = E2
            daughter(2) = Be10
            gamma_width = 0.15
            e_kinetic = 0.0
            n_extra = 2
            ll = 2
            l_p_s = 2
        CASE (2)
            daughter(1) = E2
            daughter(2) = Be10
            gamma_width = 0.15
            e_kinetic = 0.0
            n_extra = 2
            ll = 2
            l_p_s = 2
        CASE (3)
            daughter(1) = E1
            daughter(2) = Be10
            gamma_width = 0.15
            e_kinetic = 0.0
            n_extra = 2
            ll = 1
            l_p_s = 1
        END SELECT
    END SELECT frag_A_Be
!??????????????????????????????????????????????????
CASE DEFAULT
    WRITE (UNIT=*, FMT=*) "decay of unstable fragment not known"
    WRITE (UNIT=*, FMT=*) "z=",parent%z
    WRITE (UNIT=*, FMT=*) "A=",parent%A
    WRITE (UNIT=*, FMT=*) "level=",parent%level
!    CALL crash ()
    STOP
END SELECT frag_z

decay_time = exp_decay(gamma_width)
DO j=1,n_extra
    daughter(j)%id = parent%id
    daughter(j)%time = parent%time + decay_time
END DO



!I have made this very simple at the moment,
!we ignore recoil from gamma"s
!otherwise decay is isotropic for classical angular
!distributions.

daughter(1)%id = parent%id
daughter(2)%id = parent%id
daughter(3)%id = parent%id
       
IF (para%i_angle) THEN
    IF (e_kinetic == 0) THEN   ! gamma ray emission
        daughter(1)%vel = parent%vel
    ELSE
        SELECT CASE (n_extra)
        CASE (2)  ! 2body decay
            IF (para%quantum) THEN
                if(e_kinetic <= 0) THEN
                   WRITE (UNIT=*, FMT=*) "ek<0 for non_statistical",parent%Z,parent%A,&
                   parent%level,parent%ex,e_kinetic
!                   CALL crash ()
                END IF
                 CALL quantum_lp(parent,daughter(1:2),&
                 ll,l_p_s,e_kinetic)
            ELSE
                CALL RANDOM_NUMBER(ran)
                theta = ACOSD(1.0-2.0*RAN)
                CALL RANDOM_NUMBER(ran)
                phi = 360.0*RAN
                ared = REAL(daughter(1)%A*daughter(2)%A)&
                /REAL(daughter(1)%A+daughter(2)%A)
                vrel = (2.0*e_kinetic/ared)**0.5*0.9794
                v1 = vrel*ared/REAL(daughter(1)%A)
                v2 = vrel - v1
                daughter(1)%vel = vell(v1,theta,phi) + parent%vel 

                theta = 180.0-theta
                phi = MODULO(phi+180.0,360.0)

                daughter(2)%vel = vell(v2,theta,phi) + parent%vel 
                daughter(1)%ex = 0.
                daughter(2)%ex = 0.
            END IF                                        
        CASE (3)  ! 3-body decay

!           zero momentum
            p_f = 0.0

!           choose initial velocities , isotropic       
            DO i=1,3
                CALL RANDOM_NUMBER(ran)
                theta = ACOSD(1.0-2.0*RAN)
                CALL RANDOM_NUMBER(ran)
                phi = 360.0*RAN
                CALL RANDOM_NUMBER(ran)
                vv_f(i) = RAN**0.33333333/REAL(daughter(i)%A)
                v_f(i,1) = vv_f(i)*SIND(theta)*COSD(phi)
                v_f(i,2) = vv_f(i)*SIND(theta)*SIND(phi)
                v_f(i,3) = vv_f(i)*COSD(theta)
!               find total momentum     
                DO j=1,3
                    p_f(j) = p_f(j) + v_f(i,j)*REAL(daughter(i)%A)
                END DO
            END DO
!           find centre of mass velocity        
            p_f = p_f/REAL(parent%A)

!           remove centre of mass motion  and find total kinetic energy
            e = 0                   
            DO i=1,3
                vv_f(i) = 0.0
                DO j=1,3
                    v_f(i,j) = v_f(i,j) - p_f(j)
                    vv_f(i) = vv_f(i) + v_f(i,j)**2.0
                END DO
                e = e + 0.5*REAL(daughter(i)%A)*vv_f(i)**2.0*1.042
                vv_f(i) = vv_f(i)**0.5
            END DO
!           find angles and scale velocity to get desired total
!           kinetic energy      
            DO i=1,3
                theta = ACOSD(v_f(i,3)/vv_f(i))
                phi = ATAN2D(v_f(i,2),v_f(i,1))
                vv_f(i) = vv_f(i)*(e_kinetic/e)**0.5

!               add to velocity of parent       
                daughter(i)%vel = vell(vv_f(i),theta,phi) + parent%vel 
                daughter(i)%ex = 0.
           END DO
        END SELECT
    END IF
END IF

RETURN 
END SUBROUTINE non_statistical
END MODULE nonstatistical_mod