text:
#include "tmcoct.h"
#include <gsl/gsl_rng.h>

const double G_COS_0_D = 1.0 - 1.0E-14;
const double G_COS_90_D = 1.0E-7;

#define STANDARDTEST 0
/* testing program using fixed rnd seed. */

#define PARTIALREFLECTION 1
/* 1=split photon, 0=statistical reflection. */

#define COSZERO (1.0-1.0E-14)
/* cosine of about 1e-6 rad. */

#define COS90D  1.0E-6
/* cosine of about 1.57 - 1e-6 rad. */

const gsl_rng_type *T_RNG_GSL;
gsl_rng *r_RNG_GSL;

/***********************************************************
 *	A random number generator from Numerical Recipes in C.
 ****/
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC 1.0E-9

#undef MBIG
#undef MSEED
#undef MZ
#undef FAC

/***********************************************************
 *	Generate a random number between 0 and 1.  Take a
 *	number as seed the first time entering the function.
 *	The seed is limited to 1<<15.
 *	We found that when idum is too large, ran3 may return
 *	numbers beyond 0 and 1.
 ****/
double RandomNum(void) {
    static Boolean first_time = 1;
    static unsigned long int idum;    /* seed for ran3. */

    if (first_time) {
        T_RNG_GSL = gsl_rng_default;
        r_RNG_GSL = gsl_rng_alloc(T_RNG_GSL);

#if STANDARDTEST /* Use fixed seed to test the program. */
        printf("*** Using seed to test. ***\n");
    printf("*** Remove the fixed seed after the test. ***\n");
    idum =  7923;

#else
        idum = (unsigned long int) time(NULL) % (1 << 15);
        /* use 16-bit integer as the seed. */
#endif
        printf("Random seed = %lu \n", idum);

        //ran3(&idum);
        gsl_rng_set(r_RNG_GSL, idum);

        first_time = 0;
        //idum = 1;
    }
    // TODO FIXME remove 1.0 -
    return (1.0 - gsl_rng_uniform(r_RNG_GSL));
}

/***********************************************************
 *	Compute the specular reflection.
 *
 *	If the first layer is a turbid medium, use the Fresnel
 *	reflection from the boundary of the first layer as the
 *	specular reflectance.
 *
 *	If the first layer is glass, multiple reflections in
 *	the first layer is considered to get the specular
 *	reflectance.
 *
 *	The subroutine assumes the Layerspecs array is correctly
 *	initialized.
 ****/
double Rspecular(double ni, double nt) {
    double r1;
    /* direct reflections from the 1st and 2nd layers. */
    double temp;
    temp = (ni - nt) / (ni + nt);
    r1 = temp * temp;
    return (r1);
}

/***********************************************************
 *	Initialize a photon packet.
 ****/
void LaunchPhoton(SimulationStruct *in_parm, double Rspecular, TetrahedronPhotonStruct *Photon_Ptr,
                  TetrahedronPhotonStruct *PhotonCont_Ptr, TetrahedronStruct *TetrahedronRoot) {
    if (Photon_Ptr->FstBackReflectionFlag) {
        //Most often case for biased simulations
        //The photon will be re-initialized from where the bias was applied
        long double LikelihoodRatioTmp = Photon_Ptr->LikelihoodRatioAfterFstBias;
        CopyPhotonStruct(PhotonCont_Ptr, Photon_Ptr);

        //Apply one unbiased random scatter to the split photon
        Spin(regionspecs[Photon_Ptr->tetrahedron->region].g, Photon_Ptr);

        if (LikelihoodRatioTmp < 1)
            // Adjust the likelihoodRation to make the simulation unbiased
            Photon_Ptr->LikelihoodRatio = 1 - LikelihoodRatioTmp;
        else
            Photon_Ptr->LikelihoodRatio = 1;

    } else {

        Photon_Ptr->w = 1.0 - Rspecular;

        Photon_Ptr->x = in_parm->probe_x;
        Photon_Ptr->y = in_parm->probe_y;
        Photon_Ptr->z = in_parm->probe_z;

        Photon_Ptr->s = 0.0;
        Photon_Ptr->sleft = 0.0;

        Photon_Ptr->ux = Photon_Ptr->uy = 0.0;
        Photon_Ptr->uz = 1.0;
        Photon_Ptr->dead = 0.0;

        Photon_Ptr->OpticalPath = 0.0;
        Photon_Ptr->MaxDepth = 0.0;
        Photon_Ptr->LikelihoodRatio = 1.0;
        Photon_Ptr->LikelihoodRatioAfterFstBias = 1.0;

        Photon_Ptr->FstBackReflectionFlag = 0;
        Photon_Ptr->LocationFstBias = -1.0;
        Photon_Ptr->NumBackwardsSpecularReflections = 0;

        Photon_Ptr->LikelihoodRatioIncreaseFstBias = 0;

        Photon_Ptr->LikelihoodRatioIncreaseFstBias = 0;

        Photon_Ptr->NextTetrahedron = -1;

        Photon_Ptr->tetrahedron = TetrahedronRoot;
    }
}

/***********************************************************
 *	Choose (sample) a new theta angle for photon propagation
 *	according to the anisotropy.
 *
 *	If anisotropy g is 0, then
 *		cos(theta) = 2*rand-1.
 *	otherwise
 *		sample according to the Henyey-Greenstein function.
 *
 *	Returns the cosine of the polar deflection angle theta.
 ****/
double SpinTheta(double g) {
    double cost;

    if (g == 0.0)
        cost = 2.0 * RandomNum() - 1.0;
    else {
        double temp = (1.0 - g * g) / (1.0 - g + 2.0 * g * RandomNum());
        cost = (1.0 + g * g - temp * temp) / (2.0 * g);
        if (cost < -1.0) cost = -1.0;
        else if (cost > 1.0) cost = 1.0;
    }
    return (cost);
}

double SpinThetaForwardFstBias(double g) {
    double cost;

    if (g == 0.0)
        cost = RandomNum();
    else {
        double RandomNumTmp = RandomNum();
        double temp = RandomNumTmp / (1 - g) + (1 - RandomNumTmp) / sqrt(g * g + 1);
        cost = (g * g + 1 - 1. / (temp * temp)) / (2 * g);
        if (cost < -1) cost = -1;
        else if (cost > 1) cost = 1;
    }
    return (cost);
}

/***********************************************************
 *	Choose a new direction for photon propagation by
 *	sampling the polar deflection angle theta and the
 *	azimuthal angle psi.
 *
 *	Note:
 *  	theta: 0 - pi so sin(theta) is always positive
 *  	feel free to use sqrt() for cos(theta).
 *
 *  	psi:   0 - 2pi
 *  	for 0-pi  sin(psi) is +
 *  	for pi-2pi sin(psi) is -
 ****/

void Spin(double g,
          TetrahedronPhotonStruct *Photon_Ptr) {
    double cost, sint;    /* cosine and sine of the */
    /* polar deflection angle theta. */
    double cosp, sinp;    /* cosine and sine of the */
    /* azimuthal angle psi. */
    double ux = Photon_Ptr->ux;
    double uy = Photon_Ptr->uy;
    double uz = Photon_Ptr->uz;
    double psi;

    cost = SpinTheta(g);
    sint = sqrt(1.0 - cost * cost);
    /* sqrt() is faster than sin(). */

    psi = 2.0 * PI * RandomNum(); /* spin psi 0-2pi. */
    cosp = cos(psi);
    if (psi < PI)
        sinp = sqrt(1.0 - cosp * cosp);
        /* sqrt() is faster than sin(). */
    else
        sinp = -sqrt(1.0 - cosp * cosp);

    if (fabs(uz) > COSZERO) {    /* normal incident. */
        Photon_Ptr->ux = sint * cosp;
        Photon_Ptr->uy = sint * sinp;
        Photon_Ptr->uz = cost * SIGN(uz);
        /* SIGN() is faster than division. */
    } else {
        /* regular incident. */
        double temp = sqrt(1.0 - uz * uz);
        Photon_Ptr->ux = sint * (ux * uz * cosp - uy * sinp)
                         / temp + ux * cost;
        Photon_Ptr->uy = sint * (uy * uz * cosp + ux * sinp)
                         / temp + uy * cost;
        Photon_Ptr->uz = -sint * cosp * temp + uz * cost;
    }

}

void SpinBias37(double g,
                TetrahedronPhotonStruct *Photon_Ptr,
                TetrahedronPhotonStruct *PhotonCont_Ptr,
                SimulationStruct *In_Ptr)
{
    double g_squared = sq(g);

    double cost, sint;    /* cosine and sine of the */
    /* polar deflection angle theta. */
    double cosp, sinp;    /* cosine and sine of the */
    /* azimuthal angle psi. */
    double ux = Photon_Ptr->ux;
    double uy = Photon_Ptr->uy;
    double uz = Photon_Ptr->uz;
    double ux_Orig = ux;
    double uy_Orig = uy;
    double uz_Orig = uz;
    double psi;

    double BackwardBiasCoefficient = In_Ptr->BackwardBiasCoefficient;
    double costg;
    double costg1, costg2;

    double BiasCoefficientTmp = 0;

    short int ReachedTargetOpticalDepthFlag = 0;
    short int ThisIsFirstBackwardBiasFlag = 0;

    if (Photon_Ptr->z > In_Ptr->TargetDepthMin
        && Photon_Ptr->z < In_Ptr->TargetDepthMax
        && Photon_Ptr->uz > 0
        && !Photon_Ptr->FstBackReflectionFlag) {
        // The bias backwards will be applied only if the photon is going forward
        // The status of the photon prior to the bias will be saved

        CopyPhotonStruct(Photon_Ptr, PhotonCont_Ptr);
        ReachedTargetOpticalDepthFlag = 1;
        Photon_Ptr->FstBackReflectionFlag = 1; // Bias backwards only once
        Photon_Ptr->LocationFstBias = Photon_Ptr->z;
        BiasCoefficientTmp = BackwardBiasCoefficient;
        ThisIsFirstBackwardBiasFlag = 1;
    }

    /**********************************
	 ** Biased Direction towards probe
	 **********************************/
    double vx = In_Ptr->probe_x-Photon_Ptr->x;
    double vy = In_Ptr->probe_y-Photon_Ptr->y;
    double vz = In_Ptr->probe_z-Photon_Ptr->z;
    double LengthVector = sqrt(sq(vx) + sq(vy) + sq(vz));
    vx /= LengthVector;
    vy /= LengthVector;
    vz /= LengthVector;
    /*********************************/

    if ((Photon_Ptr->FstBackReflectionFlag
         || Photon_Ptr->NumBackwardsSpecularReflections > 0)
        && !ReachedTargetOpticalDepthFlag)
    {
        // It was biased at least once before, and is moving backwards
        ReachedTargetOpticalDepthFlag = 2;

        double NextStepSize = -log(In_Ptr->rndStepSizeInTissue) / regionspecs[Photon_Ptr->tetrahedron->region].muas;
        double CurrentDistanceToOrigin = sqrt(sq(Photon_Ptr->x-In_Ptr->probe_x) + sq(Photon_Ptr->y-In_Ptr->probe_y) + sq(Photon_Ptr->z-In_Ptr->probe_z));

        if (NextStepSize >= CurrentDistanceToOrigin
            && acos(-vz) <= In_Ptr->MaxCollectingAngleDeg * PI / 180) {
            ReachedTargetOpticalDepthFlag = 1;
        }

        BiasCoefficientTmp = BackwardBiasCoefficient;

    }

    int BiasFunctionRandomlySelected = 0;
    if (ReachedTargetOpticalDepthFlag) {
        // Photon reached target optical layer or may undergo an additional biased
        // scattering or unbiased scattering

        // BiasFunctionRandomlySelected=1 means use biased scattering and 2 means unbiased scattering
        if (RandomNum() <= In_Ptr->ProbabilityAdditionalBias)
            BiasFunctionRandomlySelected = 1;
        else
            BiasFunctionRandomlySelected = 2;

        if (ReachedTargetOpticalDepthFlag == 1
            || BiasFunctionRandomlySelected == 1) {
            /*************************************************************************
			** The photon is within the target depth and going forward
			** The additional biased scattering is randomly chosen
			** So the scattering is biased Henyey-Greenstein scattering
			*************************************************************************/
            cost = SpinThetaForwardFstBias(BiasCoefficientTmp);
            ux = vx;
            uy = vy;
            uz = vz;
        } else {
            /**************************************************************************************
			** The photon is within the target depth but the scattering is randomly selected is
			** unbiased scattering
			** or the photon is already going backward or it is out of target depth
			**************************************************************************************/
            cost = SpinTheta(g);
        }
    } else {
        /**************************************************************************
		**  The photon is not within the target depth or it is not going forward
		**  so do unbiased scattering
		**************************************************************************/
        cost = SpinTheta(g);
    }

    if (cost < -1)
        cost = -1;
    else if (cost > 1)
        cost = 1;

    sint = sqrt(1.0 - cost * cost);
    /* sqrt() is faster than sin(). */

    psi = 2.0 * PI * RandomNum(); /* spin psi 0-2pi. */
    cosp = cos(psi);
    if (psi < PI)
        sinp = sqrt(1.0 - cosp * cosp);
        /* sqrt() is faster than sin(). */
    else
        sinp = -sqrt(1.0 - cosp * cosp);

    if (fabs(uz) > COSZERO) {    /* normal incident. */
        Photon_Ptr->ux = sint * cosp;
        Photon_Ptr->uy = sint * sinp;
        Photon_Ptr->uz = cost * SIGN(uz);
        /* SIGN() is faster than division. */
    } else {        /* regular incident. */
        double temp = sqrt(1.0 - uz * uz);
        Photon_Ptr->ux = sint * (ux * uz * cosp - uy * sinp)
                         / temp + ux * cost;
        Photon_Ptr->uy = sint * (uy * uz * cosp + ux * sinp)
                         / temp + uy * cost;
        Photon_Ptr->uz = -sint * cosp * temp + uz * cost;
    }

    costg = ux_Orig * Photon_Ptr->ux
            + uy_Orig * Photon_Ptr->uy
            + uz_Orig * Photon_Ptr->uz;
    costg2 = costg;
    costg1 = vx * Photon_Ptr->ux
             + vy * Photon_Ptr->uy
             + vz * Photon_Ptr->uz;

    if (BiasCoefficientTmp) {
        double one_plus_a_squared = 1 + sq(BiasCoefficientTmp);
        double sqrt_one_plus_a_squared = sqrt(one_plus_a_squared);
        double LikelihoodRatioIncreaseFactor;
        if (ReachedTargetOpticalDepthFlag == 1)
            /****************************************************************************************
			 ** Likelihood for the first Bias scattering. Equation (8) of the paper:
			 ** Malektaji, Siavash, Ivan T. Lima, and Sherif S. Sherif. "Monte Carlo simulation of
			 ** optical coherence tomography for turbid media with arbitrary spatial distributions."
			 ** Journal of biomedical optics 19.4 (2014): 046001-046001.
			 ****************************************************************************************/
            LikelihoodRatioIncreaseFactor = ((1.0 - g_squared) * (sqrt_one_plus_a_squared - 1.0 + BiasCoefficientTmp) *
                                             sqrt(cube(one_plus_a_squared - 2.0 * BiasCoefficientTmp * cost))) /
                                            (2.0 * BiasCoefficientTmp * (1.0 - BiasCoefficientTmp) *
                                             sqrt_one_plus_a_squared * sqrt(cube(1.0 + g_squared - 2.0 * g * costg)));
        else {
            double cost1, cost2;
            if (BiasFunctionRandomlySelected == 1) {
                cost1 = cost;
                cost2 = costg2;
            } else {
                cost1 = costg1;
                cost2 = cost;
            }
            /*******************************************************************************************************
			 **  The likelihood ratio of additional biased scatterings, whether the biased or the unbiased
			 **  probability density function is randomly selected, is calculated according to the equation (9)
			 **  of the paper:
			 **  Malektaji, Siavash, Ivan T. Lima, and Sherif S. Sherif. "Monte Carlo simulation of
			 **  optical coherence tomography for turbid media with arbitrary spatial distributions."
			 **  Journal of biomedical optics 19.4 (2014): 046001-046001.
			 *******************************************************************************************************/
            double pdf1 = (sqrt_one_plus_a_squared * BiasCoefficientTmp * (1.0 - BiasCoefficientTmp)) /
                          ((sqrt_one_plus_a_squared - 1.0 + BiasCoefficientTmp) *
                           sqrt(cube(one_plus_a_squared - 2.0 * BiasCoefficientTmp * cost1)));

            double pdf2 = (1.0 - g_squared) / (2.0 * cube(sqrt(1.0 + g_squared - 2.0 * g * cost2)));

            LikelihoodRatioIncreaseFactor = pdf2 /
                                            (In_Ptr->ProbabilityAdditionalBias * pdf1 +
                                             (1.0 - In_Ptr->ProbabilityAdditionalBias) * pdf2);

        }

        Photon_Ptr->LikelihoodRatio *= LikelihoodRatioIncreaseFactor;
        if (ThisIsFirstBackwardBiasFlag == 1) {
            // In case there was a sure backward bias and that was the very first one
            Photon_Ptr->LikelihoodRatioAfterFstBias = Photon_Ptr->LikelihoodRatio;
            Photon_Ptr->LikelihoodRatioIncreaseFstBias = LikelihoodRatioIncreaseFactor;
        }
    }
}

/***********************************************************
 *	Move the photon s away in the current layer of medium.
 ****/
void Hop(TetrahedronPhotonStruct *Photon_Ptr) {
    double s = Photon_Ptr->s;
    Photon_Ptr->x += s * Photon_Ptr->ux;
    Photon_Ptr->y += s * Photon_Ptr->uy;
    Photon_Ptr->z += s * Photon_Ptr->uz;

    Photon_Ptr->OpticalPath += s;
    if (Photon_Ptr->MaxDepth < Photon_Ptr->z)
        Photon_Ptr->MaxDepth = Photon_Ptr->z;
}
/***********************************************************
 *	If uz != 0, return the photon step size in glass,
 *	Otherwise, return 0.
 *
 *	The step size is the distance between the current
 *	position and the boundary in the photon direction.
 *
 *	Make sure uz !=0 before calling this function.
 ****/

/***********************************************************
 *	Pick a step size for a photon packet when it is in
 *	tissue.
 *	If the member sleft is zero, make a new step size
 *	with: -log(rnd)/(mua+mus).
 *	Otherwise, pick up the leftover in sleft.
 *
 *	Layer is the index to layer.
 *	In_Ptr is the input parameters.
 ****/
void StepSizeInTissue(TetrahedronPhotonStruct *Photon_Ptr) {

    if (Photon_Ptr->sleft == 0.0) {  /* make a new step. */
        double rnd;

        do rnd = RandomNum();
        while (rnd <= 0.0);    /* avoid zero. */

        Photon_Ptr->s = -log(rnd) * regionspecs[Photon_Ptr->tetrahedron->region].rmuas;
    } else {    /* take the leftover. */
        Photon_Ptr->s = Photon_Ptr->sleft * regionspecs[Photon_Ptr->tetrahedron->region].rmuas;
        Photon_Ptr->sleft = 0.0;
    }
}

/***********************************************************
 *	Check if the step will hit the boundary.
 *	Return 1 if hit boundary.
 *	Return 0 otherwise.
 *
 * 	If the projected step hits the boundary, the members
 *	s and sleft of Photon_Ptr are updated.
 ****/
Boolean HitBoundary(TetrahedronPhotonStruct *Photon_Ptr) {
    TetrahedronStruct *tetrahedron;
    tetrahedron = Photon_Ptr->tetrahedron;

    double min_distance = 1e10;
    int index_of_tetrahedron_with_min_distance = -1;

    double cos_normals_and_photon_direction[4];
    double distance_from_face_in_photon_direction;
    double perpendicular_distance;

    for (int i = 0; i < 4; i++) {
        cos_normals_and_photon_direction[i] = (Photon_Ptr->tetrahedron->Faces[i].Nx) * Photon_Ptr->ux +
                      (Photon_Ptr->tetrahedron->Faces[i].Ny) * Photon_Ptr->uy +
                      (Photon_Ptr->tetrahedron->Faces[i].Nz) * Photon_Ptr->uz;
    }

    for (int i = 0; i < 4; i++) {
        if (cos_normals_and_photon_direction[i] < 0) {
            perpendicular_distance = ((Photon_Ptr->tetrahedron->Faces[i].Nx) * Photon_Ptr->x +
                      (Photon_Ptr->tetrahedron->Faces[i].Ny) * Photon_Ptr->y +
                      (Photon_Ptr->tetrahedron->Faces[i].Nz) * Photon_Ptr->z + Photon_Ptr->tetrahedron->Faces[i].d);
            distance_from_face_in_photon_direction = -perpendicular_distance / cos_normals_and_photon_direction[i];

            if (distance_from_face_in_photon_direction < min_distance) {
                min_distance = distance_from_face_in_photon_direction;
                index_of_tetrahedron_with_min_distance = i;
                Photon_Ptr->MinCos = cos_normals_and_photon_direction[i];
            }
        }
    }

    Boolean hit_boundary;
    if (Photon_Ptr->s > min_distance) {
        /* not horizontal & crossing. */

        Photon_Ptr->sleft = (Photon_Ptr->s - min_distance) * regionspecs[tetrahedron->region].muas;
        Photon_Ptr->s = min_distance;
        hit_boundary = 1;
        Photon_Ptr->NextTetrahedron = index_of_tetrahedron_with_min_distance;

    } else {

        hit_boundary = 0;
        Photon_Ptr->NextTetrahedron = -1;

    }
    return hit_boundary;
}

/***********************************************************
 *	Drop photon weight inside the tissue (not glass).
 *
 *  The photon is assumed not dead.
 *
 *	The weight drop is dw = w*mua/(mua+mus).
 *
 *	The dropped weight is assigned to the absorption array
 *	elements.
 ****/
void Drop(TetrahedronPhotonStruct *Photon_Ptr) {
    double dwa;    /* absorbed weight.*/
    TetrahedronStruct *tetrahedron = Photon_Ptr->tetrahedron;

    dwa = Photon_Ptr->w * regionspecs[tetrahedron->region].mua_muas;
    Photon_Ptr->w -= dwa;
}

/***********************************************************
 *	The photon weight is small, and the photon packet tries
 *	to survive a roulette.
 ****/
void Roulette(TetrahedronPhotonStruct *Photon_Ptr) {
    if (Photon_Ptr->w == 0.0)
        Photon_Ptr->dead = 1;
    else if (RandomNum() < CHANCE) /* survived the roulette.*/
        Photon_Ptr->w /= CHANCE;
    else
        Photon_Ptr->dead = 1;
}

/***********************************************************
 *	Compute the Fresnel reflectance.
 *
 *	Make sure that the cosine of the incident angle a1
 *	is positive, and the case when the angle is greater
 *	than the critical angle is ruled out.
 *
 * 	Avoid trigonometric function operations as much as
 *	possible, because they are computation-intensive.
 ****/
double
RFresnel(SimulationStruct *In_Ptr, TetrahedronPhotonStruct *Photon_Ptr, double *ux_reflected, double *uy_reflected,
         double *uz_reflected, double *ux_refracted, double *uy_refracted, double *uz_refracted) {

    double r;

    double cosa, sina, sinb, cosb;

    cosa = -(Photon_Ptr->MinCos);

    double ni = regionspecs[Photon_Ptr->tetrahedron->region].n;
    double nt;
    if (Photon_Ptr->tetrahedron->adjTetrahedrons[Photon_Ptr->NextTetrahedron] == NULL)
        nt = regionspecs[0].n;
    else {

        int next_region = Photon_Ptr->tetrahedron->adjTetrahedrons[Photon_Ptr->NextTetrahedron]->region;
        nt = regionspecs[next_region].n;
    }

    sina = sqrt(1.0 - cosa * cosa);
    sinb = ni * sina / nt;
    cosb = sqrt(1.0 - sinb * sinb);

    if (cosa > G_COS_0_D) {
        r = (ni - nt)
            / (ni + nt);
        r *= r;
        cosb = 1;
    } else if (cosa < G_COS_90_D) {
        r = 1.0;
    } else {
        sina = sqrt(1 - cosa * cosa);
        sinb = ni * sina / nt;
        if (sinb >= 1) {
            r = 1.0;
        } else {
            cosb = sqrt(1.0 - sinb * sinb);
            double cosacosb = cosa * cosb;
            double sinasinb = sina * sinb;
            double sinacosb = sina * cosb;
            double cosasinb = cosa * sinb;
            double cap = cosacosb - sinasinb;
            double cam = cosacosb + sinasinb;
            double sap = sinacosb + cosasinb;
            double sam = sinacosb - cosasinb;
            r = 0.5 * sam * sam * (cam * cam + cap * cap) / (sap * sap * cam * cam);
        }
    }

    double sbdivsa = ni / nt;
    *ux_refracted = -cosb * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Nx +
                    sbdivsa * (cosa * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Nx + Photon_Ptr->ux);
    *uy_refracted = -cosb * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Ny +
                    sbdivsa * (cosa * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Ny + Photon_Ptr->uy);
    *uz_refracted = -cosb * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Nz +
                    sbdivsa * (cosa * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Nz + Photon_Ptr->uz);

    *ux_reflected = 2 * cosa * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Nx + Photon_Ptr->ux;
    *uy_reflected = 2 * cosa * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Ny + Photon_Ptr->uy;
    *uz_reflected = 2 * cosa * Photon_Ptr->tetrahedron->Faces[Photon_Ptr->NextTetrahedron].Nz + Photon_Ptr->uz;

    return r;

}

/***********************************************************
 *	Record the photon weight exiting the first layer(uz<0),
 *	no matter whether the layer is glass or not, to the
 *	reflection array.
 *
 *	Update the photon weight as well.
 ***********************************************************/
void RecordValidationR(double Refl, SimulationStruct *In_Ptr, TetrahedronPhotonStruct *Photon_Ptr, OutStruct *Out_Ptr) {
    double x = Photon_Ptr->x;
    double y = Photon_Ptr->y;

    short SpatialFilterFlag = 0;
    if (sqrt(sq(x-In_Ptr->probe_x) + sq(y-In_Ptr->probe_y)) <
        In_Ptr->MaxCollectingRadius        /* The position of photon should be within the radius of the probe */
        && acos(-Photon_Ptr->uz) < In_Ptr->MaxCollectingAngleDeg * PI / 180) {
        if ((In_Ptr->TypeBias == 0 || In_Ptr->TypeBias == 3))
            SpatialFilterFlag = 1;

        if ((Photon_Ptr->FstBackReflectionFlag && In_Ptr->TypeBias == 37 ||
             Photon_Ptr->LocationFstBias == Photon_Ptr->MaxDepth
             && In_Ptr->TypeBias != 37 ||
             Photon_Ptr->NumBackwardsSpecularReflections > 0))
            SpatialFilterFlag = 1;
    }

    if (SpatialFilterFlag) {

        In_Ptr->NumFilteredPhotons++;
        double FilterOpticalPath = Photon_Ptr->OpticalPath;
        double FilterOpticalDepth = FilterOpticalPath / 2.0;

        unsigned int VectorPosition = (unsigned int) (
                (FilterOpticalPath / 2.0 - In_Ptr->OpticalDepthShift)
                / (In_Ptr->CoherenceLengthSource / NUM_SUBSTEPS_RESOLUTION));

        int ii;
        for (ii = VectorPosition - NUM_SUBSTEPS_RESOLUTION / 2;
             ii < VectorPosition + NUM_SUBSTEPS_RESOLUTION / 2;
             ii++) {
            if (ii >= 0 && ii < In_Ptr->NumOpticalDepthLengthSteps) {
                short FilterClassI_Flag =
                        Photon_Ptr->MaxDepth > (FilterOpticalDepth
                                                - In_Ptr->CoherenceLengthSource / 2.0);

                double tmpPhotonContribution = Photon_Ptr->w * (1.0 - Refl) * Photon_Ptr->LikelihoodRatio;
                if (FilterClassI_Flag) {
                    if (ii == VectorPosition)
                        In_Ptr->NumFilteredPhotonsClassI++;
                    Out_Ptr->ReflectanceClassI_Sum[ii] += tmpPhotonContribution;
                    if (Out_Ptr->ReflectanceClassI_Max[ii] < tmpPhotonContribution)
                        Out_Ptr->ReflectanceClassI_Max[ii] = tmpPhotonContribution;
                    Out_Ptr->ReflectanceClassI_SumSq[ii] += sq(tmpPhotonContribution);
                    Out_Ptr->NumClassI_PhotonsFilteredInRange[ii]++;

                } else {
                    if (ii == VectorPosition)
                        In_Ptr->NumFilteredPhotonsClassII++;
                    if (Out_Ptr->ReflectanceClassII_Max[ii] < tmpPhotonContribution)
                        Out_Ptr->ReflectanceClassII_Max[ii] = tmpPhotonContribution;
                    Out_Ptr->ReflectanceClassII_Sum[ii] += tmpPhotonContribution;
                    Out_Ptr->ReflectanceClassII_SumSq[ii] += sq(tmpPhotonContribution);
                    Out_Ptr->NumClassII_PhotonsFilteredInRange[ii]++;

                }
            }
        }
    }

    Photon_Ptr->w *= Refl;

}

/***********************************************************
 *	Decide whether the photon will be transmitted  or be
 *	reflected on the bottom boundary (uz>0) of the current
 *	layer.
 *
 *	If the photon is transmitted, move the photon to
 *	"layer+1". If "layer" is the last layer, record the
 *	transmitted weight as transmittance. See comments for
 *	CrossUpOrNot.
 *
 *	Update the photon parmameters.
 ****/
void CrossOrNot(SimulationStruct *In_Ptr, TetrahedronPhotonStruct *Photon_Ptr, OutStruct *Out_Ptr) {

    double uz_reflected = 0, ux_reflected = 0, uy_reflected = 0;    /* cosines of transmission alpha. */
    double uz_refracted = 0, ux_refracted = 0, uy_refracted = 0;

    double r = 0.0;                        /* reflectance */

    r = RFresnel(In_Ptr, Photon_Ptr, &ux_reflected, &uy_reflected, &uz_reflected, &ux_refracted, &uy_refracted,
                 &uz_refracted);

#if PARTIALREFLECTION
    if (Photon_Ptr->tetrahedron->adjTetrahedrons[Photon_Ptr->NextTetrahedron] == NULL && r < 1.0) { /* We use actual photon splitting if the photon is exiting the media
																							 ** We store the part the is exiting and the rest will be traced */
        Photon_Ptr->ux = ux_refracted;
        Photon_Ptr->uy = uy_refracted;
        Photon_Ptr->uz = uz_refracted;

        if (Photon_Ptr->uz < 0)
            RecordValidationR(r, In_Ptr, Photon_Ptr, Out_Ptr);
        else {
            Photon_Ptr->NumBackwardsSpecularReflections++;
            Photon_Ptr->w *= r;
        }
        Photon_Ptr->ux = ux_reflected;
        Photon_Ptr->uy = uy_reflected;
        Photon_Ptr->uz = uz_reflected;

    }
        /* If photon is not exiting the medium we use the statistical splitting */
    else if (RandomNum() > r) {/* transmitted */

        Photon_Ptr->tetrahedron = Photon_Ptr->tetrahedron->adjTetrahedrons[Photon_Ptr->NextTetrahedron];

        Photon_Ptr->ux = ux_refracted;
        Photon_Ptr->uy = uy_refracted;
        Photon_Ptr->uz = uz_refracted;

        Photon_Ptr->NextTetrahedron = -1;

    } else                        /* reflected. */
    {
        Photon_Ptr->NextTetrahedron = -1;

        if (Photon_Ptr->uz > 0)
            Photon_Ptr->NumBackwardsSpecularReflections++;

        Photon_Ptr->ux = ux_reflected;
        Photon_Ptr->uy = uy_reflected;
        Photon_Ptr->uz = uz_reflected;

    }
#else	/* use just statistical splitting */
    if(RandomNum() > r) {		/* transmitted. */

    if(Photon_Ptr->tetrahedron->adjTetrahedrons[Photon_Ptr->NextTetrahedron]==NULL) {

      Photon_Ptr->ux = ux_refracted;
      Photon_Ptr->uy = uy_refracted;
      Photon_Ptr->uz = uz_refracted;

      if (Photon_Ptr->uz<0)
            RecordValidationR(0.0, In_Ptr, Photon_Ptr, Out_Ptr);

      Photon_Ptr->dead = 1;
    }
    else {

      /* Siavash Refraction Direction Cosines */
      Photon_Ptr->ux = ux_refracted;
      Photon_Ptr->uy = uy_refracted;
      Photon_Ptr->uz = uz_refracted;
      Photon_Ptr->tetrahedron=Photon_Ptr->tetrahedron->adjTetrahedrons[Photon_Ptr->NextTetrahedron];
      Photon_Ptr->NextTetrahedron=-1;
    }
  }
  else {
    /* reflected. */

    Photon_Ptr->ux = ux_reflected;
    Photon_Ptr->uy = uy_reflected;
    Photon_Ptr->uz = uz_reflected;
    Photon_Ptr->NextTetrahedron=-1;
  }
#endif

}

/***********************************************************
 ****/
/***********************************************************
 *	Move the photon packet in glass layer.
 *	Horizontal photons are killed because they will
 *	never interact with tissue again.
 ****/
/***********************************************************
 *	Set a step size, move the photon, drop some weight,
 *	choose a new photon direction for propagation.
 *
 *	When a step size is long enough for the photon to
 *	hit an interface, this step is divided into two steps.
 *	First, move the photon to the boundary free of
 *	absorption or scattering, then decide whether the
 *	photon is reflected or transmitted.
 *	Then move the photon in the current or transmission
 *	medium with the unfinished stepsize to interaction
 *	site.  If the unfinished stepsize is still too long,
 *	repeat the above process.
 ****/
void HopDropSpinInTissue(SimulationStruct *In_Ptr,
                         TetrahedronPhotonStruct *Photon_Ptr,
                         TetrahedronPhotonStruct *PhotonCont_Ptr,
                         OutStruct *Out_Ptr) {
    StepSizeInTissue(Photon_Ptr);
    if (HitBoundary(Photon_Ptr)) {
        Hop(Photon_Ptr);    /* move to boundary plane. */
        CrossOrNot(In_Ptr, Photon_Ptr, Out_Ptr);
    } else {
        Hop(Photon_Ptr);
        Drop(Photon_Ptr);

        switch (In_Ptr->TypeBias) {
            case 0:                        /* No importance sampling */
                Spin(regionspecs[Photon_Ptr->tetrahedron->region].g, Photon_Ptr);
                break;
            case 37:                    /* BOE2 Biasing by Lima et al. */
                SpinBias37(regionspecs[Photon_Ptr->tetrahedron->region].g, Photon_Ptr, PhotonCont_Ptr, In_Ptr);
        }
    }
}

void HopDropSpin(SimulationStruct *In_Ptr, TetrahedronPhotonStruct *Photon_Ptr, TetrahedronPhotonStruct *PhotonCont_Ptr,
                 OutStruct *Out_Ptr) {
    HopDropSpinInTissue(In_Ptr, Photon_Ptr, PhotonCont_Ptr, Out_Ptr);

    if (Photon_Ptr->w < WEIGHT && !Photon_Ptr->dead)
        Roulette(Photon_Ptr);

}

void CopyPhotonStruct(TetrahedronPhotonStruct *OrigPhoton_Ptr, TetrahedronPhotonStruct *DestPhoton_Ptr) {

    DestPhoton_Ptr->x = OrigPhoton_Ptr->x;
    DestPhoton_Ptr->y = OrigPhoton_Ptr->y;
    DestPhoton_Ptr->z = OrigPhoton_Ptr->z;    /* Cartesian coordinates.[cm] */
    DestPhoton_Ptr->ux = OrigPhoton_Ptr->ux;
    DestPhoton_Ptr->uy = OrigPhoton_Ptr->uy;
    DestPhoton_Ptr->uz = OrigPhoton_Ptr->uz;/* directional cosines of a photon. */
    DestPhoton_Ptr->w = OrigPhoton_Ptr->w;            /* weight. */
    DestPhoton_Ptr->dead = OrigPhoton_Ptr->dead;        /* 1 if photon is terminated. */

    DestPhoton_Ptr->MinCos = OrigPhoton_Ptr->MinCos;
    DestPhoton_Ptr->tetrahedron = OrigPhoton_Ptr->tetrahedron;
    DestPhoton_Ptr->NextTetrahedron = OrigPhoton_Ptr->NextTetrahedron;

    DestPhoton_Ptr->s = OrigPhoton_Ptr->s;            /* current step size. [cm]. */
    DestPhoton_Ptr->sleft = OrigPhoton_Ptr->sleft;        /* step size left. dimensionless [-]. */

    DestPhoton_Ptr->OpticalPath = OrigPhoton_Ptr->OpticalPath;
    DestPhoton_Ptr->MaxDepth = OrigPhoton_Ptr->MaxDepth;
    DestPhoton_Ptr->LikelihoodRatio = OrigPhoton_Ptr->LikelihoodRatio;
    DestPhoton_Ptr->FstBackReflectionFlag = OrigPhoton_Ptr->FstBackReflectionFlag;
    DestPhoton_Ptr->LocationFstBias = OrigPhoton_Ptr->LocationFstBias;
    DestPhoton_Ptr->NumBackwardsSpecularReflections = OrigPhoton_Ptr->NumBackwardsSpecularReflections;

}


text:
/* spmatrix/test.c
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_spmatrix.h>

int status = 0;

#define BASE_GSL_COMPLEX_LONG
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX_LONG

#define BASE_GSL_COMPLEX
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX

#define BASE_GSL_COMPLEX_FLOAT
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX_FLOAT

#define BASE_LONG_DOUBLE
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_LONG_DOUBLE

#define BASE_DOUBLE
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#define BASE_ULONG
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_ULONG

#define BASE_LONG
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_LONG

#define BASE_UINT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_UINT

#define BASE_INT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_INT

#define BASE_USHORT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_USHORT

#define BASE_SHORT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_SHORT

#define BASE_UCHAR
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_UCHAR

#define BASE_CHAR
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_CHAR

int
main (void)
{
  const size_t num = 3;
  const size_t M[] = { 53, 40, 30 };
  const size_t N[] = { 107, 20, 30 };
  const double density[] = { 0.3, 0.2, 0.5 };
  gsl_rng * r = gsl_rng_alloc(gsl_rng_default);
  size_t i;

  for (i = 0; i < num; ++i)
    {
      test_all (M[i], N[i], density[i], r);
      test_float_all (M[i], N[i], density[i], r);
      test_long_double_all (M[i], N[i], density[i], r);
      test_ulong_all (M[i], N[i], density[i], r);
      test_long_all (M[i], N[i], density[i], r);
      test_uint_all (M[i], N[i], density[i], r);
      test_int_all (M[i], N[i], density[i], r);
      test_ushort_all (M[i], N[i], density[i], r);
      test_short_all (M[i], N[i], density[i], r);
      test_uchar_all (M[i], N[i], density[i], r);
      test_char_all (M[i], N[i], density[i], r);

      test_complex_all (M[i], N[i], density[i], r);
      test_complex_float_all (M[i], N[i], density[i], r);
      test_complex_long_double_all (M[i], N[i], density[i], r);
    }

  gsl_rng_free(r);

  exit (gsl_test_summary ());
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
//  This matrix class is a C++ wrapper for the GNU Scientific Library
//  Copyright (C)  ULP-IPB Strasbourg

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

#ifndef _vector_#type#_h
#define _vector_#type#_h

#ifdef __HP_aCC
#include <iostream.h>
#else
#include <iostream>
#endif

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector#typeext#.h>
#include <gsl/gsl_blas.h>
#include <gslwrap/vector_double.h>

//#define NDEBUG 0

#include <assert.h>
namespace gsl
{

#ifndef __HP_aCC
	using std::ostream;
//using std::string;
//using std::runtime_error;
#endif

class vector#typeext#_view;

class vector#typeext#
{
protected:
	gsl_vector#typeext# *gsldata;
	void free(){if(gsldata) gsl_vector#typeext#_free(gsldata);gsldata=NULL;}
	void alloc(size_t n) {gsldata=gsl_vector#typeext#_alloc(n);}
	void calloc(size_t n){gsldata=gsl_vector#typeext#_calloc(n);}
public:
	typedef #type# value_type;
	vector#typeext#() : gsldata(NULL) {;}
	vector#typeext#( const vector#typeext# &other ):gsldata(NULL) {copy(other);}
	template<class oclass>
	vector#typeext#( const oclass &other ):gsldata(NULL) {copy(other);}
	~vector#typeext#(){free();}
	vector#typeext#(const size_t& n,bool clear=true)
	{
		if(clear){this->calloc(n);}
		else     {this->alloc(n);}
	}
	vector#typeext#(const int& n,bool clear=true)
	{
		if(clear){this->calloc(n);}
		else     {this->alloc(n);}
	}

	void resize(size_t n);

	template <class oclass>
		void copy(const oclass &other)
		{
			if ( static_cast<const void *>( this ) == static_cast<const void *>( &other ) )
				return;

			if (!other.is_set())
			{
				gsldata=NULL;
				return;
			}
			resize(other.size());
			for (size_t i=0;i<size();i++)
			{
				gsl_vector#typeext#_set(gsldata, i, (#type#)other[i]);
			}
		}
	void copy(const vector#typeext#& other);
	bool is_set() const{if (gsldata) return true; else return false;}
//	void clone(vector#typeext#& other);

//	size_t size() const {if (!gsldata) {cout << "vector#typeext#::size vector not initialized" << endl; exit(-1);}return gsldata->size;}
	size_t size() const {assert (gsldata); return gsldata->size;}

	/** for interfacing with gsl c */
/*  	gsl_vector#typeext#       *gslobj()       {if (!gsldata){cout << "vector#typeext#::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return gsldata;} */
/*  	const gsl_vector#typeext# *gslobj() const {if (!gsldata){cout << "vector#typeext#::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return gsldata;} */
	gsl_vector#typeext#       *gslobj()       {assert(gsldata);return gsldata;}
	const gsl_vector#typeext# *gslobj() const {assert(gsldata);return gsldata;}

	static vector#typeext#_view create_vector_view( const gsl_vector#typeext#_view &other );

// ********Accessing vector elements

//  Unlike FORTRAN compilers, C compilers do not usually provide support for range checking of vectors and matrices (2). However, the functions gsl_vector#typeext#_get and gsl_vector#typeext#_set can perform range checking for you and report an error if you attempt to access elements outside the allowed range.

//  The functions for accessing the elements of a vector or matrix are defined in `gsl_vector#typeext#.h' and declared extern inline to eliminate function-call overhead. If necessary you can turn off range checking completely without modifying any source files by recompiling your program with the preprocessor definition GSL_RANGE_CHECK_OFF. Provided your compiler supports inline functions the effect of turning off range checking is to replace calls to gsl_vector#typeext#_get(v,i) by v->data[i*v->stride] and and calls to gsl_vector#typeext#_set(v,i,x) by v->data[i*v->stride]=x. Thus there should be no performance penalty for using the range checking functions when range checking is turned off.

//      This function returns the i-th element of a vector v. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked and 0 is returned.
	#type# get(size_t i) const {return gsl_vector#typeext#_get(gsldata,i);}

//      This function sets the value of the i-th element of a vector v to x. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked.
	void  set(size_t i,#type# x){gsl_vector#typeext#_set(gsldata,i,x);}

//      These functions return a pointer to the i-th element of a vector v. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked
	#type#       &operator[](size_t i)       { return *gsl_vector#typeext#_ptr(gsldata,i);}
	const #type# &operator[](size_t i) const { return *gsl_vector#typeext#_ptr(gsldata,i);}

	#type#       &operator()(size_t i)       { return *gsl_vector#typeext#_ptr(gsldata,i);}
	const #type# &operator()(size_t i) const { return *gsl_vector#typeext#_ptr(gsldata,i);}

//  ***** Initializing vector elements

//      This function sets all the elements of the vector v to the value x.
	void set_all(#type# x){gsl_vector#typeext#_set_all (gsldata,x);}
//      This function sets all the elements of the vector v to zero.
	void set_zero(){gsl_vector#typeext#_set_zero (gsldata);}

//      This function makes a basis vector by setting all the elements of the vector v to zero except for the i-th element which is set to one.
	int set_basis (size_t i) {return gsl_vector#typeext#_set_basis (gsldata,i);}

//  **** Reading and writing vectors

//  The library provides functions for reading and writing vectors to a file as binary data or formatted text.

//      This function writes the elements of the vector v to the stream stream in binary format. The return value is 0 for success and GSL_EFAILED if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
	int fwrite (FILE * stream) const {return gsl_vector#typeext#_fwrite (stream, gsldata);}

//      This function reads into the vector v from the open stream stream in binary format. The vector v must be preallocated with the correct length since the function uses the size of v to determine how many bytes to read. The return value is 0 for success and GSL_EFAILED if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
	int fread (FILE * stream) {return gsl_vector#typeext#_fread (stream, gsldata);}

	void load( const char *filename );
	///
	void save( const char *filename ) const;

//      This function writes the elements of the vector v line-by-line to the stream stream using the format specifier format, which should be one of the %g, %e or %f formats for floating point numbers and %d for integers. The function returns 0 for success and GSL_EFAILED if there was a problem writing to the file.
	int fprintf (FILE * stream, const char * format) const {return gsl_vector#typeext#_fprintf (stream, gsldata,format) ;}

//      This function reads formatted data from the stream stream into the vector v. The vector v must be preallocated with the correct length since the function uses the size of v to determine how many numbers to read. The function returns 0 for success and GSL_EFAILED if there was a problem reading from the file.
	int fscanf (FILE * stream)  {return gsl_vector#typeext#_fscanf (stream, gsldata); }

//  ******* Vector views

//  In addition to creating vectors from slices of blocks it is also possible to slice vectors and create vector views. For example, a subvector of another vector can be described with a view, or two views can be made which provide access to the even and odd elements of a vector.

//  A vector view is a temporary object, stored on the stack, which can be used to operate on a subset of vector elements. Vector views can be defined for both constant and non-constant vectors, using separate types that preserve constness. A vector view has the type gsl_vector#typeext#_view and a constant vector view has the type gsl_vector#typeext#_const_view. In both cases the elements of the view can be accessed as a gsl_vector#typeext# using the vector component of the view object. A pointer to a vector of type gsl_vector#typeext# * or const gsl_vector#typeext# * can be obtained by taking the address of this component with the & operator.

//      These functions return a vector view of a subvector of another vector v. The start of the new vector is offset by offset elements from the start of the original
//      vector. The new vector has n elements. Mathematically, the i-th element of the new vector v' is given by,

//      v'(i) = v->data[(offset + i)*v->stride]

//      where the index i runs from 0 to n-1.

//      The data pointer of the returned vector struct is set to null if the combined parameters (offset,n) overrun the end of the original vector.

//      The new vector is only a view of the block underlying the original vector, v. The block containing the elements of v is not owned by the new vector. When the
//      new vector goes out of scope the original vector v and its block will continue to exist. The original memory can only be deallocated by freeing the original vector.
//      Of course, the original vector should not be deallocated while the new vector is still in use.

//      The function gsl_vector#typeext#_const_subvector is equivalent to gsl_vector#typeext#_subvector but can be used for vectors which are declared const.
	vector#typeext#_view subvector (size_t offset, size_t n);
	const vector#typeext#_view subvector (size_t offset, size_t n) const;
//	vector#typeext#_const_view subvector (size_t offset, size_t n) const;

//  	class view
//  	{
//  		gsl_vector#typeext#_view *gsldata;
//  	public:
//  		view();
//  	};
//  	view subvector(size_t offset, size_t n)
//  	{
//  		return view(gsl_vector#typeext#_subvector(gsldata,offset,n);
//  	}
//  	const view subvector(size_t offset, size_t n) const
//  	{
//  		return view(gsl_vector#typeext#_const_subvector(gsldata,offset,n);
//  	}

//  Function: gsl_vector#typeext# gsl_vector#typeext#_subvector_with_stride (gsl_vector#typeext# *v, size_t offset, size_t stride, size_t n)
//  Function: gsl_vector#typeext#_const_view gsl_vector#typeext#_const_subvector_with_stride (const gsl_vector#typeext# * v, size_t offset, size_t stride, size_t n)
//      These functions return a vector view of a subvector of another vector v with an additional stride argument. The subvector is formed in the same way as for
//      gsl_vector#typeext#_subvector but the new vector has n elements with a step-size of stride from one element to the next in the original vector. Mathematically,
//      the i-th element of the new vector v' is given by,

//      v'(i) = v->data[(offset + i*stride)*v->stride]

//      where the index i runs from 0 to n-1.

//      Note that subvector views give direct access to the underlying elements of the original vector. For example, the following code will zero the even elements of the
//      vector v of length n, while leaving the odd elements untouched,

//      gsl_vector#typeext#_view v_even = gsl_vector#typeext#_subvector_with_stride (v, 0, 2, n/2);
//      gsl_vector#typeext#_set_zero (&v_even.vector);

//      A vector view can be passed to any subroutine which takes a vector argument just as a directly allocated vector would be, using &view.vector. For example, the
//      following code computes the norm of odd elements of v using the BLAS routine DNRM2,

//      gsl_vector#typeext#_view v_odd = gsl_vector#typeext#_subvector_with_stride (v, 1, 2, n/2);
//      double r = gsl_blas_dnrm2 (&v_odd.vector);

//      The function gsl_vector#typeext#_const_subvector_with_stride is equivalent to gsl_vector#typeext#_subvector_with_stride but can be used for
//      vectors which are declared const.

//  Function: gsl_vector#typeext#_view gsl_vector#typeext#_complex_real (gsl_vector#typeext#_complex *v)
//  Function: gsl_vector#typeext#_const_view gsl_vector#typeext#_complex_const_real (const gsl_vector#typeext#_complex *v)
//      These functions return a vector view of the real parts of the complex vector v.

//      The function gsl_vector#typeext#_complex_const_real is equivalent to gsl_vector#typeext#_complex_real but can be used for vectors which are declared
//      const.

//  Function: gsl_vector#typeext#_view gsl_vector#typeext#_complex_imag (gsl_vector#typeext#_complex *v)
//  Function: gsl_vector#typeext#_const_view gsl_vector#typeext#_complex_const_imag (const gsl_vector#typeext#_complex *v)
//      These functions return a vector view of the imaginary parts of the complex vector v.

//      The function gsl_vector#typeext#_complex_const_imag is equivalent to gsl_vector#typeext#_complex_imag but can be used for vectors which are declared
//      const.

//  Function: gsl_vector#typeext#_view gsl_vector#typeext#_view_array (double *base, size_t n)
//  Function: gsl_vector#typeext#_const_view gsl_vector#typeext#_const_view_array (const double *base, size_t n)
//      These functions return a vector view of an array. The start of the new vector is given by base and has n elements. Mathematically, the i-th element of the new
//      vector v' is given by,

//      v'(i) = base[i]

//      where the index i runs from 0 to n-1.

//      The array containing the elements of v is not owned by the new vector view. When the view goes out of scope the original array will continue to exist. The
//      original memory can only be deallocated by freeing the original pointer base. Of course, the original array should not be deallocated while the view is still in use.

//      The function gsl_vector#typeext#_const_view_array is equivalent to gsl_vector#typeext#_view_array but can be used for vectors which are declared const.

//  Function: gsl_vector#typeext#_view gsl_vector#typeext#_view_array_with_stride (double * base, size_t stride, size_t n)
//  Function: gsl_vector#typeext#_const_view gsl_vector#typeext#_const_view_array_with_stride (const double * base, size_t stride, size_t n)
//      These functions return a vector view of an array base with an additional stride argument. The subvector is formed in the same way as for
//      gsl_vector#typeext#_view_array but the new vector has n elements with a step-size of stride from one element to the next in the original array. Mathematically,
//      the i-th element of the new vector v' is given by,

//      v'(i) = base[i*stride]

//      where the index i runs from 0 to n-1.

//      Note that the view gives direct access to the underlying elements of the original array. A vector view can be passed to any subroutine which takes a vector
//      argument just as a directly allocated vector would be, using &view.vector.

//      The function gsl_vector#typeext#_const_view_array_with_stride is equivalent to gsl_vector#typeext#_view_array_with_stride but can be used for
//      arrays which are declared const.

//  ************* Copying vectors

//  Common operations on vectors such as addition and multiplication are available in the BLAS part of the library (see section BLAS Support). However, it is useful to have a small number of utility functions which do not require the full BLAS code. The following functions fall into this category.

//      This function copies the elements of the vector src into the vector dest.
	vector#typeext#& operator=(const vector#typeext#& other){copy(other);return (*this);}

//  Function: int gsl_vector#typeext#_swap (gsl_vector#typeext# * v, gsl_vector#typeext# * w)
//      This function exchanges the elements of the vectors v and w by copying. The two vectors must have the same length.

//  ***** Exchanging elements

//  The following function can be used to exchange, or permute, the elements of a vector.

//  Function: int gsl_vector#typeext#_swap_elements (gsl_vector#typeext# * v, size_t i, size_t j)
//      This function exchanges the i-th and j-th elements of the vector v in-place.
	int swap_elements (size_t i, size_t j) {return gsl_vector#typeext#_swap_elements (gsldata, i,j);}

//  Function: int gsl_vector#typeext#_reverse (gsl_vector#typeext# * v)
//      This function reverses the order of the elements of the vector v.
	int reverse () {return  gsl_vector#typeext#_reverse (gsldata) ;}

// ******* Vector operations

//  The following operations are only defined for real vectors.

//      This function adds the elements of vector b to the elements of vector a, a'_i = a_i + b_i. The two vectors must have the same length.
	int operator+=(const vector#typeext# &other) {return gsl_vector#typeext#_add (gsldata, other.gsldata);}

//      This function subtracts the elements of vector b from the elements of vector a, a'_i = a_i - b_i. The two vectors must have the same length.
	int operator-=(const vector#typeext# &other) {return gsl_vector#typeext#_sub (gsldata, other.gsldata);}

//  Function: int gsl_vector#typeext#_mul (gsl_vector#typeext# * a, const gsl_vector#typeext# * b)
//      This function multiplies the elements of vector a by the elements of vector b, a'_i = a_i * b_i. The two vectors must have the same length.
	int operator*=(const vector#typeext# &other) {return gsl_vector#typeext#_mul (gsldata, other.gsldata);}

//      This function divides the elements of vector a by the elements of vector b, a'_i = a_i / b_i. The two vectors must have the same length.
	int operator/=(const vector#typeext# &other) {return gsl_vector#typeext#_div (gsldata, other.gsldata);}

//      This function multiplies the elements of vector a by the constant factor x, a'_i = x a_i.
	int operator*=(#type# x) {return gsl_vector#typeext#_scale (gsldata, x);}

//  Function: int gsl_vector#typeext#_add_constant (gsl_vector#typeext# * a, const double x)
//      This function adds the constant value x to the elements of the vector a, a'_i = a_i + x.
	int operator+=(#type# x) {return gsl_vector#typeext#_add_constant (gsldata,x);}

//      This function multiplies the elements of vector a by the constant factor x, a'_i = x a_i.
	int operator/=(#type# x) {return gsl_vector#typeext#_scale (gsldata, 1/x);}

// bool operators:
	bool operator==(const vector#typeext#& other) const;
	bool operator!=(const vector#typeext#& other) const { return (!((*this)==other));}

// stream output:
//	friend ostream& operator<< ( ostream& os, const vector#typeext#& vect );
	/** returns sum of all the vector elements. */
    #type# sum() const;
	// returns sqrt(v.t*v);
    double norm2() const;

// **** Finding maximum and minimum elements of vectors

//      This function returns the maximum value in the vector v.
    double max() const{return gsl_vector#typeext#_max (gsldata) ;}

//  Function: double gsl_vector#typeext#_min (const gsl_vector#typeext# * v)
//      This function returns the minimum value in the vector v.
    double min() const{return gsl_vector#typeext#_min (gsldata) ;}

//  Function: void gsl_vector#typeext#_minmax (const gsl_vector#typeext# * v, double * min_out, double * max_out)
//      This function returns the minimum and maximum values in the vector v, storing them in min_out and max_out.

//      This function returns the index of the maximum value in the vector v. When there are several equal maximum elements then the lowest index is returned.
	size_t max_index(){return gsl_vector#typeext#_max_index (gsldata);}

//  Function: size_t gsl_vector#typeext#_min_index (const gsl_vector#typeext# * v)
//      This function returns the index of the minimum value in the vector v. When there are several equal minimum elements then the lowest index is returned.
	size_t min_index(){return gsl_vector#typeext#_min_index (gsldata);}

//  Function: void gsl_vector#typeext#_minmax_index (const gsl_vector#typeext# * v, size_t * imin, size_t * imax)
//      This function returns the indices of the minimum and maximum values in the vector v, storing them in imin and imax. When there are several equal minimum
//      or maximum elements then the lowest indices are returned.

//  Vector properties

//  Function: int gsl_vector#typeext#_isnull (const gsl_vector#typeext# * v)
//      This function returns 1 if all the elements of the vector v are zero, and 0 otherwise. };
	bool isnull(){return gsl_vector#typeext#_isnull (gsldata);}
};

// When you add create a view it will stick to its with the view until you call change_view
// ex:
// matrix_float m(5,5);
// vector_float v(5);
// // ...
// m.column(3) = v; //the 3rd column of the matrix m will equal v.
class vector#typeext#_view : public vector#typeext#
{
 public:
	vector#typeext#_view(const vector#typeext#&     other) :vector#typeext#(){init(other);}
	vector#typeext#_view(const vector#typeext#_view& other):vector#typeext#(){init(other);}
	vector#typeext#_view(const gsl_vector#typeext#& gsl_other) : vector#typeext#() {init_with_gsl_vector(gsl_other);}

	void init(const vector#typeext#& other);
	void init_with_gsl_vector(const gsl_vector#typeext#& gsl_other);
	void change_view(const vector#typeext#& other){init(other);}
 private:
};

ostream& operator<< ( ostream& os, const vector#typeext# & vect );

// vector_type<>::type is a template interface to vector_?
// it is usefull for in templated situations for getting the correct vector type
#define tmp_type_is#typeext#
#ifdef tmp_type_is
typedef vector vector_double;
template<class T>
struct vector_type  {typedef vector_double   type;};

template<class T>
struct value_type  {typedef double   type;};

#else
template<> struct vector_type<#type#> {typedef vector#typeext# type;};
#endif
#undef tmp_type_is#typeext#

}
#endif// _vector_#type#_h

ext: h

max_stars_repo_name: entn-at/GlottDNN

max_issues_repo_name: entn-at/GlottDNN

max_forks_repo_name: entn-at/GlottDNN


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"

#if (defined(SMOOTH_PHI) || defined(SMOOTH_ROTB) || defined(BSMOOTH))

/*! Structure for communication during the density computation. Holds data that is sent to other processors.
 */
static struct smoothdata_in
{
  MyDouble Pos[3];
  MyFloat Hsml;
  int NodeList[NODELISTLENGTH];
}
 *SmoothDataIn, *SmoothDataGet;

static struct smoothdata_out
{
#ifdef SMOOTH_PHI
  MyFloat SmoothPhi;
  MyFloat SmoothDivB;
#endif
#ifdef SMOOTH_ROTB
  MyFloat SmoothRotB[3];
#endif
#ifdef BSMOOTH
  MyFloat BSmooth[3];
#endif
#if defined(BLACK_HOLES)
  MyLongDouble SmoothedEntr;
#endif
  MyFloat DensityNorm;
}
 *SmoothDataResult, *SmoothDataOut;

void smoothed_values(void)
{
  int ngrp, sendTask, recvTask, place, nexport, nimport;
  int i, j, ndone, ndone_flag, dummy;
  double timeall = 0, timecomp1 = 0, timecomp2 = 0, timecommsumm1 = 0, timecommsumm2 = 0, timewait1 =
    0, timewait2 = 0;
  double timecomp, timecomm, timewait;
  double tstart, tend, t0, t1;

#ifdef BSMOOTH
  int Smooth_Flag = 0;
  double dB[3];
#endif

  /* Display information message that this step is executed on Task 0 ... */
  if(ThisTask == 0)
    {
      printf("Updating SPH interpolants for:"
#ifdef SMOOTH_PHI
	     " (divB and Phi(Dedner)"
#endif /* SMOOTH_PHI */
#ifdef SMOOTH_ROTB
	     " (rotB)"
#endif /* SMOOTH_ROTB */
#ifdef BSMOOTH
	     " (B)"
#endif /* BSMOOTH */
	     "\n");
#ifdef BSMOOTH
      printf("Flag_FullStep = %d, Main TimestepCounts = %d\n", Flag_FullStep, All.MainTimestepCounts);
#endif
    }
#ifdef BSMOOTH
  if(Flag_FullStep == 1)
    {
      if((All.MainTimestepCounts % All.BSmoothInt == 0) && (All.BSmoothInt >= 0))
	{
	  Smooth_Flag = 1;
	  if(ThisTask == 0)
	    printf("Smoothing B %d, %f\n", All.BSmoothInt, All.BSmoothFrac);
	}
      All.MainTimestepCounts++;
    }
#endif

  Ngblist = (int *) mymalloc(NumPart * sizeof(int));

  All.BunchSize =
    (int) ((All.BufferSize * 1024 * 1024) / (sizeof(struct data_index) + sizeof(struct data_nodelist) +
					     sizeof(struct smoothdata_in) + sizeof(struct smoothdata_out) +
					     sizemax(sizeof(struct smoothdata_in),
						     sizeof(struct smoothdata_out))));
  DataIndexTable = (struct data_index *) mymalloc(All.BunchSize * sizeof(struct data_index));
  DataNodeList = (struct data_nodelist *) mymalloc(All.BunchSize * sizeof(struct data_nodelist));

  CPU_Step[CPU_SMTHMISC] += measure_time();
  t0 = second();

  i = FirstActiveParticle;	/* begin with this index */

  do
    {
      for(j = 0; j < NTask; j++)
	{
	  Send_count[j] = 0;
	  Exportflag[j] = -1;
	}

      /* do local particles and prepare export list */
      tstart = second();
      for(nexport = 0; i >= 0; i = NextActiveParticle[i])
	{
	  if(density_isactive(i))
	    {
	      if(smoothed_evaluate(i, 0, &nexport, Send_count) < 0)
		break;
	    }
	}
      tend = second();
      timecomp1 += timediff(tstart, tend);

#ifdef MYSORT
      mysort_dataindex(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);
#else
      qsort(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);
#endif

      tstart = second();

      MPI_Allgather(Send_count, NTask, MPI_INT, Sendcount_matrix, NTask, MPI_INT, MPI_COMM_WORLD);

      tend = second();
      timewait1 += timediff(tstart, tend);

      for(j = 0, nimport = 0, Recv_offset[0] = 0, Send_offset[0] = 0; j < NTask; j++)
	{
	  Recv_count[j] = Sendcount_matrix[j * NTask + ThisTask];
	  nimport += Recv_count[j];

	  if(j > 0)
	    {
	      Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
	      Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
	    }
	}

      SmoothDataGet = (struct smoothdata_in *) mymalloc(nimport * sizeof(struct smoothdata_in));
      SmoothDataIn = (struct smoothdata_in *) mymalloc(nexport * sizeof(struct smoothdata_in));

      /* prepare particle data for export */
      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;

	  SmoothDataIn[j].Pos[0] = P[place].Pos[0];
	  SmoothDataIn[j].Pos[1] = P[place].Pos[1];
	  SmoothDataIn[j].Pos[2] = P[place].Pos[2];
	  SmoothDataIn[j].Hsml = PPP[place].Hsml;

	  memcpy(SmoothDataIn[j].NodeList,
		 DataNodeList[DataIndexTable[j].IndexGet].NodeList, NODELISTLENGTH * sizeof(int));
	}

      /* exchange particle data */
      tstart = second();
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;

	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* get the particles */
		  MPI_Sendrecv(&SmoothDataIn[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct smoothdata_in), MPI_BYTE,
			       recvTask, TAG_SMOOTH_A,
			       &SmoothDataGet[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct smoothdata_in), MPI_BYTE,
			       recvTask, TAG_SMOOTH_A, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		}
	    }
	}
      tend = second();
      timecommsumm1 += timediff(tstart, tend);

      myfree(SmoothDataIn);
      SmoothDataResult = (struct smoothdata_out *) mymalloc(nimport * sizeof(struct smoothdata_out));
      SmoothDataOut = (struct smoothdata_out *) mymalloc(nexport * sizeof(struct smoothdata_out));

      /* now do the particles that were sent to us */

      tstart = second();
      for(j = 0; j < nimport; j++)
	smoothed_evaluate(j, 1, &dummy, &dummy);
      tend = second();
      timecomp2 += timediff(tstart, tend);

      if(i < 0)
	ndone_flag = 1;
      else
	ndone_flag = 0;

      tstart = second();
      MPI_Allreduce(&ndone_flag, &ndone, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
      tend = second();
      timewait2 += timediff(tstart, tend);

      /* get the result */
      tstart = second();
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;
	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* send the results */
		  MPI_Sendrecv(&SmoothDataResult[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct smoothdata_out),
			       MPI_BYTE, recvTask, TAG_SMOOTH_B,
			       &SmoothDataOut[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct smoothdata_out),
			       MPI_BYTE, recvTask, TAG_SMOOTH_B, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		}
	    }
	}
      tend = second();
      timecommsumm2 += timediff(tstart, tend);

      /* add the result to the local particles */
      tstart = second();
      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;

	  if(P[place].Type == 0)
	    {
#ifdef SMOOTH_PHI
	      SphP[place].SmoothPhi += SmoothDataOut[j].SmoothPhi;
	      SphP[place].SmoothDivB += SmoothDataOut[j].SmoothDivB;
#endif /* SMOOTH_PHI */

#ifdef SMOOTH_ROTB
	      SphP[place].SmoothedRotB[0] += SmoothDataOut[j].SmoothRotB[0];
	      SphP[place].SmoothedRotB[1] += SmoothDataOut[j].SmoothRotB[1];
	      SphP[place].SmoothedRotB[2] += SmoothDataOut[j].SmoothRotB[2];
#endif /* SMOOTH_ROTB */

#ifdef BSMOOTH
	      SphP[place].BSmooth[0] += SmoothDataOut[j].BSmooth[0];
	      SphP[place].BSmooth[1] += SmoothDataOut[j].BSmooth[1];
	      SphP[place].BSmooth[2] += SmoothDataOut[j].BSmooth[2];
#endif /* BSMOOTH */

	      SphP[place].DensityNorm += SmoothDataOut[j].DensityNorm;
	    }
	}
      tend = second();
      timecomp1 += timediff(tstart, tend);

      myfree(SmoothDataOut);
      myfree(SmoothDataResult);
      myfree(SmoothDataGet);
    }
  while(ndone < NTask);

  myfree(DataNodeList);
  myfree(DataIndexTable);
  myfree(Ngblist);

  /* do final operations on results */
  tstart = second();
  for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
    {
      if(density_isactive(i))
	if(P[i].Type == 0)
	  {
#ifdef SMOOTH_PHI
	    SphP[i].SmoothPhi /= SphP[i].DensityNorm;
	    SphP[i].SmoothDivB /= SphP[i].DensityNorm;
#endif /* SMOOTH_PHI */

#ifdef SMOOTH_ROTB
	    SphP[i].SmoothedRotB[0] /= SphP[i].DensityNorm;
	    SphP[i].SmoothedRotB[1] /= SphP[i].DensityNorm;
	    SphP[i].SmoothedRotB[2] /= SphP[i].DensityNorm;
#endif /* SMOOTH_ROTB */

#ifdef BSMOOTH
	    SphP[i].BSmooth[0] /= SphP[i].DensityNorm;
	    SphP[i].BSmooth[1] /= SphP[i].DensityNorm;
	    SphP[i].BSmooth[2] /= SphP[i].DensityNorm;

	    if(Smooth_Flag == 1)
	      {
		dB[0] = All.BSmoothFrac * (SphP[i].BSmooth[0] - SphP[i].BPred[0]);
		dB[1] = All.BSmoothFrac * (SphP[i].BSmooth[1] - SphP[i].BPred[1]);
		dB[2] = All.BSmoothFrac * (SphP[i].BSmooth[2] - SphP[i].BPred[2]);
		SphP[i].BPred[0] += dB[0];
		SphP[i].BPred[1] += dB[1];
		SphP[i].BPred[2] += dB[2];
#ifndef EULERPOTENTIALS
		SphP[i].B[0] += dB[0];
		SphP[i].B[1] += dB[1];
		SphP[i].B[2] += dB[2];
#endif
	      }
#endif /* BSMOOTH */
	  }
    }
  tend = second();
  timecomp1 += timediff(tstart, tend);

  /* collect some timing information */

  t1 = WallclockTime = second();
  timeall += timediff(t0, t1);

  timecomp = timecomp1 + timecomp2;
  timewait = timewait1 + timewait2;
  timecomm = timecommsumm1 + timecommsumm2;

  CPU_Step[CPU_SMTHCOMPUTE] += timecomp;
  CPU_Step[CPU_SMTHWAIT] += timewait;
  CPU_Step[CPU_SMTHCOMM] += timecomm;
  CPU_Step[CPU_SMTHMISC] += timeall - (timecomp + timewait + timecomm);
}

/*! This function represents the core of the SPH density computation. The
*  target particle may either be local, or reside in the communication
*  buffer.
*/
int smoothed_evaluate(int target, int mode, int *nexport, int *nsend_local)
{
  int j, n, listindex = 0;
  int startnode, numngb_inbox;
  double h, h2, hinv, hinv3, hinv4;
  double wk, dwk;
  double dx, dy, dz, r, r2, u, mass_j;
  MyFloat *pos;
  double DensityNorm = 0;

#ifdef SMOOTH_PHI
  double SmoothPhi = 0.0;
  double SmoothDivB = 0.0;
#endif /* SMOOTH_PHI */

#ifdef SMOOTH_ROTB
  double smoothrotb[3];
#endif /* SMOOTH_ROTB */

#ifdef BSMOOTH
  double BSmooth[3];
#endif /* BSMOOTH */

#ifdef SMOOTH_ROTB
  smoothrotb[0] = smoothrotb[1] = smoothrotb[2] = 0;
#endif /* SMOOTH_ROTB */

#ifdef BSMOOTH
  BSmooth[0] = BSmooth[1] = BSmooth[2] = 0;
#endif /* BSMOOTH */

  if(mode == 0)
    {
      pos = P[target].Pos;
      h = PPP[target].Hsml;
    }
  else
    {
      pos = SmoothDataGet[target].Pos;
      h = SmoothDataGet[target].Hsml;
    }

  h2 = h * h;
  hinv = 1.0 / h;
#ifndef  TWODIMS
  hinv3 = hinv * hinv * hinv;
#else
  hinv3 = hinv * hinv / boxSize_Z;
#endif
  hinv4 = hinv3 * hinv;

  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = SmoothDataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
	  numngb_inbox = ngb_treefind_variable(&pos[0], h, target, &startnode, mode, nexport, nsend_local);

	  if(numngb_inbox < 0)
	    return -1;

	  for(n = 0; n < numngb_inbox; n++)
	    {
	      j = Ngblist[n];

	      dx = pos[0] - P[j].Pos[0];
	      dy = pos[1] - P[j].Pos[1];
	      dz = pos[2] - P[j].Pos[2];

#ifdef PERIODIC			/*  now find the closest image in the given box size  */
	      if(dx > boxHalf_X)
		dx -= boxSize_X;
	      if(dx < -boxHalf_X)
		dx += boxSize_X;
	      if(dy > boxHalf_Y)
		dy -= boxSize_Y;
	      if(dy < -boxHalf_Y)
		dy += boxSize_Y;
	      if(dz > boxHalf_Z)
		dz -= boxSize_Z;
	      if(dz < -boxHalf_Z)
		dz += boxSize_Z;
#endif
	      r2 = dx * dx + dy * dy + dz * dz;

	      if(r2 < h2)
		{
		  r = sqrt(r2);

		  u = r * hinv;

		  if(u < 0.5)
		    {
		      wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		      dwk = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
		    }
		  else
		    {
		      wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
		      dwk = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
		    }

		  mass_j = P[j].Mass;
		  wk /= SphP[j].d.Density;

#ifdef SMOOTH_PHI
		  SmoothPhi += mass_j * wk * SphP[j].PhiPred;
		  SmoothDivB += mass_j * wk * SphP[j].divB;
#endif /* SMOOTH_PHI */

#ifdef SMOOTH_ROTB
		  smoothrotb[0] += mass_j * wk * SphP[j].RotB[0];
		  smoothrotb[1] += mass_j * wk * SphP[j].RotB[1];
		  smoothrotb[2] += mass_j * wk * SphP[j].RotB[2];
#endif /* SMOOTH_ROTB */

#ifdef BSMOOTH
		  BSmooth[0] += mass_j * wk * SphP[j].BPred[0];
		  BSmooth[1] += mass_j * wk * SphP[j].BPred[1];
		  BSmooth[2] += mass_j * wk * SphP[j].BPred[2];
#endif /* BSMOOTH */

		  DensityNorm += mass_j * wk;
		}
	    }
	}
      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = SmoothDataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  if(mode == 0)
    {
#ifdef SMOOTH_PHI
      SphP[target].SmoothPhi = SmoothPhi;
      SphP[target].SmoothDivB = SmoothDivB;
#endif /* SMOOTH_PHI */

#ifdef SMOOTH_ROTB
      SphP[target].SmoothedRotB[0] = smoothrotb[0];
      SphP[target].SmoothedRotB[1] = smoothrotb[1];
      SphP[target].SmoothedRotB[2] = smoothrotb[2];
#endif /* SMOOTH_ROTB */

#ifdef BSMOOTH
      SphP[target].BSmooth[0] = BSmooth[0];
      SphP[target].BSmooth[1] = BSmooth[1];
      SphP[target].BSmooth[2] = BSmooth[2];
#endif /* BSMOOTH */
      SphP[target].DensityNorm = DensityNorm;
    }
  else
    {
#ifdef SMOOTH_PHI
      SmoothDataResult[target].SmoothPhi = SmoothPhi;
      SmoothDataResult[target].SmoothDivB = SmoothDivB;
#endif /* SMOOTH_PHI */

#ifdef SMOOTH_ROTB
      SmoothDataResult[target].SmoothRotB[0] = smoothrotb[0];
      SmoothDataResult[target].SmoothRotB[1] = smoothrotb[1];
      SmoothDataResult[target].SmoothRotB[2] = smoothrotb[2];
#endif /* SMOOTH_ROTB */

#ifdef BSMOOTH
      SmoothDataResult[target].BSmooth[0] = BSmooth[0];
      SmoothDataResult[target].BSmooth[1] = BSmooth[1];
      SmoothDataResult[target].BSmooth[2] = BSmooth[2];
#endif /* BSMOOTH */
      SmoothDataResult[target].DensityNorm = DensityNorm;
    }

  return 0;
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/**
 * This file is part of the "libterminal" project
 *   Copyright (c) 2019-2020 Christian Parpart <christian@parpart.family>
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once

#include <terminal/GraphicsAttributes.h>
#include <terminal/Line.h>
#include <terminal/primitives.h>

#include <crispy/algorithm.h>
#include <crispy/assert.h>
#include <crispy/ring.h>

#include <unicode/convert.h>

#include <range/v3/algorithm/copy.hpp>
#include <range/v3/iterator/insert_iterators.hpp>
#include <range/v3/view/iota.hpp>

#include <gsl/span>
#include <gsl/span_ext>

#include <algorithm>
#include <array>
#include <sstream>
#include <string>
#include <string_view>
#include <utility>

namespace terminal
{

// {{{ Margin
struct Margin
{
    struct Horizontal
    {
        ColumnOffset from;
        ColumnOffset
            to; // TODO: call it begin and end and have end point to to+1 to avoid unnecessary +1's later

        [[nodiscard]] constexpr ColumnCount length() const noexcept
        {
            return unbox<ColumnCount>(to - from) + ColumnCount(1);
        }
        [[nodiscard]] constexpr bool contains(ColumnOffset _value) const noexcept
        {
            return from <= _value && _value <= to;
        }
        [[nodiscard]] constexpr bool operator==(Horizontal rhs) const noexcept
        {
            return from == rhs.from && to == rhs.to;
        }
        [[nodiscard]] constexpr bool operator!=(Horizontal rhs) const noexcept { return !(*this == rhs); }
    };

    struct Vertical
    {
        LineOffset from;
        // TODO: call it begin and end and have end point to to+1 to avoid unnecessary +1's later
        LineOffset to;

        [[nodiscard]] constexpr LineCount length() const noexcept
        {
            return unbox<LineCount>(to - from) + LineCount(1);
        }
        [[nodiscard]] constexpr bool contains(LineOffset _value) const noexcept
        {
            return from <= _value && _value <= to;
        }
        [[nodiscard]] constexpr bool operator==(Vertical const& rhs) const noexcept
        {
            return from == rhs.from && to == rhs.to;
        }
        [[nodiscard]] constexpr bool operator!=(Vertical const& rhs) const noexcept
        {
            return !(*this == rhs);
        }
    };

    Vertical vertical {};     // top-bottom
    Horizontal horizontal {}; // left-right
};

constexpr bool operator==(Margin const& a, PageSize b) noexcept
{
    return a.horizontal.from.value == 0 && a.horizontal.to.value + 1 == b.columns.value
           && a.vertical.from.value == 0 && a.vertical.to.value + 1 == b.lines.value;
}

constexpr bool operator!=(Margin const& a, PageSize b) noexcept
{
    return !(a == b);
}
// }}}

template <typename Cell>
using Lines = crispy::ring<Line<Cell>>;

/**
 * Represents a logical grid line, i.e. a sequence lines that were written without
 * an explicit linefeed, triggering an auto-wrap.
 */
template <typename Cell>
struct LogicalLine
{
    LineOffset top {};
    LineOffset bottom {};
    std::vector<std::reference_wrapper<Line<Cell>>> lines {};

    [[nodiscard]] Line<Cell> joinWithRightTrimmed() const
    {
        // TODO: determine final line's column count and pass it to ctor.
        typename Line<Cell>::Buffer output;
        auto lineFlags = lines.front().get().flags();
        for (Line<Cell> const& line: lines)
            for (Cell const& cell: line.cells())
                output.emplace_back(cell);

        while (!output.empty() && output.back().empty())
            output.pop_back();

        return Line<Cell>(output, lineFlags);
    }

    [[nodiscard]] std::string text() const
    {
        std::string output;
        for (auto const& line: lines)
            output += line.get().toUtf8();
        return output;
    }
};

template <typename Cell>
bool operator==(LogicalLine<Cell> const& a, LogicalLine<Cell> const& b) noexcept
{
    return a.top == b.top && a.bottom == b.bottom;
}

template <typename Cell>
bool operator!=(LogicalLine<Cell> const& a, LogicalLine<Cell> const& b) noexcept
{
    return !(a == b);
}

template <typename Cell>
struct LogicalLines
{
    LineOffset topMostLine;
    LineOffset bottomMostLine;
    std::reference_wrapper<Lines<Cell>> lines;

    struct iterator // {{{
    {
        std::reference_wrapper<Lines<Cell>> lines;
        LineOffset top;
        LineOffset next; // index to next logical line's beginning
        LineOffset bottom;
        LogicalLine<Cell> current;

        iterator(std::reference_wrapper<Lines<Cell>> _lines,
                 LineOffset _top,
                 LineOffset _next,
                 LineOffset _bottom):
            lines { _lines }, top { _top }, next { _next }, bottom { _bottom }
        {
            Require(_top <= next);
            Require(next <= _bottom + 1);
            ++*this;
        }

        LogicalLine<Cell> const& operator*() const noexcept { return current; }
        LogicalLine<Cell> const* operator->() const noexcept { return &current; }

        iterator& operator++()
        {
            if (next == bottom + 1)
            {
                current.top = next;
                current.bottom = next;
                return *this;
            }

            Require(!lines.get()[unbox<int>(next)].wrapped());

            current.top = LineOffset::cast_from(next);
            current.lines.clear();
            do
                current.lines.emplace_back(lines.get()[unbox<int>(next++)]);
            while (next <= bottom && lines.get()[unbox<int>(next)].wrapped());

            current.bottom = LineOffset::cast_from(next - 1);

            return *this;
        }

        iterator& operator--()
        {
            if (next == top - 1)
            {
                current.top = top - 1;
                current.bottom = top - 1;
                return *this;
            }

            auto const bottomMost = next - 1;
            do
                --next;
            while (lines.get()[unbox<int>(next)].wrapped());
            auto const topMost = next;

            current.top = topMost;
            current.bottom = bottomMost;

            current.lines.clear();
            for (auto i = topMost; i <= bottomMost; ++i)
                current.lines.emplace_back(lines.get()[unbox<int>(i)]);

            return *this;
        }

        iterator& operator++(int)
        {
            auto c = *this;
            ++*this;
            return c;
        }
        iterator& operator--(int)
        {
            auto c = *this;
            --*this;
            return c;
        }

        bool operator==(iterator const& other) const noexcept { return current == other.current; }
        bool operator!=(iterator const& other) const noexcept { return current != other.current; }
    }; // }}}

    iterator begin() const { return iterator(lines, topMostLine, topMostLine, bottomMostLine); }
    iterator end() const { return iterator(lines, topMostLine, bottomMostLine + 1, bottomMostLine); }
};

template <typename Cell>
struct ReverseLogicalLines
{
    LineOffset topMostLine;
    LineOffset bottomMostLine;
    std::reference_wrapper<Lines<Cell>> lines;

    struct iterator // {{{
    {
        std::reference_wrapper<Lines<Cell>> lines;
        LineOffset top;
        LineOffset next; // index to next logical line's beginning
        LineOffset bottom;
        LogicalLine<Cell> current;

        iterator(std::reference_wrapper<Lines<Cell>> _lines,
                 LineOffset _top,
                 LineOffset _next,
                 LineOffset _bottom):
            lines { _lines }, top { _top }, next { _next }, bottom { _bottom }
        {
            Require(_top - 1 <= next);
            Require(next <= _bottom);
            ++*this;
        }

        LogicalLine<Cell> const& operator*() const noexcept { return current; }

        iterator& operator--()
        {
            if (next == bottom + 1)
            {
                current.top = bottom + 1;
                current.bottom = bottom + 1;
                return *this;
            }

            Require(!lines.get()[unbox<int>(next)].wrapped());

            current.top = LineOffset::cast_from(next);
            current.lines.clear();
            do
                current.lines.emplace_back(lines.get()[unbox<int>(next++)]);
            while (next <= bottom && lines.get()[unbox<int>(next)].wrapped());

            current.bottom = LineOffset::cast_from(next - 1);

            return *this;
        }

        iterator& operator++()
        {
            if (next == top - 1)
            {
                current.top = next;
                current.bottom = next;
                return *this;
            }

            auto const bottomMost = next;
            while (lines.get()[unbox<int>(next)].wrapped())
                --next;
            auto const topMost = next;
            --next; // jump to next logical line's bottom line above the current logical one

            current.top = topMost;
            current.bottom = bottomMost;

            current.lines.clear();
            for (auto i = topMost; i <= bottomMost; ++i)
                current.lines.emplace_back(lines.get()[unbox<int>(i)]);

            return *this;
        }

        iterator& operator++(int)
        {
            auto c = *this;
            ++*this;
            return c;
        }
        iterator& operator--(int)
        {
            auto c = *this;
            --*this;
            return c;
        }

        bool operator==(iterator const& other) const noexcept { return current == other.current; }
        bool operator!=(iterator const& other) const noexcept { return current != other.current; }
    }; // }}}

    iterator begin() const { return iterator(lines, topMostLine, bottomMostLine, bottomMostLine); }
    iterator end() const { return iterator(lines, topMostLine, topMostLine - 1, bottomMostLine); }
};

/**
 * Manages the screen grid buffer (main screen + scrollback history).
 *
 *
Future motivations

 *
 *

 *
manages text reflow upon resize
 *
manages underlying disk storage for very old scrollback history lines.
 *

 *
 *
Layout

 *
 *

 *      +0========================-3+   <-- scrollback top
 *      |1                        -2|
 *      |2   Scrollback history   -1|
 *      |3                         0|   <-- scrollback bottom
 *      +4-------------------------1+   <-- main page top
 *      |5                         2|
 *      |6   main page area        3|
 *      |7                         4|   <-- main page bottom
 *      +---------------------------+
 *       ^                          ^
 *       1                          pageSize.columns
 *

 */
template <typename Cell>
class Grid
{
    // TODO: Rename all "History" to "Scrollback"?
  public:
    Grid(PageSize _pageSize, bool _reflowOnResize, LineCount _maxHistoryLineCount);

    Grid(): Grid(PageSize { LineCount(25), ColumnCount(80) }, false, LineCount(0)) {}

    void reset();

    // {{{ grid global properties
    [[nodiscard]] LineCount maxHistoryLineCount() const noexcept { return maxHistoryLineCount_; }
    void setMaxHistoryLineCount(LineCount _maxHistoryLineCount);

    [[nodiscard]] LineCount totalLineCount() const noexcept { return maxHistoryLineCount_ + pageSize_.lines; }

    [[nodiscard]] LineCount historyLineCount() const noexcept
    {
        return std::min(maxHistoryLineCount_, linesUsed_ - pageSize_.lines);
    }

    [[nodiscard]] bool reflowOnResize() const noexcept { return reflowOnResize_; }
    void setReflowOnResize(bool _enabled) { reflowOnResize_ = _enabled; }

    [[nodiscard]] PageSize pageSize() const noexcept { return pageSize_; }

    /// Resizes the main page area of the grid and adapts the scrollback area's width accordingly.
    ///
    /// @param _pageSize          new size of the main page area
    /// @param _currentCursorPos  current cursor position
    /// @param _wrapPending       AutoWrap is on and a wrap is pending
    ///
    /// @returns updated cursor position.
    [[nodiscard]] CellLocation resize(PageSize _pageSize, CellLocation _currentCursorPos, bool _wrapPending);
    // }}}

    // {{{ Line API
    /// @returns reference to Line at given relative offset @p _line.
    Line<Cell>& lineAt(LineOffset _line) noexcept;
    Line<Cell> const& lineAt(LineOffset _line) const noexcept;

    gsl::span<Cell const> lineBuffer(LineOffset _line) const noexcept { return lineAt(_line).cells(); }
    gsl::span<Cell const> lineBufferRightTrimmed(LineOffset _line) const noexcept;

    [[nodiscard]] std::string lineText(LineOffset _line) const;
    [[nodiscard]] std::string lineTextTrimmed(LineOffset _line) const;
    [[nodiscard]] std::string lineText(Line<Cell> const& _line) const;

    void setLineText(LineOffset _line, std::string_view _text);

    // void resetLine(LineOffset _line, GraphicsAttributes _attribs) noexcept
    // { lineAt(_line).reset(_attribs); }

    [[nodiscard]] ColumnCount lineLength(LineOffset _line) const noexcept { return lineAt(_line).size(); }
    [[nodiscard]] bool isLineBlank(LineOffset _line) const noexcept;
    [[nodiscard]] bool isLineWrapped(LineOffset _line) const noexcept;

    [[nodiscard]] int computeLogicalLineNumberFromBottom(LineCount _n) const noexcept;

    [[nodiscard]] size_t zero_index() const noexcept { return lines_.zero_index(); }
    // }}}

    /// Gets a reference to the cell relative to screen origin (top left, 0:0).
    [[nodiscard]] Cell& useCellAt(LineOffset _line, ColumnOffset _column) noexcept;
    [[nodiscard]] Cell& at(LineOffset _line, ColumnOffset _column) noexcept;
    [[nodiscard]] Cell const& at(LineOffset _line, ColumnOffset _column) const noexcept;

    // page view API
    gsl::span<Line<Cell>> pageAtScrollOffset(ScrollOffset _scrollOffset);
    gsl::span<Line<Cell> const> pageAtScrollOffset(ScrollOffset _scrollOffset) const;
    gsl::span<Line<Cell>> mainPage();
    gsl::span<Line<Cell> const> mainPage() const;

    LogicalLines<Cell> logicalLines()
    {
        return LogicalLines<Cell> { boxed_cast<LineOffset>(-historyLineCount()),
                                    boxed_cast<LineOffset>(pageSize_.lines - 1),
                                    lines_ };
    }

    ReverseLogicalLines<Cell> logicalLinesReverse()
    {
        return ReverseLogicalLines<Cell> { boxed_cast<LineOffset>(-historyLineCount()),
                                           boxed_cast<LineOffset>(pageSize_.lines - 1),
                                           lines_ };
    }

    // {{{ buffer manipulation

    /// Completely deletes all scrollback lines.
    void clearHistory();

    /// Scrolls up by @p _n lines within the given margin.
    ///
    /// @param _n number of lines to scroll up within the given margin.
    /// @param _defaultAttributes SGR attributes the newly created grid cells will be initialized with.
    /// @param _margin the margin coordinates to perform the scrolling action into.
    LineCount scrollUp(LineCount _n, GraphicsAttributes _defaultAttributes, Margin _margin) noexcept;

    /// Scrolls up main page by @p _n lines and re-initializes grid cells with @p _defaultAttributes.
    LineCount scrollUp(LineCount _n, GraphicsAttributes _defaultAttributes = {}) noexcept;

    /// Scrolls down by @p _n lines within the given margin.
    ///
    /// @param _n number of lines to scroll down within the given margin.
    /// @param _defaultAttributes SGR attributes the newly created grid cells will be initialized with.
    /// @param _margin the margin coordinates to perform the scrolling action into.
    void scrollDown(LineCount _n, GraphicsAttributes const& _defaultAttributes, Margin const& _margin);

    // Scrolls the data within the margins to the left filling the new space on the right with empty cells.
    void scrollLeft(GraphicsAttributes _defaultAttributes, Margin _margin) noexcept;
    // }}}

    // {{{ Rendering API
    /// Renders the full screen by passing every grid cell to the callback.
    template <typename RendererT>
    void render(RendererT&& _render, ScrollOffset _scrollOffset = {}) const;

    /// Takes text-screenshot of the main page.
    [[nodiscard]] std::string renderMainPageText() const;

    /// Renders the full grid's text characters.
    ///
    /// Empty cells are represented as strings and lines split by LF.
    [[nodiscard]] std::string renderAllText() const;
    // }}}

    [[nodiscard]] constexpr LineFlags defaultLineFlags() const noexcept;

    [[nodiscard]] constexpr LineCount linesUsed() const noexcept;

    void verifyState() const;

  private:
    CellLocation growLines(LineCount _newHeight, CellLocation _cursor);
    void appendNewLines(LineCount _count, GraphicsAttributes _attr);
    void clampHistory();

    // {{{ buffer helpers
    void resizeBuffers(PageSize _newSize)
    {
        auto const newTotalLineCount = historyLineCount() + _newSize.lines;
        lines_.resize(unbox<size_t>(newTotalLineCount));
        pageSize_ = _newSize;
    }

    void rezeroBuffers() noexcept { lines_.rezero(); }

    void rotateBuffers(int offset) noexcept { lines_.rotate(offset); }

    void rotateBuffersLeft(LineCount count) noexcept { lines_.rotate_left(unbox<size_t>(count)); }

    void rotateBuffersRight(LineCount count) noexcept { lines_.rotate_right(unbox<size_t>(count)); }
    // }}}

    // private fields
    //
    PageSize pageSize_;
    bool reflowOnResize_ = false;
    LineCount maxHistoryLineCount_;

    // Number of lines is at least the sum of maxHistoryLineCount_ + pageSize_.lines,
    // because shrinking the page height does not necessarily
    // have to resize the array (as optimization).
    Lines<Cell> lines_;

    // Number of lines used in the Lines buffer.
    LineCount linesUsed_;
};

template <typename Cell>
std::ostream& dumpGrid(std::ostream& os, Grid<Cell> const& grid);

template <typename Cell>
std::string dumpGrid(Grid<Cell> const& grid);

// {{{ impl
template <typename Cell>
constexpr LineFlags Grid<Cell>::defaultLineFlags() const noexcept
{
    return reflowOnResize_ ? LineFlags::Wrappable : LineFlags::None;
}

template <typename Cell>
constexpr LineCount Grid<Cell>::linesUsed() const noexcept
{
    return linesUsed_;
}

template <typename Cell>
bool Grid<Cell>::isLineWrapped(LineOffset _line) const noexcept
{
    return _line >= -boxed_cast<LineOffset>(historyLineCount())
           && boxed_cast<LineCount>(_line) < pageSize_.lines && lineAt(_line).wrapped();
}

template <typename Cell>
template <typename RendererT>
void Grid<Cell>::render(RendererT&& _render, ScrollOffset _scrollOffset) const
{
    assert(!_scrollOffset || unbox<LineCount>(_scrollOffset) <= historyLineCount());

    auto y = LineOffset(0);
    for (int i = -*_scrollOffset, e = i + *pageSize_.lines; i != e; ++i, ++y)
    {
        auto x = ColumnOffset(0);
        Line<Cell> const& line = lines_[i];
        if (line.isTrivialBuffer())
            _render.renderTrivialLine(line.trivialBuffer(), y);
        else
        {
            _render.startLine(y);
            for (Cell const& cell: line.cells())
                _render.renderCell(cell, y, x++);
            _render.endLine();
        }
    }
    _render.finish();
}
// }}}

} // namespace terminal

// {{{ fmt formatter
namespace fmt
{

template <>
struct formatter<terminal::Margin::Horizontal>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(const terminal::Margin::Horizontal range, FormatContext& ctx)
    {
        return fmt::format_to(ctx.out(), "{}..{}", range.from, range.to);
    }
};

template <>
struct formatter<terminal::Margin::Vertical>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(const terminal::Margin::Vertical range, FormatContext& ctx)
    {
        return fmt::format_to(ctx.out(), "{}..{}", range.from, range.to);
    }
};

} // namespace fmt
// }}}

max_stars_repo_name: christianparpart/libterminal

max_issues_repo_name: christianparpart/libterminal

max_forks_repo_name: christianparpart/libterminal