# Входные параметры: широта, долгота, уровень
proc Geo:LatLon2XY { lat lon z } {
  set pi 3.1415926535897932
  set latrad [expr { $lat*$pi/180.0 }]
  set lonrad [expr { $lon*$pi/180.0 }]

  set bm0    [expr { 256 * 2**$z/2 }]   <-------------------------------------- !!!!
  set flatX  [expr { int($bm0 * (1 + $lonrad/$pi)) }]
  set flatY  [expr { int($bm0 * (1 - 0.5*log((1+sin($latrad))/(1-sin($latrad)))/$pi)) }]  
  set tilesX [expr { int($flatX/256) }]
  set tilesY [expr { int($flatY/256) }]
  set picX   [expr { $flatX % 256 }]
  set picY   [expr { $flatY % 256 }]

  # Выходные параметры: 
  # номера тайлов по X и Y
  # координаты точки в тайле
  # координаты точки в плоскости проэкции Земли
  return [list $tilesX $tilesY $picX $picY $flatX $flatY]
}

        static void  LonLat2XY (out uint tilesX
                              , out uint tilesY
                              , out uint picX
                              , out uint picY
                              , out double flatX
                              , out double flatY
                              , double lat
                              , double lon
                              , short  zoom
                              ){
           double pi = 3.1415926535897932;
           double latrad = (lat*pi)/180.0;
           l.WriteLine(IMPORTANCELEVEL.Debug,"lat  radian/input {0}/{1}",latrad, lat);
           double lonrad = (lon*pi)/180.0;

         //  double a = 6378137.0;
         //  double k = 0.0818191908426;
           double bm0  = 256 *  Math.Pow( 2.0, (double)(zoom/2));    <---------------------------------------------
           l.WriteLine(IMPORTANCELEVEL.Debug,"bm0/zoom {0}/{1}",bm0, zoom);

           flatX = (bm0 * (1+ lonrad/pi));
           flatY = (bm0 * (1- 0.5 * Math.Log(
                            (1+ Math.Sin(latrad))
                                   /
                                    (1-Math.Sin(latrad)) 
                                           )/pi));
           l.WriteLine(IMPORTANCELEVEL.Debug,"flatY {0}/{1}",flatY, coor2str(flatY));
           tilesX = (uint)(flatX)/256;
           tilesY = (uint)(flatY)/256;
           picX  =  (uint)(flatX)%256;
           picY  =  (uint)(flatY)%256;
        
        }

//           double bm0  = 256 *  Math.Pow( 2.0, (double)(zoom/2));
           double bm0  = 256 *  Math.Pow( 2.0, (double)zoom)/2.0;

// on, lat
//  , double& X_m, double& Y_m  -- 	целая часть это номер тайла, дробная часть умноженная на 
//  256 это номер пикселя в тайле
//
void tiled_map::CoordToMercator(double X, double Y, double& X_m, double& Y_m,int level)
{
	double bm = bm0(level);
    X_m = bm * (1 + X * rad() / M_PI);
    Y_m = bm * (1 - 0.5 * log((1 + sin(Y * rad())) / (1 - sin(Y * rad()))) / M_PI);
}


/* 

http://wiki.openstreetmap.org/wiki/Mercator#C

double y2lat_d(double y) { return rad2deg(2 * atan(exp(  deg2rad(y) ) ) - M_PI/2); }
double x2lon_d(double x) { return x; }
double lat2y_d(double lat) { return rad2deg(log(tan(M_PI/4+ deg2rad(lat)/2))); }
double lon2x_d(double lon) { return lon; }
 
// The following functions take or return there results in something close to meters, along the equator 
// Это не настоящие метры, а "меркаторовские". С настоящими совпадает только на экваторе, для прочих широт для получения настоящих метров надо длину отрезка домножить на косинус широты (годится для небольших отрезков, для больших будет слишком большая погрешность).
double y2lat_m(double y) { return rad2deg(2 * atan(exp( (y / earth_radius ) )) - M_PI/2); }
double x2lon_m(double x) { return rad2deg(x / earth_radius); }
double lat2y_m(double lat) { return earth_radius * log(tan(M_PI/4+ deg2rad(lat)/2)); }
double lon2x_m(double lon) { return deg2rad(lon) * earth_radius; }
*/

void tiled_map::MercatorToCoord(double X_m, double Y_m, double& X, double& Y,int level)
{
	double bm = bm0(level);
    X = M_PI * (X_m / bm - 1) / rad();
    double tmp = pow(M_E, 2 * M_PI * (1 - Y_m / bm));
    Y = asin((tmp - 1) / (tmp + 1)) / rad();
}

#
# Gilbert (Peano) lines
# Artix, 2011-01-22
# Niklaus Wirth algorithm
#

from Tkinter import *
from math import *
import sys

SIZE_X=600
SIZE_Y=500
STEP=6

class Turtle:
    x = 0
    y = 0
    a = 0
    pen = 1
    canvas = None

    def __init__(self,canvas):
        self.canvas = canvas

    def setpen(self,yes):
        self.pen = yes

    def rotate(self,g):
        self.a += g

    def move(self,sx,sy):
        self.x = sx
        self.y = sy

    def forward(self,length):
        new_x = self.x+cos(pi/180*self.a)*length
        new_y = self.y+sin(pi/180*self.a)*length
        cx = SIZE_X / 2
        cy = SIZE_Y / 2
        if self.pen:
            self.canvas.create_line(cx+self.x, cy-self.y, cx+new_x, cy-new_y,fill="green")
        self.x = new_x
        self.y = new_y

    def line(self,angle,length):
        new_x = self.x+cos(pi/180*angle)*length
        new_y = self.y+sin(pi/180*angle)*length
        cx = SIZE_X / 2
        cy = SIZE_Y / 2
        if self.pen:
            self.canvas.create_line(cx+self.x, cy-self.y, cx+new_x, cy-new_y,fill="green")
        self.x = new_x
        self.y = new_y


class Peano:

    turtle = None

    def __init__(self,turtle):
        self.turtle = turtle


    def make(self,level):
        self.turtle.move(0.5*STEP*pow(2,level)-1, -0.5*STEP*pow(2,level)-1)
        self.A(level)

    def A(self,i):
        if i>0:
            self.D(i-1)
            self.turtle.line(180,STEP)
            self.A(i-1)
            self.turtle.line(90,STEP)
            self.A(i-1)
            self.turtle.line(0,STEP)
            self.B(i-1)

    def B(self,i):
        if i>0:
            self.C(i-1)
            self.turtle.line(-90,STEP)
            self.B(i-1)
            self.turtle.line(0,STEP)
            self.B(i-1)
            self.turtle.line(90,STEP)
            self.A(i-1)

    def C(self,i):
        if i>0:
            self.B(i-1)
            self.turtle.line(0,STEP)
            self.C(i-1)
            self.turtle.line(-90,STEP)
            self.C(i-1)
            self.turtle.line(-180,STEP)
            self.D(i-1)

    def D(self,i):
        if i>0:
            self.A(i-1)
            self.turtle.line(90,STEP)
            self.D(i-1)
            self.turtle.line(180,STEP)
            self.D(i-1)
            self.turtle.line(270,STEP)
            self.C(i-1)

root = Tk()
root.title("Peano")
c = Canvas(root,width=SIZE_X,height=SIZE_Y,bg="black")
c.pack()
peano = Peano(Turtle(c))
peano.make(6)

root.mainloop()

#define C 16807
#define A 2147483647.0
double drand48();
void srand48(long int);
double yz;

void srand48(long int seed)
{
yz = (double)seed;
}

double drand48()
{
long int ki;
double uu;
ki=(C* yz)/A;
yz=C* yz-ki*A;
uu=yz/(A-1);
return uu;
}

/* --- Copyright Martin Birgmeier 1996. All rights reserved. --------------
 * File:			C.win32/extern/src/rand48.c
 * Purpose:			BSD random number generator
 * Author:			Robert Duncan, May 31 1996
 */

/************************************************************************
 *                                                                      *
 * Copyright (c) 1993 Martin Birgmeier                                  *
 * All rights reserved.                                                 *
 *                                                                      *
 * You may redistribute unmodified or modified versions of this source  *
 * code provided that the above copyright notice and this and the       *
 * following conditions are retained.                                   *
 *                                                                      *
 * This software is provided ``as is'', and comes with no warranties    *
 * of any kind. I shall in no event be liable for anything that happens *
 * to anyone/anything when using this software.                         *
 *                                                                      *
 ************************************************************************/

#include <math.h>
#include <stdlib.h>

#define RAND48_SEED_0   (0x330e)
#define RAND48_SEED_1   (0xabcd)
#define RAND48_SEED_2   (0x1234)
#define RAND48_MULT_0   (0xe66d)
#define RAND48_MULT_1   (0xdeec)
#define RAND48_MULT_2   (0x0005)
#define RAND48_ADD      (0x000b)

unsigned short _rand48_seed[3] = {
		RAND48_SEED_0,
		RAND48_SEED_1,
		RAND48_SEED_2
};
unsigned short _rand48_mult[3] = {
		RAND48_MULT_0,
		RAND48_MULT_1,
		RAND48_MULT_2
};
unsigned short _rand48_add = RAND48_ADD;

void
_dorand48(unsigned short xseed[3])
{
		unsigned long accu;
		unsigned short temp[2];

		accu = (unsigned long) _rand48_mult[0] * (unsigned long) xseed[0] +
		 (unsigned long) _rand48_add;
		temp[0] = (unsigned short) accu;        /* lower 16 bits */
		accu >>= sizeof(unsigned short) * 8;
		accu += (unsigned long) _rand48_mult[0] * (unsigned long) xseed[1] +
		 (unsigned long) _rand48_mult[1] * (unsigned long) xseed[0];
		temp[1] = (unsigned short) accu;        /* middle 16 bits */
		accu >>= sizeof(unsigned short) * 8;
		accu += _rand48_mult[0] * xseed[2] + _rand48_mult[1] * xseed[1] + _rand48_mult[2] * xseed[0];
		xseed[0] = temp[0];
		xseed[1] = temp[1];
		xseed[2] = (unsigned short) accu;
}

double
erand48(unsigned short xseed[3])
{
		_dorand48(xseed);
		return ldexp((double) xseed[0], -48) +
			   ldexp((double) xseed[1], -32) +
			   ldexp((double) xseed[2], -16);
}

double
drand48(void)
{
		return erand48(_rand48_seed);
}

long
lrand48(void)
{
		_dorand48(_rand48_seed);
		return ((long) _rand48_seed[2] << 15) + ((long) _rand48_seed[1] >> 1);
}

long
nrand48(unsigned short xseed[3])
{
		_dorand48(xseed);
		return ((long) xseed[2] << 15) + ((long) xseed[1] >> 1);
}

long
mrand48(void)
{
		_dorand48(_rand48_seed);
		return ((long) _rand48_seed[2] << 16) + (long) _rand48_seed[1];
}

long
jrand48(unsigned short xseed[3])
{
		_dorand48(xseed);
		return ((long) xseed[2] << 16) + (long) xseed[1];
}

void
srand48(long seed)
{
		_rand48_seed[0] = RAND48_SEED_0;
		_rand48_seed[1] = (unsigned short) seed;
		_rand48_seed[2] = (unsigned short) (seed >> 16);
		_rand48_mult[0] = RAND48_MULT_0;
		_rand48_mult[1] = RAND48_MULT_1;
		_rand48_mult[2] = RAND48_MULT_2;
		_rand48_add = RAND48_ADD;
}

unsigned short *
seed48(unsigned short xseed[3])
{
		static unsigned short sseed[3];

		sseed[0] = _rand48_seed[0];
		sseed[1] = _rand48_seed[1];
		sseed[2] = _rand48_seed[2];
		_rand48_seed[0] = xseed[0];
		_rand48_seed[1] = xseed[1];
		_rand48_seed[2] = xseed[2];
		_rand48_mult[0] = RAND48_MULT_0;
		_rand48_mult[1] = RAND48_MULT_1;
		_rand48_mult[2] = RAND48_MULT_2;
		_rand48_add = RAND48_ADD;
		return sseed;
}

void
lcong48(unsigned short p[7])
{
		_rand48_seed[0] = p[0];
		_rand48_seed[1] = p[1];
		_rand48_seed[2] = p[2];
		_rand48_mult[0] = p[3];
		_rand48_mult[1] = p[4];
		_rand48_mult[2] = p[5];
		_rand48_add = p[6];
}

self.turtle.line(-90,STEP)
            self.B(i-1)
            self.turtle.line(0,STEP)

/*****************************************************************
 * hilbert_i2c
 * 
 * Convert an index into a Hilbert curve to a set of coordinates.
 * Inputs:
 *  nDims:      Number of coordinate axes.
 *  nBits:      Number of bits per axis.
 *  index:      The index, contains nDims*nBits bits (so nDims*nBits must be <= 8*sizeof(bitmask_t)).
 * Outputs:
 *  coord:      The list of nDims coordinates, each with nBits bits.
 * Assumptions:
 *      nDims*nBits <= (sizeof index) * (bits_per_byte)
 */

void hilbert_i2c(unsigned nDims, unsigned nBits, bitmask_t index, bitmask_t coord[]);

/*****************************************************************
 * hilbert_c2i
 * 
 * Convert coordinates of a point on a Hilbert curve to its index.
 * Inputs:
 *  nDims:      Number of coordinates.
 *  nBits:      Number of bits/coordinate.
 *  coord:      Array of n nBits-bit coordinates.
 * Outputs:
 *  index:      Output index value.  nDims*nBits bits.
 * Assumptions:
 *      nDims*nBits <= (sizeof bitmask_t) * (bits_per_byte)
 */

bitmask_t hilbert_c2i(unsigned nDims, unsigned nBits, bitmask_t const coord[]);