Geografic-Information-System/imagery/i.atcorr/gauss.cpp

#include <cmath>
#include "common.h"
#include "Gauss.h"

const long int mu2	= 48;

float Gauss::angmu[10] = { 85.f, 80.f, 70.f, 60.f, 50.f, 40.f, 30.f, 20.f, 10.f, 0.f };
float Gauss::angphi[13] = { 0.f, 30.f, 60.f, 90.f, 120.f, 150.f, 180.f, 210.f, 240.f, 270.f, 300.f, 330.f, 360.f };

/*  preliminary computations for gauss integration */
void Gauss::init()
{
    int j;   

    /* convert angphi and angmu to radians */
    for (j = 0; j < 13; ++j) angphi[j] = (float)(angphi[j] * M_PI / 180.f);
    for (j = 0; j < 10; ++j) angmu[j] =	 (float)cos(angmu[j] * M_PI / 180.f);

    /* calculate rm & gb */
	
    float anglem[mu2];
    float weightm[mu2];
    gauss (-1.f, 1.f, anglem, weightm, mu2);

    gb[STDI(-mu)]   = 0;
    gb[STDI(0)]     = 0;
    gb[STDI(mu)]    = 0;
    rm[STDI(-mu)]   = 0;
    rm[STDI(0)]     = 0;
    rm[STDI(mu)]    = 0;
    /* do shift into rm & gb */
    for (j = -mu+1; j <= -1; ++j)
    {
	rm[-j] = anglem[mu + j - 1];
	gb[-j] = weightm[mu + j - 1];
    }

    for (j = 1; j <= mu-1; ++j)
    {
	rm[2*mu - j] = anglem[mu + j - 2];
	gb[2*mu - j] = weightm[mu + j - 2];
    }

    /* calculate rp & gp */
    gauss (0.f, (float)2 * M_PI, rp, gp, np);
}


/*	Compute for a given n, the gaussian quadrature (the n gaussian angles and the
  their respective weights). The gaussian quadrature is used in numerical integration involving the
  cosine of emergent or incident direction zenith angle. */
void Gauss::gauss (float a, float b, float *x, float *w, long int n)
{
    int m = (n + 1) / 2;
    double xm = (b + a) / 2;
    double xl = (b - a) / 2;

    for(int i = 0; i < m; i++)
    {
	double 
	    z1, 
	    pp, 
	    z = cos(M_PI * (i + 0.75) / (n + 0.5));

	do {
	    double p1 = 1;
	    double p2 = 0;

	    for(int j = 0; j < n; j++)
	    {
		double p3 = p2;
		p2 = p1;
		p1 = ((2 * j + 1) * z * p2 - j * p3) / (j+1);
	    }

	    pp = n * (z * p1 - p2) / (z * z - 1);
	    z1 = z;
	    z = z1 - p1 / pp;
	} while(fabs(z - z1) > 3e-14);

	if (fabs(z) < 3e-14) z = 0;
	x[i] = (float)(xm - xl * z);
        x[n - 1 - i] = (float)(xm + xl * z);
	w[i] = (float)(2 * xl / ((1 - z * z) * pp * pp));
	w[n - 1 - i] = w[i];
    }
}