(* Finding the minimum of the function W under the energy constraint H = E *)
(* Printing input parameters and checking their validity *)

(* Input:  x0, y0 - displacements
           omx, omy - frequencies
           h30,h21,h12,h03 - cubic anharmonicity for H
           h40,h31,h22,h13,h04 - 4th degree anharmonicity for H
           w30,w21,w12,w03 - cubic anharmonicity for W
           w40,w31,w22,w13,w04 - 4th degree anharmonicity for W
Example of input line: {x0, y0, omx, omy} = {0.1, 0.2, 0.8, 1.2};
                       {h30,h21,h12,h03}={-0.1,-0.3,0.2,0};
                       {h40,h31,h22,h13,h04}={0,0,0,0,0};
                       {w30,w21,w12,w03}={0.05,0.05,0.05,-0.05};
                       {w40,w31,w22,w13,w04}={0,0,0,0,0};
*)

(* Load a package for formatted printing of numbers *)
<<printing.txt;

fails=False;
Print["<TT>"];

Print["<P><TABLE BORDER=0 CELLSPACING=0 CELLPADDING=0><TR ALIGN=CENTER><TH COLSPAN=2>Frequencies</TH><TD ROWSPAN=3>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;</TD><TH COLSPAN=2>Displacements</TH></TR>"];
Print["<TR ALIGN=CENTER><TD><I>w</I><SUB><I>x</I></SUB></TD><TD><I>w</I><SUB><I>y</I></SUB></TD><TD><I>x</I><SUB>0</SUB></TD><TD><I>y</I><SUB>0</SUB></TD></TR>"];
Print["<TR ALIGN=CENTER><TD>&nbsp; ",printF[omx,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[omy,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[x0,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[y0,10,5]," &nbsp;</TD></TR></TABLE>"];
If[!(NumberQ[omx]&&NumberQ[omy]&&omx>0&&omy>0),
Print["<BR><FONT COLOR=\"#FF0000\">Frequencies must be positive numbers!</FONT>\n"];
fails=True];
If[!(NumberQ[x0]&&NumberQ[y0]),
Print["<BR><FONT COLOR=\"#FF0000\">Displacements must be numbers!</FONT>\n"];
fails=True];

Print["<P><TABLE BORDER=0 CELLSPACING=0 CELLPADDING=0><TR ALIGN=CENTER><TH COLSPAN=4>Cubic anharmonicity of <I>H</I></TH><TD ROWSPAN=6>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;</TD><TH COLSPAN=5>Quartic anharmonicity of <I>H</I></TH></TR>"];
Print["<TR ALIGN=CENTER><TD><I>h</I><SUB>30</SUB></TD><TD><I>h</I><SUB>21</SUB></TD><TD><I>h</I><SUB>12</SUB></TD><TD><I>h</I><SUB>03</SUB></TD><TD><I>h</I><SUB>40</SUB></TD><TD><I>h</I><SUB>31</SUB></TD><TD><I>h</I><SUB>22</SUB></TD><TD><I>h</I><SUB>13</SUB></TD><TD><I>h</I><SUB>04</SUB></TD></TR>"];
Print["<TR ALIGN=CENTER><TD>&nbsp; ",printF[h30,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[h21,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[h12,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[h03,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[h40,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[h31,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[h22,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[h13,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[h04,10,5]," &nbsp;</TD></TR>"];
Print["<TR ALIGN=CENTER><TH COLSPAN=4>Cubic anharmonicity of <I>W</I></TH><TH COLSPAN=5>Quartic anharmonicity of <I>W</I></TH></TR>"];
Print["<TR ALIGN=CENTER><TD><I>w</I><SUB>30</SUB></TD><TD><I>w</I><SUB>21</SUB></TD><TD><I>w</I><SUB>12</SUB></TD><TD><I>w</I><SUB>03</SUB></TD><TD><I>w</I><SUB>40</SUB></TD><TD><I>w</I><SUB>31</SUB></TD><TD><I>w</I><SUB>22</SUB></TD><TD><I>w</I><SUB>13</SUB></TD><TD><I>w</I><SUB>04</SUB></TD></TR>"];
Print["<TR ALIGN=CENTER><TD>&nbsp; ",printF[w30,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[w21,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[w12,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[w03,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[w40,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[w31,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[w22,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[w13,10,5]," &nbsp;</TD><TD>&nbsp; ",printF[w04,10,5]," &nbsp;</TD></TR>"];
Print["</TABLE>"];

If[!(NumberQ[h30]&&NumberQ[h21]&&NumberQ[h12]&&NumberQ[h03]),
   Print["<BR><FONT COLOR=\"#FF0000\">Parameters of cubic anharmonicity for <I>H</I> must be numbers!</FONT>\n"];
   fails=True];
If[!(NumberQ[h30]&&NumberQ[h21]&&NumberQ[h12]&&NumberQ[h03]),
   Print["<BR><FONT COLOR=\"#FF0000\">Parameters of quartic anharmonicity for <I>H</I> must be numbers!</FONT>\n"];
   fails=True];
If[!(NumberQ[w30]&&NumberQ[w21]&&NumberQ[w12]&&NumberQ[w03]),
   Print["<BR><FONT COLOR=\"#FF0000\">Parameters of cubic anharmonicity for <I>W</I> must be numbers!</FONT>\n"];
   fails=True];
If[!(NumberQ[w30]&&NumberQ[w21]&&NumberQ[w12]&&NumberQ[w03]),
   Print["<BR><FONT COLOR=\"#FF0000\">Parameters of quartic anharmonicity for <I>W</I> must be numbers!</FONT>\n"];
   fails=True];
If[fails,
   Print["<P>Computations terminated because input parameters are illegitimate.</TT><HR>"];
   Exit[]
  ];
(* Testing ends *)

If[!(omx<100&&omy<100),
Print["<BR><FONT COLOR=\"#0000FF\">Warning: one of frequencies is too large. Results may be unpredictable.</FONT>"];
];
If[!(omx>0.005&&omy>0.005),
Print["<BR><FONT COLOR=\"#0000FF\">Warning: one of frequencies is too small. Results may be unpredictable.</FONT>"];
];
If[!(Abs[x0]<100&&Abs[y0]<100),
Print["<BR><FONT COLOR=\"#0000FF\">Warning: one of displacements is too large. Results may be unpredictable.</FONT>"];
];
{h3max,h4max,w3max,w4max}={10,10,10,10};
If[!(Abs[h30]<h3max&&Abs[h21]<h3max&&Abs[h12]<h3max&&Abs[h03]<h3max),
Print["<BR><FONT COLOR=\"#0000FF\">Warning: one of parameters of cubic anharmonicity for <I>H</I> is too large. Results may be unpredictable.</FONT>"];
];
If[!(Abs[h40]<h4max&&Abs[h31]<h4max&&Abs[h22]<h4max&&Abs[h13]<h4max&&Abs[h04]<h4max),
Print["<BR><FONT COLOR=\"#0000FF\">Warning: one of parameters of quartic anharmonicity for <I>H</I> is too large. Results may be unpredictable.</FONT>"];
];
If[!(Abs[w30]<w3max&&Abs[w21]<w3max&&Abs[w12]<w3max&&Abs[w03]<w3max),
Print["<BR><FONT COLOR=\"#0000FF\">Warning: one of parameters of cubic anharmonicity for <I>W</I> is too large. Results may be unpredictable.</FONT>"];
];
If[!(Abs[w40]<w4max&&Abs[w31]<w4max&&Abs[w22]<w4max&&Abs[w13]<w4max&&Abs[w04]<w4max),
Print["<BR><FONT COLOR=\"#0000FF\">Warning: one of parameters of quartic anharmonicity for <I>W</I> is too large. Results may be unpredictable.</FONT>"];
];

ndig=128;
{x0, y0, omx, omy}=SetPrecision[{x0, y0, omx, omy},ndig];
{h30,h21,h12,h03}=SetPrecision[{h30,h21,h12,h03},ndig];
{h40,h31,h22,h13,h04}=SetPrecision[{h40,h31,h22,h13,h04},ndig];
{w30,w21,w12,w03}=SetPrecision[{w30,w21,w12,w03},ndig];
{w40,w31,w22,w13,w04}=SetPrecision[{w40,w31,w22,w13,w04},ndig];
Print["</TT>"];