(* BEZIER CURVE ALGORITHM 3.6
*
* To construct the cubic Bezier curves C0, ..., Cn-1 in
* parameter form, where Ci is represented by
*
* (xi(t),yi(t)) = ( a0(i) + a1(i)*t + a2(i)*t^2 + a3(i)*t^3,
*                   b0(i) + b1(i)*t + b2(i)*t^2 + b3(i)*t^3)
*
* for 0 <= t <= 1 as determined by the left endpoint (x(i),y(i)),
* left guidepoint (x+(i),y+(i)), right endpoint (x(i+1),y(i+1)) and
* right guidepoint (x-(i+1),y-(i+1)) for each i = 0, 1, ... , n-1;
*
* INPUT:  n, ( (x(i),y(i)), i = 0,...,n ),
*           ( (x+(i),y+(i)), i = 0,...,n-1 ),
*           ( (x-(i),y-(i)), i = 1,...,n ).
*
* OUTPUT: coefficients ( a0(i), a1(i), a2(i), a3(i),
*                       b0(i), b1(i), b2(i), b3(i), i = 0, ... , n-1 ).
*)
OK = 0;
While[OK == 0,
   FLAG=Input["This is the Bezier Curve Algorithm \n
   \n
   Choice of input method\n
   1. Input entry by entry from keyboard\n
   2. Input data from a text file\n"];
     If[FLAG == 1 || FLAG == 2,
      OK=1;
   ]; 
];
If[FLAG == 1,
   OK = 0;
   While[OK == 0,
      n = Input["Input n"];
      If[n > 0,
         OK = 1;
         X[0]=Input["Input X[0]\n"];
         Y[0]=Input["Input Y[0]\n"];
         XPL[0]=Input["Input X+[0]\n"];
         YPL[0]=Input["Input Y+[0]\n"];
         For[i=1,
            i<=n-1,
            i++,
            X[i]=Input["Input X("<>ToString[i]<>")\n"];
            Y[i]=Input["Input Y("<>ToString[i]<>")\n"];
            XMI[i-1]=Input["Input X-("<>ToString[i]<>")\n"];
            YMI[i-1]=Input["Input Y-("<>ToString[i]<>")\n"];
            XPL[i]=Input["Input X+("<>ToString[i]<>")\n"];
            YPL[i]=Input["Input Y+("<>ToString[i]<>")\n"];
         ];
         X[n]=Input["Input X("<>ToString[n]<>")\n"];
         Y[n]=Input["Input Y("<>ToString[n]<>")\n"];
         XMI[n-1]=Input["Input X-("<>ToString[n]<>")\n"];
         YMI[n-1]=Input["Input Y-("<>ToString[n]<>")\n"],
         Input["Number must be a positive integer.\n
         Enter 1 [enter] to continue\n"]
      ];   
   ];
]; 
If[FLAG == 2,
   A=InputString["Has the text file been created with \n
   the data as follows?\n
   \n
   Enter 'yes' or 'no'\n"];
   Print["The input file should look like this: \n"];
   Print["  X[0]   Y[0]   X+[0]   Y+[0]\n"];
   Print["  X[1]   Y[1]   X-[1]   Y-[1]   X+[1]   Y-[1]\n"];
   Print["  ...  \n"];
   Print[" X[n-1] Y[n-1] X-[n-1] Y-[n-1] X+[n-1] Y+[n-1] \n"];
   Print["  X[n]   Y[n]   X-[n]   Y-[n]\n"];
   Print["\n"];
   If[A == "Y" || A == "y" || 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 n\n"];
         If[n > 0,
            OK = 1;
            X[0]=Read[INP,Number];
            Y[0]=Read[INP,Number];
            XPL[0]=Read[INP,Number];
            YPL[0]=Read[INP,Number];
            For[i=1,
               i<=n-1,
               i++,
               X[i]=Read[INP,Number];
               Y[i]=Read[INP,Number];
               XMI[i-1]=Read[INP,Number];
               YMI[i-1]=Read[INP,Number];
               XPL[i]=Read[INP,Number];
               YPL[i]=Read[INP,Number];
            ];
            X[n]=Read[INP,Number];
            Y[n]=Read[INP,Number];
            XMI[n-1]=Read[INP,Number];
            YMI[n-1]=Read[INP,Number];
            Close[INP],
            Input["Number must be a positive integer\n
            \n
            Press 1 [enter] to continue\n"];
         ];
      ],
      Input["Please create the input file as indicated\n
      The program will end so the input file can be created.\n
      Press 1 [enter] to continue\n"];
      OK=0;
   ];
];   
If[OK == 1,
   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,"BEZIER CURVE ALGORITHM\n"];
   Write[OUP,"A0      A1      A2      A3     on the first line  \n"];
   Write[OUP,"B0      B1      B2      B3     on the second line  \n"];
(*STEP 1*)
   For[i=0,
      i <= n-1,
      i++,
(*STEP 2*)
      A0[i]=X[i];
      B0[i]=Y[i];
      A1[i]=3*(XPL[i]-X[i]);
      B1[i]=3*(YPL[i]-Y[i]);
      A2[i]=3*(X[i]+XMI[i]-2*XPL[i]);
      B2[i]=3*(Y[i]+YMI[i]-2*YPL[i]);
      A3[i]=X[i+1]-X[i]+3*XPL[i]-3*XMI[i];
      B3[i]=Y[i+1]-Y[i]+3*YPL[i]-3*YMI[i];
(*STEP 3*)
      Write[OUP,SetPrecision[A0[i],9],"     ",SetPrecision[A1[i],9],
"     ",SetPrecision[A2[i],9],"     ",SetPrecision[A3[i],9]];
      Write[OUP,SetPrecision[B0[i],9],"     ",SetPrecision[B1[i],9],
"     ",SetPrecision[B2[i],9],"     ",SetPrecision[B3[i],9]];
      Write[OUP,"\n"];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];
Quit[];