(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input:: *)
(* SF6 als S+6  + 6 F- in Oh 14.02.12 *)
Clear[k1,k2,k3,k4,k5,k6,sig1,sig2,sig3,sig4,sig5,c,z,R1,R2,R3,R4,R5,S2,S3,S4];
c = {k1 -> 0.88, k2 -> 2.435, k5 -> 0.88, k6 -> 2.43,sig1 -> 0.3, 
sig2 -> 0.30, sig3 -> 0.30, sig4 -> 0.30, sig5 -> 0.302}; z=16.; z1=9.0;
ad = Sqrt[3./8.]; sq2=Sqrt[2];

(* S+6 He+Ne shell *)
T = 2.25*k1/R1^2+9.*k2/R2^2 /. c;
Vee=3.0*sig1/R1+12.*sig2/R2+16/(R1+R2)+24*ad/(R1+R2) /. c;
Vne=-3.0*z/R1-8.7*z/(R1+R2);
S2 = R2*4^(1/3);

(* F- He+Ne shell *)
Tf = 2.25*k5/R5^2+9.*k6/R6^2 /. c;
Veef=3.0*sig1/R5+12*sig2/R6+16/(R5+R6)+24*ad/(R5+R6) /. c;
Vnef=-3.0*z1/R5-8.3*z1/(R5+R6) /. c;
S5 = R6*4^(1/3);

T=T+6*Tf;
Vee=Vee+6*Veef;
Vne=Vne+6*Vnef;

(* S+6/6 F- *)
Vnn = (6*z*z1 + 12*z1^2/sq2+1.5*z1^2)/(S2+S5);
Vee = Vee+(6*10*10 + 12*10*10/sq2 + 1.5*10*10)/(S2+S5);
Vne = Vne -(6*z*10+6*z1*10+24*z1*10/sq2+3*z1*10)/(S2+S5);

func = T + Vnn + Vne + Vee;

t = FindMinimum[func, {R1,0.0809742}, {R2,0.546933},
    {R5,0.1491164},{R6,1.3151171},{Method -> Automatic}, {MaxIterations -> 500}]
    
N[Vne /. c /. t[[2]],10]
N[Vnn /. c /. t[[2]],10]
N[Vee /. c /. t[[2]],10]
N[-(Vee+Vne+Vnn)/T /. c /. t[[2]],10]
N[(S2+S5) /. c /. t[[2]],10]
N[0.529177*(S2+S5) /. c /. t[[2]],10]
N[ratio=-(Vne+2*Vnn)/(Vee-Vnn) /. c /. t[[2]],8]
d=S2+S5 /.c /.t[[2]]
(* xy plane *)
plot2=Graphics[{
      Circle[{0,0},R1],Circle[{0,0},S2],
      Circle[{d,0},R5],
      Circle[{d,0},S5],
      Circle[{-d,0},S5],
      Circle[{-d,0},R5],
      Circle[{0,d},R5],
      Circle[{0,d},S5],Circle[{0,-d},R5],
      Circle[{0,-d},S5]} ] /. t[[2]] /. c

Show[plot2,{AspectRatio -> Automatic,Axes -> True,
    GridLines -> Automatic, PlotRange \[Rule] {{-6,6},{-6,6}}, 
    Frame -> True}]