(*HERMITE INTERPOLATION ALGORITHM 3.3
*
* TO OBTAIN THE COEFFICIENTS OF THE HERMITE INTERPOLATING
* POLYNOMIAL H ON THE (n+1) DISTINCT NUMBERS X(0), ..., X(n)
* FOR THE FUNCTION F:
*
* INPUT:  NUMBERS X(0), X(1), ...,X(n); VALUES F(X(0)), F(X(1)),
* 	  ..., F(X(n)) AND F'(X(0)), F'(X(1)), ..., F'(X(n)).
*
* OUTPUT: NUMBERS Q(0,0), Q(1,1), ..., Q(2n+1,2n+1) WHERE
*
*         H(X) = Q(0,0) + Q(1,1) * ( X-X(0) ) + Q(2,2) *
*                ( X - X(0) )**2 + Q(3,3) * ( X - X(0) )**2 *
*                ( X - X(1) ) + Q(4,4) * ( X - X(0) )**2 *
*                ( X - X(1) )**2 + ... + Q(2n + 1,2n + 1) *
*                ( X - X(0) )**2 * ( X - X(1) )**2 * ... *
*                ( X - X(n - 1) )**2 * (X - X(n) ).
*)
OK = 0;
While[OK == 0,
   FLAG = Input["This is Hermite interpolation\n
   \n
   Choice of input method\n
   1. Input entry by entry from keyboard\n
   2. Input data from a text file\n
   3. Generate data using a function\n"];
   If[FLAG == 1 || FLAG == 2 || FLAG == 3,
      OK = 1;
   ]; 
];
If[FLAG == 1,
   OK = 0;
   While[OK != 1,
      n = Input["Input the number of data points minus one"];
      If[n > 0,
         OK = 1;
         For[i = 0,
            i <= n,
            i++,
            X[i] = Input["Input X("<>ToString[i]<>")"];
            Q[2i,0] = Input["Input F(X("<>ToString[i]<>"))"];
	    Q[2i+1,1] = Input["Input F'(X("<>ToString[i]<>"))"];	
         ],
         Input["Number must be a positive integer.\n
         Enter 1 [enter] to continue\n"]
      ];   
   ];
];
If[FLAG == 2,
   A = InputString["Has a text file been created with the data\n
            in three columns?\n\n
            Enter 'Yes' or 'No'\n"];
   If[A == "y" || A == "Y" || A == "Yes" || A == "yes",
      NAME = InputString["Input the file name\n
               In the form - name.ext\n
               For example:   output.dta\n"]; 
      INP = OpenRead[NAME];
      OK = 0;
      While[OK == 0,
         n = Input["Input the number of data points minus one\n"];
         If[n > 0,    
            For[i = 0,
               i <= n,
               i++,
               X[i] = Read[INP,Number];
               Q[2i,0] = Read[INP,Number];
  	       Q[2i+1,1] = Read[INP,Number];
            ];
            Close[INP];
            OK = 1,
            Input["Number must be a positive integer\n
            \n
            Press 1 [enter] to continue\n"];
         ];
      ],
      Input["Please create the input file in three column form\n
      with the X , F(X) , and F'(X) values in the\n
      corresponding columns.\n
   The program will end so the input file\n
   can be created. Press 1 [enter].\n"];
      OK = 0;
   ];
]; 
If[FLAG == 3,
   TEMP = Input["Input the function F(X) in terms of x.\n
   \n
   For example: Sin[x]\n"];
   F[x_] := Evaluate[TEMP];
   DER = D[TEMP,x];
   FP[x_] := Evaluate[DER];
   OK = 0;
   While[OK == 0,
   n = Input["Input the number of data points minus one"];
      If[n > 0,
         OK = 1;
         For[i = 0,
            i <= n,
            i++,
            X[i] = Input["Input X("<>ToString[i]<>")"];
            Q[2i,0] = N[F[X[i]]];
            Q[2i+1,1] = N[FP[X[i]]];
         ];
         OK = 1,
         Input["Number must be a positive integer.\n
         Enter 1 [enter] to continue\n"];
      ];
   ];
]; 
If[OK == 1,
(*STEP 1*)
   For[i = 0,
      i <= n,
      i++,
(*STEP 2*)
      Z[2i] = X[i];
      Z[2i+1] = X[i];
      Q[2i+1,0] = Q[2i,0];
(*STEP 3*)
      If[i != 0,
         Q[2i,1] = N[(Q[2i,0]-Q[2i-1,0])/
                       (Z[2i]-Z[2i-1])];
      ];
   ];
(*STEP 4*)
   K = 2n+1;
   For[i= 2,
      i <= K,
      i++,
      For[J = 2,
         J <= i,
         J++,
         Q[i,J] = N[(Q[i,J-1]-Q[i-1,J-1])/
                      (Z[i]-Z[i-J])];
      ];
   ];
(*STEP 5*)
   FLAG = Input["Select output destination\n
                 1. Screen\n
                 2. Text file\n
                 Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
                    For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP,"HERMITE INTERPOLATING POLYNOMIAL\n"];
   Write[OUP,"Data input follows:\n"];
   Write[OUP," X   F(X)       F'(X)"];
   For[i = 0,
      i <= n,
      i++,
      Write[OUP,SetPrecision[X[i],9],"  ",SetPrecision[Q[2i,0],9],
            "  ",SetPrecision[Q[2i+1,1],9]];
   ];
   Write[OUP,"The coefficients of the Hermite Interpolation\n"];
   Write[OUP,"Polynomial\n"];
   Write[OUP,"in order of increasing exponent follow:\n"];
   For[i  = 0,
      i <= K,
      i++,
      Write[OUP,"\n",SetPrecision[Q[i,i],9]];
   ];
   A = InputString["Do you wish to evaluate this polynomial?\n\n
            Enter 'Yes' or 'No'\n"];
   If[A == "y" || A == "Y" || A == "Yes" || A == "yes",
      XXX = Input["Enter a point at which to evaluate\n"];
      S = Q[K,K]*(XXX-Z[K-1]);
      For[i = 2,
         i <= K,
         i++,
         J = K-i+1;
         S = N[(S+Q[J,J])*(XXX-Z[J-1])];
      ];
      S = S+Q[0,0];
      Write[OUP,"x-value is \n",SetPrecision[XXX,9]];
      Write[OUP,"with interpolated value \n",
            SetPrecision[S,9]];
   ];                   
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];