float2 xyzToLatlon(float3 pos) {
	pos = normalize(pos);
	return float2(
		acos(pos.y),
		atan2(pos.z, pos.x)
	);
}

float3 raySphereIntersection(float3 pos, float3 dir, float radius) {
	float radiusSq = radius * radius;
	float3 L = pos;
	float tca = dot(dir, L);

	float L_2 = dot(L, L);
	if(tca <= 0 && L_2 > radiusSq) return float3(0, 0, -1);

	float d = L_2 - tca * tca;
	if(d >= radiusSq) return float3(0, 0, -1);

	float thc = radiusSq - d;
	thc = sqrt(thc);

	return float3(tca - thc, tca + thc, 1);
}

float biasFunction(float x, float bias) {
	float k = 1.0 - bias;
	k = k * k * k;
	return (x * k) / (x * k - x + 1.0);
}

float gcDistance(float2 LL1, float2 LL2) {
	float diffLong = abs(LL1.y - LL2.y);
	float a = sin(LL1.x) * sin(LL2.x);
	float b = cos(LL1.x) * cos(LL2.x) * cos(diffLong);
	float dist = acos(a + b) / UNITY_PI;
	return dist;
}

float3 rayCylinderIntersection(float3 pos, float3 orig, float2 dir, float radius) {
	//dir = normalize(dir);
	float3 L = pos - orig;
	float tca = dir.x * L.x + dir.y * L.y;
	float L_2 = L.x * L.x + L.y * L.y;
	float radiusSq = radius * radius;
	if(tca <= 0 && L_2 >= radiusSq) return float3(0, 0, -1);
	float d = L_2 - tca * tca;
	if(d >= radiusSq) return float3(0, 0, -1);
	float thc = radiusSq - d;
	thc = sqrt(thc);
	float angle = acos(dir.x / length(dir));
	d = tca - thc;
	return float3(d, angle, 1);
}

uint RNG_at(uint globalSeedHi, uint globalSeedLo, int x, int y, int z, int iters) {
		uint k,k2,k3;
		k = (x << 4) + y * 12 + z * 2421;
		k2 = (y << 4) + x * 99 + z * 1111;
		k3 = (z << 4) + (x * 2425) - y * 24525;
		uint iseed1 = (globalSeedHi + k3 * 19438) * (globalSeedLo - k * 35741) - k2 * 333;
		uint iseed2 = (globalSeedLo + k3 * 13145) * (globalSeedHi - k * 132532) - k2 * 4623;
		uint iseed = iseed1 ^ iseed2;

		[unroll(8)]
		for(int i = 0; i < iters; i++) {
				uint x = iseed;
				iseed ^= x << 13;
				iseed ^= x >> 17;
				iseed ^= x << 5;
		}
		return iseed;
}

float3 RNG3_at(uint globalSeedHi, uint globalSeedLo, int x, int y, int z, int iters) {
		int random = RNG_at(globalSeedHi, globalSeedLo, x + 948, y + 958, z, 3) & 65535;
		float aaaa = (float)random / 65536.0;
		aaaa = -1 + 2*aaaa;
		random = RNG_at(globalSeedHi * 17, globalSeedLo * 11, x + 98458, y, z + 95871, 3) & 65535;
		float bbbb = (float)random / 65536.0;
		bbbb = -1 + 2*bbbb;
		random = RNG_at(globalSeedHi * 31, globalSeedLo * 5, x + 9068, y + 9588, z, 3) & 65535;
		float cccc = (float)random / 65536.0;
		cccc = -1 + 2*cccc;
		return float3(aaaa, bbbb, cccc);
}

float perlin(float x, float y, float z) {
		float x_int_part = floor(x);
		float y_int_part = floor(y);
		float z_int_part = floor(z);

		float3 fract_part = float3(frac(x), frac(y), frac(z));
		float3 interp = fract_part*fract_part*(3-2*fract_part);

		float2 f1 = RNG3_at(3598935778, 91588935, x_int_part, y_int_part, z_int_part, 3);
		float2 f2 = RNG3_at(3598935778, 91588935, x_int_part+1, y_int_part, z_int_part, 3);
		float2 f3 = RNG3_at(3598935778, 91588935, x_int_part, y_int_part+1, z_int_part, 3);
		float2 f4 = RNG3_at(3598935778, 91588935, x_int_part+1, y_int_part+1, z_int_part, 3);

		float fr = interp.x;
		float f5 = dot(f1, fract_part) * (1-fr) + dot(f2, fract_part-float2(1,0)) * fr;
		float f6 = dot(f3, fract_part-float2(0,1)) * (1-fr) + dot(f4, fract_part-float2(1,1)) * fr;
		fr = interp.y;
		float f = f5 * (1-fr) + f6 * fr;

		f1 = RNG3_at(3598935778, 91588935, x_int_part, y_int_part, z_int_part+1, 3);
		f2 = RNG3_at(3598935778, 91588935, x_int_part+1, y_int_part, z_int_part+1, 3);
		f3 = RNG3_at(3598935778, 91588935, x_int_part, y_int_part+1, z_int_part+1, 3);
		f4 = RNG3_at(3598935778, 91588935, x_int_part+1, y_int_part+1, z_int_part+1, 3);
		fr = interp.x;
		f5 = dot(f1, fract_part) * (1-fr) + dot(f2, fract_part-float2(1,0)) * fr;
		f6 = dot(f3, fract_part-float2(0,1)) * (1-fr) + dot(f4, fract_part-float2(1,1)) * fr;
		fr = interp.y;
		float f7 = f5 * (1-fr) + f6 * fr;

		fr = interp.z;
		return f * (1-fr) + f7 * fr;
}

float perlinOctaves(float x, float y, float z, int num) {
		float f = 0;
		float ca = 1;
		float cb = 1;
		float div = 0;
		for(int i = 0; i < num; i++) {
				f += perlin(x * ca, y * ca, z * ca) * cb;
				div += cb;
				cb *= 0.5;
				ca *= 2;
		}
		return f / div;
}
