clear; clc;

latitude  = 50.70578714046577;
longitude = 4.392066457029492;

precision = 6;
verbose   = true;

[locator, steps] = maidenhead_from_latlon(latitude, longitude, precision, verbose);

fprintf("\nRESULT\n");
fprintf("  Maidenhead locator (%d chars): %s\n", precision, locator);

function [locator, steps] = maidenhead_from_latlon(latitude, longitude, precision, verbose)
  if (nargin < 3 || isempty(precision)), precision = 6; end
  if (nargin < 4 || isempty(verbose)),   verbose   = true; end

  steps = struct();

  % 0) Basic validation
  if (!isscalar(latitude) || !isscalar(longitude))
    error("latitude and longitude must be scalar numeric values.");
  end
  if (precision != 2 && precision != 4 && precision != 6 && precision != 8)
    error("precision must be one of: 2, 4, 6, 8");
  end

  steps.input.latitude  = latitude;
  steps.input.longitude = longitude;
  steps.precision       = precision;

  % 1) Normalize / clamp
  lon_norm = mod(longitude + 180, 360) - 180;

  lat_clamp = latitude;
  if (lat_clamp >= 90),  lat_clamp =  90 - 1e-12; end
  if (lat_clamp < -90),  lat_clamp = -90;         end

  steps.normalized.latitude  = lat_clamp;
  steps.normalized.longitude = lon_norm;

  if (verbose)
    fprintf("INPUT\n");
    fprintf("  latitude  = %.12f deg\n", latitude);
    fprintf("  longitude = %.12f deg\n", longitude);

    fprintf("\nNORMALIZE / CLAMP\n");
    fprintf("  longitude normalized to [-180, 180): %.12f deg\n", lon_norm);
    fprintf("  latitude clamped to [-90, 90):      %.12f deg\n", lat_clamp);
  end

  % 2) False easting/northing
  lon_adj = lon_norm + 180;   % [0, 360)
  lat_adj = lat_clamp + 90;   % [0, 180)

  steps.adjusted.lon_adj = lon_adj;
  steps.adjusted.lat_adj = lat_adj;

  if (verbose)
    fprintf("\nFALSE EASTING / NORTHING\n");
    fprintf("  lon_adj = lon_norm + 180 = %.12f deg (range [0,360))\n", lon_adj);
    fprintf("  lat_adj = lat_clamp + 90 = %.12f deg (range [0,180))\n", lat_adj);
  end

  letters18 = "ABCDEFGHIJKLMNOPQR";         % 18 letters (A..R)
  letters24 = "ABCDEFGHIJKLMNOPQRSTUVWX";   % 24 letters (A..X)

  locator = "";

  % 3) Field (base 18)
  field_lon = floor(lon_adj / 20);   % 0..17
  field_lat = floor(lat_adj / 10);   % 0..17

  steps.field.field_lon = field_lon;
  steps.field.field_lat = field_lat;

  locator(1) = letters18(field_lon + 1);
  locator(2) = letters18(field_lat + 1);

  lon_rem_field = lon_adj - field_lon * 20;
  lat_rem_field = lat_adj - field_lat * 10;

  steps.field.lon_rem_field = lon_rem_field;
  steps.field.lat_rem_field = lat_rem_field;

  if (verbose)
    fprintf("\nFIELD (pair 1, letters A..R)\n");
    fprintf("  field_lon = floor(lon_adj / 20) = %d -> '%c'\n", field_lon, locator(1));
    fprintf("  field_lat = floor(lat_adj / 10) = %d -> '%c'\n", field_lat, locator(2));
    fprintf("  lon remainder inside field: %.12f deg\n", lon_rem_field);
    fprintf("  lat remainder inside field: %.12f deg\n", lat_rem_field);
  end

  if (precision == 2), return; end

  % 4) Square (base 10 digits)
  square_lon = floor(lon_rem_field / 2);  % 0..9
  square_lat = floor(lat_rem_field / 1);  % 0..9

  steps.square.square_lon = square_lon;
  steps.square.square_lat = square_lat;

  locator(3) = char('0' + square_lon);
  locator(4) = char('0' + square_lat);

  lon_rem_square = lon_rem_field - square_lon * 2;
  lat_rem_square = lat_rem_field - square_lat * 1;

  steps.square.lon_rem_square = lon_rem_square;
  steps.square.lat_rem_square = lat_rem_square;

  if (verbose)
    fprintf("\nSQUARE (pair 2, digits 0..9)\n");
    fprintf("  square_lon = floor(lon_rem_field / 2) = %d -> '%c'\n", square_lon, locator(3));
    fprintf("  square_lat = floor(lat_rem_field / 1) = %d -> '%c'\n", square_lat, locator(4));
    fprintf("  lon remainder inside square: %.12f deg\n", lon_rem_square);
    fprintf("  lat remainder inside square: %.12f deg\n", lat_rem_square);
  end

  if (precision == 4), return; end

  % 5) Subsquare (base 24 letters)
  lon_sub_width  = 2 / 24;
  lat_sub_height = 1 / 24;

  subs_lon = floor(lon_rem_square / lon_sub_width);    % 0..23
  subs_lat = floor(lat_rem_square / lat_sub_height);   % 0..23

  steps.subsquare.subs_lon = subs_lon;
  steps.subsquare.subs_lat = subs_lat;

  locator(5) = letters24(subs_lon + 1);
  locator(6) = letters24(subs_lat + 1);

  lon_rem_sub = lon_rem_square - subs_lon * lon_sub_width;
  lat_rem_sub = lat_rem_square - subs_lat * lat_sub_height;

  steps.subsquare.lon_sub_width_deg  = lon_sub_width;
  steps.subsquare.lat_sub_height_deg = lat_sub_height;
  steps.subsquare.lon_rem_sub = lon_rem_sub;
  steps.subsquare.lat_rem_sub = lat_rem_sub;

  if (verbose)
    fprintf("\nSUBSQUARE (pair 3, letters A..X)\n");
    fprintf("  lon_sub_width  = 2/24  = %.12f deg  (= 5 arc-min)\n", lon_sub_width);
    fprintf("  lat_sub_height = 1/24  = %.12f deg  (= 2.5 arc-min)\n", lat_sub_height);
    fprintf("  subs_lon = floor(lon_rem_square / lon_sub_width) = %d -> '%c'\n", subs_lon, locator(5));
    fprintf("  subs_lat = floor(lat_rem_square / lat_sub_height) = %d -> '%c'\n", subs_lat, locator(6));
    fprintf("  lon remainder inside subsquare: %.12f deg\n", lon_rem_sub);
    fprintf("  lat remainder inside subsquare: %.12f deg\n", lat_rem_sub);
  end

  if (precision == 6), return; end

  % 6) Extended square (base 10 digits)
  lon_ext_width  = lon_sub_width / 10;
  lat_ext_height = lat_sub_height / 10;

  ext_lon = floor(lon_rem_sub / lon_ext_width);      % 0..9
  ext_lat = floor(lat_rem_sub / lat_ext_height);     % 0..9

  steps.extended.ext_lon = ext_lon;
  steps.extended.ext_lat = ext_lat;

  locator(7) = char('0' + ext_lon);
  locator(8) = char('0' + ext_lat);

  lon_rem_ext = lon_rem_sub - ext_lon * lon_ext_width;
  lat_rem_ext = lat_rem_sub - ext_lat * lat_ext_height;

  steps.extended.lon_ext_width_deg  = lon_ext_width;
  steps.extended.lat_ext_height_deg = lat_ext_height;
  steps.extended.lon_rem_ext = lon_rem_ext;
  steps.extended.lat_rem_ext = lat_rem_ext;

  if (verbose)
    fprintf("\nEXTENDED SQUARE (pair 4, digits 0..9)\n");
    fprintf("  lon_ext_width  = (2/24)/10 = %.12f deg (= 30 arc-sec)\n", lon_ext_width);
    fprintf("  lat_ext_height = (1/24)/10 = %.12f deg (= 15 arc-sec)\n", lat_ext_height);
    fprintf("  ext_lon = floor(lon_rem_sub / lon_ext_width) = %d -> '%c'\n", ext_lon, locator(7));
    fprintf("  ext_lat = floor(lat_rem_sub / lat_ext_height) = %d -> '%c'\n", ext_lat, locator(8));
    fprintf("  lon remainder inside extended square: %.12f deg\n", lon_rem_ext);
    fprintf("  lat remainder inside extended square: %.12f deg\n", lat_rem_ext);
  end
end

INPUT
  latitude  = 50.705787140466 deg
  longitude = 4.392066457029 deg

NORMALIZE / CLAMP
  longitude normalized to [-180, 180): 4.392066457030 deg
  latitude clamped to [-90, 90):      50.705787140466 deg

FALSE EASTING / NORTHING
  lon_adj = lon_norm + 180 = 184.392066457030 deg (range [0,360))
  lat_adj = lat_clamp + 90 = 140.705787140466 deg (range [0,180))

FIELD (pair 1, letters A..R)
  field_lon = floor(lon_adj / 20) = 9 -> 'J'
  field_lat = floor(lat_adj / 10) = 14 -> 'O'
  lon remainder inside field: 4.392066457030 deg
  lat remainder inside field: 0.705787140466 deg

SQUARE (pair 2, digits 0..9)
  square_lon = floor(lon_rem_field / 2) = 2 -> '2'
  square_lat = floor(lat_rem_field / 1) = 0 -> '0'
  lon remainder inside square: 0.392066457030 deg
  lat remainder inside square: 0.705787140466 deg

SUBSQUARE (pair 3, letters A..X)
  lon_sub_width  = 2/24  = 0.083333333333 deg  (= 5 arc-min)
  lat_sub_height = 1/24  = 0.041666666667 deg  (= 2.5 arc-min)
  subs_lon = floor(lon_rem_square / lon_sub_width) = 4 -> 'E'
  subs_lat = floor(lat_rem_square / lat_sub_height) = 16 -> 'Q'
  lon remainder inside subsquare: 0.058733123696 deg
  lat remainder inside subsquare: 0.039120473799 deg

RESULT
  Maidenhead locator (6 chars): JO20EQ