(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 9.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       157,          7]
NotebookDataLength[     25909,        593]
NotebookOptionsPosition[     24015,        531]
NotebookOutlinePosition[     24527,        550]
CellTagsIndexPosition[     24484,        547]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell["\<\
(*H2O4 water.  with k,sig of CH41401.pas and minim. ident. num. results \
10.01.2013 DEVELOPMENT! *)
Clear[z,sig1,sig2,sig4,k1,k2,k4,nc,R1,R2,R3,R4,w,p,vee,vne,vnn,
a1,a2,a3,b1,b2,b3,c1,c2,c3,cs,ss,d1,d2,d3,d4,d5,pi,i,j,d,t];
z=8.0; sig1=0.43; sig2=0.318; sig3=0.447; nc=6; pi=0.0;
k1=1.0; k2=1.34; k3=1.1; k4=1.8; bohr=0.529177; rad=57.29578;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{
  3.5668390630515456`*^9, {3.5668394752510695`*^9, 3.5668394888386936`*^9}, {
   3.566840005402401*^9, 3.5668400293172426`*^9}, {3.566840906678384*^9, 
   3.566840923027212*^9}, {3.5668417775723133`*^9, 3.5668417895999346`*^9}, {
   3.566841995894697*^9, 3.5668420196535387`*^9}, {3.566842057031204*^9, 
   3.5668420654552193`*^9}, {3.5668422519691467`*^9, 3.5668422563371544`*^9}, 
   3.5668422887072115`*^9, {3.566842337504097*^9, 3.566842360326937*^9}, {
   3.5668424117758274`*^9, 3.5668424218690453`*^9}, {3.566842627885008*^9, 
   3.566842684060706*^9}, {3.566842771030859*^9, 3.5668427769432697`*^9}, {
   3.5668428297025623`*^9, 3.5668428906206694`*^9}, {3.5668429241763277`*^9, 
   3.566842930135538*^9}, {3.5668971967082343`*^9, 3.5668971977222357`*^9}, {
   3.566897262025549*^9, 3.5668973340352755`*^9}, {3.5668973674349337`*^9, 
   3.5668973682305355`*^9}}],

Cell["\<\
(* C He shell, i=j each *)
Ekin = 2.25*k1/R1^2;
vee=3.0*sig1/R1;
vne=-3.0*z/R1;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.5668390682151546`*^9}}],

Cell["\<\
(* bonding pairs *)
Ekin = Ekin + 2.25*(2*k2/R2^2+k3/R3^2+k4/R4^2);
vee=vee+3.0*(2*sig2/R2+sig3/R3);
oc={-2,-2,-2,-2,-1,-1};
ch={z,1,1,0,0,0};
rr={R1,R2,R2,R3,R4,R4};
(*w=ArcCos[-1/3]/2; a=4*(R1+R2)/Sqrt[6];*)
cs=Cos[w/2];  ss=Sin[w/2];
(* d1=R1+R3;  d1=(R1+R2)*cs; *)  (* R1+R2+p is radius of outer sphere *) \
d1=0.0;
d1=(R1+R2+p)*cs; d2=(R1+R2+p)*ss; d3=(R1+R2)*cs;  d4=(R1+R2)*ss;  d5=R1+R3; \
d6=R4;
xn={0,-d1,-d1,0,0,0};
yn={0,d2,-d2,0,0,0};
zn={0,0,0,0,0,0};        
xc={0,-d3,-d3,d5,0,0};
yc={0,d4,-d4,0,0,0};
zc={0,0,0,0,d6,-d6};
vne=vne-2*(3-(p/R2)^2)/R2;
vnn = 0.0;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839078760773*^9}, {
  3.566839318626795*^9, 3.5668393238528037`*^9}, {3.566839506201524*^9, 
  3.5668395108503323`*^9}, {3.566839541207986*^9, 3.566839568664034*^9}, {
  3.566839684650237*^9, 3.56683970887708*^9}, {3.5668397799352045`*^9, 
  3.5668398201988754`*^9}, {3.566839877872177*^9, 3.5668398909606*^9}, {
  3.5668404106911125`*^9, 3.566840417539525*^9}, {3.5668404815464373`*^9, 
  3.566840691460406*^9}, {3.566840732519678*^9, 3.5668408608923035`*^9}, {
  3.566840930000425*^9, 3.5668409412168446`*^9}, {3.566841374686006*^9, 
  3.566841375388007*^9}, {3.5668421835006266`*^9, 3.5668421845770283`*^9}, {
  3.5668971767713995`*^9, 3.566897177676201*^9}}],

Cell["\<\
(*ww*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vee = vee + \
oc[[i]]*oc[[j]]/Sqrt[(xc[[i]]-xc[[j]])^2+(yc[[i]]-yc[[j]])^2+(zc[[i]]-zc[[j]])\
^2]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839083924382*^9}, {
  3.5668410500738354`*^9, 3.5668411059531336`*^9}, {3.5668412393245683`*^9, 
  3.5668412393713684`*^9}, {3.5668412896190567`*^9, 3.5668413118646955`*^9}, {
  3.566841703815384*^9, 3.5668417040025845`*^9}}],

Cell["\<\
(*nn*)
For[i = 1, i < 3, i++,
  For[j = i+1, j < 4, j++,
    vnn = vnn + \
ch[[i]]*ch[[j]]/Sqrt[(xn[[i]]-xn[[j]])^2+(yn[[i]]-yn[[j]])^2+(zn[[i]]-zn[[j]])\
^2]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.5668390890567913`*^9}, {
  3.566841524165468*^9, 3.566841530218279*^9}}],

Cell["\<\
(*nw*)
For[i = 1, i < nc+1, i++,
  For[j = 1, j < 4, j++,
    If[i != j,
      vne = vne + \
oc[[i]]*ch[[j]]/Sqrt[(xc[[i]]-xn[[j]])^2+(yc[[i]]-yn[[j]])^2+(zc[[i]]-zn[[j]])\
^2]]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.5668390979332066`*^9}, 
   3.5668412394337683`*^9, {3.566841289743857*^9, 3.566841311927096*^9}, 
   3.5668415580955276`*^9, {3.5668416231788425`*^9, 3.5668416346760626`*^9}}],

Cell[CellGroupData[{

Cell["\<\
Epot=vne+vee+vnn;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.19},{R2,1.06},{R3,1.1},{R4,0.78},{w,1.9},{p,0.56},{\
Method -> \"Newton\"}, {MaxIterations -> 500}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839104968819*^9}, {
   3.5668402974973135`*^9, 3.5668403782742558`*^9}, {3.5668414049812593`*^9, 
   3.5668414375229163`*^9}, 3.566841756605877*^9, {3.5668429499163733`*^9, 
   3.566842978277223*^9}, {3.566897388869372*^9, 3.566897408837407*^9}}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "73.04019471733528`"}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"R1", "\[Rule]", "0.19412445818729052`"}], ",", 
     RowBox[{"R2", "\[Rule]", "1.0638215912071483`"}], ",", 
     RowBox[{"R3", "\[Rule]", "1.117861147006825`"}], ",", 
     RowBox[{"R4", "\[Rule]", "0.8198885237720577`"}], ",", 
     RowBox[{"w", "\[Rule]", "1.915058475490735`"}], ",", 
     RowBox[{"p", "\[Rule]", "0.5621994829075847`"}]}], "}"}]}], 
  "}"}]], "Output",
 CellChangeTimes->{
  3.566841481015793*^9, 3.5668415755831585`*^9, 3.56684165016689*^9, 
   3.566841800098753*^9, 3.566841957081829*^9, 3.5668420294659557`*^9, 
   3.5668420712740297`*^9, 3.5668421450777593`*^9, 3.5668421925798426`*^9, 
   3.5668422607363625`*^9, 3.5668422938240204`*^9, {3.566842342714506*^9, 
   3.5668423649289455`*^9}, 3.566842426985854*^9, 3.5668427114231544`*^9, 
   3.566842786334486*^9, {3.5668428367849746`*^9, 3.566842859717015*^9}, 
   3.5668428959714785`*^9, 3.5668429348623466`*^9, 3.566842985234835*^9, {
   3.566897186973817*^9, 3.566897202683045*^9}, {3.566897267111158*^9, 
   3.5668973386216836`*^9}, 3.5668973723801427`*^9, 3.5668974129090137`*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
(* N[func, 14] *)
vne /. t[[2]]
vee /. t[[2]]
vnn /. t[[2]]
-Epot/Ekin /. t[[2]]
(R1+R2+p)*bohr /. t[[2]]
w*rad /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839112846833*^9}, {
  3.566842810779729*^9, 3.566842811091729*^9}}],

Cell[BoxData[
 RowBox[{"-", "193.10666484210472`"}]], "Output",
 CellChangeTimes->{
  3.566841481124993*^9, 3.566841575676759*^9, 3.5668416502760897`*^9, 
   3.566841800114353*^9, 3.566841957097429*^9, 3.566842029481556*^9, 
   3.56684207128963*^9, 3.566842145093359*^9, 3.566842192595443*^9, 
   3.5668422607363625`*^9, 3.5668422938396206`*^9, {3.5668423427301064`*^9, 
   3.5668423649289455`*^9}, 3.5668424270014544`*^9, 3.566842711438754*^9, 
   3.566842786334486*^9, {3.5668428367849746`*^9, 3.5668428597326145`*^9}, 
   3.5668428959714785`*^9, 3.566842934877947*^9, 3.566842985234835*^9, {
   3.566897186989417*^9, 3.5668972026986446`*^9}, {3.566897267111158*^9, 
   3.5668973386372833`*^9}, 3.566897372395743*^9, 3.566897412924614*^9}],

Cell[BoxData["37.899852799878495`"], "Output",
 CellChangeTimes->{
  3.566841481124993*^9, 3.566841575676759*^9, 3.5668416502760897`*^9, 
   3.566841800114353*^9, 3.566841957097429*^9, 3.566842029481556*^9, 
   3.56684207128963*^9, 3.566842145093359*^9, 3.566842192595443*^9, 
   3.5668422607363625`*^9, 3.5668422938396206`*^9, {3.5668423427301064`*^9, 
   3.5668423649289455`*^9}, 3.5668424270014544`*^9, 3.566842711438754*^9, 
   3.566842786334486*^9, {3.5668428367849746`*^9, 3.5668428597326145`*^9}, 
   3.5668428959714785`*^9, 3.566842934877947*^9, 3.566842985234835*^9, {
   3.566897186989417*^9, 3.5668972026986446`*^9}, {3.566897267111158*^9, 
   3.5668973386372833`*^9}, 3.566897372395743*^9, 3.566897412924614*^9}],

Cell[BoxData["9.126422607554025`"], "Output",
 CellChangeTimes->{
  3.566841481124993*^9, 3.566841575676759*^9, 3.5668416502760897`*^9, 
   3.566841800114353*^9, 3.566841957097429*^9, 3.566842029481556*^9, 
   3.56684207128963*^9, 3.566842145093359*^9, 3.566842192595443*^9, 
   3.5668422607363625`*^9, 3.5668422938396206`*^9, {3.5668423427301064`*^9, 
   3.5668423649289455`*^9}, 3.5668424270014544`*^9, 3.566842711438754*^9, 
   3.566842786334486*^9, {3.5668428367849746`*^9, 3.5668428597326145`*^9}, 
   3.5668428959714785`*^9, 3.566842934877947*^9, 3.566842985234835*^9, {
   3.566897186989417*^9, 3.5668972026986446`*^9}, {3.566897267111158*^9, 
   3.5668973386372833`*^9}, 3.566897372395743*^9, 3.566897412924614*^9}],

Cell[BoxData["1.9999999999999782`"], "Output",
 CellChangeTimes->{
  3.566841481124993*^9, 3.566841575676759*^9, 3.5668416502760897`*^9, 
   3.566841800114353*^9, 3.566841957097429*^9, 3.566842029481556*^9, 
   3.56684207128963*^9, 3.566842145093359*^9, 3.566842192595443*^9, 
   3.5668422607363625`*^9, 3.5668422938396206`*^9, {3.5668423427301064`*^9, 
   3.5668423649289455`*^9}, 3.5668424270014544`*^9, 3.566842711438754*^9, 
   3.566842786334486*^9, {3.5668428367849746`*^9, 3.5668428597326145`*^9}, 
   3.5668428959714785`*^9, 3.566842934877947*^9, 3.566842985234835*^9, {
   3.566897186989417*^9, 3.5668972026986446`*^9}, {3.566897267111158*^9, 
   3.5668973386372833`*^9}, 3.566897372395743*^9, 3.566897412924614*^9}],

Cell[BoxData["0.9631791523469879`"], "Output",
 CellChangeTimes->{
  3.566841481124993*^9, 3.566841575676759*^9, 3.5668416502760897`*^9, 
   3.566841800114353*^9, 3.566841957097429*^9, 3.566842029481556*^9, 
   3.56684207128963*^9, 3.566842145093359*^9, 3.566842192595443*^9, 
   3.5668422607363625`*^9, 3.5668422938396206`*^9, {3.5668423427301064`*^9, 
   3.5668423649289455`*^9}, 3.5668424270014544`*^9, 3.566842711438754*^9, 
   3.566842786334486*^9, {3.5668428367849746`*^9, 3.5668428597326145`*^9}, 
   3.5668428959714785`*^9, 3.566842934877947*^9, 3.566842985234835*^9, {
   3.566897186989417*^9, 3.5668972026986446`*^9}, {3.566897267111158*^9, 
   3.5668973386372833`*^9}, 3.566897372395743*^9, 3.566897412924614*^9}],

Cell[BoxData["109.72476909885255`"], "Output",
 CellChangeTimes->{
  3.566841481124993*^9, 3.566841575676759*^9, 3.5668416502760897`*^9, 
   3.566841800114353*^9, 3.566841957097429*^9, 3.566842029481556*^9, 
   3.56684207128963*^9, 3.566842145093359*^9, 3.566842192595443*^9, 
   3.5668422607363625`*^9, 3.5668422938396206`*^9, {3.5668423427301064`*^9, 
   3.5668423649289455`*^9}, 3.5668424270014544`*^9, 3.566842711438754*^9, 
   3.566842786334486*^9, {3.5668428367849746`*^9, 3.5668428597326145`*^9}, 
   3.5668428959714785`*^9, 3.566842934877947*^9, 3.566842985234835*^9, {
   3.566897186989417*^9, 3.5668972026986446`*^9}, {3.566897267111158*^9, 
   3.5668973386372833`*^9}, 3.566897372395743*^9, 3.566897412924614*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
(* Hvnfx: sum of H2-nuclei forces *)
Hvnfx = 0.0;
  For[j = 1, j < 6, j++,
    If[j != 2,
    Hvnfx = Hvnfx - \
(xn[[j]]-xn[[2]])*ch[[3]]*ch[[j]]/((xn[[2]]-xn[[j]])^2+(yn[[2]]-yn[[j]])^2+(\
zn[[2]]-zn[[j]])^2)^(3/2)]]

(* Hvefx: sum of H2 - clouds forces *)
Hvefx = 0.0;
  For[j = 1, j < nc+1, j++,
    If[j != 2,
    Hvefx = Hvefx- \
(xc[[j]]-xn[[2]])*oc[[j]]*ch[[2]]/((xc[[j]]-xn[[2]])^2+(yc[[j]]-yn[[2]])^2+(\
zc[[j]]-zn[[2]])^2)^(3/2)]]
Hforcetotx = Hvnfx+Hvefx/. t[[2]]
Hforcennx=Hvnfx /. t[[2]]
HForcenex=Hvefx /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839137994077*^9}}],

Cell[BoxData[
 RowBox[{"-", "0.5220699732010926`"}]], "Output",
 CellChangeTimes->{
  3.566841481717794*^9, 3.56684157626956*^9, 3.5668416508376913`*^9, 
   3.5668418001299534`*^9, 3.566841957113029*^9, 3.566842029497156*^9, 
   3.5668420713052297`*^9, 3.566842145108959*^9, 3.5668421926110425`*^9, 
   3.5668422607675624`*^9, 3.5668422938552203`*^9, {3.5668423427457066`*^9, 
   3.566842364960145*^9}, 3.5668424270170546`*^9, 3.5668427114543543`*^9, 
   3.566842786365686*^9, {3.5668428368161745`*^9, 3.5668428597482147`*^9}, 
   3.566842896002679*^9, 3.566842934893547*^9, 3.566842985266035*^9, {
   3.5668971870050173`*^9, 3.566897202714245*^9}, {3.566897267126758*^9, 
   3.5668973386528835`*^9}, 3.5668973724113426`*^9, 3.566897412940214*^9}],

Cell[BoxData[
 RowBox[{"-", "1.38980854519985`"}]], "Output",
 CellChangeTimes->{
  3.566841481717794*^9, 3.56684157626956*^9, 3.5668416508376913`*^9, 
   3.5668418001299534`*^9, 3.566841957113029*^9, 3.566842029497156*^9, 
   3.5668420713052297`*^9, 3.566842145108959*^9, 3.5668421926110425`*^9, 
   3.5668422607675624`*^9, 3.5668422938552203`*^9, {3.5668423427457066`*^9, 
   3.566842364960145*^9}, 3.5668424270170546`*^9, 3.5668427114543543`*^9, 
   3.566842786365686*^9, {3.5668428368161745`*^9, 3.5668428597482147`*^9}, 
   3.566842896002679*^9, 3.566842934893547*^9, 3.566842985266035*^9, {
   3.5668971870050173`*^9, 3.566897202714245*^9}, {3.566897267126758*^9, 
   3.5668973386528835`*^9}, 3.5668973724113426`*^9, 3.566897412940214*^9}],

Cell[BoxData["0.8677385719987575`"], "Output",
 CellChangeTimes->{
  3.566841481717794*^9, 3.56684157626956*^9, 3.5668416508376913`*^9, 
   3.5668418001299534`*^9, 3.566841957113029*^9, 3.566842029497156*^9, 
   3.5668420713052297`*^9, 3.566842145108959*^9, 3.5668421926110425`*^9, 
   3.5668422607675624`*^9, 3.5668422938552203`*^9, {3.5668423427457066`*^9, 
   3.566842364960145*^9}, 3.5668424270170546`*^9, 3.5668427114543543`*^9, 
   3.566842786365686*^9, {3.5668428368161745`*^9, 3.5668428597482147`*^9}, 
   3.566842896002679*^9, 3.566842934893547*^9, 3.566842985266035*^9, {
   3.5668971870050173`*^9, 3.566897202714245*^9}, {3.566897267126758*^9, 
   3.5668973386528835`*^9}, 3.5668973724113426`*^9, 3.566897412940214*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
(* Hvnfy: sum of H2-nuclei forces *)
Hvnfy = 0.0;
  For[j = 1, j < 6, j++,
    If[j != 2,
    Hvnfy = Hvnfy - \
(yn[[j]]-yn[[2]])*ch[[2]]*ch[[j]]/((xn[[2]]-xn[[j]])^2+(yn[[2]]-yn[[j]])^2+(\
zn[[2]]-zn[[j]])^2)^(3/2)]]

(* Hvefy: sum of H2 - clouds forces *)
Hvefy = 0.0;
  For[j = 1, j < nc+1, j++,
    If[j != 2,
    Hvefy = Hvefy - \
(yc[[j]]-yn[[2]])*oc[[j]]*ch[[2]]/((xc[[j]]-xn[[2]])^2+(yc[[j]]-yn[[2]])^2+(\
zc[[j]]-zn[[2]])^2)^(3/2)]]
Hforcetoty = Hvnfy+Hvefy /. t[[2]]
Hforcenny=Hvnfy /. t[[2]]
HForceney=Hvefy /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839148118495*^9}}],

Cell[BoxData["0.774612011809959`"], "Output",
 CellChangeTimes->{
  3.5668414820609946`*^9, 3.5668415766283607`*^9, 3.5668416512276917`*^9, 
   3.566841800145553*^9, 3.5668419571286287`*^9, 3.5668420295283566`*^9, 
   3.5668420713208294`*^9, 3.5668421451245594`*^9, 3.5668421926266427`*^9, 
   3.566842260783162*^9, 3.56684229387082*^9, {3.5668423427613063`*^9, 
   3.566842364975745*^9}, 3.5668424270326543`*^9, 3.5668427114699545`*^9, 
   3.566842786365686*^9, {3.5668428368161745`*^9, 3.566842859763815*^9}, 
   3.5668428960182786`*^9, 3.566842934909147*^9, 3.5668429852816353`*^9, {
   3.5668971870206175`*^9, 3.566897202729845*^9}, {3.566897267142358*^9, 
   3.5668973386684837`*^9}, 3.566897372426943*^9, 3.566897412955814*^9}],

Cell[BoxData["2.087580469045268`"], "Output",
 CellChangeTimes->{
  3.5668414820609946`*^9, 3.5668415766283607`*^9, 3.5668416512276917`*^9, 
   3.566841800145553*^9, 3.5668419571286287`*^9, 3.5668420295283566`*^9, 
   3.5668420713208294`*^9, 3.5668421451245594`*^9, 3.5668421926266427`*^9, 
   3.566842260783162*^9, 3.56684229387082*^9, {3.5668423427613063`*^9, 
   3.566842364975745*^9}, 3.5668424270326543`*^9, 3.5668427114699545`*^9, 
   3.566842786365686*^9, {3.5668428368161745`*^9, 3.566842859763815*^9}, 
   3.5668428960182786`*^9, 3.566842934909147*^9, 3.5668429852816353`*^9, {
   3.5668971870206175`*^9, 3.566897202729845*^9}, {3.566897267142358*^9, 
   3.5668973386684837`*^9}, 3.566897372426943*^9, 3.566897412955814*^9}],

Cell[BoxData[
 RowBox[{"-", "1.3129684572353089`"}]], "Output",
 CellChangeTimes->{
  3.5668414820609946`*^9, 3.5668415766283607`*^9, 3.5668416512276917`*^9, 
   3.566841800145553*^9, 3.5668419571286287`*^9, 3.5668420295283566`*^9, 
   3.5668420713208294`*^9, 3.5668421451245594`*^9, 3.5668421926266427`*^9, 
   3.566842260783162*^9, 3.56684229387082*^9, {3.5668423427613063`*^9, 
   3.566842364975745*^9}, 3.5668424270326543`*^9, 3.5668427114699545`*^9, 
   3.566842786365686*^9, {3.5668428368161745`*^9, 3.566842859763815*^9}, 
   3.5668428960182786`*^9, 3.566842934909147*^9, 3.5668429852816353`*^9, {
   3.5668971870206175`*^9, 3.566897202729845*^9}, {3.566897267142358*^9, 
   3.5668973386684837`*^9}, 3.566897372426943*^9, 3.566897412955814*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
(* Wvnfy: sum of w2-nuclei forces *)
Wvnfy = 0.0;
  For[j = 1, j < 6, j++,
    Wvnfy = Wvnfy - \
(yn[[j]]-zc[[2]])*oc[[2]]*ch[[j]]/((xc[[2]]-xn[[j]])^2+(yc[[2]]-yn[[j]])^2+(\
zc[[2]]-zn[[j]])^2)^(3/2)]

(* Wvefy: sum of W2 - clouds forces *)
Wvefy = 0.0;
  For[j = 1, j < nc+1, j++,
    If[j != 2,
    Wvefy = Wvefy - \
(yc[[j]]-zc[[2]])*oc[[j]]*oc[[2]]/((xc[[j]]-xc[[2]])^2+(yc[[j]]-yc[[2]])^2+(\
zc[[j]]-zc[[2]])^2)^(3/2)]]
Wforcetoty = Wvnfy+Wvefy /. t[[2]]
Wforcenny=Wvnfy /. t[[2]]
WForceney=Wvefy /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.5668391587889137`*^9}}],

Cell[BoxData["17.043532511185177`"], "Output",
 CellChangeTimes->{
  3.566841482419795*^9, 3.5668415769715614`*^9, 3.5668416515708923`*^9, 
   3.566841800176753*^9, 3.566841957144229*^9, 3.5668420295439563`*^9, 
   3.5668420713364296`*^9, 3.566842145140159*^9, 3.566842192642243*^9, 
   3.5668422607987623`*^9, 3.5668422939020205`*^9, {3.5668423427769065`*^9, 
   3.5668423649913454`*^9}, 3.5668424270482545`*^9, 3.5668427114855547`*^9, 
   3.5668427863812857`*^9, {3.5668428368317747`*^9, 3.566842859779415*^9}, 
   3.5668428960338783`*^9, 3.5668429349247465`*^9, 3.5668429853128357`*^9, {
   3.566897187036217*^9, 3.566897202745445*^9}, {3.566897267157958*^9, 
   3.566897338684084*^9}, 3.566897372442543*^9, 3.566897412971414*^9}],

Cell[BoxData["16.571054653852084`"], "Output",
 CellChangeTimes->{
  3.566841482419795*^9, 3.5668415769715614`*^9, 3.5668416515708923`*^9, 
   3.566841800176753*^9, 3.566841957144229*^9, 3.5668420295439563`*^9, 
   3.5668420713364296`*^9, 3.566842145140159*^9, 3.566842192642243*^9, 
   3.5668422607987623`*^9, 3.5668422939020205`*^9, {3.5668423427769065`*^9, 
   3.5668423649913454`*^9}, 3.5668424270482545`*^9, 3.5668427114855547`*^9, 
   3.5668427863812857`*^9, {3.5668428368317747`*^9, 3.566842859779415*^9}, 
   3.5668428960338783`*^9, 3.5668429349247465`*^9, 3.5668429853128357`*^9, {
   3.566897187036217*^9, 3.566897202745445*^9}, {3.566897267157958*^9, 
   3.566897338684084*^9}, 3.566897372442543*^9, 3.566897412971414*^9}],

Cell[BoxData["0.472477857333093`"], "Output",
 CellChangeTimes->{
  3.566841482419795*^9, 3.5668415769715614`*^9, 3.5668416515708923`*^9, 
   3.566841800176753*^9, 3.566841957144229*^9, 3.5668420295439563`*^9, 
   3.5668420713364296`*^9, 3.566842145140159*^9, 3.566842192642243*^9, 
   3.5668422607987623`*^9, 3.5668422939020205`*^9, {3.5668423427769065`*^9, 
   3.5668423649913454`*^9}, 3.5668424270482545`*^9, 3.5668427114855547`*^9, 
   3.5668427863812857`*^9, {3.5668428368317747`*^9, 3.566842859779415*^9}, 
   3.5668428960338783`*^9, 3.5668429349247465`*^9, 3.5668429853128357`*^9, {
   3.566897187036217*^9, 3.566897202745445*^9}, {3.566897267157958*^9, 
   3.566897338684084*^9}, 3.566897372442543*^9, 3.566897412971414*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
plot1=Graphics[{Circle[{xc[[1]],yc[[1]]},R1],      \
Circle[{xc[[2]],yc[[2]]},R2],Circle[{xc[[3]],yc[[3]]},R2],Circle[{xc[[4]],yc[[\
4]]},R3],Circle[{xc[[5]],yc[[5]]},R4],Disk[{xn[[1]],yn[[1]]},0.08], \
Disk[{xn[[2]],yn[[2]]},0.08], Disk[{xn[[3]],yn[[3]]},0.08]} ] /. t[[2]];

Show[plot1,{AspectRatio \[Rule] Automatic,Axes -> True,GridLines -> \
Automatic, PlotRange \[Rule] {{-3,3},{-3,3}}, Frame -> True}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668391703953342`*^9, 3.5668391835773573`*^9}, {
  3.5668412395117683`*^9, 3.566841239776969*^9}, {3.5668412898218565`*^9, 
  3.566841290087057*^9}, {3.56684186652367*^9, 3.566841909423745*^9}}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0, 0}, 0.19412445818729052`], 
   CircleBox[{-0.724001113542698, 1.0287131042112805`}, 1.0638215912071483`], 
   CircleBox[{-0.724001113542698, -1.0287131042112805`}, 1.0638215912071483`],
    CircleBox[{1.3119856051941154`, 0}, 1.117861147006825], 
   CircleBox[{0, 0}, 0.8198885237720577], DiskBox[{0, 0}, 0.08], 
   DiskBox[{-1.0475706750943732`, 1.4884641209788483`}, 0.08], 
   DiskBox[{-1.0475706750943732`, -1.4884641209788483`}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-3, 3}, {-3, 3}}]], "Output",
 CellChangeTimes->{
  3.566841482794196*^9, 3.566841577345962*^9, 3.5668416519608927`*^9, 
   3.5668418001923532`*^9, 3.566841957159829*^9, 3.566842029559556*^9, 
   3.56684207135203*^9, 3.566842145155759*^9, 3.5668421926578426`*^9, 
   3.5668422608143625`*^9, 3.5668422939176207`*^9, {3.566842342792506*^9, 
   3.5668423650069456`*^9}, 3.566842427063854*^9, 3.5668427115011544`*^9, 
   3.566842786396886*^9, {3.566842836847375*^9, 3.566842859810615*^9}, 
   3.5668428960494785`*^9, 3.5668429349403467`*^9, 3.566842985328436*^9, {
   3.566897187051817*^9, 3.5668972027610445`*^9}, {3.566897267173558*^9, 
   3.5668973386996837`*^9}, 3.5668973724581428`*^9, 3.566897412987014*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
plot2=Graphics[{Circle[{xc[[1]],zc[[1]]},R1],      \
Circle[{xc[[2]],zc[[2]]},R2],Circle[{xc[[3]],zc[[3]]},R2],Circle[{xc[[4]],zc[[\
4]]},R3],Circle[{xc[[5]],zc[[5]]},R4],Circle[{xc[[6]],zc[[6]]},R4],Disk[{xn[[\
1]],zn[[1]]},0.08], Disk[{xn[[2]],zn[[2]]},0.08], \
Disk[{xn[[3]],zn[[3]]},0.08]} ] /. t[[2]];

Show[plot2,{AspectRatio \[Rule] Automatic,Axes -> True,GridLines -> \
Automatic, PlotRange \[Rule] {{-3,3},{-3,3}}, Frame -> True}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668391703953342`*^9, 3.5668391881949654`*^9}, {
  3.5668412398393693`*^9, 3.5668412401201696`*^9}, {3.566841311459095*^9, 
  3.5668413117242956`*^9}, {3.5668419284557786`*^9, 3.5668419511070185`*^9}, {
  3.566842124953724*^9, 3.566842139087349*^9}}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0, 0}, 0.19412445818729052`], 
   CircleBox[{-0.724001113542698, 0}, 1.0638215912071483`], 
   CircleBox[{-0.724001113542698, 0}, 1.0638215912071483`], 
   CircleBox[{1.3119856051941154`, 0}, 1.117861147006825], 
   CircleBox[{0, 0.8198885237720577}, 0.8198885237720577], 
   CircleBox[{0, -0.8198885237720577}, 0.8198885237720577], 
   DiskBox[{0, 0}, 0.08], DiskBox[{-1.0475706750943732`, 0}, 0.08], 
   DiskBox[{-1.0475706750943732`, 0}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-3, 3}, {-3, 3}}]], "Output",
 CellChangeTimes->{
  3.566841482981396*^9, 3.566841577501962*^9, 3.5668416521168933`*^9, 
   3.5668418002079535`*^9, 3.5668419571754293`*^9, 3.566842029575156*^9, 
   3.5668420713832293`*^9, 3.566842145171359*^9, 3.566842192673443*^9, 
   3.5668422608299627`*^9, 3.566842293933221*^9, {3.5668423428081064`*^9, 
   3.5668423650225453`*^9}, 3.5668424270794544`*^9, 3.566842711516754*^9, 
   3.5668427864280863`*^9, {3.566842836878575*^9, 3.566842859826215*^9}, 
   3.566842896080679*^9, 3.566842934955947*^9, 3.5668429853440356`*^9, {
   3.566897187067417*^9, 3.5668972027766447`*^9}, {3.5668972672047577`*^9, 
   3.5668973387152834`*^9}, 3.5668973724893427`*^9, 3.566897413002614*^9}]
}, Open  ]],

Cell["", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668391703953342`*^9, 3.566839177774147*^9}}]
},
WindowToolbars->"EditBar",
WindowSize->{993, 964},
WindowMargins->{{Automatic, 354}, {3, Automatic}},
PrintingCopies->1,
PrintingPageRange->{Automatic, Automatic},
PrivateNotebookOptions->{"VersionedStylesheet"->{"Default.nb"[8.] -> False}},
FrontEndVersion->"9.0 for Microsoft Windows (64-bit) (November 20, 2012)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[557, 20, 1284, 21, 99, "Input"],
Cell[1844, 43, 199, 7, 82, "Input"],
Cell[2046, 52, 1338, 32, 320, "Input"],
Cell[3387, 86, 489, 12, 82, "Input"],
Cell[3879, 100, 329, 10, 82, "Input"],
Cell[4211, 112, 451, 12, 99, "Input"],
Cell[CellGroupData[{
Cell[4687, 128, 506, 10, 82, "Input"],
Cell[5196, 140, 1190, 22, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[6423, 167, 288, 11, 133, "Input"],
Cell[6714, 180, 741, 11, 31, "Output"],
Cell[7458, 193, 724, 10, 31, "Output"],
Cell[8185, 205, 723, 10, 31, "Output"],
Cell[8911, 217, 724, 10, 31, "Output"],
Cell[9638, 229, 724, 10, 31, "Output"],
Cell[10365, 241, 724, 10, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[11126, 256, 644, 21, 286, "Input"],
Cell[11773, 279, 747, 11, 31, "Output"],
Cell[12523, 292, 745, 11, 31, "Output"],
Cell[13271, 305, 730, 10, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[14038, 320, 646, 21, 286, "Input"],
Cell[14687, 343, 733, 10, 31, "Output"],
Cell[15423, 355, 733, 10, 31, "Output"],
Cell[16159, 367, 751, 11, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[16947, 383, 632, 20, 269, "Input"],
Cell[17582, 405, 733, 10, 31, "Output"],
Cell[18318, 417, 733, 10, 31, "Output"],
Cell[19054, 429, 732, 10, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[19823, 444, 675, 12, 116, "Input"],
Cell[20501, 458, 1287, 22, 379, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[21825, 485, 762, 14, 133, "Input"],
Cell[22590, 501, 1296, 23, 379, "Output"]
}, Open  ]],
Cell[23901, 527, 110, 2, 31, "Input"]
}
]
*)

(* End of internal cache information *)