{-------------------------------------------------------------------------
 Module:    OSM Tile Coordinate Math (Shared)

 Comment:   Pure functions for converting between geographic coordinates
            (latitude/longitude), Mercator Y, and OSM tile indices.
            Shared between TGlobe4, TGlobe5, and other consumers.
-------------------------------------------------------------------------}
{$IFDEF FPC}{$MODE DELPHI}{$ENDIF}
{$WEAKPACKAGEUNIT ON}
unit uCICOSMTileMath;

interface

uses
  Math, uCICOSMTileCommon;

function  OSMMercatorY(LatDeg: Double): Double;
function  OSMMercatorYInv(Y: Double): Double;
procedure OSMLatLonToTile(Lat, Lon: Double; Z: Integer; out TX, TY: Integer);
procedure OSMTileToLatLon(Z, TX, TY: Integer;
  out LatT, LatB, LonL, LonR: Double);
function  OSMScaleFactorToZoom(ScaleFactor: Extended;
  GlobeUnitsPerDeg: Integer; MaxZoom: Integer): Integer;

implementation

function OSMMercatorY(LatDeg: Double): Double;
var
  LatRad: Double;
begin
  LatRad := LatDeg * (Pi / 180.0);
  Result := Ln(Tan(Pi / 4.0 + LatRad / 2.0));
end;

function OSMMercatorYInv(Y: Double): Double;
begin
  Result := (2.0 * ArcTan(Exp(Y)) - Pi / 2.0) * (180.0 / Pi);
end;

procedure OSMLatLonToTile(Lat, Lon: Double; Z: Integer; out TX, TY: Integer);
var
  N    : Double;
  LatR : Double;
  MaxN : Integer;
begin
  if Lat > OSM_LAT_LIMIT then
    Lat := OSM_LAT_LIMIT
  else if Lat < -OSM_LAT_LIMIT then
    Lat := -OSM_LAT_LIMIT;

  N    := Power(2.0, Z);
  LatR := Lat * (Pi / 180.0);

  TX := Trunc((Lon + 180.0) / 360.0 * N);
  TY := Trunc((1.0 - Ln(Tan(LatR) + 1.0 / Cos(LatR)) / Pi) / 2.0 * N);

  MaxN := Round(N) - 1;

  if TX < 0 then TX := 0
  else if TX > MaxN then TX := MaxN;

  if TY < 0 then TY := 0
  else if TY > MaxN then TY := MaxN;
end;

procedure OSMTileToLatLon(Z, TX, TY: Integer;
  out LatT, LatB, LonL, LonR: Double);
var
  N: Double;
begin
  N    := Power(2.0, Z);
  LonL := TX / N * 360.0 - 180.0;
  LonR := (TX + 1) / N * 360.0 - 180.0;
  LatT := OSMMercatorYInv(Pi - 2.0 * Pi * TY / N);
  LatB := OSMMercatorYInv(Pi - 2.0 * Pi * (TY + 1) / N);
end;

function OSMScaleFactorToZoom(ScaleFactor: Extended;
  GlobeUnitsPerDeg: Integer; MaxZoom: Integer): Integer;
var
  Z: Integer;
begin
  if ScaleFactor <= 0 then
  begin
    Result := 0;
    Exit;
  end;
  Z := Round(Log2(ScaleFactor * 360.0 * GlobeUnitsPerDeg / OSM_TILE_SIZE));
  if Z < 0 then
    Z := 0
  else if Z > MaxZoom then
    Z := MaxZoom;
  Result := Z;
end;

end.