(***********************  filename: nbsolveAposBR0.m  ************************)

(*****************************************************************************)
(*                                                                           *)
(*  This Mathematica program takes from data.m and the user inputs           *)
(*  compatible with A > 0 and computes in terms of the input parameters      *)
(*  H0, Q0, and Hmax the parameters A, B, R0 that satisfy the equations      *)
(*                                                                           *)
(*                          2                                                *)
(*                        H0    A   3 k    B                                 *)
(*                        --- = - - --- - --- ,                              *)
(*                         2    3     2     6                                *)
(*                        c         R0    R0                                 *)
(*                                                                           *)
(*                          2             3/2                                *)
(*                      Hmax    A   2 (-k)                                   *)
(*                      ----- = - + --------- ,                              *)
(*                        2     3       1/2                                  *)
(*                       c           3 B                                     *)
(*                                                                           *)
(*                                  c  (k R0  + B)                           *)
(*  and                    Q0 = 1 + -------------- ,                         *)
(*                                       2   6                               *)
(*                                     H0  R0                                *)
(*                                                                           *)
(*  and computes from them other parameters of interest.                     *) 
(*                                                                           *)
(*****************************************************************************)

$MaxPrecision = Infinity
$MaxExtraPrecision = Infinity

USimplify[x_] := Simplify[x, Assumptions -> {Meter > 0, Centimeter > 0,
                                             Second > 0, Year > 0}]

(*
t0 := present epoch in seconds

H0 := H(t0)
Q0 := Q(t0)
R0 := R(t0)

 X := H0/c
 Y := Hmax/c
 Z := Q0

k0 := k/R0^2
B0 := B/R0^6

X^2 = H0^2/c^2
    = A/3 - k/R0^2 - B/(3 R0^6)
    = A/3 - k0 - B0/3

Y^2 = Hmax^2/c^2
    = A/3 + 2 (-k)^(3/2)/(3 Sqrt[B])
    = A/3 + 2 (-k0)^(3/2)/(3 Sqrt[B0])

  Z = Q0
    = 1 + k/((H0^2/c^2) R0^2) + B/((H0^2/c^2) R0^6)
    = 1 + k0/X^2 + B0/X^2
*)

Clear[X, Y, Z, A, B, k0, R0, B0]

X = H0sec/c; Y = Hmaxsec/c; Z = Q0;

R0 = Sqrt[k/k0]

B0 = X^2 (Z - 1) - k0

A = 3 X^2 + 3 k0 + B0

B = B0 R0^6

(* The formula below for k0 was obtained from solving the X^2, Y^2, and Z    *)
(* equations in an ancillary Mathematica session.                            *)

k0 = USimplify[(3 (X^2 (Z + 2) - 3 Y^2)(X^2 (Z - 2) + Y^2) +
               Sqrt[3] Sqrt[(X^2 (Z + 2) - 3 Y^2)^3 (X^2 (3 Z - 2) - Y^2)])/
               (24 (X^2 - Y^2))]

RHmax = (B/(-k))^(1/4)

Q[R_] := 1 + c^2 (k R^4 + B)/(H[R]^2 R^6)

H[R_] := c Sqrt[(A R^6 - 3 k R^4 - B)/(3 R^6)]


Print[StringForm[" "]]

Print[StringForm["`1` = ``,     `3` = ``,     \!\(R\_0\) = ``",
                  iA, N[N[A,1000]], iB, N[N[B,1000]], N[N[R0,1000]] ]]

Print[StringForm[" "]]

Print[StringForm["\!\(c\^2\) `1` = `` = ``,     \!\(c\^2\) `4` = ``",
                  iA, N[N[c^2 A,1000]],
                         N[N[c^2 A (31557600 Second/Year)^2,1000]],
                  iB, N[N[c^2 B,1000]] ]]

lhs1 = N[N[USimplify[A Centimeter^2],1000]] 

if1 := If[lhs1 > 0, (* then1 *)

(* The formula for Rminsqr was obtained by solving the equation              *)
(*                                                                           *)
(*            (1/c^2) H^2(R) = (A R^6 - 3 k R^4 - B)/(3 R^6) = 0             *)
(*                                                                           *)
(* for R^2 in an ancillary Mathematica session.                              *)

Rminsqr = If[USimplify[(A^2 B + 4 k^3) Centimeter^6] >= 0,
             USimplify[(k/A) (1 + 2^(1/3) k/
                       (A^2 B + 2 k^3 + Sqrt[A^2 B (A^2 B + 4 k^3)])^(1/3) +
                       (A^2 B + 2 k^3 + Sqrt[A^2 B (A^2 B + 4 k^3)])^(1/3)/
                       (2^(1/3) k))],
             (k/A)(1 - 2 Cos[theta/3])];

theta = USimplify[Pi/2 - ArcTan[(A^2 B + 2 k^3)/Sqrt[-A^2 B (A^2 B + 4 k^3)]]];

Rmin = Sqrt[Rminsqr];

Print[StringForm[" "]];

Print[StringForm["\!\(R\_min\) = ``,     \!\(R\_\!\(H\_max\)\) = ``",
                  N[N[Rmin,1000]]  ,     N[N[RHmax,1000]]            ]];

(* The formula for RQminsqr was obtained by solving the equation             *)
(*                                                                           *)
(*       dQ/dR = -2 c^4 (A k R^6 + 3 A B R^2 - 4 B k)/(R^6 H^4(R)) = 0       *)
(*                                                                           *)
(* for R^2 in an ancillary Mathematica session.                              *)

RQminsqr = USimplify[(B/A)^(1/3)((2 + Sqrt[4 + A^2 B/k^3])^(1/3) +
                                 (2 - Sqrt[4 + A^2 B/k^3])^(1/3))];

RQmin = Sqrt[RQminsqr]; Qmin = Q[RQmin];

Print[StringForm[" "]];

Print[StringForm["\!\(R\_\!\(Q\_min\)\) = ``,     \!\(Q\_min\) = ``",
	            N[N[RQmin,1000]]          ,     N[N[Qmin,1000]]    ]];

Print[StringForm[" "]];

Print[StringForm[
"\!\(R\_\!\(H\_max\)\)/\!\(R\_min\) = ``,    \
 \!\(R\_\!\(Q\_min\)\)/\!\(R\_min\) = ``,     \!\(R\_0\)/\!\(R\_min\) = ``",
 N[RHmax/Rmin,6]                        ,    \
 N[RQmin/Rmin,6]                        ,     N[R0/Rmin,6]                  ]];

Print[StringForm[" "]];

Print[StringForm[StringReplace["`1`/\!\(R\_min\%2\) = 
                                curvature of space at time of bounce = ``",
                               Whitespace -> " "], ik, N[N[k/Rmin^2,1000]] ]];

Print[StringForm[" "]];

Print[StringForm[StringReplace["\!\(k\_0\) = `1`/\!\(R\_0\%2\)
                                = curvature of space at present epoch = ``",
                               Whitespace -> " "], ik, N[N[k0,1000]]        ]];


Print[StringForm[" "]];

newtimeused = TimeUsed[] - totaltimeused;

totaltimeused = totaltimeused + newtimeused;

Print[StringForm["Time used this run  = `` = ``",
                  newtimeused seconds, (newtimeused/60) minutes]];

Print[StringForm["Total time used = `` = ``",
                  totaltimeused seconds, (totaltimeused/60) minutes]];

Print[StringForm["Memory in use = ``",
                  N[MemoryInUse[]/2^20,7] megabytes]]

(* end of then1 *),

(* else1 *)
   Throw[Print[StringForm[" "]]; \
         Print[StringForm["Condition `1` > 0 not satisfied:  `2` = ``.",
                           iA, iA, lhs1 ]]; \

         Print[StringForm["Try increasing Q0." ]];

         Q0 = Input["Enter new Q0:  "];

         Print[StringForm[" "]];

         Print[StringForm["\!\(Q\_0\) = ``",
                           N[Q0]            ]];

         <<nbsolveAposBR0.m ] ]

Catch[if1]

(**********************  End of file nbsolveAposBR0.m  ***********************)
