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

#define NEQU 2
#define NSTEPS 2000

int error;

typedef long double (*function)(long double, long double *);
function Vector[21];

long double TNTargetF1(long double T, long double *V)
/*--------------------------------------------------*/
/*- This is the first differential equation.       -*/
/*-                                                -*/
/*-  dx[1]                                         -*/
/*-  -----  = TNTargetF1(T, X[1], X[2], ... X[m])  -*/
/*-    dt                                          -*/
/*-                                                -*/
/*- The vector V is defined:                       -*/
/*-      T    = T                                  -*/
/*-      V[1] = X[1]                               -*/
/*-      V[2] = X[2]                               -*/
/*-           .                                    -*/
/*-           .                                    -*/
/*-           .                                    -*/
/*-      V[m] = X[m]                               -*/
/*-                                                -*/
/*- where m is the number of coupled equations.    -*/
/*--------------------------------------------------*/

{
  return V[2];
} /* function TNTargetF1 */

long double TNTargetF2(long double T, long double *V)
/*--------------------------------------------------*/
/*- This is the 2nd differential equation.         -*/
/*-                                                -*/
/*-  dx[2]                                         -*/
/*-  -----  = TNTargetF1(T, X[1], X[2], ... X[m])  -*/
/*-    dt                                          -*/
/*-                                                -*/
/*- The vector V is defined:                       -*/
/*-      V[1] = X[1]                               -*/
/*-      V[2] = X[2]                               -*/
/*-           .                                    -*/
/*-           .                                    -*/
/*-           .                                    -*/
/*-      V[m] = X[m]                               -*/
/*-                                                -*/
/*- where m is the number of coupled equations.    -*/
/*--------------------------------------------------*/
{
  return -sinl(V[1]);
} /* function TNTargetF2 */

long double TNTargetF3(long double T, long double *V) { return 0; }
long double TNTargetF4(long double T, long double *V) { return 0; }
long double TNTargetF5(long double T, long double *V) { return 0; }
long double TNTargetF6(long double T, long double *V) { return 0; }
long double TNTargetF7(long double T, long double *V) { return 0; }
long double TNTargetF8(long double T, long double *V) { return 0; }
long double TNTargetF9(long double T, long double *V) { return 0; }
long double TNTargetF10(long double T, long double *V) { return 0; }
long double TNTargetF11(long double T, long double *V) { return 0; }
long double TNTargetF12(long double T, long double *V) { return 0; }
long double TNTargetF13(long double T, long double *V) { return 0; }
long double TNTargetF14(long double T, long double *V) { return 0; }
long double TNTargetF15(long double T, long double *V) { return 0; }
long double TNTargetF16(long double T, long double *V) { return 0; }
long double TNTargetF17(long double T, long double *V) { return 0; }
long double TNTargetF18(long double T, long double *V) { return 0; }
long double TNTargetF19(long double T, long double *V) { return 0; }
long double TNTargetF20(long double T, long double *V) { return 0; }

void AssignFunc() {
  Vector[0] = NULL;
  Vector[1] = &TNTargetF1;
  Vector[2] = &TNTargetF2;
  Vector[3] = &TNTargetF3;
  Vector[4] = &TNTargetF4;
  Vector[5] = &TNTargetF5;
  Vector[6] = &TNTargetF6;
  Vector[7] = &TNTargetF7;
  Vector[8] = &TNTargetF8;
  Vector[9] = &TNTargetF9;
  Vector[10] = &TNTargetF10;
  Vector[11] = &TNTargetF11;
  Vector[12] = &TNTargetF12;
  Vector[13] = &TNTargetF13;
  Vector[14] = &TNTargetF14;
  Vector[15] = &TNTargetF15;
  Vector[16] = &TNTargetF16;
  Vector[17] = &TNTargetF17;
  Vector[18] = &TNTargetF18;
  Vector[19] = &TNTargetF19;
  Vector[20] = &TNTargetF20;
}

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-    Input:  LowerLimit, UpperLimit, XInitial, NumReturn, NumIntervals     -*/
/*-   Output:  TValues, XValues, Error                                       -*/
/*-                                                                          -*/
/*-  Purpose:  This unit solves an initial value problem (a first order      -*/
/*-            differential equation with initial condition specified).      -*/
/*-            Given a function of the form                                  -*/
/*-                       dx/dt = TNTargetF(T, X)                            -*/
/*-            and an initial condition                                      -*/
/*-                       X[LowerLimit] = const                              -*/
/*-            the Adams-Bashforth/Adams-Moulton four step predictor/        -*/
/*-            corrector method is used to solve for X in  the interval      -*/
/*-            [LowerLimit, UpperLimit].  The interval is divided up into    -*/
/*-            NumIntervals equal sub-intervals.                             -*/
/*-                                                                          -*/
/*-  User-defined Functions: TNTargetF(T, X : real) : real;                  -*/
/*-                                                                          -*/
/*-      User-defined Types: TNvector = array[0..TNArraySize] of real;       -*/
/*-                                                                          -*/
/*-    Global Variables:  LowerLimit   : real;    Lower limit of integration -*/
/*-                       UpperLimit   : real;    Upper limit of integration -*/
/*-                       XInitial     : real;    Value of X at LowerLimit   -*/
/*-                       NumReturn    : integer; Number of values to return -*/
/*-                       NumIntervals : integer; Number of subintervals     -*/
/*-                       TValues      : TNvector; Value of T between limits -*/
/*-                       XValues      : TNvector; Value of X at TValues     -*/
/*-                       Error        : byte;    Flags if something goes    -*/
/*-                                               wrong                      -*/
/*-                                                                          -*/
/*-               Notes:  NumIntervals should be an integer                  -*/
/*-                       multiple of NumReturn if you want                  -*/
/*-                       evenly spaced values returned.                     -*/
/*-                                                                          -*/
/*-              Errors:  0: No errors                                       -*/
/*-                       1: NumReturn <= 0                                  -*/
/*-                       2: NumIntervals < NumReturn                        -*/
/*-                       3: LowerLimit = UpperLimit                         -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void GetFirstFourValues(long double Spacing,
			      long double *T,
			      long double *X,
			      long double **Values,
			      long double (*DerivValues)[NEQU+1][6]);  // Runge-Kutta method;
void GetFirstFourValues38(long double Spacing,
			      long double *T,
			      long double *X,
			      long double **Values,
			      long double (*DerivValues)[NEQU+1][6]);  // Runge-Kutta method;
void GetFirstFourValuesGill(long double Spacing,
			      long double *T,
			      long double *X,
			      long double **Values,
			      long double (*DerivValues)[NEQU+1][6]);  // Runge-Kutta method;


void Adams(long double LowerLimit,
		long double UpperLimit,
		long double *XInitial,
		int NumReturn,
		int NumIntervals,
		unsigned char *Error)
{
  long double Spacing;                 /* Size of each subinterval */
  long double T;                       /* Value of T */
  long double X[NEQU+1], NewX[NEQU+1];     /* Iteration variables */
  int Index,FuncIndex;            /* Counters */
  long double DerivValues[NEQU+1][6];    /* Derivative at last 4 X values */
  long double Values[NEQU+1][6];
  int i_print;                    /* A counter for printing output */

  *Error = 0;
  if ( NumReturn <= 0 )
    *Error = 1;
  if ( NumIntervals < NumReturn )
    *Error = 2;
  if ( LowerLimit == UpperLimit )
    *Error = 3;
  if ( *Error == 0 )
  {
    AssignFunc();         /* Assignment of function pointers to the array */
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
    for (Index = 0; Index <= 5;Index++) {
      DerivValues[FuncIndex][Index] = 0;
    }
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
    {
      DerivValues[FuncIndex][1] = Vector[FuncIndex](LowerLimit, XInitial);
    }
    Spacing = (UpperLimit - LowerLimit) / NumIntervals;
    T = LowerLimit;
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++) {
      X[FuncIndex] = XInitial[FuncIndex];
    }
    settiming();
    printf("%e  ",T);
    for(i_print=1;i_print<=NEQU;i_print++) printf("%e  ",X[i_print]);
    printf("\n");
    GetFirstFourValues(Spacing, &T, X, Values, &DerivValues);
    for (Index = 5; Index <= NumIntervals; Index++)
    {
      {
	for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
	{
	  /* Apply predictor */
	  NewX[FuncIndex] = X[FuncIndex] + Spacing / 720 *
			     (1901 * DerivValues[FuncIndex][5] -
			      2774 * DerivValues[FuncIndex][4] +
			      2616 * DerivValues[FuncIndex][3] -
			      1274 * DerivValues[FuncIndex][2] +
			       251 * DerivValues[FuncIndex][1]);
	}
	T += Spacing;
	for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
	{
	  DerivValues[FuncIndex][1] = DerivValues[FuncIndex][2];
	  DerivValues[FuncIndex][2] = DerivValues[FuncIndex][3];
	  DerivValues[FuncIndex][3] = DerivValues[FuncIndex][4];
	  DerivValues[FuncIndex][4] = DerivValues[FuncIndex][5];
	  DerivValues[FuncIndex][5] = Vector[FuncIndex](T, NewX);
	}
	for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
	{
	  /* Apply corrector */
	  NewX[FuncIndex] = X[FuncIndex] + Spacing / 720 *
			    (251 * DerivValues[FuncIndex][5] +
			     646 * DerivValues[FuncIndex][4] -
			     264 * DerivValues[FuncIndex][3] +
			     106 * DerivValues[FuncIndex][2] -
			      19 * DerivValues[FuncIndex][1]);
	  X[FuncIndex] = NewX[FuncIndex];
	}
	T -= Spacing;
      }
      T += Spacing;
 
//      printf("%e  ",T);
//      for(i_print=1;i_print<=NEQU;i_print++) printf("%e  ",X[i_print]);
//      printf("\n");
    }
  }
  return X[2];
} /* procedure Adams */

void GetFirstFourValues(long double Spacing,
			      long double *T,
			      long double *X,
			      long double **Values,
			      long double (*DerivValues)[NEQU+1][6])

/*---------------------------------------------------------------*/
/*- Input: Spacing, T, X, Values, DerivValues                   -*/
/*- Output: T, X, Values, DerivValues                           -*/
/*-                                                             -*/
/*- This procedure uses the Runge-Kutta one step method to      -*/
/*- solve for the first three (T + i-Spacing, 1<=i<=3) Values.  -*/
/*- The derivative at these points are also calculated          -*/
/*- (DerivValues).                                              -*/
/*---------------------------------------------------------------*/

{
  int Index,FuncIndex;

  long double HalfSpacing;
  long double F1[NEQU+1], F2[NEQU+1], F3[NEQU+1], F4[NEQU+1], Hx[NEQU+1];
  int i_print;

  HalfSpacing = Spacing / 2;

  for (Index = 1; Index < 5; Index++)
  {
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F1[FuncIndex] = Spacing * Vector[FuncIndex](*T, X);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] + F1[FuncIndex]/2;
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F2[FuncIndex] = Spacing * Vector[FuncIndex](*T + HalfSpacing, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] + F2[FuncIndex]/2;
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F3[FuncIndex] = Spacing * Vector[FuncIndex](*T + HalfSpacing, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] + F3[FuncIndex];
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F4[FuncIndex] = Spacing * Vector[FuncIndex](*T + Spacing, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      X[FuncIndex] = X[FuncIndex] +
      (   F1[FuncIndex] + 2 * F2[FuncIndex] +
      2 * F3[FuncIndex] +     F4[FuncIndex]) / 6;
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      (*DerivValues)[FuncIndex][Index + 1] = Vector[FuncIndex](*T, X);
    *T += Spacing;
    timing((int)Spacing);
//    printf("%e  ",*T);
//    for(i_print=1;i_print<=NEQU;i_print++) printf("%e  ",X[i_print]);
//    printf("\n");
  }
} /* procedure GetFirstFourValues */


void GetFirstFourValues38(long double Spacing,
			      long double *T,
			      long double *X,
			      long double **Values,
			      long double (*DerivValues)[NEQU+1][6])

/*---------------------------------------------------------------*/
/*- Input: Spacing, T, X, Values, DerivValues                   -*/
/*- Output: T, X, Values, DerivValues                           -*/
/*-                                                             -*/
/*- This procedure uses the Runge-Kutta one step method to      -*/
/*- solve for the first three (T + i-Spacing, 1<=i<=3) Values.  -*/
/*- The derivative at these points are also calculated          -*/
/*- (DerivValues).                                              -*/
/*---------------------------------------------------------------*/

{
  int Index,FuncIndex;

  long double HalfSpacing;
  long double F1[NEQU+1], F2[NEQU+1], F3[NEQU+1], F4[NEQU+1], Hx[NEQU+1];
  int i_print;

  HalfSpacing = Spacing / 2;

  for (Index = 1; Index < 5; Index++)
  {
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F1[FuncIndex] = Spacing * Vector[FuncIndex](*T, X);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] + F1[FuncIndex]/3;
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F2[FuncIndex] = Spacing * Vector[FuncIndex](*T + Spacing/3, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] - F1[FuncIndex]/3 + F2[FuncIndex];
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F3[FuncIndex] = Spacing * Vector[FuncIndex](*T + 2*Spacing/3, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] + F1[FuncIndex] - F2[FuncIndex] + F3[FuncIndex];
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F4[FuncIndex] = Spacing * Vector[FuncIndex](*T + Spacing, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      X[FuncIndex] = X[FuncIndex] +
      (F1[FuncIndex] + 3 * F2[FuncIndex] + 3 * F3[FuncIndex] + F4[FuncIndex]) / 8;
    *T += Spacing;
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      (*DerivValues)[FuncIndex][Index + 1] = Vector[FuncIndex](*T, X);
//    printf("%e  ",*T);
//    for(i_print=1;i_print<=NEQU;i_print++) printf("%e  ",X[i_print]);
//    printf("\n");
  }
} /* procedure GetFirstFourValues */

void GetFirstFourValuesGill(long double Spacing,
			      long double *T,
			      long double *X,
			      long double **Values,
			      long double (*DerivValues)[NEQU+1][6])

/*---------------------------------------------------------------*/
/*- Input: Spacing, T, X, Values, DerivValues                   -*/
/*- Output: T, X, Values, DerivValues                           -*/
/*-                                                             -*/
/*- This procedure uses the Runge-Kutta one step method to      -*/
/*- solve for the first three (T + i-Spacing, 1<=i<=3) Values.  -*/
/*- The derivative at these points are also calculated          -*/
/*- (DerivValues).                                              -*/
/*---------------------------------------------------------------*/

{
  int Index,FuncIndex;

  long double HalfSpacing;
  long double F1[NEQU+1], F2[NEQU+1], F3[NEQU+1], F4[NEQU+1], Hx[NEQU+1];
  int i_print;

  HalfSpacing = Spacing / 2;

  for (Index = 1; Index < 5; Index++)
  {
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F1[FuncIndex] = Spacing * Vector[FuncIndex](*T, X);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] + F1[FuncIndex]/2;
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F2[FuncIndex] = Spacing * Vector[FuncIndex](*T + HalfSpacing, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] + 1/2*(sqrt(2)-1)*F1[FuncIndex] + (1-1/sqrt(2))*F2[FuncIndex];
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F3[FuncIndex] = Spacing * Vector[FuncIndex](*T + HalfSpacing, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++)
      Hx[FuncIndex] = X[FuncIndex] - 1/sqrt(2)*F2[FuncIndex] + (1+1/sqrt(2))*F3[FuncIndex];
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      F4[FuncIndex] = Spacing * Vector[FuncIndex](*T + Spacing, Hx);
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      X[FuncIndex] = X[FuncIndex] +
      (F1[FuncIndex] + (2-sqrt(2))*F2[FuncIndex] +  (2+sqrt(2))*F3[FuncIndex] + F4[FuncIndex]) / 6;
    *T += Spacing;
    for (FuncIndex = 1; FuncIndex <= NEQU; FuncIndex++ )
      (*DerivValues)[FuncIndex][Index + 1] = Vector[FuncIndex](*T, X);
//    printf("%e  ",*T);
//    for(i_print=1;i_print<=NEQU;i_print++) printf("%e  ",X[i_print]);
//    printf("\n");
  }
} /* procedure GetFirstFourValues */

long double fkt(int n)      /*  faktorial  */
{
  long double m = 1;
  if (n == 0) return 1;
  if (n < 0) {
    printf("Factorial is not possible to compute due to minus.\n");
    exit(0);
  }
  while (n > 1) {
    m *= n;
    n--;
  }
  return m;
}

long double an(double a, int n)
{
  int i;
  long double b;

  b = a;

  if (n == 0) return 1;
  n--;
  while (n--) { a *= b+1; b++; }

  return a;
}

double F(double a, double b, double c, double x)
{
  int n = 0;

  long double vh,v = 0;

  for (;;) {
   vh = v;
   v += an(a,n)*an(b,n)*pow(x,n)/fkt(n)/an(c,n);
   n++;
   if (fabs(vh - v) < 1e-7 ) break;
  }
  return v;
}

long double getK(double m)
{
  double fi,l,g,k;

  if (m <= .965925)
    return F(.5,.5,1,m*m)/4;
  else
    printf("m is too big. It is impossible to compute K\n");
  } 
}

double JacobiCN(double z, double m)
{
  long double pocatek[NEQU+1]= {0,0,1};
   
  return Adams(0,z,pocatek,NSTEPS,NSTEPS,&error);
}

void main() {
  int error = 0;
  long double pocatek[NEQU+1]= {0,0,1};

  printf("\n\n\n\n\n\n\n\n\n\n\n\n\n");
  Adams(0,1000*M_PI,pocatek,20,20,&error);
}
