// Set OVERLAPPING to 1 in order to fill in the edges of the box, with
// c endcaps and b overlapping a.  With the proper choice of parameters,
// (setting rim slds to sld, core sld to solvent, rim thickness to thickness
// and subtracting 2*thickness from length, this should match the hollow
// rectangular prism.)  Set it to 0 for the documented behaviour.
#define OVERLAPPING 0
static double
form_volume(double length_a, double length_b, double length_c,
    double thick_rim_a, double thick_rim_b, double thick_rim_c)
{
    return
#if OVERLAPPING
        // Hollow rectangular prism only includes the volume of the shell
        // so uncomment the next line when comparing.  Solid rectangular
        // prism, or parallelepiped want filled cores, so comment when
        // comparing.
        //-length_a * length_b * length_c +
        (length_a + 2.0*thick_rim_a) *
        (length_b + 2.0*thick_rim_b) *
        (length_c + 2.0*thick_rim_c);
#else
        length_a * length_b * length_c +
        2.0 * thick_rim_a * length_b * length_c +
        2.0 * length_a * thick_rim_b * length_c +
        2.0 * length_a * length_b * thick_rim_c;
#endif
}

static double
Iq(double q,
    double core_sld,
    double arim_sld,
    double brim_sld,
    double crim_sld,
    double solvent_sld,
    double length_a,
    double length_b,
    double length_c,
    double thick_rim_a,
    double thick_rim_b,
    double thick_rim_c)
{
    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)

    const double mu = 0.5 * q * length_b;

    // Scale sides by B
    const double a_over_b = length_a / length_b;
    const double c_over_b = length_c / length_b;

    double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b;
    double tB_over_b = 1+ 2.0*thick_rim_b/length_b;
    double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b;

    double Vin = length_a * length_b * length_c;
#if OVERLAPPING
    const double capA_area = length_b*length_c;
    const double capB_area = (length_a+2.*thick_rim_a)*length_c;
    const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b);
#else
    const double capA_area = length_b*length_c;
    const double capB_area = length_a*length_c;
    const double capC_area = length_a*length_b;
#endif
    const double Va = length_a * capA_area;
    const double Vb = length_b * capB_area;
    const double Vc = length_c * capC_area;
    const double Vat = Va + 2.0 * thick_rim_a * capA_area;
    const double Vbt = Vb + 2.0 * thick_rim_b * capB_area;
    const double Vct = Vc + 2.0 * thick_rim_c * capC_area;

    // Scale factors (note that drC is not used later)
    const double dr0 = (core_sld-solvent_sld);
    const double drA = (arim_sld-solvent_sld);
    const double drB = (brim_sld-solvent_sld);
    const double drC = (crim_sld-solvent_sld);

    // outer integral (with gauss points), integration limits = 0, 1
    double outer_sum = 0; //initialize integral
    for( int i=0; i<GAUSS_N; i++) {
        double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
        double mu_proj = mu * sqrt(1.0-sigma*sigma);

        // inner integral (with gauss points), integration limits = 0, pi/2
        const double siC = sas_sinx_x(mu * sigma * c_over_b);
        const double siCt = sas_sinx_x(mu * sigma * tC_over_b);
        double inner_sum = 0.0;
        for(int j=0; j<GAUSS_N; j++) {
            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
            double sin_uu, cos_uu;
            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
            const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b);
            const double siB = sas_sinx_x(mu_proj * cos_uu );
            const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b);
            const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b);

#if OVERLAPPING
            const double f = dr0*Vin*siA*siB*siC
                + drA*(Vat*siAt-Va*siA)*siB*siC
                + drB*siAt*(Vbt*siBt-Vb*siB)*siC
                + drC*siAt*siBt*(Vct*siCt-Vc*siC);
#else
            const double f = dr0*Vin*siA*siB*siC
                + drA*(Vat*siAt-Va*siA)*siB*siC
                + drB*siA*(Vbt*siBt-Vb*siB)*siC
                + drC*siA*siB*(Vct*siCt-Vc*siC);
#endif

            inner_sum += GAUSS_W[j] * f * f;
        }
        inner_sum *= 0.5;
        // now sum up the outer integral
        outer_sum += GAUSS_W[i] * inner_sum;
    }
    outer_sum *= 0.5;

    //convert from [1e-12 A-1] to [cm-1]
    return 1.0e-4 * outer_sum;
}

static double
Iqxy(double qa, double qb, double qc,
    double core_sld,
    double arim_sld,
    double brim_sld,
    double crim_sld,
    double solvent_sld,
    double length_a,
    double length_b,
    double length_c,
    double thick_rim_a,
    double thick_rim_b,
    double thick_rim_c)
{
    // cspkernel in csparallelepiped recoded here
    const double dr0 = core_sld-solvent_sld;
    const double drA = arim_sld-solvent_sld;
    const double drB = brim_sld-solvent_sld;
    const double drC = crim_sld-solvent_sld;

    double Vin = length_a * length_b * length_c;
#if OVERLAPPING
    const double capA_area = length_b*length_c;
    const double capB_area = (length_a+2.*thick_rim_a)*length_c;
    const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b);
#else
    const double capA_area = length_b*length_c;
    const double capB_area = length_a*length_c;
    const double capC_area = length_a*length_b;
#endif
    const double Va = length_a * capA_area;
    const double Vb = length_b * capB_area;
    const double Vc = length_c * capC_area;
    const double Vat = Va + 2.0 * thick_rim_a * capA_area;
    const double Vbt = Vb + 2.0 * thick_rim_b * capB_area;
    const double Vct = Vc + 2.0 * thick_rim_c * capC_area;

    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
    // the scaling by B.
    const double tA = length_a + 2.0*thick_rim_a;
    const double tB = length_b + 2.0*thick_rim_b;
    const double tC = length_c + 2.0*thick_rim_c;
    const double siA = sas_sinx_x(0.5*length_a*qa);
    const double siB = sas_sinx_x(0.5*length_b*qb);
    const double siC = sas_sinx_x(0.5*length_c*qc);
    const double siAt = sas_sinx_x(0.5*tA*qa);
    const double siBt = sas_sinx_x(0.5*tB*qb);
    const double siCt = sas_sinx_x(0.5*tC*qc);

#if OVERLAPPING
    const double f = dr0*Vin*siA*siB*siC
        + drA*(Vat*siAt-Va*siA)*siB*siC
        + drB*siAt*(Vbt*siBt-Vb*siB)*siC
        + drC*siAt*siBt*(Vct*siCt-Vc*siC);
#else
    const double f = dr0*Vin*siA*siB*siC
        + drA*(Vat*siAt-Va*siA)*siB*siC
        + drB*siA*(Vbt*siBt-Vb*siB)*siC
        + drC*siA*siB*(Vct*siCt-Vc*siC);
#endif

    return 1.0e-4 * f * f;
}