(* CLAMPED CUBIC SPLINE ALGORITHM 3.5
*
* 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)) = f'(x(0)) and S'(x(n)) = f'(x(n)):
*
* 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;
*          FP0 = f'(x(0)); FPN = f'(x(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 Clamped 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] = Evaluate[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 with x-values been created?\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];
      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] = Evaluate[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,
   FP0 = Input["Input F'(X(0))\n"];
   FPN = Input["Input F'(X(n))\n"];
];
If[OK == 1,
(*STEP 1*)
   M = n-1;
   For[i = 0,
      i <= M,
      i++,
      H[i] = X[i+1]-X[i];
   ];
(*STEP 2, use XA instead of ALPHA*)
   XA[0] = 3.0*(A[1]-A[0])/H[0]-3.0*FP0;
   XA[n] = 3.0*FPN-3.0*(A[n]-A[n-1])/H[n-1];
(*STEP 3*)
   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 4*)
(* steps 4,5,6, and part of 7 solve the tridiagonal system using
   Algorithm 6.7 *)
(* use XL, XU, XZ in place of L, MU, and Z resp. *)
   XL[0] = 2.0*H[0];
   XU[0] = 0.5;
   XZ[0] = XA[0]/XL[0];
(*STEP 5*)
   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 6*)
   XL[n] = H[n-1]*(2.0-XU[n-1]);
   XZ[n] = (XA[n]-H[n-1]*XZ[n-1])/XL[n];
   c[n] = XZ[n];
(*STEP 7*)
   For[i = 1,
      i <= n,
      i++,
      J = n-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 8*)
   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,"CLAMPED 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[];