(* NATURAL CUBIC SPLINE ALGORITHM 3.4
*                                     
*  To construct the cubic spline interpolant S for the function f,
*  defined at the numbers x(0) < x(1) < ... < x(n), satisfying
*  S''(x(0)) = S''(x(n)) = 0:
*
*  INPUT:   n; x(0), x(1), ..., x(N); either generate A(I) = f(x(I))
*           for I = 0, 1, ..., n or input A(I) for I = 0, 1, ..., n.
*
*  OUTPUT:  A(J), B(J), C(J), D(J) for J = 0, 1, ..., n - 1.
*
*  NOTE:    S(x) = A(J) + B(J) * ( x - x(J) ) + C(J) * ( x - x(J) )**2 +
*          D(J) * ( x - x(J) )**3 for x(J) <= x < x(J + 1)
*)
Input["This is the natural cubic spline interpolation.\n
\n
   Be ready to choose type of input method\n
\n
   Press 1 [enter] to continue\n"];
OK = 0;
While[OK == 0,
   FLAG=Input["   1. Input entry by entry from keyboard\n
   2. Input data from a text file\n
   3. Generate data using a function F with nodes entered from keyboard\n
   4. Generate data using a function F with nodes entered from a text file\n"];
   If[FLAG >= 1 && FLAG <=4,
      OK=1;
   ];
];
If[FLAG==1,
   OK = 0;
   While[OK == 0,
      n=Input["Input n\n"];
      If[n>0,
         OK = 1;
         For[i=0,
            i<=n,
            i++,
            X[i]=Input["Input X("<>ToString[i]<>")\n"];
            A[i]=Input["Input F(X("<>ToString[i]<>"))\n"];
         ],
         Input["Number must be a positive integer\n
         \n
         Press 1 [enter] to continue\n"];
      ];
   ];
];
If[FLAG == 2,
   AA=InputString["Has a text file been created with the data\n
            in two columns?\n\n
            Enter 'Yes' or 'No'\n"];
   If[AA == "y" || AA == "Y" || AA == "Yes" || AA == "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 n\n"];
         If[n > 0,    
            For[i = 0,
               i <= n,
               i++,
               X[i]=Read[INP,Number];
               A[i]=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 two column form\n
       with the X and F(X) values in\n 
      the corresponding columns.\n
   The program will end so the input file\n
   can be created.\n
      Press 1 [enter] to continue\n"];
      OK=0;
   ];
];
If[FLAG == 3,
   TEMP=Input["Input the function F(X) in terms of x.\n
   \n
   For example: Cos[x]\n"];
   F[x_] := Evaluate[TEMP];
   OK = 0;
   While[OK==0,
   n = Input["Input n"];
      If[n> 0,
         OK = 1;
         For[i = 0,
            i <= n,
            i++,
            X[i]=Input["Input X("<>ToString[i]<>")"];
            A[i]=F[X[i]];
         ];
         OK=1,
         Input["Number must be a positive integer.\n
         Enter 1 [enter] to continue\n"];
      ];
   ];
];  
If[FLAG == 4,
   AA=InputString["Has the text file\n
        with x-values been created?\n
   Enter 'Yes' or 'No'\n"];
   If[AA == "y" || AA == "Y" || 
           AA == "Yes" || AA == "yes",
      NAME = InputString["Input the file name\n
               In the form - name.ext\n
               For example:   output.dta\n"]; 
      INP=OpenRead[NAME];
      TEMP=Input["Input the function F(X) in terms of x.\n
      \n
      For example: Cos[x]\n"];
      F[x_] := Evaluate[TEMP];
      OK=0;
      While[OK == 0,
         n = Input["Input n\n"];
         If[n > 0, 
            OK=1;   
            For[i = 0,
               i <= n,
               i++,
               X[i]=Read[INP,Number];
               A[i]=F[X[i]];
            ];
            Close[INP],
            Input["Number must be a positive integer\n
            \n
            Press 1 [enter] to continue\n"];
         ];
      ],
      Input["The program will end so the input file\n
      can be created.\n
      Press 1 [enter] to continue\n"];
      OK=0;
   ];
];
If[OK == 1,
   M = n-1;
(*STEP 1*)
   For[i=0,
      i <= M,
      i++,
      H[i]=X[i+1]-X[i];
   ];
(*STEP 2, use XA in place of ALPHA*)
   For[i=1,
      i <= M,
      i++,
      XA[i]=3.0*(A[i+1]*H[i-1]-A[i]*(X[i+1]-X[i-1])+A[i-1]*H[i])/(H[i]*H[i-1]);
   ];
(*STEP 3*)
(* steps 3,4,5 and part of 6 solve the tridiagonal system using
   Algorithm 6.7 use XL, XU, XZ in place of L, MU, Z resp.*)
   XL[0]=1;
   XU[0]=0;
   XZ[0]=0;
(*STEP 4*)
   For[i=1,
      i <= M,
      i++,
      XL[i]=2.0*(X[i+1]-X[i-1])-H[i-1]*XU[i-1];
      XU[i]=H[i]/XL[i];
      XZ[i]=(XA[i]-H[i-1]*XZ[i-1])/XL[i];
   ];
(*STEP 5*)
   XL[n]=1;
   XZ[n]=0;
   c[n]=XZ[n];
(*STEP 6*)
   For[i=0,
      i <= M,
      i++,
      J=M-i;
      c[J]=XZ[J]-XU[J]*c[J+1];
      B[J]=(A[J+1]-A[J])/H[J]-H[J]*(c[J+1]+2.0*c[J])/3.0;
      d[J]=(c[J+1]-c[J])/(3.0*H[J]);
   ];
(*STEP 7*)
   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,"NATURAL CUBIC SPLINE INTERPOLATION\n"];
   Write[OUP,"The numbers X(0), ..., X(n) are:\n"];
   For[i=0,
      i<=n,
      i++,
      Write[OUP," \n",SetPrecision[X[i],9]];
   ];
   Write[OUP,"\n\nThe coefficients of the spline on the\n"];
   Write[OUP,"subintervals are:\n"];
   Write[OUP,"for i=0, ..., n-1\n"];
   Write[OUP,"   A(I)    B(I)   C(I)   D(I)  \n"];
   For[i=0,
      i <= M,
      i++,
      Write[OUP,SetPrecision[A[i],9]," ",SetPrecision[B[i],9]
         ," ",SetPrecision[c[i],9]," ",SetPrecision[d[i],9]];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];
Quit[];