/*
 * Berechnung von Uebergangsboegen nach Dr.-Ing. Peter Schuhr
 *
 * Aus:
 *
 * Anwendung der ROMBERG-Integration zur Uebergangsbogenberechnung aus
 * Kruemmungsbildern.
 *
 * Dr.-Ing. Peter Schuhr, Schienen der Welt, 17(1986), Januar-Heft,
 * S. 17-21
 *
 * Originales BASIC-Programm in C neu geschrieben:
 * Joerg Wunsch, 2007-11-19
 */

/*
 * Copyright (c) 2007 Joerg Wunsch
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE DEVELOPERS ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE DEVELOPERS BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include <assert.h>
#include <math.h>
#include <stdio.h>

enum kruemmungslinie
{
  KEINE_AUSWAHL,
  KLOTHOIDE,
  PARABEL,
  KLEIN,
  BLOSS,
  LETZTE_AUSWAHL
};

/*
 * Eingangsparameter
 */
static enum kruemmungslinie Kz;	/* Kennzahl fuer Kruemmungslinie */
static double Ra, Re;		/* Anfangs- und Endradius */
static double Lz;		/* Laenge bis Z */
static double Le;		/* Gesamtlaenge des Uebergangsbogens */

static double Ka, Ke;		/* Kruemmung bei Ra/Re */

static double const eps = 1e-6;
static double const rho = 200 / M_PI;

/*
 * Berechnung des T-Wertes in Abhaengigkeit der gewaehlten
 * Kruemmungslinie fuer den Uebergangsbogen.
 *
 * Anders als im Original von SCHUHR wird die Ableitung des Sinus oder
 * Kosinus vom Aurufer durchgefuehrt, statt sie in einer eher
 * umstaendlich formulierten Fallunterscheidung innerhalb von FNT()
 * selbst auszufuehren.
 */
static double
T(double L1)
{
  switch (Kz)
    {
    case KLOTHOIDE:
      return Ka * L1 +
      (Ke - Ka) / Le * (L1 * L1) / 2.0;

    case PARABEL:
      if (L1 > Le / 2)
	return Le * (Ka + Ke) / 2.0 -
	(Ke * (Le - L1) -
	 2.0 / 3.0 * (Ke - Ka) / Le / Le * pow(Le - L1, 3.0));
      else
	return Ka * L1 +
	2.0 / 3.0 * (Ke - Ka) / Le / Le * pow(L1, 3.0);

    case KLEIN:
      return Ka * L1 +
      (Ke - Ka) * ((L1 * L1) / 2.0 / Le +
		   Le / 4.0 / M_PI / M_PI *
		   (cos(2.0 * M_PI * L1 / Le) - 1.0));

    case BLOSS:
      return Ka * L1 +
      (Ke - Ka) / (Le * Le) * pow(L1, 3.0) -
      (Ke - Ka) / 2.0 / pow(Le, 3.0) * pow(L1, 4.0);

    default:
      return NAN;
    }
}

/*
 * Berechnung der Kruemmung im Punkt, Fallunterscheidung wie bei T().
 *
 * Diese Berechnung ist im originalen BASIC-Programm nicht vorhanden,
 * die entsprechenden Formeln sind aber im Text des Artikels erwaehnt.
 */
static double
K(double L1)
{
  switch (Kz)
    {
    case KLOTHOIDE:
      return Ka +
      (Ke - Ka) * L1 / Le;

    case PARABEL:
      if (L1 > Le / 2)
	return Ke -
	2.0 * (Ke - Ka) * (1.0 - L1 / Le) * (1.0 - L1 / Le);
      else
	return Ke +
	2.0 * (Ke - Ka) * (L1 / Le) * (L1 / Le);

    case KLEIN:
      return Ka +
      (Ke - Ka) * (L1 / Le - sin(2.0 * M_PI * L1 / Le) / (2 * M_PI));

    case BLOSS:
      return Ka +
      (Ke - Ka) * (3.0 * L1 * L1 / (Le * Le) -
		   2.0 * pow(L1, 3.0) / pow(Le, 3.0));

    default:
      return NAN;
    }
}


/*
 * Integration nach ROMBERG fuer die Berechnung der kartesischen
 * Koordinaten.
 *
 * Im Gegensatz zum Original "Romb" als Funktion ausgefuehrt, die ihr
 * Ergebnis als return-Wert abliefert.
 */
static double
romberg(double a, double b, double (*transfer)(double))
{
  double t[31];
  double r, v;
  double h, h1, h2;
  unsigned int i, k, m, n;

  h = b - a;
  r = v = 0.0;

  if (h == 0.0)
    return 0.0;

  t[0] = h * (transfer(T(a)) + transfer(T(b))) / 2;

  for (n = 1, k = 1;;)
    {
      for (;; k++)
	{
	  t[k] = 0.0;
	  h /= 2;
	  h1 = t[0];
	  t[0] = 0.0;

	  for (i = 1; i <= n; i++)
	    t[0] += transfer(T(a + (2 * i - 1) * h));

	  t[0] *= h;
	  t[0] += h1 / 2;

	  n *= 2;

	  for (i = 1, m = 1; i <= k; i++)
	    {
	      m *= 4;
	      assert(i < sizeof(t) / sizeof(t[0]));
	      h2 = t[i];
	      t[i] = t[i - 1] + (t[i - 1] - h1) / (m - 1);
	      h1 = h2;
	    }

	  assert(k < sizeof(t) / sizeof(t[0]));
	  if (fabs(t[k] - v) < eps)
	    return t[k];

	  v = t[k];
	}
    }
}

int
main(void)
{
  double k, t, x, y;
  int i;
  char line[100];

  for (;;)
    {
      printf("Kz    Ra     Re     Lz     Le: ");
      fflush(stdout);
      if (fgets(line, sizeof line, stdin) == NULL)
	{
	  fprintf(stderr, "Read error\n");
	  return 1;
	}
      if (sscanf(line, "%i %lf %lf %lf %lf", &i, &Ra, &Re, &Lz, &Le) != 5)
	{
	  fprintf(stderr, "Parse error\n");
	  return 1;
	}
      if (i > KEINE_AUSWAHL && i < LETZTE_AUSWAHL)
	{
	  Kz = (enum kruemmungslinie)i;
	}
      else
	{
	  fprintf(stderr, "Ungueltige Auswahl fuer Kz\n");
	  return 1;
	}
      Ka = (Ra == 0.0)? 0.0: (1 / Ra);
      Ke = (Re == 0.0)? 0.0: (1 / Re);
      y = romberg(0.0, Lz, sin);
      x = romberg(0.0, Lz, cos);
      t = T(Lz) * rho;		/* Winkelausgabe in gon */
      k = K(Lz) * 1000;		/* Kruemmungsausgabe in 1000/m */
      printf("t = %12.6f, y = %12.6f, x = %12.6f, k = %12.6f\n",
	     t, y, x, k);
    }
}