#! /client/bin/perl

  $B = 48;    # 48 Grad nördlicher Breite
  $M = 6370;  # Maßstab in mm

$S = &sind($B);
$SS = $S**2;
$C = &cosd($B);
$CS = $C * $S;

  $l = -2.89;  # 2.89 Grad westlich der Kartenmitte
  $b = 2.12;   # 2.12 Grad nördlich der Kartenmitte

$t = &cosd($B+$b);
if (abs($l) < 2.0) {
  $u = &cosrd($l);
} else {
  $u = 1 - &cosd($l);
};
$v = &sind($b);
$w = &cosd($b);

$x = ($SS*$u - 1) * $v - $CS * $w * $u;
$y = $t * &sind($l);
$z = $C * $t * &cosd($l) + $S * &sind($B+$b);

$xy = $x**2 + $y**2;
if ($xy < 1e-7) {
  $xyz = $xy;
  $xyzc = 0;
} elsif ($xy < 1) {
  if ($xy < 0.01) {
    $xyz = $xy + $xy**2/4 + $xy**3/8;
  } else {
    $xyz = (2 - 2 * sqrt(1-$xy));
  };
  $xyzc = sqrt($xyz/$xy) - 1;
} else {
  $xyz = 2;
  $xyzc = 0.4142135624;
};

$d1 = $M * sqrt($xyz);
$xmap = -$M * ($x + $xyzc * $x);
$ymap =  $M * ($y + $xyzc * $y);

print "l=$l  b=$b  d1=$d1  x=$x  xmap=$xmap  y=$y  ymap=$ymap\n\n";

exit;

sub sind {
  # sin x;  x in degrees

  return (sin(0.01745329252 * $_[0]));
};


sub cosd {
  # cos x;  x in degrees

  return (cos(0.01745329252 * $_[0]));
};


sub cosrd {
  # 1 - cos x;  x in degrees, for small x

  local ($q) = $_[0]**2;

  return ((-3.866323852e-09 * $q + 0.0001523087099) * $q);
};


sub corrsin {
  # (arcsin sqrt x / sqrt x) - 1;  for small x

  return (((0.04464285714 * $_[0] + 0.075) * $_[0] + 1/6) * $_[0]);
};


sub arcsin {
  # arcsin x;  for x up to 0.7
  local ($x) = $_[0];

  local ($y) = $x + $x * &corrsin($x**2);

  return ($y - (sin($y) - $x) / cos($y));
};


sub arccos {
  # arccos x;  for x up to 0.7
  local ($x) = $_[0];

  local ($y) = 1.570796327 - $x - $x * &corrsin($x**2);

  return ($y + (cos($y) - $x) / sin($y));
};