AppendTo[$Path," "];

(*   Before beginning evaluation of the notebook for the first time, replace the     *)
(*   blank space between the " marks above with the path to CosmicEvolution:         *)
(*                                                                                   *)
(*             Windows example:  "C:\\Folder\\Subfolder\\...\\CosmicEvolution"       *)
(*                                                                                   *)
(*   Unix, Linux, OS X example:  "Home Directory/Subdirectory/.../CosmicEvolution"   *)
(*                                                                                   *)
(*                When done replacing, evaluate this cell.                           *)

Print[StandardForm[StringForm["For a sample run with A > 0, A = 0, or A < 0, push the appropriate button:"]]]
GridBox[Transpose[{{ButtonBox["Atype = Apos", Active -> True]},
                   {ButtonBox["Atype = Azero", Active -> True]},
                   {ButtonBox["Atype = Aneg", Active -> True]}}]] // DisplayForm
Print[StandardForm[StringForm["Now  a) choose Evaluate Notebook from the Kernel -> Evaluation menu,"]]]
Print[StandardForm[StringForm["or   b) evaluate individual inputs manually."]]]

<< "nbdata.m"

<< "nbsolveABR0.m"

<< "nbProcEqS.m"

<< "nbIterateS.m"

(*  If S0 = Smax, the program jumps down to the first "Switch[S0, Smax, ...]" cell.  *)

Switch[S0, 1, Print[StringForm[" (searching for times of \!\(100\^th\) doubling, \!\(200\^th\)  doubling, and present epoch)"]], _, Null] 
Switch[S0, 1, troots = {FindRoot[Log[2, Ssol[t]] == 100, {t, 10^(-15)}],
                        FindRoot[Log[2, Ssol[t]] == 200, {t, 10^(-15)}],
                        FindRoot[Rsol[t] == N[R0, 6], {t, targettime/2}]}, _, Null]

Switch[S0, 1, t100x2 = t /. troots[[1]]; t200x2 = t /. troots[[2]]; t0 = t /. troots[[3]]; 
   Print[StringForm["Time of \!\(100\^th\) doubling = `` seconds = `` years",
                     t100x2, EngineeringForm[Convert[t100x2 Second, Year]/Year]]];
   Print[StringForm["Time of \!\(200\^th\) doubling = `` seconds = `` years",
                     t200x2, EngineeringForm[Convert[t200x2 Second, Year]/Year]]];
   Print[StringForm["Time of present epoch = `` seconds = `` years", t0,
                     EngineeringForm[Convert[t0 Second, Year]/Year]]], _, Null]

Switch[S0, 1, H100x2 = Hsol[t100x2]; H200x2 = Hsol[t200x2]; 
   Print[StringForm["`1` at \!\(100\^th\) doubling = `` = ``",
                     iH, H100x2, Convert[H100x2, Kilo Meter/(Mega Parsec Second)]]];
   Print[StringForm["`1` at \!\(200\^th\) doubling = `` = ``",
                     iH, H200x2, Convert[H200x2, Kilo Meter/(Mega Parsec Second)]]];
   Print[StringForm["`1` at present epoch = `` = ``",
                     iH, Hsol[t0], Convert[Hsol[t0], Kilo Meter/(Mega Parsec Second)]]], _, Null]

Switch[S0, 1, R100x2 = Rsol[t100x2]; R200x2 = Rsol[t200x2]; 
   Print[StringForm["`1` at \!\(100\^th\) doubling = ``", iR, R100x2]];
   Print[StringForm["`1` at \!\(200\^th\) doubling = ``", iR, R200x2]];
   Print[StringForm["`1` at present epoch = ``;     `3` at present epoch = `4` at present epoch/\!\(`5`\_min\) = ``", iR, Rsol[t0], iS, iR, iR, Rsol[t0]/N[Rmin,1000]]], _, Null]

Switch[S0, 1, 
   Print[StringForm["(`1` at \!\(100\^th\) doubling)/(`2` at \!\(H\_max\)) = ``",
                      iR, iR,  R100x2/N[RHmax, 1000]]];
   Print[StringForm["(`1` at \!\(200\^th\) doubling)/(`2` at \!\(H\_max\)) = ``",
                      iR, iR,  R200x2/N[RHmax, 1000]]];
   Print[StringForm["(`1` at present epoch)/(`2` at \!\(H\_max\)) = ``",
                      iR, iR,  Rsol[t0]/N[RHmax, 1000]]], _, Null]

Switch[S0, 1, Q100x2 = Qsol[t100x2]; Q200x2 = Qsol[t200x2]; 
   Print[StringForm["`1` at \!\(100\^th\) doubling = `2`(``) = ``",
                     iQ, iQ, N[t100x2], Q100x2]];
   Print[StringForm["`1` at \!\(200\^th\) doubling = `2`(``) = ``",
                     iQ, iQ, N[t200x2], Q200x2]];
 If[Min[Q100x2, Q200x2] >= 0,
   Print[StringForm[
   "Coasting (`1` ≈ 0):    From   `2` = `` seconds    to    `4` = `` seconds = `` years,  `7` < ``.",
              iQ, it, NumberForm[t100x2,3],
              it, NumberForm[t200x2,3], EngineeringForm[N[Convert[t200x2 Second,Year]/Year,6],3],
              iQ, NumberForm[Max[Q100x2, Q200x2],1]]],
   Print[StringForm["Decelerating (`1` < 0).", iQ]]], _, Null]

Switch[S0, 1, If[Atype == Apos,
                 Print[StringForm["Onward and upward to exponential acceleration (`1` ∼ 1):", iQ]],
                 Null, If[Min[Qsol[t0/2], Qsol[t0], Qsol[2 t0], Qsol[targettime]] >= 0,
                          Print[StringForm["More coasting:"]],
                          Print[StringForm["Coasting and decelerating (`1` < 0):",iQ]] ]];
   Print[StringForm["           `1`(\!\(t\_0\)/2) = `2`(``) = ``",
                                iQ, iQ, t0/2, Qsol[t0/2]]]; 
   Print[StringForm["              `1`(\!\(t\_0\)) = `2`(``) = ``",
                                   iQ, iQ, t0, Qsol[t0]]];
   Print[StringForm["           `1`(2 \!\(t\_0\)) = `2`(``) = ``",
                                iQ, iQ, 2 t0, Qsol[2 t0]]];
   Print[StringForm["  `1`(targettime) = `2`(``) = ``",
                       iQ, iQ, N[targettime, 6], NumberForm[N[Qsol[targettime],1000],6]]], _, Null]

SetOptions[ParametricPlot, PlotStyle -> Thickness[.003]];
SetOptions[ParametricPlot, PlotPoints -> 1000];
$TextStyle = {"Courier New", FontSize -> 18};

sgnS[a_, b_] := Sign[Ssol[b] - Ssol[a]];
interpolt[lambda_, a_, b_] := (1 - lambda) Log[2, 1 + a] + lambda Log[2, 1 + b];
interpolS[lambda_, a_, b_] := If[sgnS[a, b] == 1,
                                 (1 - lambda) Log[2, Ssol[a]] + lambda Log[2, Ssol[b]],
                                 (1 - lambda) Log[2, Ssol[b]] + lambda Log[2, Ssol[a]]];

Switch[S0, 1,
    g0 = ParametricPlot[{Log[2, 1 + t], Log[2, Ssol[t]]}, {t, 0, 10^(-13)}, 
                        AxesOrigin -> {0, 70},
                        AxesLabel -> {"\!\(log\_2\)(1 + t)", "\!\(log\_2\)(S(t))"}, 
                        ImageSize -> {1050, 650},
                        PlotRange -> {{0, interpolt[105/100, 0, 10^(-13)]},
                                      {70, interpolS[105/100, 0, 10^(-13)]}},
                        DisplayFunction -> Identity]; 
   g00 = Graphics[{Thickness[.004], Line[{{2 10^(-16), 0},                   (* Fills spurious *)
                                          {2 10^(-16), 91}}]}];              (* gap in graph.  *)
    P1 = Graphics[Point[{Log[2, 1 + 10^(-13)], Log[2, Ssol[10^(-13)]]}]];
  T01a = Graphics[Text[StyleForm[
                      "\!\(P\_1\):(\!\(log\_2\)(1 + \!\(10\^-13\)), \!\(log\_2\)(S(\!\(10\^-13\))))",
                                 FontWeight -> Bold],
                       Offset[{-170, 20}, {Log[2, 1 + 10^(-13)], Log[2, Ssol[10^(-13)]]}]]]; 
  T01b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 10^(-13)],6],4],        
                                                           NumberForm[N[Log[2,Ssol[10^(-13)]],6],4]],
                                 FontWeight -> Bold],
                       Offset[{-112, -25}, {Log[2,1 + 10^(-13)], Log[2,Ssol[10^(-13)]]}]]], _, Null];

Switch[S0, 1, Show[g0, g00, P1, T01a, T01b, DisplayFunction -> $DisplayFunction], _, Null]

Switch[S0, 1,
    g1 = ParametricPlot[{Log[2, 1 + t], Log[2, Ssol[t]]}, {t, 10^(-13), 1}, 
                        AxesOrigin -> {0, interpolS[-2/10, 10^(-13), 1]},
                        AxesLabel -> {"\!\(log\_2\)(1 + t)", "\!\(log\_2\)(S(t))"}, 
                        ImageSize -> {1050, 650},
                        PlotRange -> {{0, interpolt[105/100, 10^(-13), 1]},
                                      {interpolS[-2/10, 10^(-13), 1],
                                       interpolS[12/10, 10^(-13), 1]}},
                        DisplayFunction -> Identity];                     
   g10 = Graphics[{Thickness[.004], Line[{{.001, Log[2, Ssol[10^(-13)]]},       (* Fills spurious *)
                                          {.001, 134}}]}];                      (* gap in graph.  *) 
  T11a = Graphics[Text[StyleForm[
                      "\!\(P\_1\):(\!\(log\_2\)(1 + \!\(10\^-13\)), \!\(log\_2\)(S(\!\(10\^-13\))))",
                                 FontWeight -> Bold],
                       Offset[{195, 0}, {Log[2, 1 + 10^(-13)], Log[2, Ssol[10^(-13)]]}]]]; 
  T11b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 10^(-13)],6],4],        
                                                           NumberForm[N[Log[2,Ssol[10^(-13)]],6],4]],
                                 FontWeight -> Bold],
                       Offset[{160, -25}, {Log[2,1 + 10^(-13)], Log[2,Ssol[10^(-13)]]}]]];
    P2 = Graphics[Point[{Log[2, 1 + 1], Log[2, Ssol[1]]}]];
  T12a = Graphics[Text[StyleForm["\!\(P\_2\):(\!\(log\_2\)(1 + 1), \!\(log\_2\)(S(1))",
                                 FontWeight -> Bold],
                       Offset[{-132, 20}, {Log[2, 1 + 1], Log[2, Ssol[1]]}]]];
  T12b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 1],6],4],        
                                                           NumberForm[N[Log[2,Ssol[1]],6],4]],
                                 FontWeight -> Bold], 
                       Offset[{-83, -25}, {Log[2, 1 + 1], Log[2, Ssol[1]]}]]], _, Null];   

Switch[S0, 1, Show[g1, g10, P1, T11a, T11b, P2, T12a, T12b, DisplayFunction -> $DisplayFunction], _, Null]

Switch[S0, 1,
    g2 = ParametricPlot[{Log[2, 1 + t], Log[2, Ssol[t]]}, {t, 1, 10^3}, 
                        AxesOrigin -> {interpolt[-5/100, 1, 10^3],
                                       interpolS[-2/10, 1, 10^3]},
                        AxesLabel -> {"\!\(log\_2\)(1 + t)", "\!\(log\_2\)(S(t))"}, 
                        ImageSize -> {1050, 650},
                        PlotRange -> {{interpolt[-5/100, 1, 10^3],
                                       interpolt[105/100, 1, 10^3]},
                                      {interpolS[-2/10, 1, 10^3],
                                       interpolS[12/10, 1, 10^3]}},
                        DisplayFunction -> Identity]; 
  T21a = Graphics[Text[StyleForm["\!\(P\_2\):(\!\(log\_2\)(1 + 1), \!\(log\_2\)(S(1))",
                                 FontWeight -> Bold],
                       Offset[{132, -20}, {Log[2, 1 + 1], Log[2, Ssol[1]]}]]];
  T21b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 1],6],4],        
                                                           NumberForm[N[Log[2,Ssol[1]],6],4]],
                                 FontWeight -> Bold],                          
                       Offset[{105, -45}, {Log[2, 1 + 1], Log[2, Ssol[1]]}]]]; 
    P3 = Graphics[Point[{Log[2, 1 + 10^3], Log[2, Ssol[10^3]]}]];
  T22a = Graphics[Text[StyleForm["\!\(P\_3\):(\!\(log\_2\)(1 + 10^3), \!\(log\_2\)(S(10^3))",
                                 FontWeight -> Bold],
                       Offset[{-162, 45}, {Log[2, 1 + 10^3], Log[2, Ssol[10^3]]}]]];
  T22b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 10^3],6],4],        
                                                           NumberForm[N[Log[2,Ssol[10^3]],6],4]],
                                 FontWeight -> Bold], 
                       Offset[{-80, 20}, {Log[2, 1 + 10^3], Log[2, Ssol[10^3]]}]]], _, Null]; 

Switch[S0, 1, Show[g2, P2, T21a, T21b, P3, T22a, T22b, DisplayFunction -> $DisplayFunction], _, Null]

Switch[S0, 1,
    g3 = ParametricPlot[{Log[2, 1 + t], Log[2, Ssol[t]]}, {t, 10^3, 10^10}, 
                        AxesOrigin -> {interpolt[-5/100, 10^3, 10^10],
                                       interpolS[-2/10, 10^3, 10^10]},
                        AxesLabel -> {"\!\(log\_2\)(1 + t)", "\!\(log\_2\)(S(t))"}, 
                        ImageSize -> {1050, 650},
                        PlotRange -> {{interpolt[-5/100, 10^3, 10^10],
                                       interpolt[105/100, 10^3, 10^10]},
                                      {interpolS[-2/10, 10^3, 10^10],
                                       interpolS[12/10, 10^3, 10^10]}},
                        DisplayFunction -> Identity];
  T31a = Graphics[Text[StyleForm["\!\(P\_3\):(\!\(log\_2\)(1 + 10^3), \!\(log\_2\)(S(10^3))",
                                 FontWeight -> Bold],
                       Offset[{165, -20}, {Log[2, 1 + 10^3], Log[2, Ssol[10^3]]}]]];
  T31b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 10^3],6],4],        
                                                           NumberForm[N[Log[2,Ssol[10^3]],6],4]],
                                 FontWeight -> Bold], 
                       Offset[{107, -45}, {Log[2, 1 + 10^3], Log[2, Ssol[10^3]]}]]]; 
    P4 = Graphics[Point[{Log[2, 1 + 10^10], Log[2, Ssol[10^10]]}]];
  T32a = Graphics[Text[StyleForm["\!\(P\_4\):(\!\(log\_2\)(1 + 10^10), \!\(log\_2\)(S(10^10))",
                                 FontWeight -> Bold],
                       Offset[{-173, 45}, {Log[2, 1 + 10^10], Log[2, Ssol[10^10]]}]]];
  T32b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 10^10],6],4],        
                                                           NumberForm[N[Log[2,Ssol[10^10]],6],4]],
                                 FontWeight -> Bold], 
                       Offset[{-80, 20}, {Log[2, 1 + 10^10], Log[2, Ssol[10^10]]}]]], _, Null]; 

Switch[S0, 1, Show[g3, P3, T31a, T31b, P4, T32a, T32b, DisplayFunction -> $DisplayFunction], _, Null]

Switch[S0, 1,
    g4 = ParametricPlot[{Log[2, 1 + t], Log[2, Ssol[t]]}, {t, 10^10, t0}, 
                        AxesOrigin -> {interpolt[-5/100, 10^10, t0],
                                       interpolS[-2/10, 10^10, t0]},
                        AxesLabel -> {"\!\(log\_2\)(1 + t)", "\!\(log\_2\)(S(t))"}, 
                        ImageSize -> {1050, 650},
                        PlotRange -> {{interpolt[-5/100, 10^10, t0],
                                       interpolt[105/100, 10^10, t0]},
                                      {interpolS[-2/10, 10^10, t0],
                                       interpolS[12/10, 10^10, t0]}}, 
                        DisplayFunction -> Identity];
  T41a = Graphics[Text[StyleForm["\!\(P\_4\):(\!\(log\_2\)(1 + 10^10), \!\(log\_2\)(S(10^10))",
                                 FontWeight -> Bold],
                       Offset[{175, -20}, {Log[2, 1 + 10^10], Log[2, Ssol[10^10]]}]]];
  T41b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 10^10],6],4],        
                                                           NumberForm[N[Log[2,Ssol[10^10]],6],4]],
                                 FontWeight -> Bold], 
                       Offset[{103, -45}, {Log[2, 1 + 10^10], Log[2, Ssol[10^10]]}]]]; 
    P5 = Graphics[Point[{Log[2, 1 + t0], Log[2, Ssol[t0]]}]];
  T42a = Graphics[Text[StyleForm[
                             "\!\(P\_5\):(\!\(log\_2\)(1 + \!\(t\_0\)), \!\(log\_2\)(S(\!\(t\_0\)))",
                                 FontWeight -> Bold],
                       Offset[{-138, 45}, {Log[2, 1 + t0], Log[2, Ssol[t0]]}]]]; 
  T42b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + t0],6],4],        
                                                           NumberForm[N[Log[2,Ssol[t0]],6],4]],
                                 FontWeight -> Bold], 
                       Offset[{-80, 20}, {Log[2, 1 + t0], Log[2, Ssol[t0]]}]]], _, Null]; 

Switch[S0, 1, Show[g4, P4, T41a, T41b, P5, T42a, T42b, DisplayFunction -> $DisplayFunction], _, Null]

Switch[S0, 1,
    g5 = ParametricPlot[{Log[2, 1 + t], Log[2, Ssol[t]]}, {t, t0, 2 t0}, 
                        AxesOrigin -> {interpolt[-5/100, t0, 2 t0],
                                       interpolS[-2/10, t0, 2 t0]},
                        AxesLabel -> {"\!\(Log\_2\)(1 + t)", "\!\(Log\_2\)(S(t))"}, 
                        ImageSize -> {1050, 650},
                        PlotRange -> {{interpolt[-5/100, t0, 2 t0],
                                       interpolt[105/100, t0, 2 t0]},
                                      {interpolS[-2/10, t0, 2 t0],
                                       interpolS[12/10, t0, 2 t0]}}, 
                        DisplayFunction -> Identity];
  T51a = Graphics[Text[StyleForm[
                             "\!\(P\_5\):(\!\(log\_2\)(1 + \!\(t\_0\)), \!\(log\_2\)(S(\!\(t\_0\)))",
                                 FontWeight -> Bold],
                       Offset[{140, -20}, {Log[2, 1 + t0], Log[2, Ssol[t0]]}]]];
  T51b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + t0],6],4],        
                                                           NumberForm[N[Log[2,Ssol[t0]],6],4]],
                                 FontWeight -> Bold], 
                       Offset[{103, -45}, {Log[2, 1 + t0], Log[2, Ssol[t0]]}]]];
   P6 = Graphics[Point[{Log[2, 1 + 2 t0], Log[2, Ssol[2 t0]]}]];
 T52a = Graphics[Text[StyleForm[
                           "\!\(P\_6\):(\!\(log\_2\)(1 + 2\!\(t\_0\)), \!\(log\_2\)(S(2\!\(t\_0\)))",
                                FontWeight -> Bold],
                       Offset[{-148, 45}, {Log[2, 1 + 2 t0], Log[2, Ssol[2 t0]]}]]];
 T52b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 2 t0],6],4],        
                                                          NumberForm[N[Log[2,Ssol[2 t0]],6],4]],
                                FontWeight -> Bold], 
                      Offset[{-80, 20}, {Log[2, 1 + 2 t0], Log[2, Ssol[2 t0]]}]]], _, Null]; 

Switch[S0, 1, Show[g5, P5, T51a, T51b, P6, T52a, T52b, DisplayFunction -> $DisplayFunction], _, Null]

Switch[S0, 1,
    g6 = ParametricPlot[{Log[2, 1 + t], Log[2, Ssol[t]]}, {t, 2 t0, targettime}, 
                        AxesOrigin -> {interpolt[-5/100, 2 t0, targettime],
                                       interpolS[-2/10, 2 t0, targettime]},
                        AxesLabel -> {"\!\(Log\_2\)(1 + t)", "\!\(Log\_2\)(S(t))"}, 
                        ImageSize -> {1050, 650},
                        PlotRange -> {{interpolt[-5/100, 2 t0, targettime],
                                       interpolt[105/100, 2 t0, targettime]},
                                      {interpolS[-2/10, 2 t0, targettime],
                                       interpolS[12/10, 2 t0, targettime]}}, 
                        DisplayFunction -> Identity];
  T61a = Graphics[Text[StyleForm[
                           "\!\(P\_6\):(\!\(log\_2\)(1 + 2\!\(t\_0\)), \!\(log\_2\)(S(\!\(2t\_0\)))",
                                 FontWeight -> Bold],
                       Offset[{150, -20 + 65 (1 - sgnS[2 t0, targettime])/2},
                              {Log[2, 1 + 2 t0], Log[2, Ssol[2 t0]]}]]];
  T61b = Graphics[Text[StyleForm[StringForm[" = (``, ``)", NumberForm[N[Log[2,1 + 2 t0],6],4],        
                                                           NumberForm[N[Log[2,Ssol[2 t0]],6],4]],
                                 FontWeight -> Bold], 
                       Offset[{99, -45 + 65 (1 - sgnS[2 t0, targettime])/2},
                              {Log[2, 1 + 2 t0], Log[2, Ssol[2 t0]]}]]], _, Null]; 

Switch[S0, 1, Show[g6, P6, T61a, T61b, DisplayFunction -> $DisplayFunction], _, Null]

Switch[S0, 1,
  g000 = Graphics[{Thickness[.004], Line[{{2 10^(-16), 0},                  (* Fills spurious *)
                                          {2 10^(-16), 140}}]}];            (* gaps in graph. *)
    T1 = Graphics[Text[StyleForm["\!\(P\_1\)", FontWeight -> Bold],
                       Offset[{17, 0}, {Log[2, 1 + 10^(-13)], Log[2, Ssol[10^(-13)]]}]]];
    T2 = Graphics[Text[StyleForm["\!\(P\_2\)", FontWeight -> Bold],
                       Offset[{13, -13}, {Log[2, 1 + 1], Log[2, Ssol[1]]}]]];
    T3 = Graphics[Text[StyleForm["\!\(P\_3\)", FontWeight -> Bold],
                       Offset[{13, -13}, {Log[2, 1 + 10^3], Log[2, Ssol[10^3]]}]]];
    T4 = Graphics[Text[StyleForm["\!\(P\_4\)", FontWeight -> Bold],
                       Offset[{13, -13}, {Log[2, 1 + 10^10], Log[2, Ssol[10^10]]}]]];
    T5 = Graphics[Text[StyleForm["\!\(P\_5\)", FontWeight -> Bold],
                       Offset[{5, -15}, {Log[2, 1 + t0], Log[2, Ssol[t0]]}]]];
    T6 = Graphics[Text[StyleForm["\!\(P\_6\)", FontWeight -> Bold],
                       Offset[{-10, 15}, {Log[2, 1 + 2 t0], Log[2, Ssol[2 t0]]}]]];
    Show[g0, g000, g1, g2, g3, g4, g5, g6, P1, T1, P2, T2, P3, T3, P4, T4, P5, T5, P6, T6,
                   AxesOrigin -> {0, 0}, 
                   PlotRange -> {{0, (105/100)Log[2, 1 + targettime]},
                                 {0, (12/10)Log[2, Ssol[targettime]]}},
                   DisplayFunction -> $DisplayFunction], _, Null]

(*  If the graphs extend beyond the righthand margin, try lowering the
     magnification factor in  the Format -> Magnification menu.               *)

Switch[S0, 1,
   Print[StringForm["            \!\(log\_2\)(1 + ``) = \!\(log\_2\)(1 + \!\(t\_0\)) = ``", 
                     NumberForm[N[t0,6],3], Log[2, 1 + t0]]];
   Print[StringForm["                          
                     \!\(log\_2\)(`1`(\!\(t\_0\))) = ``", iS, Log[2, Ssol[t0]]]];
   Print[StringForm["         \!\(log\_2\)(1 + ``) = \!\(log\_2\)(1 + 2 \!\(t\_0\)) = ``",
                     NumberForm[N[2 t0,6],3], Log[2, 1 + 2 t0]]];
   Print[StringForm["                       
                     \!\(log\_2\)(`1`(2 \!\(t\_0\))) = ``", iS, Log[2, Ssol[2 t0]]]];
   Print[StringForm["\!\(log\_2\)(1 + ``) = \!\(log\_2\)(1 + targettime) = ``",
                     NumberForm[N[targettime,6],3], N[1 + targettime]]];
   Print[StringForm["              
                     \!\(log\_2\)(`1`(targettime)) = ``",
                     iS, N[Log[2, Ssol[targettime]]]]], _, Null]

(*  At this point you can choose to
      a) re-evaluate the "nbIterateS.m" cell to extend the interval of integration,
      b) re-evaluate other input cells,
      c) create and evaluate new input cells,
      d) delete the "Atype = ..." cell that appears after the "buttons" cell at the beginning,
          and then
            i) re-evaluate the notebook selecting anew the parameters Atype, Q0, S0, and targettime, or
           ii) modify the inputs H0 and Hmax in nbdata.m and re-evaluate the notebook,
      e) close the notebook,  or
      f) do something else.                                                                                                                       *)

Switch[S0, Smax, Print[StringForm[" (searching for time of bounce)"]], _, Null]
Switch[S0, Smax, tbounce = FindRoot[Ssol[t] == 1, {t, targettime/2},
                                    AccuracyGoal -> 50, PrecisionGoal -> 50,
                                    WorkingPrecision -> 100, MaxIterations -> 1000], _, Null]

Switch[S0, Smax, t1 = t /. tbounce[[1]];
   Print[StringForm["\!\(t\_1\) = ``   (seconds)", t1]];
   Print[StringForm["   = ``   (years)", EngineeringForm[Convert[t1 Second, Year]/Year]]];
   Print[StringForm["   = time between bounce and maximum expansion"]];  
   Print[StringForm["   = time between maximum expansion and next bounce"]];
   Print[StringForm["`` < \!\(t\_1\) < ``   (seconds);     `` < \!\(t\_1\) < ``   (years)",
                     N[Floor[t1 (10^6/targettime)] (targettime/10^6)],
                     N[Ceiling[t1 (10^6/targettime)] (targettime/10^6)],
                     EngineeringForm[Convert[N[Floor[t1 (10^6/targettime)] (targettime/10^6)]*
                                               Second, Year]/Year,7],
                     EngineeringForm[Convert[N[Ceiling[t1 (10^6/targettime)] (targettime/10^6)]*
                                               Second, Year]/Year,7] ]];
   Print[StringForm[""]];
   Print[StringForm["`1`(``) = ``", iS, N[Floor[t1 (10^6/targettime)] (targettime/10^6)],
                                        N[Ssol[N[Floor[t1 (10^6/targettime)] (targettime/10^6)]]]]];
   Print[StringForm["`1`(\!\(t\_1\)) = ``   (≈ 1)", iS,Ssol[t1] ]];
   Print[StringForm["`1`(``) = ``    (`4` < 0 because integration couldn't track sharp turnaround at bounce)",            iS,       N[Ceiling[t1 (10^6/targettime)] (targettime/10^6)],
                               N[Ssol[N[Ceiling[t1 (10^6/targettime)] (targettime/10^6)]]], iS]];
   Print[StringForm[""]];
   Print[StringForm["`1`(``) = `` = ``",
                     iH, N[Floor[t1 (10^6/targettime)] (targettime/10^6)],
                         N[Hsol[N[Floor[t1 (10^6/targettime)] (targettime/10^6)]]],
                         N[Convert[Hsol[N[Floor[t1 (10^6/targettime)] (targettime/10^6)]],
                                   (Kilo*Meter)/(Second*(Mega*Parsec))]] ]];
   Print[StringForm["`1`(``) = `` = ``",
                     iH, N[Ceiling[t1 (10^6/targettime)] (targettime/10^6)],
                         N[Hsol[N[Ceiling[t1 (10^6/targettime)] (targettime/10^6)]]],
                         N[Convert[Hsol[N[Ceiling[t1 (10^6/targettime)] (targettime/10^6)]],
                                   (Kilo*Meter)/(Second*(Mega*Parsec))]]]];
   Print[StringForm[""]];
   Print[StringForm[StringReplace["time between successive bounces = 2 \!\(t\_1\) = ``   (seconds) = ``   (years)",                    Whitespace -> " "], NumberForm[2 t1,3],
                                                      EngineeringForm[Convert[2 t1 Second,
                                                                              Year]/Year,3]]],
_, Null]

(*  At this point you can choose to
      a) re-evaluate the "nbIterateS.m" cell to extend the interval of integration,
      b) re-evaluate other input cells,
      c) create and evaluate new input cells,
      d) delete the "Atype = ..." cell that appears after the "buttons" cell at the beginning,
          and then
            i) re-evaluate the notebook selecting anew the parameters Atype, Q0, S0, and targettime, or
           ii) modify the inputs H0 and Hmax in nbdata.m and re-evaluate the notebook,
      e) close the notebook,  or
      f) do something else.                                                                                                                       *)

Print[StringForm["$Path = ``", $Path]]
savepath = $Path;
Print[StringForm["If the $Path variable as shown has unwanted or repeated appendages, you can remove the last n of them with the command $Path = Drop[$Path,-n] .  If you happen to drop more than intended, you can recover with the command $Path = savepath ."]]