#cython: wraparound=False
#cython: boundscheck=False
#cython: cdivision=True
#cython: nonecheck=False
#cython: profile=False
#cython: infer_types=False
include 'clpt_nonlinear_header.pxi'

from .clpt_commons_bc1 cimport cfwx, cfwt, cfN


cdef int NL_kinematics=0 # to use cfstrain_donnell in cfN


def calc_k0L_attempt(np.ndarray[cDOUBLE, ndim=1] coeffs,
             double alpharad, double r2, double L, double tLA,
             np.ndarray[cDOUBLE, ndim=2] F,
             int m1, int m2, int n2,
             int freqx, int freqt, int num_cores, str method,
             np.ndarray[cDOUBLE, ndim=1] c0, int m0, int n0):
    cdef double sina, cosa, xa, xb, ta, tb
    cdef int c, row, col
    cdef int i1, k1, i2, j2, k2, l2
    cdef np.ndarray[cINT, ndim=1] rows, cols
    cdef np.ndarray[cDOUBLE, ndim=1] k0Lv

    cdef double sini1x, cosi1x, cosk1x
    cdef double sini2x, cosi2x, sink2x, cosk2x, sinj2t, cosj2t, sinl2t, cosl2t

    cdef double r, x, t, wx, wt, w0x, w0t
    cdef double tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6
    cdef double tmp7, tmp8, tmp9, tmp10, tmp11
    cdef int i, j, nx, nt
    cdef double A11, A12, A16, A22, A26, A66
    cdef double B11, B12, B16, B22, B26, B66

    cdef double **weightx
    cdef double **weightt
    cdef double **xs
    cdef double **ts

    freqx = 8
    freqt = 8

    fdim = 3*m1 + 3*2*m2*n2 + 3*m1**2 + 2*6*m1*m2*n2 + 12*m2**2*n2**2

    rows = np.zeros((fdim,), dtype=INT)
    cols = np.zeros((fdim,), dtype=INT)
    k0Lv = np.zeros((fdim,), dtype=DOUBLE)

    sina = sin(alpharad)
    cosa = cos(alpharad)

    A11 = F[0, 0]
    A12 = F[0, 1]
    A16 = F[0, 2]
    A22 = F[1, 1]
    A26 = F[1, 2]
    A66 = F[2, 2]
    B11 = F[0, 3]
    B12 = F[0, 4]
    B16 = F[0, 5]
    B22 = F[1, 4]
    B26 = F[1, 5]
    B66 = F[2, 5]

    xa = 0.
    xb = L
    ta = 0.
    tb = 2*pi

    nx = max(m1-i0, m2-i0)
    nt = n2-j0

    with nogil:
        weightx = <double **>malloc(nx * sizeof(double *))
        weightt = <double **>malloc(nt * sizeof(double *))
        xs = <double **>malloc(nx * sizeof(double *))
        ts = <double **>malloc(nt * sizeof(double *))
        for i in range(nx):
            weightx[i] = <double *>malloc(freqx*(i+1) * sizeof(double))
            xs[i] = <double *>malloc(freqx*(i+1) * sizeof(double))
            trapz_wp(freqx*(i+1), xa, xb, weightx[i], xs[i])

        for j in range(nt):
            weightt[j] = <double *>malloc(freqt*(j+1) * sizeof(double))
            ts[j] = <double *>malloc(freqt*(j+1) * sizeof(double))
            trapz_wp(freqt*(j+1), ta, tb, weightt[j], ts[j])

        c = -1

        for k1 in range(i0, m1+i0):
            col = (k1-i0)*num1 + num0

            nx = k1-i0
            nt = 0

            # k0L_01

            tmp0 = 0
            tmp1 = 0
            tmp2 = 0

            for i in range(freqx*(nx+1)):
                x = xs[nx][i]
                r = r2 + sina*x

                cosk1x = cos(pi*k1*x/L)

                for j in range(freqt*(nt+1)):
                    t = ts[nt][j]

                    cfwx(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wx)
                    cfwt(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wt)
                    cfw0x(&x, &t, 1, &c0[0], L, m0, n0, &w0x, funcnum)
                    cfw0t(&x, &t, 1, &c0[0], L, m0, n0, &w0t, funcnum)

                    tmp0 += weightx[nx][i]*weightt[nt][j]*(-pi*cosk1x*k1*(A16*r*(w0t + wt) - A26*sina*(L - x)*(w0t + wt) + r*(w0x + wx)*(A11*r + A12*sina*(-L + x)))/(L**2*cosa*r))
                    tmp1 += weightx[nx][i]*weightt[nt][j]*(-pi*cosk1x*k1*r2*(r + sina*(L - x))*(A16*r*(w0x + wx) + A66*(w0t + wt))/(L**2*r))
                    tmp2 += weightx[nx][i]*weightt[nt][j]*(-pi*cosk1x*k1*((L - x)*(A16*r*(w0x + wx) + A66*(w0t + wt))*sin(t - tLA) + (A16*r*(w0t + wt) - A26*sina*(L - x)*(w0t + wt) + r*(w0x + wx)*(A11*r + A12*sina*(-L + x)))*cos(t - tLA))/(L**2*cosa*r))
            c += 1
            rows[c] = 0
            cols[c] = col+2
            k0Lv[c] = tmp0
            c += 1
            rows[c] = 1
            cols[c] = col+2
            k0Lv[c] = tmp1
            c += 1
            rows[c] = 2
            cols[c] = col+2
            k0Lv[c] = tmp2

        for k2 in range(i0, m2+i0):
            nx = k2-i0
            for l2 in range(j0, n2+j0):
                col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1

                nt = l2-j0

                # k0L_02

                tmp0 = 0
                tmp1 = 0
                tmp2 = 0
                tmp3 = 0
                tmp4 = 0
                tmp5 = 0

                for i in range(freqx*(nx+1)):
                    x = xs[nx][i]
                    r = r2 + sina*x

                    sink2x = sin(pi*k2*x/L)
                    cosk2x = cos(pi*k2*x/L)

                    for j in range(freqt*(nt+1)):
                        t = ts[nt][j]

                        sinl2t = sin(l2*t)
                        cosl2t = cos(l2*t)

                        cfwx(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wx)
                        cfwt(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wt)
                        cfw0x(&x, &t, 1, &c0[0], L, m0, n0, &w0x, funcnum)
                        cfw0t(&x, &t, 1, &c0[0], L, m0, n0, &w0t, funcnum)

                        tmp0 += weightx[nx][i]*weightt[nt][j]*(-(L*cosl2t*l2*sink2x*(w0t + wt)*(A12*r + A22*sina*(-L + x)) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(A11*r + A12*sina*(-L + x)) + r*(A16*r + A26*sina*(-L + x))*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt)))/(L**2*cosa*r**2))
                        tmp1 += weightx[nx][i]*weightt[nt][j]*((L*l2*sink2x*sinl2t*(w0t + wt)*(A12*r + A22*sina*(-L + x)) - pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(A11*r + A12*sina*(-L + x)) + r*(A16*r + A26*sina*(-L + x))*(L*l2*sink2x*sinl2t*(w0x + wx) - pi*cosk2x*cosl2t*k2*(w0t + wt)))/(L**2*cosa*r**2))
                        tmp2 += weightx[nx][i]*weightt[nt][j]*(-r2*(r + sina*(L - x))*(L*cosl2t*l2*sink2x*(A26*(w0t + wt) + A66*r*(w0x + wx)) + pi*cosk2x*k2*r*sinl2t*(A16*r*(w0x + wx) + A66*(w0t + wt)))/(L**2*r**2))
                        tmp3 += weightx[nx][i]*weightt[nt][j]*(r2*(r + sina*(L - x))*(L*l2*sink2x*sinl2t*(A26*(w0t + wt) + A66*r*(w0x + wx)) - pi*cosk2x*cosl2t*k2*r*(A16*r*(w0x + wx) + A66*(w0t + wt)))/(L**2*r**2))
                        tmp4 += weightx[nx][i]*weightt[nt][j]*(-(L*cosl2t*l2*sink2x*(w0t + wt)*(A26*(L - x)*sin(t - tLA) + (A12*r + A22*sina*(-L + x))*cos(t - tLA)) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(A16*(L - x)*sin(t - tLA) + (A11*r + A12*sina*(-L + x))*cos(t - tLA)) + r*(A66*(L - x)*sin(t - tLA) + (A16*r + A26*sina*(-L + x))*cos(t - tLA))*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt)))/(L**2*cosa*r**2))
                        tmp5 += weightx[nx][i]*weightt[nt][j]*((L*l2*sink2x*sinl2t*(w0t + wt)*(A26*(L - x)*sin(t - tLA) + (A12*r + A22*sina*(-L + x))*cos(t - tLA)) - pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(A16*(L - x)*sin(t - tLA) + (A11*r + A12*sina*(-L + x))*cos(t - tLA)) + r*(A66*(L - x)*sin(t - tLA) + (A16*r + A26*sina*(-L + x))*cos(t - tLA))*(L*l2*sink2x*sinl2t*(w0x + wx) - pi*cosk2x*cosl2t*k2*(w0t + wt)))/(L**2*cosa*r**2))

                c += 1
                rows[c] = 0
                cols[c] = col+4
                k0Lv[c] = tmp0
                c += 1
                rows[c] = 0
                cols[c] = col+5
                k0Lv[c] = tmp1
                c += 1
                rows[c] = 1
                cols[c] = col+4
                k0Lv[c] = tmp2
                c += 1
                rows[c] = 1
                cols[c] = col+5
                k0Lv[c] = tmp3
                c += 1
                rows[c] = 2
                cols[c] = col+4
                k0Lv[c] = tmp4
                c += 1
                rows[c] = 2
                cols[c] = col+5
                k0Lv[c] = tmp5

        for i1 in range(i0, m1+i0):
            row = (i1-i0)*num1 + num0
            for k1 in range(i0, m1+i0):
                col = (k1-i0)*num1 + num0

                if k1 > i1:
                    nx = k1-i0
                else:
                    nx = i1-i0
                nt = 0

                # k0L_11

                tmp0 = 0
                tmp1 = 0
                tmp2 = 0

                for i in range(freqx*(nx+1)):
                    x = xs[nx][i]
                    r = r2 + sina*x

                    sini1x = sin(pi*i1*x/L)
                    cosi1x = cos(pi*i1*x/L)
                    cosk1x = cos(pi*k1*x/L)

                    for j in range(freqt*(nt+1)):
                        t = ts[nt][j]

                        cfwx(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wx)
                        cfwt(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wt)
                        cfw0x(&x, &t, 1, &c0[0], L, m0, n0, &w0x, funcnum)
                        cfw0t(&x, &t, 1, &c0[0], L, m0, n0, &w0t, funcnum)

                        tmp0 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(L*sina*sini1x*(A12*r*(w0x + wx) + A26*(w0t + wt)) + pi*cosi1x*i1*r*(A11*r*(w0x + wx) + A16*(w0t + wt)))/(L**2*r))
                        tmp1 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(A16*r*(w0x + wx) + A66*(w0t + wt))*(-L*sina*sini1x + pi*cosi1x*i1*r)/(L**2*r))
                        tmp2 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(-pi*L*cosi1x*i1*sina*(B12*r*(w0x + wx) + B26*(w0t + wt)) + sini1x*(A26*L**2*cosa*(w0t + wt) + r*(pi**2*B16*i1**2*(w0t + wt) + (w0x + wx)*(A12*L**2*cosa + pi**2*B11*i1**2*r))))/(L**3*r))

                c += 1
                rows[c] = row+0
                cols[c] = col+2
                k0Lv[c] = tmp0
                c += 1
                rows[c] = row+1
                cols[c] = col+2
                k0Lv[c] = tmp1
                c += 1
                rows[c] = row+2
                cols[c] = col+2
                k0Lv[c] = tmp2

            for k2 in range(i0, m2+i0):
                for l2 in range(j0, n2+j0):
                    col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1

                    if k2 > i1:
                        nx = k2-i0
                    else:
                        nx = i1-i0
                    nt = l2-j0

                    # k0L_12

                    tmp0 = 0
                    tmp1 = 0
                    tmp2 = 0
                    tmp3 = 0
                    tmp4 = 0
                    tmp5 = 0

                    for i in range(freqx*(nx+1)):
                        x = xs[nx][i]
                        r = r2 + sina*x

                        sini1x = sin(pi*i1*x/L)
                        cosi1x = cos(pi*i1*x/L)
                        sink2x = sin(pi*k2*x/L)
                        cosk2x = cos(pi*k2*x/L)

                        for j in range(freqt*(nt+1)):
                            t = ts[nt][j]

                            sinl2t = sin(l2*t)
                            cosl2t = cos(l2*t)

                            cfwx(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wx)
                            cfwt(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wt)
                            cfw0x(&x, &t, 1, &c0[0], L, m0, n0, &w0x, funcnum)
                            cfw0t(&x, &t, 1, &c0[0], L, m0, n0, &w0t, funcnum)

                            tmp0 += weightx[nx][i]*weightt[nt][j]*((L*cosl2t*l2*sink2x*(w0t + wt)*(pi*A12*cosi1x*i1*r + A22*L*sina*sini1x) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(pi*A11*cosi1x*i1*r + A12*L*sina*sini1x) + r*(pi*A16*cosi1x*i1*r + A26*L*sina*sini1x)*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt)))/(L**2*r**2))
                            tmp1 += weightx[nx][i]*weightt[nt][j]*((-L*l2*sink2x*sinl2t*(w0t + wt)*(pi*A12*cosi1x*i1*r + A22*L*sina*sini1x) + pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(pi*A11*cosi1x*i1*r + A12*L*sina*sini1x) + r*(pi*A16*cosi1x*i1*r + A26*L*sina*sini1x)*(-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt)))/(L**2*r**2))
                            tmp2 += weightx[nx][i]*weightt[nt][j]*((-L*sina*sini1x + pi*cosi1x*i1*r)*(L*cosl2t*l2*sink2x*(A26*(w0t + wt) + A66*r*(w0x + wx)) + pi*cosk2x*k2*r*sinl2t*(A16*r*(w0x + wx) + A66*(w0t + wt)))/(L**2*r**2))
                            tmp3 += weightx[nx][i]*weightt[nt][j]*((-L*sina*sini1x + pi*cosi1x*i1*r)*(-L*l2*sink2x*sinl2t*(A26*(w0t + wt) + A66*r*(w0x + wx)) + pi*cosk2x*cosl2t*k2*r*(A16*r*(w0x + wx) + A66*(w0t + wt)))/(L**2*r**2))
                            tmp4 += weightx[nx][i]*weightt[nt][j]*((L*cosl2t*l2*sink2x*(w0t + wt)*(-pi*B22*L*cosi1x*i1*sina + sini1x*(A22*L**2*cosa + pi**2*B12*i1**2*r)) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(-pi*B12*L*cosi1x*i1*sina + sini1x*(A12*L**2*cosa + pi**2*B11*i1**2*r)) + r*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt))*(-pi*B26*L*cosi1x*i1*sina + sini1x*(A26*L**2*cosa + pi**2*B16*i1**2*r)))/(L**3*r**2))
                            tmp5 += weightx[nx][i]*weightt[nt][j]*((-L*l2*sink2x*sinl2t*(w0t + wt)*(-pi*B22*L*cosi1x*i1*sina + sini1x*(A22*L**2*cosa + pi**2*B12*i1**2*r)) + pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(-pi*B12*L*cosi1x*i1*sina + sini1x*(A12*L**2*cosa + pi**2*B11*i1**2*r)) + r*(-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt))*(-pi*B26*L*cosi1x*i1*sina + sini1x*(A26*L**2*cosa + pi**2*B16*i1**2*r)))/(L**3*r**2))

                    c += 1
                    rows[c] = row+0
                    cols[c] = col+4
                    k0Lv[c] = tmp0
                    c += 1
                    rows[c] = row+0
                    cols[c] = col+5
                    k0Lv[c] = tmp1
                    c += 1
                    rows[c] = row+1
                    cols[c] = col+4
                    k0Lv[c] = tmp2
                    c += 1
                    rows[c] = row+1
                    cols[c] = col+5
                    k0Lv[c] = tmp3
                    c += 1
                    rows[c] = row+2
                    cols[c] = col+4
                    k0Lv[c] = tmp4
                    c += 1
                    rows[c] = row+2
                    cols[c] = col+5
                    k0Lv[c] = tmp5

        for i2 in range(i0, m2+i0):
            for j2 in range(j0, n2+j0):
                row = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1
                for k1 in range(i0, m1+i0):
                    col = (k1-i0)*num1 + num0

                    if i2 > k1:
                        nx = i2-i0
                    else:
                        nx = k1-i0
                    nt = j2-j0

                    # k0L_21

                    tmp0 = 0
                    tmp1 = 0
                    tmp2 = 0
                    tmp3 = 0
                    tmp4 = 0
                    tmp5 = 0

                    for i in range(freqx*(nx+1)):
                        x = xs[nx][i]
                        r = r2 + sina*x

                        sini2x = sin(pi*i2*x/L)
                        cosi2x = cos(pi*i2*x/L)
                        cosk1x = cos(pi*k1*x/L)

                        for j in range(freqt*(nt+1)):
                            t = ts[nt][j]

                            sinj2t = sin(j2*t)
                            cosj2t = cos(j2*t)

                            cfwx(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wx)
                            cfwt(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wt)
                            cfw0x(&x, &t, 1, &c0[0], L, m0, n0, &w0x, funcnum)
                            cfw0t(&x, &t, 1, &c0[0], L, m0, n0, &w0t, funcnum)

                            tmp0 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(L*sini2x*(cosj2t*j2*(A16*r*(w0x + wx) + A66*(w0t + wt)) + sina*sinj2t*(A12*r*(w0x + wx) + A26*(w0t + wt))) + pi*cosi2x*i2*r*sinj2t*(A11*r*(w0x + wx) + A16*(w0t + wt)))/(L**2*r))
                            tmp1 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(-L*j2*sini2x*sinj2t*(A16*r*(w0x + wx) + A66*(w0t + wt)) + cosj2t*(L*sina*sini2x*(A12*r*(w0x + wx) + A26*(w0t + wt)) + pi*cosi2x*i2*r*(A11*r*(w0x + wx) + A16*(w0t + wt))))/(L**2*r))
                            tmp2 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(L*sini2x*(cosj2t*j2*(A12*r*(w0x + wx) + A26*(w0t + wt)) - sina*sinj2t*(A16*r*(w0x + wx) + A66*(w0t + wt))) + pi*cosi2x*i2*r*sinj2t*(A16*r*(w0x + wx) + A66*(w0t + wt)))/(L**2*r))
                            tmp3 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(-L*j2*sini2x*sinj2t*(A12*r*(w0x + wx) + A26*(w0t + wt)) + cosj2t*(A16*r*(w0x + wx) + A66*(w0t + wt))*(-L*sina*sini2x + pi*cosi2x*i2*r))/(L**2*r))
                            tmp4 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(L*cosj2t*j2*(B16*r*(w0x + wx) + B66*(w0t + wt))*(L*sina*sini2x - 2*pi*cosi2x*i2*r) + sinj2t*(-pi*L*cosi2x*i2*r*sina*(B12*r*(w0x + wx) + B26*(w0t + wt)) + sini2x*(B26*L**2*j2**2*(w0t + wt) + r*(A26*L**2*cosa*(w0t + wt) + pi**2*B16*i2**2*r*(w0t + wt) + (w0x + wx)*(A12*L**2*cosa*r + pi**2*B11*i2**2*r**2 + B12*L**2*j2**2)))))/(L**3*r**2))
                            tmp5 += weightx[nx][i]*weightt[nt][j]*(pi*cosk1x*k1*(-L*j2*sinj2t*(B16*r*(w0x + wx) + B66*(w0t + wt))*(L*sina*sini2x - 2*pi*cosi2x*i2*r) + cosj2t*(-pi*L*cosi2x*i2*r*sina*(B12*r*(w0x + wx) + B26*(w0t + wt)) + sini2x*(B26*L**2*j2**2*(w0t + wt) + r*(A26*L**2*cosa*(w0t + wt) + pi**2*B16*i2**2*r*(w0t + wt) + (w0x + wx)*(A12*L**2*cosa*r + pi**2*B11*i2**2*r**2 + B12*L**2*j2**2)))))/(L**3*r**2))

                    c += 1
                    rows[c] = row+0
                    cols[c] = col+2
                    k0Lv[c] = tmp0
                    c += 1
                    rows[c] = row+1
                    cols[c] = col+2
                    k0Lv[c] = tmp1
                    c += 1
                    rows[c] = row+2
                    cols[c] = col+2
                    k0Lv[c] = tmp2
                    c += 1
                    rows[c] = row+3
                    cols[c] = col+2
                    k0Lv[c] = tmp3
                    c += 1
                    rows[c] = row+4
                    cols[c] = col+2
                    k0Lv[c] = tmp4
                    c += 1
                    rows[c] = row+5
                    cols[c] = col+2
                    k0Lv[c] = tmp5

                for k2 in range(i0, m2+i0):
                    for l2 in range(j0, n2+j0):
                        col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1

                        if i2 > k2:
                            nx = i2-i0
                        else:
                            nx = k2-i0

                        if j2 > l2:
                            nt = j2-j0
                        else:
                            nt = l2-j0

                        # k0L_22

                        tmp0 = 0
                        tmp1 = 0
                        tmp2 = 0
                        tmp3 = 0
                        tmp4 = 0
                        tmp5 = 0
                        tmp6 = 0
                        tmp7 = 0
                        tmp8 = 0
                        tmp9 = 0
                        tmp10 = 0
                        tmp11 = 0

                        for i in range(freqx*(nx+1)):
                            x = xs[nx][i]
                            r = r2 + sina*x

                            sini2x = sin(pi*i2*x/L)
                            cosi2x = cos(pi*i2*x/L)
                            sink2x = sin(pi*k2*x/L)
                            cosk2x = cos(pi*k2*x/L)

                            for j in range(freqt*(nt+1)):
                                t = ts[nt][j]

                                sinj2t = sin(j2*t)
                                cosj2t = cos(j2*t)
                                sinl2t = sin(l2*t)
                                cosl2t = cos(l2*t)

                                cfwx(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wx)
                                cfwt(&coeffs[0], m1, m2, n2, &x, &t, 1, L, &wt)
                                cfw0x(&x, &t, 1, &c0[0], L, m0, n0, &w0x, funcnum)
                                cfw0t(&x, &t, 1, &c0[0], L, m0, n0, &w0t, funcnum)

                                tmp0 += weightx[nx][i]*weightt[nt][j]*((L*cosl2t*l2*sink2x*(w0t + wt)*(pi*A12*cosi2x*i2*r*sinj2t + L*sini2x*(A22*sina*sinj2t + A26*cosj2t*j2)) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(pi*A11*cosi2x*i2*r*sinj2t + L*sini2x*(A12*sina*sinj2t + A16*cosj2t*j2)) + r*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt))*(pi*A16*cosi2x*i2*r*sinj2t + L*sini2x*(A26*sina*sinj2t + A66*cosj2t*j2)))/(L**2*r**2))
                                tmp1 += weightx[nx][i]*weightt[nt][j]*((-L*l2*sink2x*sinl2t*(w0t + wt)*(pi*A12*cosi2x*i2*r*sinj2t + L*sini2x*(A22*sina*sinj2t + A26*cosj2t*j2)) + pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(pi*A11*cosi2x*i2*r*sinj2t + L*sini2x*(A12*sina*sinj2t + A16*cosj2t*j2)) + r*(-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt))*(pi*A16*cosi2x*i2*r*sinj2t + L*sini2x*(A26*sina*sinj2t + A66*cosj2t*j2)))/(L**2*r**2))
                                tmp2 += weightx[nx][i]*weightt[nt][j]*((L*cosl2t*l2*sink2x*(w0t + wt)*(-A26*L*j2*sini2x*sinj2t + cosj2t*(pi*A12*cosi2x*i2*r + A22*L*sina*sini2x)) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(-A16*L*j2*sini2x*sinj2t + cosj2t*(pi*A11*cosi2x*i2*r + A12*L*sina*sini2x)) + r*(-A66*L*j2*sini2x*sinj2t + cosj2t*(pi*A16*cosi2x*i2*r + A26*L*sina*sini2x))*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt)))/(L**2*r**2))
                                tmp3 += weightx[nx][i]*weightt[nt][j]*((-L*l2*sink2x*sinl2t*(w0t + wt)*(-A26*L*j2*sini2x*sinj2t + cosj2t*(pi*A12*cosi2x*i2*r + A22*L*sina*sini2x)) + pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(-A16*L*j2*sini2x*sinj2t + cosj2t*(pi*A11*cosi2x*i2*r + A12*L*sina*sini2x)) + r*(-A66*L*j2*sini2x*sinj2t + cosj2t*(pi*A16*cosi2x*i2*r + A26*L*sina*sini2x))*(-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt)))/(L**2*r**2))
                                tmp4 += weightx[nx][i]*weightt[nt][j]*((L*cosl2t*l2*sink2x*(w0t + wt)*(pi*A26*cosi2x*i2*r*sinj2t + L*sini2x*(A22*cosj2t*j2 - A26*sina*sinj2t)) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(pi*A16*cosi2x*i2*r*sinj2t + L*sini2x*(A12*cosj2t*j2 - A16*sina*sinj2t)) + r*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt))*(pi*A66*cosi2x*i2*r*sinj2t + L*sini2x*(A26*cosj2t*j2 - A66*sina*sinj2t)))/(L**2*r**2))
                                tmp5 += weightx[nx][i]*weightt[nt][j]*((-L*l2*sink2x*sinl2t*(w0t + wt)*(pi*A26*cosi2x*i2*r*sinj2t + L*sini2x*(A22*cosj2t*j2 - A26*sina*sinj2t)) + pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(pi*A16*cosi2x*i2*r*sinj2t + L*sini2x*(A12*cosj2t*j2 - A16*sina*sinj2t)) + r*(-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt))*(pi*A66*cosi2x*i2*r*sinj2t + L*sini2x*(A26*cosj2t*j2 - A66*sina*sinj2t)))/(L**2*r**2))
                                tmp6 += weightx[nx][i]*weightt[nt][j]*((L*cosl2t*l2*sink2x*(w0t + wt)*(-A22*L*j2*sini2x*sinj2t + A26*cosj2t*(-L*sina*sini2x + pi*cosi2x*i2*r)) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(-A12*L*j2*sini2x*sinj2t + A16*cosj2t*(-L*sina*sini2x + pi*cosi2x*i2*r)) + r*(-A26*L*j2*sini2x*sinj2t + A66*cosj2t*(-L*sina*sini2x + pi*cosi2x*i2*r))*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt)))/(L**2*r**2))
                                tmp7 += weightx[nx][i]*weightt[nt][j]*((L*l2*sink2x*sinl2t*(w0t + wt)*(A22*L*j2*sini2x*sinj2t + cosj2t*(A26*L*sina*sini2x - pi*A26*cosi2x*i2*r)) + pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(-A12*L*j2*sini2x*sinj2t + A16*cosj2t*(-L*sina*sini2x + pi*cosi2x*i2*r)) + r*(A26*L*j2*sini2x*sinj2t + cosj2t*(A66*L*sina*sini2x - pi*A66*cosi2x*i2*r))*(L*l2*sink2x*sinl2t*(w0x + wx) - pi*cosk2x*cosl2t*k2*(w0t + wt)))/(L**2*r**2))
                                tmp8 += weightx[nx][i]*weightt[nt][j]*((L*cosl2t*l2*sink2x*(w0t + wt)*(B26*L*cosj2t*j2*(L*sina*sini2x - 2*pi*cosi2x*i2*r) + sinj2t*(-pi*B22*L*cosi2x*i2*r*sina + sini2x*(A22*L**2*cosa*r + pi**2*B12*i2**2*r**2 + B22*L**2*j2**2))) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(B16*L*cosj2t*j2*(L*sina*sini2x - 2*pi*cosi2x*i2*r) + sinj2t*(-pi*B12*L*cosi2x*i2*r*sina + sini2x*(A12*L**2*cosa*r + pi**2*B11*i2**2*r**2 + B12*L**2*j2**2))) + r*(B66*L*cosj2t*j2*(L*sina*sini2x - 2*pi*cosi2x*i2*r) + sinj2t*(-pi*B26*L*cosi2x*i2*r*sina + sini2x*(A26*L**2*cosa*r + pi**2*B16*i2**2*r**2 + B26*L**2*j2**2)))*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt)))/(L**3*r**3))
                                tmp9 += weightx[nx][i]*weightt[nt][j]*((-L*l2*sink2x*sinl2t*(w0t + wt)*(B26*L*cosj2t*j2*(L*sina*sini2x - 2*pi*cosi2x*i2*r) + sinj2t*(-pi*B22*L*cosi2x*i2*r*sina + sini2x*(A22*L**2*cosa*r + pi**2*B12*i2**2*r**2 + B22*L**2*j2**2))) + pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(B16*L*cosj2t*j2*(L*sina*sini2x - 2*pi*cosi2x*i2*r) + sinj2t*(-pi*B12*L*cosi2x*i2*r*sina + sini2x*(A12*L**2*cosa*r + pi**2*B11*i2**2*r**2 + B12*L**2*j2**2))) + r*(B66*L*cosj2t*j2*(L*sina*sini2x - 2*pi*cosi2x*i2*r) + sinj2t*(-pi*B26*L*cosi2x*i2*r*sina + sini2x*(A26*L**2*cosa*r + pi**2*B16*i2**2*r**2 + B26*L**2*j2**2)))*(-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt)))/(L**3*r**3))
                                tmp10 += weightx[nx][i]*weightt[nt][j]*((L*cosl2t*l2*sink2x*(w0t + wt)*(B26*L*j2*sinj2t*(-L*sina*sini2x + 2*pi*cosi2x*i2*r) + cosj2t*(-pi*B22*L*cosi2x*i2*r*sina + sini2x*(A22*L**2*cosa*r + pi**2*B12*i2**2*r**2 + B22*L**2*j2**2))) + pi*cosk2x*k2*r**2*sinl2t*(w0x + wx)*(B16*L*j2*sinj2t*(-L*sina*sini2x + 2*pi*cosi2x*i2*r) + cosj2t*(-pi*B12*L*cosi2x*i2*r*sina + sini2x*(A12*L**2*cosa*r + pi**2*B11*i2**2*r**2 + B12*L**2*j2**2))) + r*(B66*L*j2*sinj2t*(-L*sina*sini2x + 2*pi*cosi2x*i2*r) + cosj2t*(-pi*B26*L*cosi2x*i2*r*sina + sini2x*(A26*L**2*cosa*r + pi**2*B16*i2**2*r**2 + B26*L**2*j2**2)))*(L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt)))/(L**3*r**3))
                                tmp11 += weightx[nx][i]*weightt[nt][j]*((-L*l2*sink2x*sinl2t*(w0t + wt)*(B26*L*j2*sinj2t*(-L*sina*sini2x + 2*pi*cosi2x*i2*r) + cosj2t*(-pi*B22*L*cosi2x*i2*r*sina + sini2x*(A22*L**2*cosa*r + pi**2*B12*i2**2*r**2 + B22*L**2*j2**2))) + pi*cosk2x*cosl2t*k2*r**2*(w0x + wx)*(B16*L*j2*sinj2t*(-L*sina*sini2x + 2*pi*cosi2x*i2*r) + cosj2t*(-pi*B12*L*cosi2x*i2*r*sina + sini2x*(A12*L**2*cosa*r + pi**2*B11*i2**2*r**2 + B12*L**2*j2**2))) + r*(B66*L*j2*sinj2t*(-L*sina*sini2x + 2*pi*cosi2x*i2*r) + cosj2t*(-pi*B26*L*cosi2x*i2*r*sina + sini2x*(A26*L**2*cosa*r + pi**2*B16*i2**2*r**2 + B26*L**2*j2**2)))*(-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt)))/(L**3*r**3))

                        c += 1
                        rows[c] = row+0
                        cols[c] = col+4
                        k0Lv[c] = tmp0
                        c += 1
                        rows[c] = row+0
                        cols[c] = col+5
                        k0Lv[c] = tmp1
                        c += 1
                        rows[c] = row+1
                        cols[c] = col+4
                        k0Lv[c] = tmp2
                        c += 1
                        rows[c] = row+1
                        cols[c] = col+5
                        k0Lv[c] = tmp3
                        c += 1
                        rows[c] = row+2
                        cols[c] = col+4
                        k0Lv[c] = tmp4
                        c += 1
                        rows[c] = row+2
                        cols[c] = col+5
                        k0Lv[c] = tmp5
                        c += 1
                        rows[c] = row+3
                        cols[c] = col+4
                        k0Lv[c] = tmp6
                        c += 1
                        rows[c] = row+3
                        cols[c] = col+5
                        k0Lv[c] = tmp7
                        c += 1
                        rows[c] = row+4
                        cols[c] = col+4
                        k0Lv[c] = tmp8
                        c += 1
                        rows[c] = row+4
                        cols[c] = col+5
                        k0Lv[c] = tmp9
                        c += 1
                        rows[c] = row+5
                        cols[c] = col+4
                        k0Lv[c] = tmp10
                        c += 1
                        rows[c] = row+5
                        cols[c] = col+5
                        k0Lv[c] = tmp11

    nx = max(m1-i0, m2-i0)
    nt = n2-j0

    with nogil:
        for i in range(nx):
            free(weightx[i])
            free(xs[i])
        for j in range(nt):
            free(weightt[j])
            free(ts[j])
    free(weightx)
    free(xs)
    free(weightt)
    free(ts)

    size = num0 + num1*m1 + num2*m2*n2

    k0L = coo_matrix((k0Lv, (rows, cols)), shape=(size, size))

    return k0L


def calc_k0L(np.ndarray[cDOUBLE, ndim=1] coeffs,
             double alpharad, double r2, double L, double tLA,
             np.ndarray[cDOUBLE, ndim=2] F,
             int m1, int m2, int n2,
             int nx, int nt, int num_cores, str method,
             np.ndarray[cDOUBLE, ndim=1] c0, int m0, int n0):
    cdef double sina, cosa, xa, xb, ta, tb
    cdef int c, row, col
    cdef int i1, k1, i2, j2, k2, l2
    cdef np.ndarray[cINT, ndim=1] rows, cols
    cdef np.ndarray[cDOUBLE, ndim=1] k0Lv

    cdef cc_attributes args

    fdim = 3*m1 + 3*2*m2*n2 + 3*m1**2 + 2*6*m1*m2*n2 + 12*m2**2*n2**2

    k0Lv = np.zeros((fdim,), dtype=DOUBLE)
    rows = np.zeros((fdim,), dtype=INT)
    cols = np.zeros((fdim,), dtype=INT)

    sina = sin(alpharad)
    cosa = cos(alpharad)

    args.sina = &sina
    args.cosa = &cosa
    args.tLA = &tLA
    args.r2 = &r2
    args.L = &L
    args.F = &F[0,0]
    args.m1 = &m1
    args.m2 = &m2
    args.n2 = &n2
    args.coeffs = &coeffs[0]
    args.c0 = &c0[0]
    args.m0 = &m0
    args.n0 = &n0

    xa = 0.
    xb = L
    ta = 0.
    tb = 2*pi

    # numerical integration
    integratev(<void *>cfk0L, fdim, k0Lv, xa, xb, nx, ta, tb, nt,
               &args, num_cores, method)

    c = -1

    for k1 in range(i0, m1+i0):
        col = (k1-i0)*num1 + num0
        # k0L_01
        c += 1
        rows[c] = 0
        cols[c] = col+2
        c += 1
        rows[c] = 1
        cols[c] = col+2
        c += 1
        rows[c] = 2
        cols[c] = col+2

    for k2 in range(i0, m2+i0):
        for l2 in range(j0, n2+j0):
            col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1
            # k0L_02
            c += 1
            rows[c] = 0
            cols[c] = col+4
            c += 1
            rows[c] = 0
            cols[c] = col+5
            c += 1
            rows[c] = 1
            cols[c] = col+4
            c += 1
            rows[c] = 1
            cols[c] = col+5
            c += 1
            rows[c] = 2
            cols[c] = col+4
            c += 1
            rows[c] = 2
            cols[c] = col+5

    for i1 in range(i0, m1+i0):
        row = (i1-i0)*num1 + num0
        for k1 in range(i0, m1+i0):
            col = (k1-i0)*num1 + num0
            # k0L_11
            c += 1
            rows[c] = row+0
            cols[c] = col+2
            c += 1
            rows[c] = row+1
            cols[c] = col+2
            c += 1
            rows[c] = row+2
            cols[c] = col+2

        for k2 in range(i0, m2+i0):
            for l2 in range(j0, n2+j0):
                col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1
                # k0L_12
                c += 1
                rows[c] = row+0
                cols[c] = col+4
                c += 1
                rows[c] = row+0
                cols[c] = col+5
                c += 1
                rows[c] = row+1
                cols[c] = col+4
                c += 1
                rows[c] = row+1
                cols[c] = col+5
                c += 1
                rows[c] = row+2
                cols[c] = col+4
                c += 1
                rows[c] = row+2
                cols[c] = col+5

    for i2 in range(i0, m2+i0):
        for j2 in range(j0, n2+j0):
            row = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1
            for k1 in range(i0, m1+i0):
                col = (k1-i0)*num1 + num0
                # k0L_21
                c += 1
                rows[c] = row+0
                cols[c] = col+2
                c += 1
                rows[c] = row+1
                cols[c] = col+2
                c += 1
                rows[c] = row+2
                cols[c] = col+2
                c += 1
                rows[c] = row+3
                cols[c] = col+2
                c += 1
                rows[c] = row+4
                cols[c] = col+2
                c += 1
                rows[c] = row+5
                cols[c] = col+2

            for k2 in range(i0, m2+i0):
                for l2 in range(j0, n2+j0):
                    col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1

                    # k0L_22
                    c += 1
                    rows[c] = row+0
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+0
                    cols[c] = col+5
                    c += 1
                    rows[c] = row+1
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+1
                    cols[c] = col+5
                    c += 1
                    rows[c] = row+2
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+2
                    cols[c] = col+5
                    c += 1
                    rows[c] = row+3
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+3
                    cols[c] = col+5
                    c += 1
                    rows[c] = row+4
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+4
                    cols[c] = col+5
                    c += 1
                    rows[c] = row+5
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+5
                    cols[c] = col+5

    size = num0 + num1*m1 + num2*m2*n2

    k0L = coo_matrix((k0Lv, (rows, cols)), shape=(size, size))

    return k0L


cdef void cfk0L(int npts, double *xs, double *ts, double *out,
                double *alphas, double *betas, void *args) nogil:
    cdef int i1, k1, i2, j2, k2, l2
    cdef int c, i, pos

    cdef double A11, A12, A16, A22, A26, A66
    cdef double B11, B12, B16, B22, B26, B66

    cdef double p00, p01, p02, p10, p11, p12
    cdef double p20, p21, p22, p30, p31, p32
    cdef double p40, p41, p42, p50, p51, p52
    cdef double q02, q04, q05, q14, q15, q22, q24, q25

    cdef double r, x, t, alpha, beta

    cdef double *F
    cdef double *coeffs
    cdef double *c0
    cdef double sina, cosa, tLA, r2, L
    cdef int m0, n0, m1, m2, n2
    cdef double wx, wt, w0x, w0t

    cdef cc_attributes *args_in=<cc_attributes *>args

    sina = args_in.sina[0]
    cosa = args_in.cosa[0]
    tLA = args_in.tLA[0]
    r2 = args_in.r2[0]
    L = args_in.L[0]
    F = args_in.F
    m1 = args_in.m1[0]
    m2 = args_in.m2[0]
    n2 = args_in.n2[0]
    coeffs = args_in.coeffs
    c0 = args_in.c0
    m0 = args_in.m0[0]
    n0 = args_in.n0[0]

    A11 = F[0]
    A12 = F[1]
    A16 = F[2]
    A22 = F[7]
    A26 = F[8]
    A66 = F[14]
    B11 = F[3]
    B12 = F[4]
    B16 = F[5]
    B22 = F[10]
    B26 = F[11]
    B66 = F[17]

    cdef double sini1x, cosi1x, cosk1x, sini2x, cosi2x, sink2x, cosk2x
    cdef double sinl2t, cosl2t, sinj2t, cosj2t
    cdef double *vsini1x = <double *>malloc(m1 * sizeof(double))
    cdef double *vcosi1x = <double *>malloc(m1 * sizeof(double))
    cdef double *vsini2x = <double *>malloc(m2 * sizeof(double))
    cdef double *vcosi2x = <double *>malloc(m2 * sizeof(double))
    cdef double *vsinj2t = <double *>malloc(n2 * sizeof(double))
    cdef double *vcosj2t = <double *>malloc(n2 * sizeof(double))
    cdef double *k0Lq_1_q02 = <double *>malloc(m1 * sizeof(double))
    cdef double *k0Lq_1_q22 = <double *>malloc(m1 * sizeof(double))
    cdef double *k0Lq_2_q04 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *k0Lq_2_q05 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *k0Lq_2_q14 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *k0Lq_2_q15 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *k0Lq_2_q24 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *k0Lq_2_q25 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *wxs = <double *>malloc(npts * sizeof(double))
    cdef double *wts = <double *>malloc(npts * sizeof(double))
    cdef double *w0xs = <double *>malloc(npts * sizeof(double))
    cdef double *w0ts = <double *>malloc(npts * sizeof(double))

    cfwx(coeffs, m1, m2, n2, xs, ts, npts, L, wxs)
    cfwt(coeffs, m1, m2, n2, xs, ts, npts, L, wts)
    cfw0x(xs, ts, npts, c0, L, m0, n0, w0xs, funcnum)
    cfw0t(xs, ts, npts, c0, L, m0, n0, w0ts, funcnum)

    for i in range(npts):
        x = xs[i]
        t = ts[i]
        wx = wxs[i]
        wt = wts[i]
        w0x = w0xs[i]
        w0t = w0ts[i]
        alpha = alphas[i]
        beta = betas[i]

        for i1 in range(i0, m1+i0):
            vsini1x[i1-i0] = sin(pi*i1*x/L)
            vcosi1x[i1-i0] = cos(pi*i1*x/L)

        for i2 in range(i0, m2+i0):
            vsini2x[i2-i0] = sin(pi*i2*x/L)
            vcosi2x[i2-i0] = cos(pi*i2*x/L)

        for j2 in range(j0, n2+j0):
            vsinj2t[j2-j0] = sin(j2*t)
            vcosj2t[j2-j0] = cos(j2*t)

        r = r2 + sina*x

        c = -1

        # p_0
        p00 = (-A11*r + A12*sina*(L - x))/(L*cosa)
        p01 = (-A12*r + A22*sina*(L - x))/(L*cosa)
        p02 = (-A16*r + A26*sina*(L - x))/(L*cosa)
        p10 = -A16*r2*(r + sina*(L - x))/L
        p11 = -A26*r2*(r + sina*(L - x))/L
        p12 = -A66*r2*(r + sina*(L - x))/L
        p20 = (A16*(-L + x)*sin(t - tLA) - (A11*r + A12*sina*(-L + x))*cos(t - tLA))/(L*cosa)
        p21 = (A26*(-L + x)*sin(t - tLA) - (A12*r + A22*sina*(-L + x))*cos(t - tLA))/(L*cosa)
        p22 = (A66*(-L + x)*sin(t - tLA) - (A16*r + A26*sina*(-L + x))*cos(t - tLA))/(L*cosa)

        for k1 in range(i0, m1+i0):
            cosk1x = vcosi1x[k1-i0]
            # q_1
            q02 = pi*cosk1x*k1*(w0x + wx)/L
            q22 = pi*cosk1x*k1*(w0t + wt)/(L*r)

            # k0L_01
            c += 1
            out[c] = beta*out[c] + alpha*(p00*q02 + p02*q22)
            c += 1
            out[c] = beta*out[c] + alpha*(p10*q02 + p12*q22)
            c += 1
            out[c] = beta*out[c] + alpha*(p20*q02 + p22*q22)

            # create buffer q_1
            k0Lq_1_q02[k1-i0] = q02
            k0Lq_1_q22[k1-i0] = q22

        for k2 in range(i0, m2+i0):
            sink2x = vsini2x[k2-i0]
            cosk2x = vcosi2x[k2-i0]
            for l2 in range(j0, n2+j0):
                sinl2t = vsinj2t[l2-j0]
                cosl2t = vcosj2t[l2-j0]
                # q_2
                q04 = pi*cosk2x*k2*sinl2t*(w0x + wx)/L
                q05 = pi*cosk2x*cosl2t*k2*(w0x + wx)/L
                q14 = cosl2t*l2*sink2x*(w0t + wt)/(r*r)
                q15 = -l2*sink2x*sinl2t*(w0t + wt)/(r*r)
                q24 = (L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt))/(L*r)
                q25 = (-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt))/(L*r)

                # k0L_02
                c += 1
                out[c] = beta*out[c] + alpha*(p00*q04 + p01*q14 + p02*q24)
                c += 1
                out[c] = beta*out[c] + alpha*(p00*q05 + p01*q15 + p02*q25)
                c += 1
                out[c] = beta*out[c] + alpha*(p10*q04 + p11*q14 + p12*q24)
                c += 1
                out[c] = beta*out[c] + alpha*(p10*q05 + p11*q15 + p12*q25)
                c += 1
                out[c] = beta*out[c] + alpha*(p20*q04 + p21*q14 + p22*q24)
                c += 1
                out[c] = beta*out[c] + alpha*(p20*q05 + p21*q15 + p22*q25)

                # create buffer q_2
                pos = (k2-i0)*n2 + (l2-j0)
                k0Lq_2_q04[pos] = q04
                k0Lq_2_q05[pos] = q05
                k0Lq_2_q14[pos] = q14
                k0Lq_2_q15[pos] = q15
                k0Lq_2_q24[pos] = q24
                k0Lq_2_q25[pos] = q25

        for i1 in range(i0, m1+i0):
            sini1x = vsini1x[i1-i0]
            cosi1x = vcosi1x[i1-i0]
            # p_1
            p00 = pi*A11*i1*r*cosi1x/L + A12*sina*sini1x
            p01 = pi*A12*i1*r*cosi1x/L + A22*sina*sini1x
            p02 = pi*A16*i1*r*cosi1x/L + A26*sina*sini1x
            p10 = -A16*sina*sini1x + pi*A16*i1*r*cosi1x/L
            p11 = -A26*sina*sini1x + pi*A26*i1*r*cosi1x/L
            p12 = -A66*sina*sini1x + pi*A66*i1*r*cosi1x/L
            p20 = (-pi*B12*L*i1*sina*cosi1x + (A12*(L*L)*cosa + (pi*pi)*B11*(i1*i1)*r)*sini1x)/(L*L)
            p21 = (-pi*B22*L*i1*sina*cosi1x + (A22*(L*L)*cosa + (pi*pi)*B12*(i1*i1)*r)*sini1x)/(L*L)
            p22 = (-pi*B26*L*i1*sina*cosi1x + (A26*(L*L)*cosa + (pi*pi)*B16*(i1*i1)*r)*sini1x)/(L*L)

            for k1 in range(i0, m1+i0):
                # access buffer q_1
                q02 = k0Lq_1_q02[k1-i0]
                q22 = k0Lq_1_q22[k1-i0]

                # k0L_11
                c += 1
                out[c] = beta*out[c] + alpha*(p00*q02 + p02*q22)
                c += 1
                out[c] = beta*out[c] + alpha*(p10*q02 + p12*q22)
                c += 1
                out[c] = beta*out[c] + alpha*(p20*q02 + p22*q22)

            for k2 in range(i0, m2+i0):
                for l2 in range(j0, n2+j0):
                    # access buffer q_2
                    pos = (k2-i0)*n2 + (l2-j0)
                    q04 = k0Lq_2_q04[pos]
                    q05 = k0Lq_2_q05[pos]
                    q14 = k0Lq_2_q14[pos]
                    q15 = k0Lq_2_q15[pos]
                    q24 = k0Lq_2_q24[pos]
                    q25 = k0Lq_2_q25[pos]

                    # k0L_12
                    c += 1
                    out[c] = beta*out[c] + alpha*(p00*q04 + p01*q14 + p02*q24)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p00*q05 + p01*q15 + p02*q25)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p10*q04 + p11*q14 + p12*q24)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p10*q05 + p11*q15 + p12*q25)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p20*q04 + p21*q14 + p22*q24)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p20*q05 + p21*q15 + p22*q25)

        for i2 in range(i0, m2+i0):
            sini2x = vsini2x[i2-i0]
            cosi2x = vcosi2x[i2-i0]
            for j2 in range(j0, n2+j0):
                sinj2t = vsinj2t[j2-j0]
                cosj2t = vcosj2t[j2-j0]
                # p_2
                p00 = pi*A11*i2*r*sinj2t*cosi2x/L + (A12*sina*sinj2t + A16*j2*cosj2t)*sini2x
                p01 = pi*A12*i2*r*sinj2t*cosi2x/L + (A22*sina*sinj2t + A26*j2*cosj2t)*sini2x
                p02 = pi*A16*i2*r*sinj2t*cosi2x/L + (A26*sina*sinj2t + A66*j2*cosj2t)*sini2x
                p10 = pi*A11*i2*r*cosj2t*cosi2x/L + (A12*sina*cosj2t - A16*j2*sinj2t)*sini2x
                p11 = pi*A12*i2*r*cosj2t*cosi2x/L + (A22*sina*cosj2t - A26*j2*sinj2t)*sini2x
                p12 = pi*A16*i2*r*cosj2t*cosi2x/L + (A26*sina*cosj2t - A66*j2*sinj2t)*sini2x
                p20 = pi*A16*i2*r*sinj2t*cosi2x/L + (A12*j2*cosj2t - A16*sina*sinj2t)*sini2x
                p21 = pi*A26*i2*r*sinj2t*cosi2x/L + (A22*j2*cosj2t - A26*sina*sinj2t)*sini2x
                p22 = pi*A66*i2*r*sinj2t*cosi2x/L + (A26*j2*cosj2t - A66*sina*sinj2t)*sini2x
                p30 = pi*A16*i2*r*cosj2t*cosi2x/L - (A12*j2*sinj2t + A16*sina*cosj2t)*sini2x
                p31 = pi*A26*i2*r*cosj2t*cosi2x/L - (A22*j2*sinj2t + A26*sina*cosj2t)*sini2x
                p32 = pi*A66*i2*r*cosj2t*cosi2x/L - (A26*j2*sinj2t + A66*sina*cosj2t)*sini2x
                p40 = (-pi*L*i2*r*(B12*sina*sinj2t + 2*B16*j2*cosj2t)*cosi2x + (B16*(L*L)*j2*sina*cosj2t + (B12*(L*L)*(j2*j2) + r*(A12*(L*L)*cosa + (pi*pi)*B11*(i2*i2)*r))*sinj2t)*sini2x)/((L*L)*r)
                p41 = (-pi*L*i2*r*(B22*sina*sinj2t + 2*B26*j2*cosj2t)*cosi2x + (B26*(L*L)*j2*sina*cosj2t + (B22*(L*L)*(j2*j2) + r*(A22*(L*L)*cosa + (pi*pi)*B12*(i2*i2)*r))*sinj2t)*sini2x)/((L*L)*r)
                p42 = (-pi*L*i2*r*(B26*sina*sinj2t + 2*B66*j2*cosj2t)*cosi2x + (B66*(L*L)*j2*sina*cosj2t + (B26*(L*L)*(j2*j2) + r*(A26*(L*L)*cosa + (pi*pi)*B16*(i2*i2)*r))*sinj2t)*sini2x)/((L*L)*r)
                p50 = (pi*L*i2*r*(-B12*sina*cosj2t + 2*B16*j2*sinj2t)*cosi2x + (-B16*(L*L)*j2*sina*sinj2t + (B12*(L*L)*(j2*j2) + r*(A12*(L*L)*cosa + (pi*pi)*B11*(i2*i2)*r))*cosj2t)*sini2x)/((L*L)*r)
                p51 = (pi*L*i2*r*(-B22*sina*cosj2t + 2*B26*j2*sinj2t)*cosi2x + (-B26*(L*L)*j2*sina*sinj2t + (B22*(L*L)*(j2*j2) + r*(A22*(L*L)*cosa + (pi*pi)*B12*(i2*i2)*r))*cosj2t)*sini2x)/((L*L)*r)
                p52 = (pi*L*i2*r*(-B26*sina*cosj2t + 2*B66*j2*sinj2t)*cosi2x + (-B66*(L*L)*j2*sina*sinj2t + (B26*(L*L)*(j2*j2) + r*(A26*(L*L)*cosa + (pi*pi)*B16*(i2*i2)*r))*cosj2t)*sini2x)/((L*L)*r)

                for k1 in range(i0, m1+i0):
                    # access buffer q_1
                    q02 = k0Lq_1_q02[k1-i0]
                    q22 = k0Lq_1_q22[k1-i0]

                    # k0L_21
                    c += 1
                    out[c] = beta*out[c] + alpha*(p00*q02 + p02*q22)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p10*q02 + p12*q22)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p20*q02 + p22*q22)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p30*q02 + p32*q22)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p40*q02 + p42*q22)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p50*q02 + p52*q22)

                for k2 in range(i0, m2+i0):
                    for l2 in range(j0, n2+j0):
                        # access buffer q_2
                        pos = (k2-i0)*n2 + (l2-j0)
                        q04 = k0Lq_2_q04[pos]
                        q05 = k0Lq_2_q05[pos]
                        q14 = k0Lq_2_q14[pos]
                        q15 = k0Lq_2_q15[pos]
                        q24 = k0Lq_2_q24[pos]
                        q25 = k0Lq_2_q25[pos]
                        # k0L_22
                        c += 1
                        out[c] = beta*out[c] + alpha*(p00*q04 + p01*q14 + p02*q24)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p00*q05 + p01*q15 + p02*q25)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p10*q04 + p11*q14 + p12*q24)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p10*q05 + p11*q15 + p12*q25)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p20*q04 + p21*q14 + p22*q24)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p20*q05 + p21*q15 + p22*q25)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p30*q04 + p31*q14 + p32*q24)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p30*q05 + p31*q15 + p32*q25)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p40*q04 + p41*q14 + p42*q24)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p40*q05 + p41*q15 + p42*q25)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p50*q04 + p51*q14 + p52*q24)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p50*q05 + p51*q15 + p52*q25)

    free(wxs)
    free(wts)
    free(w0xs)
    free(w0ts)

    free(vsini1x)
    free(vcosi1x)
    free(vsini2x)
    free(vcosi2x)
    free(vsinj2t)
    free(vcosj2t)
    free(k0Lq_1_q02)
    free(k0Lq_1_q22)
    free(k0Lq_2_q04)
    free(k0Lq_2_q05)
    free(k0Lq_2_q14)
    free(k0Lq_2_q15)
    free(k0Lq_2_q24)
    free(k0Lq_2_q25)


def calc_kG(np.ndarray[cDOUBLE, ndim=1] coeffs,
            double alpharad, double r2, double L, double tLA,
            np.ndarray[cDOUBLE, ndim=2] F,
            int m1, int m2, int n2,
            int nx, int nt, int num_cores, str method,
            np.ndarray[cDOUBLE, ndim=1] c0, int m0, int n0):
    cdef double sina, cosa, xa, xb, ta, tb
    cdef int c, row, col
    cdef int i1, k1, i2, j2, k2, l2
    cdef int size

    cdef np.ndarray[cINT, ndim=1] rows, cols
    cdef np.ndarray[cDOUBLE, ndim=1] kGv

    cdef cc_attributes args

    fdim = 1*m1**2 + 2*m1*m2*n2 + 4*m2**2*n2**2

    rows = np.zeros((fdim,), dtype=INT)
    cols = np.zeros((fdim,), dtype=INT)
    kGv = np.zeros((fdim,), dtype=DOUBLE)

    sina = sin(alpharad)
    cosa = cos(alpharad)

    args.sina = &sina
    args.cosa = &cosa
    args.tLA = &tLA
    args.r2 = &r2
    args.L = &L
    args.F = &F[0,0]
    args.m1 = &m1
    args.m2 = &m2
    args.n2 = &n2
    args.coeffs = &coeffs[0]
    args.c0 = &c0[0]
    args.m0 = &m0
    args.n0 = &n0

    xa = 0.
    xb = L
    ta = 0.
    tb = 2*pi

    # numerical integration
    integratev(<void *>cfkG, fdim, kGv, xa, xb, nx, ta, tb, nt,
               &args, num_cores, method)

    c = -1

    for i1 in range(i0, m1+i0):
        row = (i1-i0)*num1 + num0
        for k1 in range(i0, m1+i0):
            col = (k1-i0)*num1 + num0

            #NOTE symmetry
            if row > col:
                continue

            # kG_11
            c += 1
            rows[c] = row+2
            cols[c] = col+2

        for k2 in range(i0, m2+i0):
            for l2 in range(j0, n2+j0):
                col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1
                # kG_12
                c += 1
                rows[c] = row+2
                cols[c] = col+4
                c += 1
                rows[c] = row+2
                cols[c] = col+5

    for i2 in range(i0, m2+i0):
        for j2 in range(j0, n2+j0):
            row = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1
            for k2 in range(i0, m2+i0):
                for l2 in range(j0, n2+j0):
                    col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1

                    #NOTE symmetry
                    if row > col:
                        continue

                    # kG_22
                    c += 1
                    rows[c] = row+4
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+4
                    cols[c] = col+5
                    c += 1
                    rows[c] = row+5
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+5
                    cols[c] = col+5

    size = num0 + num1*m1 + num2*m2*n2

    kG = coo_matrix((kGv, (rows, cols)), shape=(size, size))

    return kG


cdef void cfkG(int npts, double *xs, double *ts, double *out,
               double *alphas, double *betas, void *args) nogil:
    cdef int i1, k1, i2, j2, k2, l2

    cdef double q02, q04, q05, q14, q15
    cdef double p20, p21, p40, p41, p50, p51

    cdef double r, x, t, alpha, beta
    cdef int c, i, pos, row, col

    cdef double *F
    cdef double *coeffs
    cdef double *c0
    cdef double r2, L, sina, cosa, tLA
    cdef int m0, n0, m1, m2, n2
    cdef double Nxx, Ntt, Nxt

    cdef cc_attributes *args_in=<cc_attributes *>args

    sina = args_in.sina[0]
    cosa = args_in.cosa[0]
    tLA = args_in.tLA[0]
    r2 = args_in.r2[0]
    L = args_in.L[0]
    F = args_in.F
    m1 = args_in.m1[0]
    m2 = args_in.m2[0]
    n2 = args_in.n2[0]
    coeffs = args_in.coeffs
    c0 = args_in.c0
    m0 = args_in.m0[0]
    n0 = args_in.n0[0]

    cdef double cosi1x, cosk1x, sini2x, cosi2x, sink2x, cosk2x
    cdef double sinl2t, cosl2t, sinj2t, cosj2t
    cdef double *vcosi1x = <double *>malloc(m1 * sizeof(double))
    cdef double *vsini2x = <double *>malloc(m2 * sizeof(double))
    cdef double *vcosi2x = <double *>malloc(m2 * sizeof(double))
    cdef double *vsinj2t = <double *>malloc(n2 * sizeof(double))
    cdef double *vcosj2t = <double *>malloc(n2 * sizeof(double))
    cdef double *kGq_2_q04 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *kGq_2_q05 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *kGq_2_q14 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *kGq_2_q15 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *Ns = <double *>malloc(e_num*npts * sizeof(double))

    cfN(coeffs, sina, cosa, tLA, xs, ts, npts, r2, L, F, m1, m2, n2,
        c0, m0, n0, funcnum, Ns, NL_kinematics)

    for i in range(npts):
        x = xs[i]
        t = ts[i]
        Nxx = Ns[e_num*i + 0]
        Ntt = Ns[e_num*i + 1]
        Nxt = Ns[e_num*i + 2]
        alpha = alphas[i]
        beta = betas[i]

        for i1 in range(i0, m1+i0):
            vcosi1x[i1-i0] = cos(pi*i1*x/L)

        for i2 in range(i0, m2+i0):
            vsini2x[i2-i0] = sin(pi*i2*x/L)
            vcosi2x[i2-i0] = cos(pi*i2*x/L)

        for j2 in range(j0, n2+j0):
            vsinj2t[j2-j0] = sin(j2*t)
            vcosj2t[j2-j0] = cos(j2*t)

        r = r2 + sina*x

        c = -1

        for k2 in range(i0, m2+i0):
            sink2x = vsini2x[k2-i0]
            cosk2x = vcosi2x[k2-i0]
            for l2 in range(j0, n2+j0):
                sinl2t = vsinj2t[l2-j0]
                cosl2t = vcosj2t[l2-j0]
                # q_2
                q04 = pi*k2*sinl2t*cosk2x/L
                q05 = pi*k2*cosl2t*cosk2x/L
                q14 = l2*sink2x*cosl2t/r
                q15 = -l2*sinl2t*sink2x/r
                # create buffer q_2
                pos = (k2-i0)*n2 + (l2-j0)
                kGq_2_q04[pos] = q04
                kGq_2_q05[pos] = q05
                kGq_2_q14[pos] = q14
                kGq_2_q15[pos] = q15

        for i1 in range(i0, m1+i0):
            row = (i1-i0)*num1 + num0

            cosi1x = vcosi1x[i1-i0]
            # p_1
            p20 = pi*Nxx*i1*r*cosi1x/L
            p21 = pi*Nxt*i1*r*cosi1x/L

            for k1 in range(i0, m1+i0):
                col = (k1-i0)*num1 + num0

                #NOTE symmetry
                if row > col:
                    continue

                cosk1x = vcosi1x[k1-i0]
                # q_1
                q02 = pi*k1*cosk1x/L

                # kG_11
                c += 1
                out[c] = beta*out[c] + alpha*(p20*q02)

            for k2 in range(i0, m2+i0):
                for l2 in range(j0, n2+j0):
                    # access buffer q_2
                    pos = (k2-i0)*n2 + (l2-j0)
                    q04 = kGq_2_q04[pos]
                    q05 = kGq_2_q05[pos]
                    q14 = kGq_2_q14[pos]
                    q15 = kGq_2_q15[pos]
                    # kG_12
                    c += 1
                    out[c] = beta*out[c] + alpha*(p20*q04 + p21*q14)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p20*q05 + p21*q15)

        for i2 in range(i0, m2+i0):
            sini2x = vsini2x[i2-i0]
            cosi2x = vcosi2x[i2-i0]
            for j2 in range(j0, n2+j0):
                row = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1

                sinj2t = vsinj2t[j2-j0]
                cosj2t = vcosj2t[j2-j0]
                # p_2
                p40 = Nxt*j2*sini2x*cosj2t + pi*Nxx*i2*r*sinj2t*cosi2x/L
                p41 = Ntt*j2*sini2x*cosj2t + pi*Nxt*i2*r*sinj2t*cosi2x/L
                p50 = -Nxt*j2*sinj2t*sini2x + pi*Nxx*i2*r*cosj2t*cosi2x/L
                p51 = -Ntt*j2*sinj2t*sini2x + pi*Nxt*i2*r*cosj2t*cosi2x/L

                for k2 in range(i0, m2+i0):
                    for l2 in range(j0, n2+j0):
                        col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1

                        #NOTE symmetry
                        if row > col:
                            continue

                        # access buffer q_2
                        pos = (k2-i0)*n2 + (l2-j0)
                        q04 = kGq_2_q04[pos]
                        q05 = kGq_2_q05[pos]
                        q14 = kGq_2_q14[pos]
                        q15 = kGq_2_q15[pos]
                        # kG_22
                        c += 1
                        out[c] = beta*out[c] + alpha*(p40*q04 + p41*q14)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p40*q05 + p41*q15)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p50*q04 + p51*q14)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p50*q05 + p51*q15)

    free(Ns)

    free(vcosi1x)
    free(vsini2x)
    free(vcosi2x)
    free(vsinj2t)
    free(vcosj2t)
    free(kGq_2_q04)
    free(kGq_2_q05)
    free(kGq_2_q14)
    free(kGq_2_q15)


def calc_kLL(np.ndarray[cDOUBLE, ndim=1] coeffs,
             double alpharad, double r2, double L, double tLA,
             np.ndarray[cDOUBLE, ndim=2] F,
             int m1, int m2, int n2,
             int nx, int nt, int num_cores, str method,
             np.ndarray[cDOUBLE, ndim=1] c0, int m0, int n0):
    cdef double sina, cosa, xa, xb, ta, tb
    cdef int c, row, col
    cdef int i1, k1, i2, j2, k2, l2
    cdef int size

    cdef np.ndarray[cINT, ndim=1] rows, cols
    cdef np.ndarray[cDOUBLE, ndim=1] kLLv

    cdef cc_attributes args

    fdim = 1*m1**2 + 2*m1*m2*n2 + 4*m2**2*n2**2

    rows = np.zeros((fdim,), dtype=INT)
    cols = np.zeros((fdim,), dtype=INT)
    kLLv = np.zeros((fdim,), dtype=DOUBLE)

    sina = sin(alpharad)
    cosa = cos(alpharad)

    args.sina = &sina
    args.cosa = &cosa
    args.tLA = &tLA
    args.r2 = &r2
    args.L = &L
    args.F = &F[0,0]
    args.m1 = &m1
    args.m2 = &m2
    args.n2 = &n2
    args.coeffs = &coeffs[0]
    args.c0 = &c0[0]
    args.m0 = &m0
    args.n0 = &n0

    xa = 0.
    xb = L
    ta = 0.
    tb = 2*pi

    # numerical integration
    integratev(<void *>cfkLL, fdim, kLLv, xa, xb, nx, ta, tb, nt,
               &args, num_cores, method)

    c = -1

    for i1 in range(i0, m1+i0):
        row = (i1-i0)*num1 + num0
        for k1 in range(i0, m1+i0):
            col = (k1-i0)*num1 + num0

            #NOTE symmetry
            if row > col:
                continue

            # kLL_11
            c += 1
            rows[c] = row+2
            cols[c] = col+2

        for k2 in range(i0, m2+i0):
            for l2 in range(j0, n2+j0):
                col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1
                # kLL_12
                c += 1
                rows[c] = row+2
                cols[c] = col+4
                c += 1
                rows[c] = row+2
                cols[c] = col+5

    for i2 in range(i0, m2+i0):
        for j2 in range(j0, n2+j0):
            row = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1
            for k2 in range(i0, m2+i0):
                for l2 in range(j0, n2+j0):
                    col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1

                    #NOTE symmetry
                    if row > col:
                        continue

                    # kLL_22
                    c += 1
                    rows[c] = row+4
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+4
                    cols[c] = col+5
                    c += 1
                    rows[c] = row+5
                    cols[c] = col+4
                    c += 1
                    rows[c] = row+5
                    cols[c] = col+5

    size = num0 + num1*m1 + num2*m2*n2

    kLL = coo_matrix((kLLv, (rows, cols)), shape=(size, size))

    return kLL


cdef void cfkLL(int npts, double *xs, double *ts, double *out,
                double *alphas, double *betas, void *args) nogil:
    cdef int i1, k1, i2, j2, k2, l2
    cdef int c, i, pos, row, col

    cdef double A11, A12, A16, A22, A26, A66

    cdef double p20, p21, p22, p40, p41, p42
    cdef double p50, p51, p52
    cdef double q02, q04, q05, q14, q15, q22, q24, q25

    cdef double r, x, t, alpha, beta

    cdef double *F
    cdef double *coeffs
    cdef double *c0
    cdef double sina, r2, L
    cdef int m0, n0, m1, m2, n2
    cdef double wx, wt, w0x, w0t

    cdef cc_attributes *args_in=<cc_attributes *>args
    sina = args_in.sina[0]
    r2 = args_in.r2[0]
    L = args_in.L[0]
    F = args_in.F
    m1 = args_in.m1[0]
    m2 = args_in.m2[0]
    n2 = args_in.n2[0]
    coeffs = args_in.coeffs
    c0 = args_in.c0
    m0 = args_in.m0[0]
    n0 = args_in.n0[0]

    A11 = F[0]
    A12 = F[1]
    A16 = F[2]
    A22 = F[7]
    A26 = F[8]
    A66 = F[14]

    cdef double sini1x, cosi1x, cosk1x, sini2x, cosi2x, sink2x, cosk2x
    cdef double sinl2t, cosl2t, sinj2t, cosj2t
    cdef double *vcosi1x = <double *>malloc(m1 * sizeof(double))
    cdef double *vsini2x = <double *>malloc(m2 * sizeof(double))
    cdef double *vcosi2x = <double *>malloc(m2 * sizeof(double))
    cdef double *vsinj2t = <double *>malloc(n2 * sizeof(double))
    cdef double *vcosj2t = <double *>malloc(n2 * sizeof(double))
    cdef double *kLLq_2_q04 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *kLLq_2_q05 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *kLLq_2_q14 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *kLLq_2_q15 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *kLLq_2_q24 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *kLLq_2_q25 = <double *>malloc(m2*n2 * sizeof(double))
    cdef double *wxs = <double *>malloc(npts * sizeof(double))
    cdef double *wts = <double *>malloc(npts * sizeof(double))
    cdef double *w0xs = <double *>malloc(npts * sizeof(double))
    cdef double *w0ts = <double *>malloc(npts * sizeof(double))

    cfwx(coeffs, m1, m2, n2, xs, ts, npts, L, wxs)
    cfwt(coeffs, m1, m2, n2, xs, ts, npts, L, wts)
    cfw0x(xs, ts, npts, c0, L, m0, n0, w0xs, funcnum)
    cfw0t(xs, ts, npts, c0, L, m0, n0, w0ts, funcnum)

    for i in range(npts):
        x = xs[i]
        t = ts[i]
        wx = wxs[i]
        wt = wts[i]
        w0x = w0xs[i]
        w0t = w0ts[i]
        alpha = alphas[i]
        beta = betas[i]

        for i1 in range(i0, m1+i0):
            vcosi1x[i1-i0] = cos(pi*i1*x/L)

        for i2 in range(i0, m2+i0):
            vsini2x[i2-i0] = sin(pi*i2*x/L)
            vcosi2x[i2-i0] = cos(pi*i2*x/L)

        for j2 in range(j0, n2+j0):
            vsinj2t[j2-j0] = sin(j2*t)
            vcosj2t[j2-j0] = cos(j2*t)

        r = r2 + sina*x

        c = -1

        for k2 in range(i0, m2+i0):
            sink2x = vsini2x[k2-i0]
            cosk2x = vcosi2x[k2-i0]
            for l2 in range(j0, n2+j0):
                sinl2t = vsinj2t[l2-j0]
                cosl2t = vcosj2t[l2-j0]
                # q_2
                q04 = pi*cosk2x*k2*sinl2t*(w0x + wx)/L
                q05 = pi*cosk2x*cosl2t*k2*(w0x + wx)/L
                q14 = cosl2t*l2*sink2x*(w0t + wt)/(r*r)
                q15 = -l2*sink2x*sinl2t*(w0t + wt)/(r*r)
                q24 = (L*cosl2t*l2*sink2x*(w0x + wx) + pi*cosk2x*k2*sinl2t*(w0t + wt))/(L*r)
                q25 = (-L*l2*sink2x*sinl2t*(w0x + wx) + pi*cosk2x*cosl2t*k2*(w0t + wt))/(L*r)
                # create buffer q_2
                pos = (k2-i0)*n2 + (l2-j0)
                kLLq_2_q04[pos] = q04
                kLLq_2_q05[pos] = q05
                kLLq_2_q14[pos] = q14
                kLLq_2_q15[pos] = q15
                kLLq_2_q24[pos] = q24
                kLLq_2_q25[pos] = q25

        for i1 in range(i0, m1+i0):
            row = (i1-i0)*num1 + num0

            cosi1x = vcosi1x[i1-i0]
            # p_1
            p20 = pi*cosi1x*i1*(A11*r*(w0x + wx) + A16*(w0t + wt))/L
            p21 = pi*cosi1x*i1*(A12*r*(w0x + wx) + A26*(w0t + wt))/L
            p22 = pi*cosi1x*i1*(A16*r*(w0x + wx) + A66*(w0t + wt))/L

            for k1 in range(i0, m1+i0):
                col = (k1-i0)*num1 + num0

                #NOTE symmetry
                if row > col:
                    continue

                cosk1x = vcosi1x[k1-i0]
                # q_1
                q02 = pi*cosk1x*k1*(w0x + wx)/L
                q22 = pi*cosk1x*k1*(w0t + wt)/(L*r)
                # kLL_11
                c += 1
                out[c] = beta*out[c] + alpha*(p20*q02 + p22*q22)

            for k2 in range(i0, m2+i0):
                for l2 in range(j0, n2+j0):
                    # access buffer q_2
                    pos = (k2-i0)*n2 + (l2-j0)
                    q04 = kLLq_2_q04[pos]
                    q05 = kLLq_2_q05[pos]
                    q14 = kLLq_2_q14[pos]
                    q15 = kLLq_2_q15[pos]
                    q24 = kLLq_2_q24[pos]
                    q25 = kLLq_2_q25[pos]
                    # kLL_12
                    c += 1
                    out[c] = beta*out[c] + alpha*(p20*q04 + p21*q14 + p22*q24)
                    c += 1
                    out[c] = beta*out[c] + alpha*(p20*q05 + p21*q15 + p22*q25)

        for i2 in range(i0, m2+i0):
            sini2x = vsini2x[i2-i0]
            cosi2x = vcosi2x[i2-i0]
            for j2 in range(j0, n2+j0):
                row = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1

                sinj2t = vsinj2t[j2-j0]
                cosj2t = vcosj2t[j2-j0]
                # p_2
                p40 = (L*cosj2t*j2*sini2x*(A12*(w0t + wt) + A16*r*(w0x + wx)) + pi*cosi2x*i2*r*sinj2t*(A11*r*(w0x + wx) + A16*(w0t + wt)))/(L*r)
                p41 = (L*cosj2t*j2*sini2x*(A22*(w0t + wt) + A26*r*(w0x + wx)) + pi*cosi2x*i2*r*sinj2t*(A12*r*(w0x + wx) + A26*(w0t + wt)))/(L*r)
                p42 = (L*cosj2t*j2*sini2x*(A26*(w0t + wt) + A66*r*(w0x + wx)) + pi*cosi2x*i2*r*sinj2t*(A16*r*(w0x + wx) + A66*(w0t + wt)))/(L*r)
                p50 = (-L*j2*sini2x*sinj2t*(A12*(w0t + wt) + A16*r*(w0x + wx)) + pi*cosi2x*cosj2t*i2*r*(A11*r*(w0x + wx) + A16*(w0t + wt)))/(L*r)
                p51 = (-L*j2*sini2x*sinj2t*(A22*(w0t + wt) + A26*r*(w0x + wx)) + pi*cosi2x*cosj2t*i2*r*(A12*r*(w0x + wx) + A26*(w0t + wt)))/(L*r)
                p52 = (-L*j2*sini2x*sinj2t*(A26*(w0t + wt) + A66*r*(w0x + wx)) + pi*cosi2x*cosj2t*i2*r*(A16*r*(w0x + wx) + A66*(w0t + wt)))/(L*r)
                for k2 in range(i0, m2+i0):
                    for l2 in range(j0, n2+j0):
                        col = (k2-i0)*num2 + (l2-j0)*num2*m2 + num0 + num1*m1

                        #NOTE symmetry
                        if row > col:
                            continue

                        # access buffer q_2
                        pos = (k2-i0)*n2 + (l2-j0)
                        q04 = kLLq_2_q04[pos]
                        q05 = kLLq_2_q05[pos]
                        q14 = kLLq_2_q14[pos]
                        q15 = kLLq_2_q15[pos]
                        q24 = kLLq_2_q24[pos]
                        q25 = kLLq_2_q25[pos]
                        # kLL_22
                        c += 1
                        out[c] = beta*out[c] + alpha*(p40*q04 + p41*q14 + p42*q24)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p40*q05 + p41*q15 + p42*q25)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p50*q04 + p51*q14 + p52*q24)
                        c += 1
                        out[c] = beta*out[c] + alpha*(p50*q05 + p51*q15 + p52*q25)

    free(wxs)
    free(wts)
    free(w0xs)
    free(w0ts)

    free(vcosi1x)
    free(vsini2x)
    free(vcosi2x)
    free(vsinj2t)
    free(vcosj2t)
    free(kLLq_2_q04)
    free(kLLq_2_q05)
    free(kLLq_2_q14)
    free(kLLq_2_q15)
    free(kLLq_2_q24)
    free(kLLq_2_q25)


def calc_fint_0L_L0_LL(np.ndarray[cDOUBLE, ndim=1] coeffs,
              double alpharad, double r2, double L, double tLA,
              np.ndarray[cDOUBLE, ndim=2] F,
              int m1, int m2, int n2,
              int nx, int nt, int num_cores, str method,
              np.ndarray[cDOUBLE, ndim=1] c0, int m0, int n0):
    cdef cc_attributes args
    cdef double sina, cosa

    fdim = num0 + num1*m1 + num2*m2*n2
    fint = np.zeros((fdim,), dtype=DOUBLE)

    sina = sin(alpharad)
    cosa = cos(alpharad)

    args.sina = &sina
    args.cosa = &cosa
    args.tLA = &tLA
    args.r2 = &r2
    args.L = &L
    args.F = &F[0,0]
    args.m1 = &m1
    args.m2 = &m2
    args.n2 = &n2
    args.coeffs = &coeffs[0]
    args.c0 = &c0[0]
    args.m0 = &m0
    args.n0 = &n0

    xa = 0.
    xb = L
    ta = 0.
    tb = 2*pi

    # numerical integration
    integratev(<void *>cffint, fdim, fint, xa, xb, nx, ta, tb, nt,
               &args, num_cores, method)

    return fint


cdef void cffint(int npts, double *xs, double *ts, double *fint,
                 double *alphas, double *betas, void *args) nogil:
    cdef double A11, A12, A16, A22, A26, A66
    cdef double B11, B12, B16, B22, B26, B66
    cdef double D11, D12, D16, D22, D26, D66
    cdef double Nxx0, Ntt0, Nxt0, Mxx0, Mtt0, Mxt0
    cdef double NxxL, NttL, NxtL, MxxL, MttL, MxtL
    cdef double exx0, ett0, gxt0, kxx0, ktt0, kxt0
    cdef double exxL, ettL, gxtL, kxxL, kttL, kxtL
    cdef double x, t, wx, wt, w0x, w0t, w0
    cdef double sini1x, cosi1x, sini2x, dsini2x, cosi2x, sinj2t, cosj2t

    cdef double alpha, beta

    cdef double *F
    cdef double *c
    cdef double *c0
    cdef double sina, cosa, tLA, r, r2, L
    cdef int m0, n0, m1, m2, n2
    cdef int i1, i2, j2, i, col

    cdef cc_attributes *args_in=<cc_attributes *>args

    sina = args_in.sina[0]
    cosa = args_in.cosa[0]
    tLA = args_in.tLA[0]
    r2 = args_in.r2[0]
    L = args_in.L[0]
    F = args_in.F
    m1 = args_in.m1[0]
    m2 = args_in.m2[0]
    n2 = args_in.n2[0]
    c = args_in.coeffs
    c0 = args_in.c0
    m0 = args_in.m0[0]
    n0 = args_in.n0[0]

    A11 = F[0]
    A12 = F[1]
    A16 = F[2]
    A22 = F[7]
    A26 = F[8]
    A66 = F[14]
    B11 = F[3]
    B12 = F[4]
    B16 = F[5]
    B22 = F[10]
    B26 = F[11]
    B66 = F[17]
    D11 = F[18]
    D12 = F[19]
    D16 = F[20]
    D22 = F[25]
    D26 = F[26]
    D66 = F[32]

    cdef double *w0xs = <double *>malloc(npts * sizeof(double))
    cdef double *w0ts = <double *>malloc(npts * sizeof(double))
    cdef double *vsini1x = <double *>malloc(m1 * sizeof(double))
    cdef double *vcosi1x = <double *>malloc(m1 * sizeof(double))
    cdef double *vsini2x = <double *>malloc(m2 * sizeof(double))
    cdef double *vcosi2x = <double *>malloc(m2 * sizeof(double))
    cdef double *vsinj2t = <double *>malloc(n2 * sizeof(double))
    cdef double *vcosj2t = <double *>malloc(n2 * sizeof(double))

    cfw0x(xs, ts, npts, c0, L, m0, n0, w0xs, funcnum)
    cfw0t(xs, ts, npts, c0, L, m0, n0, w0ts, funcnum)

    for i in range(npts):
        x = xs[i]
        t = ts[i]
        w0x = w0xs[i]
        w0t = w0ts[i]
        alpha = alphas[i]
        beta = betas[i]
        r = r2 + x*sina

        for i1 in range(i0, m1+i0):
            vsini1x[i1-i0] = sin(pi*i1*x/L)
            vcosi1x[i1-i0] = cos(pi*i1*x/L)

        for i2 in range(i0, m2+i0):
            vsini2x[i2-i0] = sin(pi*i2*x/L)
            vcosi2x[i2-i0] = cos(pi*i2*x/L)

        for j2 in range(j0, n2+j0):
            vsinj2t[j2-j0] = sin(j2*t)
            vcosj2t[j2-j0] = cos(j2*t)

        # wx and wt
        wx = 0.
        wt = 0.
        for i1 in range(i0, m1+i0):
            col = (i1-i0)*num1 + num0
            wx += c[col+2]*i1*pi/L*vcosi1x[i1-i0]
        for j2 in range(j0, n2+j0):
            sinj2t = vsinj2t[j2-j0]
            cosj2t = vcosj2t[j2-j0]
            for i2 in range(i0, m2+i0):
                col = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1
                dsini2x = i2*pi/L*vcosi2x[i2-i0]
                sini2x = vsini2x[i2-i0]
                wx += c[col+4]*dsini2x*sinj2t
                wx += c[col+5]*dsini2x*cosj2t

                wt += c[col+4]*sini2x*(j2*cosj2t)
                wt += c[col+5]*sini2x*(-j2*sinj2t)

        # strains
        exx0 = 0.
        ett0 = 0.
        gxt0 = 0.
        kxx0 = 0.
        ktt0 = 0.
        kxt0 = 0.

        exxL = 0.
        ettL = 0.
        gxtL = 0.
        kxxL = 0.
        kttL = 0.
        kxtL = 0.

        #TODO if castro=1
        w0 = 0.

        exx0 = (-c[0]/(L*cosa)
                -c[2]*cos(t - tLA)/(L*cosa))
        exxL = 0.5*castro*w0x*w0x

        ett0 = (c[0]*sina*(L - x)/(L*cosa*r)
                + c[2]*sina*(L - x)*cos(t - tLA)/(L*cosa*r))
        ettL = 0.5*castro*(2*cosa*r*w0 + w0t*w0t)/(r*r)

        gxt0 = (-c[1]*r2*(r + sina*(L - x))/(L*r)
               + c[2]*(-L + x)*sin(t - tLA)/(L*cosa*r))
        gxtL = castro*w0t*w0x/r

        for i1 in range(i0, m1+i0):
            sini1x = vsini1x[i1-i0]
            cosi1x = vcosi1x[i1-i0]
            col = (i1-i0)*num1 + num0

            exx0 += pi*c[col+0]*cosi1x*i1/L
            exxL += pi*c[col+2]*cosi1x*i1*(w0x + 0.5*wx)/L

            ett0 += (c[col+0]*sina*sini1x/r
                     + c[col+2]*cosa*sini1x/r)

            gxt0 += c[col+1]*(-sina*sini1x/r + pi*cosi1x*i1/L)
            gxtL += pi*c[col+2]*cosi1x*i1*(w0t + 0.5*wt)/(L*r)

            kxx0 += (pi*pi)*c[col+2]*(i1*i1)*sini1x/(L*L)

            ktt0 += -pi*c[col+2]*cosi1x*i1*sina/(L*r)

        for j2 in range(j0, n2+j0):
            sinj2t = vsinj2t[j2-j0]
            cosj2t = vcosj2t[j2-j0]
            for i2 in range(i0, m2+i0):
                sini2x = vsini2x[i2-i0]
                cosi2x = vcosi2x[i2-i0]
                col = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1

                exx0 += (pi*c[col+0]*cosi2x*i2*sinj2t/L
                         + pi*c[col+1]*cosi2x*cosj2t*i2/L)
                exxL += (pi*c[col+4]*cosi2x*i2*sinj2t*(w0x + 0.5*wx)/L
                         + pi*c[col+5]*cosi2x*cosj2t*i2*(w0x + 0.5*wx)/L)

                ett0 += (c[col+0]*sina*sini2x*sinj2t/r
                         + c[col+1]*cosj2t*sina*sini2x/r
                         + c[col+2]*cosj2t*j2*sini2x/r
                         -c[col+3]*j2*sini2x*sinj2t/r
                         + 0.5*c[col+4]*sini2x*2*cosa*r*sinj2t/(r*r)
                         + 0.5*c[col+5]*sini2x*2*cosa*cosj2t*r/(r*r))
                ettL += (c[col+4]*sini2x*cosj2t*j2*(w0t + 0.5*wt)/(r*r)
                         - c[col+5]*sini2x*j2*sinj2t*(w0t + 0.5*wt)/(r*r))

                gxt0 += (c[col+0]*cosj2t*j2*sini2x/r
                        -c[col+1]*j2*sini2x*sinj2t/r
                        + c[col+2]*sinj2t*(-L*sina*sini2x + pi*cosi2x*i2*r)/(L*r)
                        + c[col+3]*cosj2t*(-L*sina*sini2x + pi*cosi2x*i2*r)/(L*r))
                gxtL += (c[col+4]*(L*cosj2t*j2*sini2x*(w0x + 0.5*wx) + pi*cosi2x*i2*sinj2t*(w0t + 0.5*wt))/(L*r)
                        + c[col+5]*(-L*j2*sini2x*sinj2t*(w0x + 0.5*wx) + pi*cosi2x*cosj2t*i2*(w0t + 0.5*wt))/(L*r))

                kxx0 += ((pi*pi)*c[col+4]*(i2*i2)*sini2x*sinj2t/(L*L)
                         + (pi*pi)*c[col+5]*cosj2t*(i2*i2)*sini2x/(L*L))

                ktt0 += (c[col+4]*sinj2t*(L*(j2*j2)*sini2x - pi*cosi2x*i2*r*sina)/(L*(r*r))
                         + c[col+5]*cosj2t*(L*(j2*j2)*sini2x - pi*cosi2x*i2*r*sina)/(L*(r*r)))

                kxt0 += (c[col+4]*cosj2t*j2*(L*sina*sini2x - 2*pi*cosi2x*i2*r)/(L*(r*r))
                         + c[col+5]*j2*sinj2t*(-L*sina*sini2x + 2*pi*cosi2x*i2*r)/(L*(r*r)))

        # stresses
        Nxx0 = A11*exx0 + A12*ett0 + A16*gxt0 + B11*kxx0 + B12*ktt0 + B16*kxt0
        Ntt0 = A12*exx0 + A22*ett0 + A26*gxt0 + B12*kxx0 + B22*ktt0 + B26*kxt0
        Nxt0 = A16*exx0 + A26*ett0 + A66*gxt0 + B16*kxx0 + B26*ktt0 + B66*kxt0
        Mxx0 = B11*exx0 + B12*ett0 + B16*gxt0 + D11*kxx0 + D12*ktt0 + D16*kxt0
        Mtt0 = B12*exx0 + B22*ett0 + B26*gxt0 + D12*kxx0 + D22*ktt0 + D26*kxt0
        Mxt0 = B16*exx0 + B26*ett0 + B66*gxt0 + D16*kxx0 + D26*ktt0 + D66*kxt0

        NxxL = A11*exxL + A12*ettL + A16*gxtL + B11*kxxL + B12*kttL + B16*kxtL
        NttL = A12*exxL + A22*ettL + A26*gxtL + B12*kxxL + B22*kttL + B26*kxtL
        NxtL = A16*exxL + A26*ettL + A66*gxtL + B16*kxxL + B26*kttL + B66*kxtL
        MxxL = B11*exxL + B12*ettL + B16*gxtL + D11*kxxL + D12*kttL + D16*kxtL
        MttL = B12*exxL + B22*ettL + B26*gxtL + D12*kxxL + D22*kttL + D26*kxtL
        MxtL = B16*exxL + B26*ettL + B66*gxtL + D16*kxxL + D26*kttL + D66*kxtL

        fint[0] = beta*(fint[0]) + alpha*((NttL*sina*(L - x) - NxxL*r)/(L*cosa))
        fint[1] = beta*(fint[1]) + alpha*(-NxtL*r2*(r + sina*(L - x))/L)
        fint[2] = beta*(fint[2]) + alpha*((NxtL*(-L + x)*sin(t - tLA) - (NttL*sina*(-L + x) + NxxL*r)*cos(t - tLA))/(L*cosa))

        for i1 in range(i0, m1+i0):
            col = (i1-i0)*num1 + num0
            sini1x = vsini1x[i1-i0]
            cosi1x = vcosi1x[i1-i0]
            fint[col+0] = beta*(fint[col+0]) + alpha*(NttL*sina*sini1x + pi*NxxL*cosi1x*i1*r/L)
            fint[col+1] = beta*(fint[col+1]) + alpha*(-NxtL*sina*sini1x + pi*NxtL*cosi1x*i1*r/L)
            fint[col+2] = beta*(fint[col+2]) + alpha*((-pi*L*cosi1x*i1*(MttL*sina - Nxt0*(w0t + wt) - NxtL*w0t - NxtL*wt - Nxx0*r*w0x - Nxx0*r*wx - NxxL*r*w0x - NxxL*r*wx) + sini1x*((L*L)*NttL*cosa + (pi*pi)*MxxL*(i1*i1)*r))/(L*L))

        for j2 in range(j0, n2+j0):
            sinj2t = vsinj2t[j2-j0]
            cosj2t = vcosj2t[j2-j0]
            for i2 in range(i0, m2+i0):
                sini2x = vsini2x[i2-i0]
                cosi2x = vcosi2x[i2-i0]
                col = (i2-i0)*num2 + (j2-j0)*num2*m2 + num0 + num1*m1
                fint[col+0] = beta*(fint[col+0]) + alpha*((L*sini2x*(NttL*sina*sinj2t + NxtL*cosj2t*j2) + pi*NxxL*cosi2x*i2*r*sinj2t)/L)
                fint[col+1] = beta*(fint[col+1]) + alpha*(-NxtL*j2*sini2x*sinj2t + cosj2t*(NttL*sina*sini2x + pi*NxxL*cosi2x*i2*r/L))
                fint[col+2] = beta*(fint[col+2]) + alpha*((L*sini2x*(NttL*cosj2t*j2 - NxtL*sina*sinj2t) + pi*NxtL*cosi2x*i2*r*sinj2t)/L)
                fint[col+3] = beta*(fint[col+3]) + alpha*( -NttL*j2*sini2x*sinj2t + NxtL*cosj2t*(-L*sina*sini2x + pi*cosi2x*i2*r)/L)
                fint[col+4] = beta*(fint[col+4]) + alpha*((L*cosj2t*j2*(L*sini2x*(MxtL*sina + Ntt0*(w0t + wt) + NttL*w0t + NttL*wt + Nxt0*r*w0x + Nxt0*r*wx + NxtL*r*w0x + NxtL*r*wx) - 2*pi*MxtL*cosi2x*i2*r) + sinj2t*(pi*L*cosi2x*i2*r*(-MttL*sina + Nxt0*(w0t + wt) + NxtL*w0t + NxtL*wt + Nxx0*r*w0x + Nxx0*r*wx + NxxL*r*w0x + NxxL*r*wx) + sini2x*((L*L)*MttL*(j2*j2) + r*((L*L)*NttL*cosa + (pi*pi)*MxxL*(i2*i2)*r))))/((L*L)*r))
                fint[col+5] = beta*(fint[col+5]) + alpha*((-L*j2*sinj2t*(L*sini2x*(MxtL*sina + Ntt0*(w0t + wt) + NttL*w0t + NttL*wt + Nxt0*r*w0x + Nxt0*r*wx + NxtL*r*w0x + NxtL*r*wx) - 2*pi*MxtL*cosi2x*i2*r) + cosj2t*(pi*L*cosi2x*i2*r*(-MttL*sina + Nxt0*(w0t + wt) + NxtL*w0t + NxtL*wt + Nxx0*r*w0x + Nxx0*r*wx + NxxL*r*w0x + NxxL*r*wx) + sini2x*((L*L)*MttL*(j2*j2) + r*((L*L)*NttL*cosa + (pi*pi)*MxxL*(i2*i2)*r))))/((L*L)*r))

    free(w0xs)
    free(w0ts)
    free(vsini1x)
    free(vcosi1x)
    free(vsini2x)
    free(vcosi2x)
    free(vsinj2t)
    free(vcosj2t)
