static double[] cos = new double[]{
            1,0.99992,0.99968,0.99928,0.99873,0.99801,0.99714,0.9961,0.99491,0.99356,0.99205,0.99038,0.98856,0.98658,0.98444,0.98215,0.97969,0.97709,0.97432,0.9714,0.96833,0.9651,0.96172,0.95819,0.9545,0.95067,
            0.94668,0.94254,0.93825,0.93381,0.92922,0.92448,0.9196,0.91457,0.90939,0.90407,0.89861,0.893,0.88725,0.88136,0.87533,0.86916,0.86285,0.85641,0.84983,0.84311,0.83626,0.82927,0.82216,0.81491,0.80753,
            0.80003,0.79239,0.78463,0.77675,0.76874,0.76061,0.75236,0.74399,0.7355,0.72689,0.71817,0.70933,0.70038,0.69132,0.68215,0.67287,0.66349,0.65399,0.6444,0.6347,0.6249,0.615,0.605,0.59491,0.58472,
            0.57444,0.56406,0.5536,0.54305,0.53241,0.52169,0.51089,0.5,0.48903,0.47799,0.46687,0.45568,0.44441,0.43307,0.42167,0.41019,0.39865,0.38705,0.37539,0.36366,0.35188,0.34004,0.32815,0.31621,0.30421,
            0.29217,0.28008,0.26795,0.25577,0.24355,0.2313,0.21901,0.20668,0.19432,0.18193,0.1695,0.15706,0.14459,0.13209,0.11957,0.10704,0.09449,0.08192,0.06934,0.05675,0.04414,0.03154,0.01892,0.00631,-0.00631,
            -0.01892,-0.03154,-0.04414,-0.05675,-0.06934,-0.08192,-0.09449,-0.10704,-0.11957,-0.13209,-0.14459,-0.15706,-0.1695,-0.18193,-0.19432,-0.20668,-0.21901,-0.2313,-0.24355,-0.25577,-0.26795,-0.28008,-0.29217,-0.30421,-0.31621,
            -0.32815,-0.34004,-0.35188,-0.36366,-0.37539,-0.38705,-0.39865,-0.41019,-0.42167,-0.43307,-0.44441,-0.45568,-0.46687,-0.47799,-0.48903,-0.5,-0.51089,-0.52169,-0.53241,-0.54305,-0.5536,-0.56406,-0.57444,-0.58472,-0.59491,
            -0.605,-0.615,-0.6249,-0.6347,-0.6444,-0.65399,-0.66349,-0.67287,-0.68215,-0.69132,-0.70038,-0.70933,-0.71817,-0.72689,-0.7355,-0.74399,-0.75236,-0.76061,-0.76874,-0.77675,-0.78463,-0.79239,-0.80003,-0.80753,-0.81491,
            -0.82216,-0.82927,-0.83626,-0.84311,-0.84983,-0.85641,-0.86285,-0.86916,-0.87533,-0.88136,-0.88725,-0.893,-0.89861,-0.90407,-0.90939,-0.91457,-0.9196,-0.92448,-0.92922,-0.93381,-0.93825,-0.94254,-0.94668,-0.95067,-0.9545,
            -0.95819,-0.96172,-0.9651,-0.96833,-0.9714,-0.97432,-0.97709,-0.97969,-0.98215,-0.98444,-0.98658,-0.98856,-0.99038,-0.99205,-0.99356,-0.99491,-0.9961,-0.99714,-0.99801,-0.99873,-0.99928,-0.99968,-0.99992,-1};


        public static byte[, ,] GetChunk(int worldChunkX, int worldChunkZ) {
            double
                min = 0.0,
                max = 1.0,
                fx, fy, fz, fxinv, fyinv, xinterp1, xinterp2, xinterp3, xinterp4, yinterp1, yinterp2, incdiv;
            int
                maxOctaves = 4,
                seed = 62476,
                worldStartX = worldChunkX * 16,
                worldStartZ = worldChunkZ * 16,
                zoneStartX = 0,
                zoneStartZ = 0,
                count = 0,
                x, y, z, z0, z1, xBound, yBound,  depth, noisex, noisey, noisez, xSize, ySize, xySize, n_x, n_y, n1_x, n1_y, m_x, m_y, n1_z, zoneDiffX, zoneDiffZ, m_z;

            //double[] yGuide = new double[] { 1.0, -1, -1, -1, 1.0 };
            //double[] yGuide = new double[] { 0.0, -1, -1, 0.2, 1.0 };
            //double[] yGuide = new double[] { 1.0, -1, -1, -1, 0.0, -1, -1, -1, 1.0 };
            double[] yGuide = new double[] { 1.0, -1, 0.25, -1, 1.0 };
            //double[] yGuide = new double[] { 0.0, -1, 0.2, -1, 0.0 };
            //double[] yGuide = new double[] { 1.0, -1, -1, 0.4, -1, 0.3, -1, -1, 1.0 };

            //double[] cos = new double[250];
            //for (int i = 0; i <= cos.GetUpperBound(0); i++)
            //    cos[i] = Math.Cos(((double)i / (double)cos.GetUpperBound(0)) * Math.PI);


            double[] noise;
            double[] map = new double[16 * 512 * 16];
            byte[, ,] ret = new byte[16, 256, 16];

            int frequency = yGuide.GetUpperBound(0) / 2,
                increment = 256 / frequency;
            double amplitude = (max - min) / 2.0,
                midpoint = min + amplitude,
                noiseMod = midpoint;

            bool first = true, hasGrass = false; ;

            Random R;

            //TimeSpan timeNoise = new TimeSpan(0);

            while (increment > 1 && count < maxOctaves) {

                //calculate start coords for the zone covering the requested chunk at the current iteration
                zoneStartX = (int)Math.Floor((double)worldStartX / (double)increment) * increment;
                zoneStartZ = (int)Math.Floor((double)worldStartZ / (double)increment) * increment;


                xSize = (int)Math.Ceiling(16.0 / increment) + 1;
                ySize = (frequency * 2) + 1;

                xBound = xSize - 1;
                yBound = ySize - 1;
                xySize = (xSize * ySize);
                noise = new double[xSize * xySize];
                m_y = increment * 721;
                //generate noise map
                for (x = 0; x <= xBound; x++) {
                    n_x = x * xySize;
                    m_x = (zoneStartX + (x * increment)) * 13758;
                    for (z = 0; z <= xBound; z++) {
                        m_z = (zoneStartZ + (z * increment)) * 6124773;
                        for (y = 0; y <= yBound; y++) {
                            n_y = y * xSize;
                            R = new Random(seed + m_x  + (y * m_y) + m_z);
                            noise[n_x + n_y + z] = (R.NextDouble() * amplitude * (R.Next(2) == 0 ? -1 : 1)) + noiseMod;
                        }
                    }
                }

                //apply specified y guide rules
                if (first) {
                    for (y = 0; y <= yBound; y++) {
                        n_y = y * xSize;
                        if (yGuide[y] >= 0) {
                            for (x = 0; x <= xBound; x++) {
                                n_x = x * xySize;
                                m_z = n_x + n_y;
                                for (z = 0; z <= xBound; z++) {
                                    noise[m_z + z] = yGuide[y];
                                }
                            }
                        }
                    }
                    first = false;
                }
                zoneDiffX = worldStartX - zoneStartX;
                zoneDiffZ = worldStartZ - zoneStartZ;
                incdiv = (double)cos.Length / increment;
                for (x = 0; x < 16; x++) {
                    noisex = x / increment;
                    fx = (1.0 - cos[(int)(((double)((x % increment) + zoneDiffX)) * incdiv)]) * 0.5;
                    fxinv = 1.0 - fx;
                    n_x = noisex * xySize;
                    n1_x = n_x + xySize;
                    m_x = x * 8192;
                    for (y = 0; y < 256; y ++) {
                        noisey = y*2 / increment;
                        fy = (1.0 - cos[(int)(((double)(y*2 % increment)) * incdiv)]) * 0.5;
                        fyinv = 1.0 - fy;
                        n_y = noisey * xSize;
                        n1_y = n_y + xSize;
                        m_z = m_x + y*32;
                        for (noisez = 0; noisez < xBound; noisez++) {
                            n1_z = noisez + 1;
                            xinterp1 = noise[n_x + n_y + noisez] * fxinv + noise[n1_x + n_y + noisez] * fx;
                            xinterp2 = noise[n_x + n1_y + noisez] * fxinv + noise[n1_x + n1_y + noisez] * fx;
                            xinterp3 = noise[n_x + n_y + n1_z] * fxinv + noise[n1_x + n_y + n1_z] * fx;
                            xinterp4 = noise[n_x + n1_y + n1_z] * fxinv + noise[n1_x + n1_y + n1_z] * fx;

                            yinterp1 = xinterp1 * fyinv + xinterp2 * fy;
                            yinterp2 = xinterp3 * fyinv + xinterp4 * fy;

                            z0 = noisez * increment;
                            z1 = Math.Min(16, z0 + increment);

                            for (z = z0; z < z1; z++) {
                                fz = (1.0 - cos[(int)((double)((z % increment) + zoneDiffZ) * incdiv)]) * 0.5;
                                map[m_z + z] += yinterp1 * (1.0 - fz) + yinterp2 * fz;
                            }
                        }
                    }
                }

                amplitude *= 0.5; //persistence assumed to be 1.0
                count++;
                frequency *= 2;
                increment /= 2;
                noiseMod = 0;
            }


            //determine block type based on threshhold and depth
            for (x = 0; x < 16; x++) {
                m_x = x * 8192;
                for (z = 0; z < 16; z++) {
                    m_z = m_x + z;
                    depth = 0;
                    hasGrass = false;
                    for (y = 255; y >= 0; y--) {
                        m_y = (y * 2) * 16;
                        if (map[m_z + m_y] > 0.25) { //threshhold
                            ret[x, y, z] = 0;
                            depth = 0;
                        }
                        else {
                            if (++depth == 1 && !hasGrass) {
                                ret[x, y, z] = 2;
                                hasGrass = true;
                            }
                            else if (depth < 5)
                                ret[x, y, z] = 3;
                            else
                                ret[x, y, z] = 1;

                        }

                    }
                }
            }

            return ret;
        }