#include <math.h>

double pi = 3.141592653589793;

double GAUSS_N = 76;

double GAUSS_Z[] = {-0.99950595,-0.99739779,-0.99360877,-0.98814445,-0.98101394,-0.97222923
,-0.96180513,-0.94975921,-0.93611178,-0.92088586,-0.90410712,-0.88580385
,-0.86600691,-0.84474969,-0.82206804,-0.79800019,-0.77258673,-0.74587051
,-0.71789659,-0.68871214,-0.65836635,-0.62691042,-0.59439737,-0.56088203
,-0.52642092,-0.49107214,-0.45489531,-0.41795142,-0.38030277,-0.34201284
,-0.3031462,-0.26376839,-0.2239458,-0.18374559,-0.14323555,-0.10248398
,-0.06155959,-0.0205314,0.0205314,0.06155959,0.10248398,0.14323555
,0.18374559,0.2239458,0.26376839,0.3031462,0.34201284,0.38030277
,0.41795142,0.45489531,0.49107214,0.52642092,0.56088203,0.59439737
,0.62691042,0.65836635,0.68871214,0.71789659,0.74587051,0.77258673
,0.79800019,0.82206804,0.84474969,0.86600691,0.88580385,0.90410712
,0.92088586,0.93611178,0.94975921,0.96180513,0.97222923,0.98101394
,0.98814445,0.99360877,0.99739779,0.99950595};

double GAUSS_W[] = {0.00126779,0.0029491,0.00462794,0.00629918,0.00795985,0.00960711
,0.01123817,0.01285028,0.01444073,0.01600683,0.01754594,0.01905546
,0.02053285,0.02197561,0.02338133,0.02474761,0.02607216,0.02735276
,0.02858722,0.02977349,0.03090955,0.03199348,0.03302347,0.03399778
,0.03491476,0.03577286,0.03657064,0.03730676,0.03797996,0.03858913
,0.03913322,0.03961133,0.04002265,0.04036647,0.04064223,0.04084946
,0.04098781,0.04105704,0.04105704,0.04098781,0.04084946,0.04064223
,0.04036647,0.04002265,0.03961133,0.03913322,0.03858913,0.03797996
,0.03730676,0.03657064,0.03577286,0.03491476,0.03399778,0.03302347
,0.03199348,0.03090955,0.02977349,0.02858722,0.02735276,0.02607216
,0.02474761,0.02338133,0.02197561,0.02053285,0.01905546,0.01754594
,0.01600683,0.01444073,0.01285028,0.01123817,0.00960711,0.00795985
,0.00629918,0.00462794,0.0029491,0.00126779};

double fq (double qab, double qc, double radius, double length)
{
    return (sas_2J1x_x(qab * radius) * sas_sinx_x(qc * 0.5 * length));
}


double form_volume (double radius, double length)
{
    return (pi * square (radius)  * length);
}


double orient_avg_1D (double q, double radius, double length)
{
    double zm,zb,total,theta,sin_theta,cos_theta,form;
    int n,i;
    double *GAUSS_Z, *GAUSS_W;

    zm = pi / 4;
    zb = pi / 4;
    total = 0.0;
    for (n=0 ; n < len(GAUSS_N) ; n++) {
        theta = GAUSS_Z[i] * zm + zb;
        sin_theta = sin(theta);
        cos_theta = cos(theta);
        form = fq(q * sin_theta, q * cos_theta, radius, length);
        total += GAUSS_W[i] * form * form * sin_theta;
    }
    return (total * zm);
}


double Iqxy (double qab, double qc, double sld, double solvent_sld, double radius, double length)
{
    double s,form;

    s = sld - solvent_sld * form_volume(radius, length);
    form = fq(qab, qc, radius, length);
    return (0.0001 * square(s * form));
}


double Iq (double q, double sld, double solvent_sld, double radius, double length)
{
    double s;

    s = sld - solvent_sld * form_volume(radius, length);
    return (0.0001 * s * s * orient_avg_1D(q, radius, length));
}