unit uHCGCI;

interface

uses
  Classes, SysUtils, Math, {MathLib, Math387,} uHCTypeDefs, uHCDefines;

//var
function atan2(y: extended; x: extended): extended;
procedure HCNewPos(var Lon, Lat: integer; const heading, distance: integer;
  const err_no: short);
procedure HCSetGroupSlivers(pGroup: PHCGroup; pHeading: long);
function HCSliverToDeg(pSlivers: long): short;
function HCSliverToDegF(pSlivers: long): double; experimental;
function HCDegSliverLonF(const pLon: double): long;
function HCDegSliverLatF(const pLat: double): long;
function SliverToDegLatF(const pSlivers: long): double;
function SliverToDegLonF(const pSlivers: long): double;
function HCDegToSliver(pDegrees: short): long;
procedure HCGetRangeBearing(const pLat1: double; const pLon1: double;
  const pLat2: double; const pLon2: double; out pRange: double; out pBearing: double);
procedure HCFindRangeBearing(const pSLat: long; const pSLon: long;
  const pTLat: long; const pTLon: long; out pRange: long; out pBearing: long);
procedure Polar2DestLatLon(const pStartLatDeg: double; const pStartLonDeg: double;
  const pHdgDeg: double; const pDistNM: double; out pDestLatDeg: double;
  out pDestLonDeg: double); experimental;
procedure VincentyDirect(const pStartLatDeg: double; const pStartLonDeg: double;
  const pDistNM: double; const pStartHdgDeg: double; out pDestLatDeg: double;
  out pDestLonDeg: double; out pDestHdgDeg: double; out pIterations: long); experimental;
function HCInterceptAngleDeg(const pSLat: double; const pSLon: double;
  const pTLat: double; const pTLon: double; const pSSpeed: double;
  const pTSpeed: double; const pTHeading: double; out pITime: long;
  out pIDistance: double; out pIAngle: double): short; experimental;

implementation

(* see FPATAN2.TM (here pruned) :
  Subject: Re: How to duplicate C function ATAN2 in Delphi?
  From: Terje Mathisen <Terje.Mathisen@hda.hydro.com>
  Date: 18 May 1995 19:25:10 GMT

  I think you should seriously consider using the FPATAN
  instruction for this!

  This x87 opcode implements an IEEE-compliant ATAN2() function,
  with full extended precision, and the hardware will handle all
  the special cases for you.

  If you have numeric exceptions enabled, and input bogus values,
  the x87 chip will raise the appropriate signal, without the need
  for upfront testing of parameters.

  A BP/TP/Delphi-compatible version would look like this: *)
function atan2(y: extended; x: extended): extended;
  assembler;
{$IFDEF fpc}
{$ASMMODE intel}
{$ENDIF}
asm
         FLD     [y] ;
         FLD     [x] ;
         FPATAN
end;

{*** Calculate new position after moving distance in the direction heading ***}
{***  place result back into Lon and Lat.  Uses Great Circle routes        ***}
procedure HCNewPos(var Lon, Lat: integer; const heading, distance: integer;
  const err_no: short);
var
  dLon, dLat, dheading, ddistance: double;
  dFinalLon, dFinalLat: double;
  EPS, deltalon: double;
begin
  dLon := Lon / SLIVER1 * DR;
  dLat := (90 - (Lat / SLIVER1)) * DR;
  dheading := heading / SLIVER1 * DR; // ALE 20090914 this looks wrong
  ddistance := distance / NM1 * DR / 60;
  EPS := 0.00000000005;

  dFinalLat := ArcSin(Sin(dLat) * Cos(ddistance) + Cos(dLat) *
    Sin(ddistance) * Cos(dheading));

  if (Abs(cos(dFinalLat)) < EPS) then
  begin
    dFinalLon := 0;
  end
  else
  begin
    deltalon := atan2(Sin(dheading) * Sin(ddistance) * Cos(dLat),
      Cos(ddistance) - Sin(dLat) * Sin(dFinalLat));
    dFinalLon := fmod(dLon + deltalon, 2 * PI);
  end;

  dFinalLat := dFinalLat / DR;
  dFinalLon := dFinalLon / DR;

  Lon := Trunc(dFinalLon * SLIVER1);
  Lat := Trunc((90 - dFinalLat) * SLIVER1);

end;

{*** end HCNewPos ***}

{*** Set heading of a group and its slivers ***}
procedure HCSetGroupSlivers(pGroup: PHCGroup; pHeading: long);
var
  lUnit: PHCUnit;
begin
  pGroup^.slivers := pHeading;
  pGroup^.course := HCSliverToDeg(pHeading);
  assert((pGroup^.course >= 0) and (pGroup^.course < 360));
  lUnit := pGroup^.Units;
  while (lUnit <> nil) do
  begin
    lUnit^.course := pGroup^.course;
    lUnit := lUnit^.Next;
  end;
  //UpdateActive_mini_report(G, (UnitPtr)-1L);
end;

{*** end HCSetGroupSlivers ***}

{*** Convert from slivers to degrees ***}
function HCSliverToDeg(pSlivers: long): short;
var
  lDeg: long;
begin
  lDeg := (pSlivers + SLIVER05) div SLIVER1;

  if (lDeg = 360) then
    lDeg := 0;

  assert((lDeg >= 0) and (lDeg < 360));

  Result := short(lDeg);
end;

{*** end HCSliverToDeg ***}

{*** Convert from slivers to degrees ***}
function HCSliverToDegF(pSlivers: long): double; experimental;
begin
  result := pSlivers / SLIVER1;

  if (result >= 360) then
    result := result - 360
  else if (result < 0) then
    result := result + 360;

  assert((result >= 0) and (result < 360));
end;

{*** end HCSliverToDegF ***}

{*** Convert from slivers to degrees for latitude ***}
function SliverToDegLatF(const pSlivers: long): double;
var
  lDeg: double;
begin
  lDeg := 90 - pSlivers / SLIVER1;

  assert((lDeg >= -90) and (lDeg <= 90));

  Result := lDeg;
end;
{*** end SliverToDegLatF ***}

{*** Convert from slivers to degrees for longitude ***}
function SliverToDegLonF(const pSlivers: long): double;
var
  lDeg: double;
begin
  lDeg := pSlivers / SLIVER1;

  //assert((lDeg >= -90) and (lDeg <= 90));

  Result := lDeg;
end;
{*** end SliverToDegLonF ***}

{*** Convert from degrees to slivers for longitude ***}
function HCDegSliverLonF(const pLon: double): long;
begin
  result := Round(((pLon)*SLIVER15 + ((pLon)/3)) / 15);
end;
{*** end HCDegSliverLonF ***}

{*** Convert from degrees to slivers for latitude ***}
function HCDegSliverLatF(const pLat: double): long;
begin
  result := Round(SLIVER90 - (((pLat)*SLIVER15+((pLat)/3))/15));
end;
{*** end HCDegSliverLatF ***}

{*** Convert from degrees to slivers ***}
function HCDegToSliver(pDegrees: short): long;
var
  lSlivers: long;
begin
  lSlivers := (pDegrees) * SLIVER1;

  Result := lSlivers;
end;

{*** end HCDegToSliver ***}

{*** Find range and bearing between points 1 and 2 ***}
procedure HCGetRangeBearing(const pLat1: double; const pLon1: double;
  const pLat2: double; const pLon2: double; out pRange: double; out pBearing: double);
const
  PI = 3.1415926535897932384626433832795;
  DR = 3.1415926535897932384626433832795 / 180; // degrees to radians
  RD = 180 / 3.1415926535897932384626433832795; // radians to degrees
var
  lEarthRadius: double;
  lLat1r: double;
  lLon1r: double;
  lLat2r: double;
  lLon2r: double; // radians
  lA: double;
  lB: double;
  lC: double;
begin
  // ALE 20140509 lat and lon values are in degrees
  //  bearing is in degrees, range is in nm
  lEarthRadius := 3443; // in nautical miles

  lLat1r := pLat1 * DR;
  lLon1r := pLon1 * DR;
  lLat2r := pLat2 * DR;
  lLon2r := pLon2 * DR;

  // ALE 20140509 figure the range
  //if (pRange > 0) then
  begin
    lA := power(sin((lLat2r - lLat1r) / 2.0), 2);
    lB := power(sin((lLon2r - lLon1r) / 2.0), 2);
    lC := sqrt(lA + cos(lLat2r) * cos(lLat1r) * lB);

    pRange := 2 * arcsin(lC) * lEarthRadius;
  end;

  //if (pBearing > 0) then
  begin
    // ALE 20140509 figure the bearing
    lA := cos(lLat2r) * sin(lLon2r - lLon1r);
    lB := cos(lLat1r) * sin(lLat2r) - sin(lLat1r) * cos(lLat2r) * cos(lLon2r - lLon1r);

    pBearing := atan2(lA, lB) * RD;
    pBearing := fmod((pBearing + 360), 360);
  end;
end;

{*** end HCGetRangeBearing ***}

{*** Find range and bearing between source and destination ***}
procedure HCFindRangeBearing(const pSLat: long; const pSLon: long;
  const pTLat: long; const pTLon: long; out pRange: long; out pBearing: long);
var
  lDSLat: double;
  lDSLon: double;
  lDTLat: double;
  lDTLon: double;
  lDrange: double;
  lDbearing: double;
begin

  // ALE 20080419 skip the whole ball of wax when S and T positions are the same
  if ((pSLat = pTLat) and (pSLon = pTLon)) then
  begin
    pRange := 0;
    pBearing := 0;
  end
  else
  begin
    lDSLat := 90 - pSLat / SLIVER1;
    lDSLon := pSLon / SLIVER1;
    lDTLat := 90 - pTLat / SLIVER1;
    lDTLon := pTLon / SLIVER1;

    HCGetRangeBearing(lDSLat, lDSLon, lDTLat, lDTLon, lDrange, lDbearing);

    pRange := Round(NM1 * lDrange);
    lDbearing := SLIVER1 * lDbearing;
    pBearing := Round(lDbearing);
  end;
end;

{*** end HCFindRangeBearing ***}

{*** Convert from altitude in meters to HC altitude band ***}
function AltitudeToAltBand(const pAltitude: HCALTITUDE_M): HCALTITUDE;
begin
  if (pAltitude <= ALT_TOODEEP) then
    Result := TooDeep
  else if (pAltitude <= ALT_VDEEP) then
    Result := Vdeep
  else if (pAltitude <= ALT_DEEP) then
    Result := Deep
  else if (pAltitude <= ALT_INTERMED) then
    Result := Intermed
  else if (pAltitude <= ALT_SHALLOW) then
    Result := Shallow
  else if (pAltitude <= ALT_SNORKEL) then
    Result := Snorkel
  else if (pAltitude <= ALT_SEA_LEVEL) then
    Result := SeaLevel
  else if (pAltitude <= ALT_SCOPE_ALT) then
    Result := ScopeAlt
  else if (pAltitude <= ALT_VLOW) then
    Result := Vlow
  else if (pAltitude <= ALT_LOW) then
    Result := Low
  else if (pAltitude <= ALT_MEDIUM) then
    Result := Medium
  else if (pAltitude <= ALT_HIGH) then
    Result := High
  else if (pAltitude <= ALT_VHIGH) then
    Result := Vhigh
  else
    Result := SubOrbit;
end;

{*** end HCAltitudeToAltBand ***}

{*** ALE 20180222 given a starting lat/lon, a heading and distance ***}
{***  calculate the destination lat/lon                            ***}
{***  lat/lon are in degrees                                       ***}
procedure Polar2DestLatLon(const pStartLatDeg: double; const pStartLonDeg: double;
  const pHdgDeg: double; const pDistNM: double; out pDestLatDeg: double;
  out pDestLonDeg: double); experimental;
var
  lDestHdgDeg: double;
  lIter: long;
begin
  VincentyDirect(pStartLatDeg, pStartLonDeg, pDistNM, pHdgDeg,
    pDestLatDeg, pDestLonDeg, lDestHdgDeg, lIter);
end;

{**
 * Vincenty direct calculation.
 *
 * @private
 * @param   [number] distance - Distance along bearing in metres.
 * @param   [number] initialBearing - Initial bearing in degrees from north.
 * @returns [Object] Object including point (destination point), finalBearing.
 * @throws  [Error]  If formula failed to converge.
 * See: https://github.com/chrisveness/geodesy
 *}
procedure VincentyDirect(const pStartLatDeg: double; const pStartLonDeg: double;
  const pDistNM: double; const pStartHdgDeg: double; out pDestLatDeg: double;
  out pDestLonDeg: double; out pDestHdgDeg: double; out pIterations: long); experimental;
const
  PI = 3.1415926535897932384626433832795;
  DR = 3.1415926535897932384626433832795 / 180; // degrees to radians
  RD = 180 / 3.1415926535897932384626433832795; // radians to degrees
  WGS84a = 6378137;
  WGS84b = 6356752.314245;
  WGS84f = 1 / 298.257223563;
var
  lIter: integer;
  lDistM: double;
  lStartLatRad, lStartLonRad: double; // ALE 20180303 radians
  lHdgRad: double; // ALE 20180302 radians
  lElA, lElB, lElF: double; // ALE 20180306 Ellipsoid characteristic
  lSinHdg, lCosHdg: double;
  lTanU1, lCosU1, lSinU1: double;
  lSigma1: double;
  lSinAlpha: double;
  lCosSqAlpha: double;
  lUSq: double;
  lA, lB: double;
  lSigma: double;
  lCos2SigmaM: double;
  lSinSigma, lCosSigma: double;
  lDeltaSigma: double;
  lSigmaPrime: double;

  lX: double;
  lLamda: double;
  lC: double;
  lL: double;
  lDestLatRad, lDestLonRad: double; // ALE 20180303 radians
  lDestHdgRad: double;
begin
  lDistM := pDistNM * 1852;
  //var φ1 = this.lat.toRadians(), λ1 = this.lon.toRadians();
  lStartLatRad := DR * pStartLatDeg;
  lStartLonRad := DR * pStartLonDeg;

  //var α1 = initialBearing.toRadians();
  lHdgRad := DR * pStartHdgDeg;
  //var s = distance;

  //var a = this.datum.ellipsoid.a, b = this.datum.ellipsoid.b, f = this.datum.ellipsoid.f;
  lElA := WGS84a;
  lElB := WGS84b;
  lElF := WGS84f;

  //var sinα1 = Math.sin(α1);
  //var cosα1 = Math.cos(α1);
  lSinHdg := sin(lHdgRad);
  lCosHdg := cos(lHdgRad);

  //var tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
  lTanU1 := (1 - lElF) * tan(lStartLatRad);
  lCosU1 := 1 / sqrt(1 + lTanU1 * lTanU1);
  lSinU1 := lTanU1 * lCosU1;

  //var σ1 = Math.atan2(tanU1, cosα1);
  lSigma1 := atan2(lTanU1, lCosHdg);
  //var sinα = cosU1 * sinα1;
  lSinAlpha := lCosU1 * lSinHdg;
  //var cosSqα = 1 - sinα*sinα;
  lCosSqAlpha := 1 - lSinAlpha * lSinAlpha;
  //var uSq = cosSqα * (a*a - b*b) / (b*b);
  lUSq := lCosSqAlpha * (lElA * lElA - lElB * lElB) / (lElB * lElB);
  //var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  lA := 1 + lUSq / 16384 * (4096 + lUSq * (-768 + lUSq * (320 - 175 * lUSq)));
  //var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
  lB := lUSq / 1024 * (256 + lUSq * (-128 + lUSq * (74 - 47 * lUSq)));

  //var cos2σM, sinσ, cosσ, Δσ;

  //var σ = s / (b*A), σʹ, iterations = 0;
  lSigma := lDistM / (lElB * lA);
  lIter := 0;
  repeat
    //cos2σM = Math.cos(2*σ1 + σ);
    lCos2SigmaM := cos(2 * lSigma1 + lSigma);
    //sinσ = Math.sin(σ);
    lSinSigma := sin(lSigma);
    //cosσ = Math.cos(σ);
    lCosSigma := cos(lSigma);
    //Δσ = B*sinσ*(cos2σM+B/4*(cosσ*(-1+2*cos2σM*cos2σM)-
    //    B/6*cos2σM*(-3+4*sinσ*sinσ)*(-3+4*cos2σM*cos2σM)));
    lDeltaSigma := lB * lSinSigma * (lCos2SigmaM + lB / 4 *
      (lCosSigma * (-1 + 2 * lCos2SigmaM * lCos2SigmaM) - lB /
      6 * lCos2SigmaM * (-3 + 4 * lSinSigma * lSinSigma) *
      (-3 + 4 * lCos2SigmaM * lCos2SigmaM)));
    //σʹ = σ;
    lSigmaPrime := lSigma;
    //σ = s / (b*A) + Δσ;
    lSigma := lDistM / (lElB * lA) + lDeltaSigma;

    Inc(lIter);
  until ((abs(lSigma - lSigmaPrime) <= 1e-12) or (lIter >= 100));

  if (lIter >= 100) then
    raise Exception.Create('Vincenty didn''t converge'); // not possible!

  //var x = sinU1*sinσ - cosU1*cosσ*cosα1;
  lX := lSinU1 * lSinSigma - lCosU1 * lCosSigma * lCosHdg;
  //var φ2 = Math.atan2(sinU1*cosσ + cosU1*sinσ*cosα1, (1-f)*Math.sqrt(sinα*sinα + x*x));
  lDestLatRad := atan2(lSinU1 * lCosSigma + lCosU1 * lSinSigma *
    lCosHdg, (1 - lElF) * sqrt(lSinAlpha * lSinAlpha + lX * lX));
  //var λ = Math.atan2(sinσ*sinα1, cosU1*cosσ - sinU1*sinσ*cosα1);
  lLamda := atan2(lSinSigma * lSinHdg, lCosU1 * lCosSigma - lSinU1 *
    lSinSigma * lCosHdg);
  //var C = f/16*cosSqα*(4+f*(4-3*cosSqα));
  lC := lElF / 16 * lCosSqAlpha * (4 + lElF * (4 - 3 * lCosSqAlpha));
  //var L = λ - (1-C) * f * sinα *
  //    (σ + C*sinσ*(cos2σM+C*cosσ*(-1+2*cos2σM*cos2σM)));
  lL := lLamda - (1 - lC) * lElF * lSinAlpha * (lSigma + lC *
    lSinSigma * (lCos2SigmaM + lC * lCosSigma * (-1 + 2 * lCos2SigmaM * lCos2SigmaM)));
  //var λ2 = (λ1+L+3*Math.PI)%(2*Math.PI) - Math.PI;  // normalise to -180..+180
  lDestLonRad := fmod((lStartLonRad + lL + 3 * PI), 2 * PI) - PI;
  // normalise to -180..+180

  //var α2 = Math.atan2(sinα, -x);
  lDestHdgRad := atan2(lSinAlpha, -lX);
  //α2 = (α2 + 2*Math.PI) % (2*Math.PI); // normalise to 0..360
  lDestHdgRad := fmod((lDestHdgRad + 2 * PI), (2 * PI)); // normalise to 0..360

  pDestLatDeg := lDestLatRad * RD;
  pDestLonDeg := lDestLonRad * RD;
  pDestHdgDeg := lDestHdgRad * RD;
  pIterations := lIter;
end;

{*============================== InterceptAngle =============================*
 *--------------------------------------------------------------------------
 * Purpose:           given a source object, this routine calculates a bearing
                      to intercept a target.
 * Called By:
 * Calls:             LongXLongDivLong,FindRangeBearing,isin,iasinl
 * Parameters:        long sLat,sLon - source position in slivers
                      long tLat,tLon - target position in slivers
                      long sSpeed    - source speed in kts (formerly du's/hr)
                      long tSpeed    - target speed in kts (formerly du's/hr)
                      long tHeading  - target heading in slivers
                      long *time     - returns the time to intecept
                      long *distance - returns nm to intercept
 * Global Inputs:
 * Global Outputs:
 * Returns:           the intercept bearing
 * Logic Description: Uses the law of sines to determine the intercept angle.
     Does this by setting up a triangle between the source, target, and the
     imaginary intercept point.  We call FindRangeBearing to determine one
     side of the triangle and the bearing couple with the targets heading
     gives us an angle.  Since the other two sides of the triangle are related
     based upon the ratio of the sSpeed and tSpeed we can then apply the law
     of sines to determine another angle.  This is then coverted to a bearing
     for the source object.  To determine the time to intecept we then apply
     the law of sines to get the side of the triangle that applies to the
     source leg and this gives us the distance the source must travel then
     we divide by sSpeed to give time.
 * Revision History:  4-15-89, MWJ, created.
 *--------------------------------------------------------------------------*
 *============================ InterceptAngle ==============================*

 * Bret's quick Trig Refresher:
 *
 *                  B
 *                 /\
 *                /  \
 *               /    \
 *              /      \
 *             /        \
 *            /          \
 *         c /            \ a
 *          /              \
 *         /                \
 *        /                  \
 *       /                    \
 *      /                      \
 *     /                        \
 *    /                          \
 *  A ---------------------------- C
 *                  b

 * The Law of Sines states that:
 *
 *  sin(A)    sin(B)   sin(C)
 *  ------  = ------ = ------
 *     a         b        c
 *
 * To solve the intercept, we solve this triangle in two parts. First we
 * have to determine the angles.
 * a = sSpeed
 * b = tSpeed
 * c = range between S and T
 *
 * A = difference between T's heading and the bearing between S and T
 * B = asin(b/a * sin(a))
 * C = 180 - A - B
 *
 * Note: for a and b we don't know the actual values.  We know that when we
 * have intercepted, the target as moved from A a distance of
 * targetspeed*intercept_time and the source has moved from B a distance
 * of sspeed*intercept_time.  We don't know intercept time, but we don't
 * actually need it, because the equation for B  contains the term b/a,
 * and so the intercept_time cancels out and we are left with sspeed/tspeed,
 * and we do know these.
 *
 * now that we know angle C, we can compute the actual distance that
 * the source travels (a) using
 * a = c*sin(A)/sin(C).  Since c is a real distance, a is as well.
 * The bearing to intercept is:
 *
 *
 *}
function HCInterceptAngleDeg(const pSLat: double; const pSLon: double;
  const pTLat: double; const pTLon: double; const pSSpeed: double;
  const pTSpeed: double; const pTHeading: double; out pITime: long;
  out pIDistance: double; out pIAngle: double): short; experimental;
var
	lAngleA, a, lAngleB, b, lAngleC, c: double;
	Range, Bearing: double;
	dBearing, daBearing, dtHeading: double;
	dist, angle: double;
        dtime: double;
begin
  if (pSSpeed <= 0) then
  begin
    Result := CFALSE;
    exit;
  end;

  HCGetRangeBearing(pTLat, pTLon, pSLat, pSLon, Range, Bearing);
  dBearing := Bearing;
  daBearing := dBearing;

  dtHeading := pTHeading;

  {*
   * note that we don't actually know a and b here.
   * They are actually a=sSpeed*time and b=tSpeed*time, where time is not known.
   * Because we end up with a/b [which is (sSpeed*time)/(tSpeed*time)] in the equation
   * to compute B below, the time terms cancel and we end up with sSpeed/tSpeed anyhow.  So
   * we just set a=sSpeed and b=tSpeed, becuase it enables us to find B.
   *
   * We will really compute a later, when when know all the angles.
   *
   *}

  a := pSSpeed;
  b := pTSpeed;
  c := Range; // in NM

  // find the angles
  lAngleA := dBearing - dtHeading;

  assert(lAngleA < 360);
  if (lAngleA < 0) then
    lAngleA := lAngleA + 360;
  if (lAngleA > 180) then
    lAngleA := 360 - lAngleA;
  if (lAngleA > 180) then
    lAngleA := 360 - lAngleA;

  assert((lAngleA >=0) and (lAngleA <=180));
  assert(a <> 0);

  // for great circle routes, the opposite path is not always
  // just path + 180
  HCGetRangeBearing(pSLat, pSLon, pTLat, pTLon, Range, Bearing);
  dBearing := Bearing;

  lAngleB := arcsin(b/a * sin(DR*lAngleA))/DR;

  if (IsInfinite(lAngleB) or IsNan(lAngleB)) then
  begin
    Result := CFALSE; // Impossible triangle - no intercept possible.
    exit;
  end;

  lAngleC := 180 - lAngleA - lAngleB;

  if ((lAngleC < 0) or (lAngleC >= 180)) then
  begin
    Result := CFALSE;
    exit;
  end;

  // find the distance, which is a
  dist := c*sin(lAngleA * DR)/sin(lAngleC * DR);

  if ((dist > MAX_DIST_NM) or IsInfinite(dist) or IsNan(dist)) then
  begin
    Result := CFALSE; // Impossible triangle - no intercept possible.
    exit;
  end;

  pIDistance := dist;


  dtime := 3600*dist/pSSpeed;

  assert(dtime >= 0);

  pITime := round(dtime + 0.5); // the 3600 is to get it into seconds - all speeds are in nm/hour

  // now compute the bearing.
  if (pTSpeed = 0) then// ALE 20140504 we were having strange effects when target wasn't moving, this tries to help
  begin
    angle := dBearing;
  end
  else
  begin
    dtHeading := dtHeading - daBearing;
    if (abs(dtHeading) > 180) then
      dtHeading := -dtHeading;

    if (dtHeading < 0) then
      angle := dBearing + lAngleB
    else
      angle := dBearing - lAngleB;
  end;
  if (angle >= 360) then
    angle := angle - 360;
  if (angle < 0) then
    angle := angle + 360;

  assert((angle >=0) and (angle <= 360));

  pIAngle := angle;

  Result := CTRUE;
end;

end.