/* This is a library of commonly used mathmatical equations and idea that
   I have come across*/


/* snoise(x) is taking noise from = 0 - 1 to 1 - -1 */
float snoise (point p)
{
	 float value = (2 * noise(p) - 1);
	
return value;

}

/* applies noise to the x,y,z values */
vector vsnoise(point p)
{
	vector value = (2 * (vector noise(p))-1);
	
return value;

}

/* Basic fractional Brownian motion */
float fBm(point p; float octaves, lacunarity, gain)
{
	float value = 0;
	float amp = 1;
	point pp = p;
	float i;
	
	for(i = 0; i < octaves; i +=1)
	{
		value += snoise(pp) * pow(lacunarity, -amp * i);
		amp *= gain; pp *= lacunarity;
	}
	
	return value;
}

/* Basic fractional Brownian motion vector */
vector vfBm(point p; float octaves, lacunarity, gain)
{
	vector value = 0;
	float amp = 1;
	point pp = p;
	float i;
	
	for(i = 0; i < octaves; i +=1)
	{
		value += vsnoise(pp) * pow(lacunarity, -amp * i);
		amp *= gain; pp *= lacunarity;
	}
	
	return value;
}

/* fractional Brownian motion with offset */
float mfBm(point p; float octaves, lacunarity, gain,offset)
{
	float value = 1;
	float amp = 1;
	point pp = p;
	float i;
	
	for(i = 0; i < octaves; i +=1)
	{
		value *= ((noise(pp)) * offset) * pow(lacunarity, -amp * i);
		amp *= gain; pp *= lacunarity;
	}
	
	return value;
}

/* turbulance fractional Brownian motion */
float turb(point p; float octaves, lacunarity, gain)
{
	float value = 0;
	float amp = 1;
	point pp = p;
	float i;
	
	for(i = 0; i < octaves; i +=1)
	{
		value += abs(snoise(pp)) * pow(lacunarity, -amp * i);
		amp *= gain; pp *= lacunarity;
	}
	
	return value;
}


/* Voronoi cell noise (a.k.a. Worley noise) -- 3-D, 1-feature version. */
void
voronoi_f1_3d (point P; float jitter; output float f1; output point pos1;)
{
    point thiscell = point (floor(xcomp(P))+0.5, floor(ycomp(P))+0.5,
			    floor(zcomp(P))+0.5);
    f1 = 1000;
    uniform float i, j, k;
    for (i = -1;  i <= 1;  i += 1) {
        for (j = -1;  j <= 1;  j += 1) {
            for (k = -1;  k <= 1;  k += 1) {
		point testcell = thiscell + vector(i,j,k);
                point pos = testcell + 
		            jitter * (vector cellnoise (testcell) - 0.5);
		vector offset = pos - P;
                float dist = offset . offset; /* actually dist^2 */
                if (dist < f1) {
                    f1 = dist;  pos1 = pos;
                }
            }
	}
    }
    f1 = sqrt(f1);
}


/* Voronoi cell noise (a.k.a. Worley noise) -- 3-D, 2-feature version. */
void
voronoi_f1f2_3d (point P; float jitter; output float f1; output point pos1; output float f2; output point pos2;float seed;)
{
    point thiscell = point (floor(xcomp(P))+0.5, floor(ycomp(P))+0.5,
			    floor(zcomp(P))+0.5);
    f1 = f2 = 1000;
    uniform float i, j, k;
    for (i = -1;  i <= 1;  i += 1) {
        for (j = -1;  j <= 1;  j += 1) {
            for (k = -1;  k <= 1;  k += 1) {
		point testcell = thiscell + vector(i,j,k);
                point pos = testcell + 
		            jitter * (vector cellnoise ((testcell),seed) - 0.5);
		vector offset = pos - P;
                float dist = offset . offset;
                if (dist < f1) {
                    f2 = f1;  pos2 = pos1;
                    f1 = dist;  pos1 = pos;
                } else if (dist < f2) {
                    f2 = dist;  pos2 = pos;
		}
            }
	}
    }
    f1 = sqrt(f1);  f2 = sqrt(f2);
}


/* Voronoi cell noise (a.k.a. Worley noise) -- 2-D, 1-feature version. */
void
voronoi_f1_2d (float ss, tt; float jitter; output float f1; output float spos1, tpos1;)
{
    float sthiscell = floor(ss)+0.5, tthiscell = floor(tt)+0.5;
    f1 = 1000;
    uniform float i, j;
    for (i = -1;  i <= 1;  i += 1) {
	float stestcell = sthiscell + i;
        for (j = -1;  j <= 1;  j += 1) {
	    float ttestcell = tthiscell + j;
	    float spos = stestcell +
		     jitter * (float cellnoise(stestcell, ttestcell) - 0.5);
	    float tpos = ttestcell +
		 jitter * (float cellnoise(stestcell+23, ttestcell-87) - 0.5);
	    float soffset = spos - ss;
	    float toffset = tpos - tt;
	    float dist = soffset*soffset + toffset*toffset;
	    if (dist < f1) {
		f1 = dist;
		spos1 = spos;  tpos1 = tpos;
	    }
	}
    }
    f1 = sqrt(f1);
}


/* Voronoi cell noise (a.k.a. Worley noise) -- 2-D, 2-feature version. */
void
voronoi_f1f2_2d (float ss, tt; float jitter; output float f1; output float spos1, tpos1; output float f2; output float spos2, tpos2;)
{
    float sthiscell = floor(ss)+0.5, tthiscell = floor(tt)+0.5;
    f1 = f2 = 1000;
    uniform float i, j;
    for (i = -1;  i <= 1;  i += 1) {
	float stestcell = sthiscell + i;
        for (j = -1;  j <= 1;  j += 1) {
	    float ttestcell = tthiscell + j;
	    float spos = stestcell +
		   jitter * (cellnoise(stestcell, ttestcell) - 0.5);
	    float tpos = ttestcell +
		   jitter * (cellnoise(stestcell+23, ttestcell-87) - 0.5);
	    float soffset = spos - ss;
	    float toffset = tpos - tt;
	    float dist = soffset*soffset + toffset*toffset;
	    if (dist < f1) {
		f2 = f1;  spos2 = spos1;  tpos2 = tpos1;
		f1 = dist;  spos1 = spos;  tpos1 = tpos;
	    } else if (dist < f2) {
		f2 = dist;
		spos2 = spos;  tpos2 = tpos;
	    }
	}
    }
    f1 = sqrt(f1);  f2 = sqrt(f2);
}