/* gplan.c
   Routines to chew through tables of perturbations.  */

#include "plantbl.h"
#if __STDC__
extern int epsiln ( double J );
double cos (double);
double sin (double);
double floor (double);
double fabs (double);
void mean_elements (double);
static int sscc (int, double, int);
double g1plan (double, struct plantbl *);
int g2plan (double, struct plantbl *, double *);
int g3plan (double, struct plantbl *, double *, int);
#else
int epsiln();
double cos (), sin (), floor (), fabs ();
static int sscc ();
#endif
#define mods3600(x) ((x) - 1.296e6 * floor ((x)/1.296e6))

static double J2000 = 2451545.0;
#if 0
/* Conversion factors between degrees and radians */
static double DTR = 1.7453292519943295769e-2;
static double RTD = 5.7295779513082320877e1;
static double RTS = 2.0626480624709635516e5;	/* arc seconds per radian */
#endif
static double STR = 4.8481368110953599359e-6;	/* radians per arc second */

/* From Simon et al (1994)  */
static double freqs[] =
{
/* Arc sec per 10000 Julian years.  */
  53810162868.8982,
  21066413643.3548,
  12959774228.3429,
  6890507749.3988,
  1092566037.7991,
  439960985.5372,
  154248119.3933,
  78655032.0744,
  52272245.1795
};

static double phases[] =
{
/* Arc sec.  */
  252.25090552 * 3600.,
  181.97980085 * 3600.,
  100.46645683 * 3600.,
  355.43299958 * 3600.,
  34.35151874 * 3600.,
  50.07744430 * 3600.,
  314.05500511 * 3600.,
  304.34866548 * 3600.,
  860492.1546,
};

static double ss[NARGS][31];
static double cc[NARGS][31];


int 
gplan (J, plan, pobj)
     double J;
     struct plantbl *plan;
     double pobj[];
{
  register double su, cu, sv, cv, T;
  double t, sl, sb, sr;
  int i, j, k, m, n, k1, ip, np, nt;
#ifdef _MSC_VER
  char FAR *p;
  double FAR *pl;
  double FAR *pb;
  double FAR *pr;
#else
  /* On some systems such as Silicon Graphics, char is unsigned
     by default.  */
#ifdef vax
  /* VAX CC rejects "signed char."  */
  char *p;
#else
#ifdef __STDC__
  signed char *p;
#else
  char *p;
#endif
#endif
  double *pl, *pb, *pr;
#endif

  T = (J - J2000) / plan->timescale;
  n = plan->maxargs;
  /* Calculate sin( i*MM ), etc. for needed multiple angles.  */
  for (i = 0; i < n; i++)
    {
      if ((j = plan->max_harmonic[i]) > 0)
	{
	  sr = (mods3600 (freqs[i] * T) + phases[i]) * STR;
	  sscc (i, sr, j);
	}
    }

  /* Point to start of table of arguments. */
  p = plan->arg_tbl;
  /* Point to tabulated cosine and sine amplitudes.  */
#ifdef _MSC_VER
  pl = (double FAR *) plan->lon_tbl;
  pb = (double FAR *) plan->lat_tbl;
  pr = (double FAR *) plan->rad_tbl;
#else
  pl = (double *) plan->lon_tbl;
  pb = (double *) plan->lat_tbl;
  pr = (double *) plan->rad_tbl;
#endif
  sl = 0.0;
  sb = 0.0;
  sr = 0.0;

  for (;;)
    {
      /* argument of sine and cosine */
      /* Number of periodic arguments. */
      np = *p++;
      if (np < 0)
	break;
      if (np == 0)
	{			/* It is a polynomial term.  */
	  nt = *p++;
	  /* Longitude polynomial. */
	  cu = *pl++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pl++;
	    }
	  sl +=  mods3600 (cu);
	  /* Latitude polynomial. */
	  cu = *pb++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pb++;
	    }
	  sb += cu;
	  /* Radius polynomial. */
	  cu = *pr++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pr++;
	    }
	  sr += cu;
	  continue;
	}
      k1 = 0;
      cv = 0.0;
      sv = 0.0;
      for (ip = 0; ip < np; ip++)
	{
	  /* What harmonic.  */
	  j = *p++;
	  /* Which planet.  */
	  m = *p++ - 1;
	  if (j)
	    {
	      k = j;
	      if (j < 0)
		k = -k;
	      k -= 1;
	      su = ss[m][k];	/* sin(k*angle) */
	      if (j < 0)
		su = -su;
	      cu = cc[m][k];
	      if (k1 == 0)
		{		/* set first angle */
		  sv = su;
		  cv = cu;
		  k1 = 1;
		}
	      else
		{		/* combine angles */
		  t = su * cv + cu * sv;
		  cv = cu * cv - su * sv;
		  sv = t;
		}
	    }
	}
      /* Highest power of T.  */
      nt = *p++;
      /* Longitude. */
      cu = *pl++;
      su = *pl++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pl++;
	  su = su * T + *pl++;
	}
      sl += cu * cv + su * sv;
      /* Latitiude. */
      cu = *pb++;
      su = *pb++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pb++;
	  su = su * T + *pb++;
	}
      sb += cu * cv + su * sv;
      /* Radius. */
      cu = *pr++;
      su = *pr++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pr++;
	  su = su * T + *pr++;
	}
      sr += cu * cv + su * sv;
    }
  pobj[0] = STR * sl;
  pobj[1] = STR * sb;
  pobj[2] = STR * plan->distance * sr + plan->distance;
  return (0);
}


/* Prepare lookup table of sin and cos ( i*Lj )
 * for required multiple angles
 */
static int 
sscc (k, arg, n)
     int k;
     double arg;
     int n;
{
  double cu, su, cv, sv, s;
  int i;

  su = sin (arg);
  cu = cos (arg);
  ss[k][0] = su;		/* sin(L) */
  cc[k][0] = cu;		/* cos(L) */
  sv = 2.0 * su * cu;
  cv = cu * cu - su * su;
  ss[k][1] = sv;		/* sin(2L) */
  cc[k][1] = cv;
  for (i = 2; i < n; i++)
    {
      s = su * cv + cu * sv;
      cv = cu * cv - su * sv;
      sv = s;
      ss[k][i] = sv;		/* sin( i+1 L ) */
      cc[k][i] = cv;
    }
  return (0);
}

/* Compute mean elements at Julian date J.  */
static double Args[NARGS];
static double LP_equinox;
static double NF_arcsec;
static double Ea_arcsec;
static double pA_precession;

void
mean_elements (J)
     double J;
{
  double x, T, T2;

  /* Time variables.  T is in Julian centuries.  */
  T = (J - 2451545.0) / 36525.0;
  T2 = T*T;

  /* Mean longitudes of planets (Simon et al, 1994)
     .047" subtracted from constant term for offset to DE403 origin. */

  /* Mercury */
  x = mods3600( 538101628.6889819 * T + 908103.213 );
  x += (6.39e-6 * T
	 - 0.0192789) * T2;
  Args[0] = STR * x;

  /* Venus */
  x = mods3600( 210664136.4335482 * T + 655127.236 );
  x += (-6.27e-6  * T
	 + 0.0059381) * T2;
  Args[1] = STR * x;

  /* Earth  */
  x = mods3600( 129597742.283429 * T + 361679.198 );
  x += (-5.23e-6 * T
	 - 2.04411e-2 ) * T2;
  Ea_arcsec = x;
  Args[2] = STR * x;

  /* Mars */
  x = mods3600(  68905077.493988 * T +  1279558.751 );
  x += (-1.043e-5 * T
	 + 0.0094264) * T2;
  Args[3] = STR * x;

  /* Jupiter */
  x = mods3600( 10925660.377991 * T + 123665.420 );
  x += ((((-3.4e-10 * T
	    + 5.91e-8) * T
	   + 4.667e-6) * T
	  + 5.706e-5) * T
         - 3.060378e-1)*T2;
  Args[4] = STR * x;

   /* Saturn */
  x = mods3600( 4399609.855372 * T + 180278.752 );
  x += (((( 8.3e-10 * T
	  - 1.452e-7) * T
	  - 1.1484e-5) * T
	   - 1.6618e-4) * T
	 + 7.561614E-1)*T2;
  Args[5] = STR * x;

  /* Uranus */
  x = mods3600( 1542481.193933 * T + 1130597.971 )
       + (0.00002156*T - 0.0175083)*T2;
  Args[6] = STR * x;

  /* Neptune */
  x = mods3600( 786550.320744 * T + 1095655.149 )
       + (-0.00000895*T + 0.0021103)*T2;
  Args[7] = STR * x;

  /* Copied from cmoon.c, DE404 version.  */
  /* Mean elongation of moon = D */
  x = mods3600( 1.6029616009939659e+09 * T + 1.0722612202445078e+06 );
  x += (((((-3.207663637426e-013 * T
	       + 2.555243317839e-011) * T
	      + 2.560078201452e-009) * T
	     - 3.702060118571e-005) * T
            + 6.9492746836058421e-03) * T /* D, t^3 */
           - 6.7352202374457519e+00) * T2; /* D, t^2 */
  Args[9] = STR * x;

  /* Mean distance of moon from its ascending node = F */
  x = mods3600( 1.7395272628437717e+09 * T + 3.3577951412884740e+05 );
  x += ((((( 4.474984866301e-013 * T
	       + 4.189032191814e-011) * T
	       - 2.790392351314e-009) * T
	      - 2.165750777942e-006) * T
	      - 7.5311878482337989e-04) * T /* F, t^3 */
	     - 1.3117809789650071e+01) * T2; /* F, t^2 */
  NF_arcsec = x;
  Args[10] = STR * x;

/* Mean anomaly of sun = l' (J. Laskar) */
  x = mods3600(1.2959658102304320e+08 * T + 1.2871027407441526e+06);
  x += ((((((((
	       1.62e-20 * T
	       - 1.0390e-17 ) * T
	      - 3.83508e-15 ) * T
	     + 4.237343e-13 ) * T
	    + 8.8555011e-11 ) * T
	   - 4.77258489e-8 ) * T
	  - 1.1297037031e-5 ) * T
	 + 8.7473717367324703e-05) * T
	- 5.5281306421783094e-01) * T2;
  Args[11] = STR * x;

  /* Mean anomaly of moon = l */
  x = mods3600( 1.7179159228846793e+09 * T + 4.8586817465825332e+05 );
  x += (((((-1.755312760154e-012 * T
		+ 3.452144225877e-011) * T
		- 2.506365935364e-008) * T
	       - 2.536291235258e-004) * T
              + 5.2099641302735818e-02) * T /* l, t^3 */
             + 3.1501359071894147e+01) * T2; /* l, t^2 */
  Args[12] = STR * x;

  /* Mean longitude of moon, re mean ecliptic and equinox of date = L  */
  x = mods3600( 1.7325643720442266e+09 * T + 7.8593980921052420e+05);
  x += ((((( 7.200592540556e-014 * T
	     + 2.235210987108e-010) * T
	    - 1.024222633731e-008) * T
	   - 6.073960534117e-005) * T
	 + 6.9017248528380490e-03) * T /* L, t^3 */
	- 5.6550460027471399e+00) * T2; /* L, t^2 */
  LP_equinox = x;
  Args[13] = STR * x;

  /* Precession of the equinox  */
 x = ((((((((( -8.66e-20*T - 4.759e-17)*T
           + 2.424e-15)*T
           + 1.3095e-12)*T
           + 1.7451e-10)*T
           - 1.8055e-8)*T
           - 0.0000235316)*T
           + 0.000076)*T
           + 1.105414)*T
           + 5028.791959)*T;
  /* Moon's longitude re fixed J2000 equinox.  */
 /*
   Args[13] -= x;
 */
   pA_precession = STR * x;

  /* Free librations.  */
  /* longitudinal libration 2.891725 years */
  x = mods3600( 4.48175409e7 * T + 8.060457e5 );
  Args[14] = STR * x;
  /* libration P, 24.2 years */
  x = mods3600(  5.36486787e6 * T - 391702.8 );
  Args[15] = STR * x;

#if 0
  Args[16] = 0.0;
#endif
  /* libration W, 74.7 years. */
 x = mods3600( 1.73573e6 * T );
 Args[17] = STR * x;
}


/* Generic program to accumulate sum of trigonometric series
   in three variables (e.g., longitude, latitude, radius)
   of the same list of arguments.  */

int
g3plan (J, plan, pobj, objnum)
     double J;
     struct plantbl *plan;
     double pobj[];
     int objnum;
{
  int i, j, k, m, n, k1, ip, np, nt;
#ifdef _MSC_VER
  char FAR *p;
  long FAR *pl;
  long FAR *pb;
  long FAR *pr;
#else
  /* On some systems such as Silicon Graphics, char is unsigned
     by default.  */
#ifdef vax
  /* VAX CC rejects "signed char."  */
  char *p;
#else
#ifdef __STDC__
  signed char *p;
#else
  char *p;
#endif
#endif
  long *pl, *pb, *pr;
#endif
  double su, cu, sv, cv;
  double T, t, sl, sb, sr;

  mean_elements (J);
#if 0
  /* For librations, moon's longitude is sidereal.  */
  if (flag)
    Args[13] -= pA_precession;
#endif
  T = (J - J2000) / plan->timescale;
  n = plan->maxargs;
  /* Calculate sin( i*MM ), etc. for needed multiple angles.  */
  for (i = 0; i < n; i++)
    {
      if ((j = plan->max_harmonic[i]) > 0)
	{
	  sscc (i, Args[i], j);
	}
    }

  /* Point to start of table of arguments. */
  p = plan->arg_tbl;
  /* Point to tabulated cosine and sine amplitudes.  */
#ifdef _MSC_VER
  pl = (long FAR *) plan->lon_tbl;
  pb = (long FAR *) plan->lat_tbl;
  pr = (long FAR *) plan->rad_tbl;
#else
  pl = (long *) plan->lon_tbl;
  pb = (long *) plan->lat_tbl;
  pr = (long *) plan->rad_tbl;
#endif
  sl = 0.0;
  sb = 0.0;
  sr = 0.0;

  for (;;)
    {
      /* argument of sine and cosine */
      /* Number of periodic arguments. */
      np = *p++;
      if (np < 0)
	break;
      if (np == 0)
	{			/* It is a polynomial term.  */
	  nt = *p++;
	  /* "Longitude" polynomial (phi). */
	  cu = *pl++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pl++;
	    }
	  /*	  sl +=  mods3600 (cu); */
	  sl += cu;
	  /* "Latitude" polynomial (theta). */
	  cu = *pb++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pb++;
	    }
	  sb += cu;
	  /* Radius polynomial (psi). */
	  cu = *pr++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pr++;
	    }
	  sr += cu;
	  continue;
	}
      k1 = 0;
      cv = 0.0;
      sv = 0.0;
      for (ip = 0; ip < np; ip++)
	{
	  /* What harmonic.  */
	  j = *p++;
	  /* Which planet.  */
	  m = *p++ - 1;
	  if (j)
	    {
/*	      k = abs (j); */
		if (j < 0)
			k = -j;
		else
			k = j;
	      k -= 1;
	      su = ss[m][k];	/* sin(k*angle) */
	      if (j < 0)
		su = -su;
	      cu = cc[m][k];
	      if (k1 == 0)
		{		/* set first angle */
		  sv = su;
		  cv = cu;
		  k1 = 1;
		}
	      else
		{		/* combine angles */	
		  t = su * cv + cu * sv;
		  cv = cu * cv - su * sv;
		  sv = t;
		}
	    }
	}
      /* Highest power of T.  */
      nt = *p++;
      /* Longitude. */
      cu = *pl++;
      su = *pl++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pl++;
	  su = su * T + *pl++;
	}
      sl += cu * cv + su * sv;
      /* Latitiude. */
      cu = *pb++;
      su = *pb++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pb++;
	  su = su * T + *pb++;
	}
      sb += cu * cv + su * sv;
      /* Radius. */
      cu = *pr++;
      su = *pr++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pr++;
	  su = su * T + *pr++;
	}
      sr += cu * cv + su * sv;
    }
  t = plan->trunclvl;
  pobj[0] = Args[objnum - 1] + STR * t * sl;
  pobj[1] = STR * t * sb;
  pobj[2] = plan->distance * (1.0 + STR * t * sr);
  return (0);
}


/* Generic program to accumulate sum of trigonometric series
   in two variables (e.g., longitude, radius)
   of the same list of arguments.  */
int 
g2plan (J, plan, pobj)
     double J;
     struct plantbl *plan;
     double pobj[];
{
  int i, j, k, m, n, k1, ip, np, nt;
#if _MSC_VER
  char FAR *p;
  long FAR *pl;
  long FAR *pr;
#else
  /* On some systems such as Silicon Graphics, char is unsigned
     by default.  */
#ifdef vax
  /* VAX CC rejects "signed char."  */
  char *p;
#else
#ifdef __STDC__
  signed char *p;
#else
  char *p;
#endif
#endif
  long *pl, *pr;
#endif
  double su, cu, sv, cv;
  double T, t, sl, sr;

  mean_elements (J);
#if 0
  /* For librations, moon's longitude is sidereal.  */
  if (flag)
    Args[13] -= pA_precession;
#endif
  T = (J - J2000) / plan->timescale;
  n = plan->maxargs;
  /* Calculate sin( i*MM ), etc. for needed multiple angles.  */
  for (i = 0; i < n; i++)
    {
      if ((j = plan->max_harmonic[i]) > 0)
	{
	  sscc (i, Args[i], j);
	}
    }

  /* Point to start of table of arguments. */
  p = plan->arg_tbl;
  /* Point to tabulated cosine and sine amplitudes.  */
#ifdef _MSC_VER
  pl = (long FAR *) plan->lon_tbl;
  pr = (long FAR *) plan->rad_tbl;
#else
  pl = (long *) plan->lon_tbl;
  pr = (long *) plan->rad_tbl;
#endif
  sl = 0.0;
  sr = 0.0;

  for (;;)
    {
      /* argument of sine and cosine */
      /* Number of periodic arguments. */
      np = *p++;
      if (np < 0)
	break;
      if (np == 0)
	{			/* It is a polynomial term.  */
	  nt = *p++;
	  /* Longitude polynomial. */
	  cu = *pl++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pl++;
	    }
	  /*	  sl +=  mods3600 (cu); */
	  sl += cu;
	  /* Radius polynomial. */
	  cu = *pr++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pr++;
	    }
	  sr += cu;
	  continue;
	}
      k1 = 0;
      cv = 0.0;
      sv = 0.0;
      for (ip = 0; ip < np; ip++)
	{
	  /* What harmonic.  */
	  j = *p++;
	  /* Which planet.  */
	  m = *p++ - 1;
	  if (j)
	    {
/*	      k = abs (j); */
		if (j < 0)
			k = -j;
		else
			k = j;
	      k -= 1;
	      su = ss[m][k];	/* sin(k*angle) */
	      if (j < 0)
		su = -su;
	      cu = cc[m][k];
	      if (k1 == 0)
		{		/* set first angle */
		  sv = su;
		  cv = cu;
		  k1 = 1;
		}
	      else
		{		/* combine angles */
		  t = su * cv + cu * sv;
		  cv = cu * cv - su * sv;
		  sv = t;
		}
	    }
	}
      /* Highest power of T.  */
      nt = *p++;
      /* Longitude. */
      cu = *pl++;
      su = *pl++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pl++;
	  su = su * T + *pl++;
	}
      sl += cu * cv + su * sv;
      /* Radius. */
      cu = *pr++;
      su = *pr++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pr++;
	  su = su * T + *pr++;
	}
      sr += cu * cv + su * sv;
    }
  t = plan->trunclvl;
  pobj[0] = t * sl;
  pobj[2] = t * sr;
  return (0);
}



/* Generic program to accumulate sum of trigonometric series
   in one variable.  */

double
g1plan (J, plan)
     double J;
     struct plantbl *plan;
{
  int i, j, k, m, k1, ip, np, nt;
#ifdef _MSC_VER
  char FAR *p;
  long FAR *pl;
#else
  /* On some systems such as Silicon Graphics, char is unsigned
     by default.  */
#ifdef vax
  /* VAX CC rejects "signed char."  */
  char *p;
#else
#ifdef __STDC__
  signed char *p;
#else
  char *p;
#endif
#endif
  long *pl;
#endif
  double su, cu, sv, cv;
  double T, t, sl;

  T = (J - J2000) / plan->timescale;
  mean_elements (J);
  /* Calculate sin( i*MM ), etc. for needed multiple angles.  */
  for (i = 0; i < NARGS; i++)
    {
      if ((j = plan->max_harmonic[i]) > 0)
	{
	  sscc (i, Args[i], j);
	}
    }

  /* Point to start of table of arguments. */
  p = plan->arg_tbl;
  /* Point to tabulated cosine and sine amplitudes.  */
#ifdef _MSC_VER
  pl = (long FAR *) plan->lon_tbl;
#else
  pl = (long *) plan->lon_tbl;
#endif
  sl = 0.0;

  for (;;)
    {
      /* argument of sine and cosine */
      /* Number of periodic arguments. */
      np = *p++;
      if (np < 0)
	break;
      if (np == 0)
	{			/* It is a polynomial term.  */
	  nt = *p++;
	  cu = *pl++;
	  for (ip = 0; ip < nt; ip++)
	    {
	      cu = cu * T + *pl++;
	    }
	  /*	  sl +=  mods3600 (cu); */
	  sl += cu;
	  continue;
	}
      k1 = 0;
      cv = 0.0;
      sv = 0.0;
      for (ip = 0; ip < np; ip++)
	{
	  /* What harmonic.  */
	  j = *p++;
	  /* Which planet.  */
	  m = *p++ - 1;
	  if (j)
	    {
/*	      k = abs (j); */
		if (j < 0)
			k = -j;
		else
			k = j;
	      k -= 1;
	      su = ss[m][k];	/* sin(k*angle) */
	      if (j < 0)
		su = -su;
	      cu = cc[m][k];
	      if (k1 == 0)
		{		/* set first angle */
		  sv = su;
		  cv = cu;
		  k1 = 1;
		}
	      else
		{		/* combine angles */
		  t = su * cv + cu * sv;
		  cv = cu * cv - su * sv;
		  sv = t;
		}
	    }
	}
      /* Highest power of T.  */
      nt = *p++;
      /* Cosine and sine coefficients.  */
      cu = *pl++;
      su = *pl++;
      for (ip = 0; ip < nt; ip++)
	{
	  cu = cu * T + *pl++;
	  su = su * T + *pl++;
	}
      sl += cu * cv + su * sv;
    }
  return (plan->trunclvl * sl);
}


/* Compute geocentric moon.  */

extern struct plantbl moonlr, moonlat;
extern double coseps, sineps;

int 
gmoon (J, rect, pol)
     double J;
     double rect[], pol[];
{
  double x, cosB, sinB, cosL, sinL;

  g2plan (J, &moonlr, pol);
  x = pol[0];
  x += LP_equinox;
  if (x < -6.48e5)
    x += 1.296e6;
  if (x > 6.48e5)
    x -= 1.296e6;
  pol[0] = STR * x;
  x = g1plan (J, &moonlat);
  pol[1] = STR * x;
  x = (1.0 + STR * pol[2]) * moonlr.distance;
  pol[2] = x;
  /* Convert ecliptic polar to equatorial rectangular coordinates.  */
  epsiln(J);
  cosB = cos(pol[1]);
  sinB = sin(pol[1]);
  cosL = cos(pol[0]);
  sinL = sin(pol[0]);
  rect[0] = cosB * cosL * x;
  rect[1] = (coseps * cosB * sinL - sineps * sinB) * x;
  rect[2] = (sineps * cosB * sinL + coseps * sinB) * x;
  return 0;
}